- Main Page
- Classes
- Files

- File List

# C:/Users/Andreas Schiffler/Desktop/SVN/NewGamePhysics/NewGamePhysics/Physics/GravityEarthGfcModel.cs

    00001 ﻿// --------------------------------------------------------
    00002 // GravityEarthGfcModel.cs
    00003 // (A. Schiffler, 2009)
    00004 // --------------------------------------------------------
    00005
    00006 namespace NewGamePhysics.Physics
    00007 {
    00008     using System;
    00009     using System.Collections.Generic;
    00010     using System.IO;
    00011
    00012     using NewGamePhysics.Mathematics;
    00013     using NewGamePhysics.Utilities;
    00014
    00019     public class GravityEarthGfcModel
    00020     {
    00024         public const double Gal = 0.01;
    00025
    00029         public const double G = 6.67428E-11;
    00030
    00035         public const double GM = 0.3986004415E+15;
    00036
    00041         public const double R = 0.6378136300E+07;
    00042
    00047         public const double omega = 7.292115E-05;
    00048
    00053         private double[,] C_lm;
    00054
    00059         private double[,] S_lm;
    00060
    00064         private int modelOrder;
    00065
    00069         private int calcOrder;
    00070
    00074         public GravityEarthGfcModel()
    00075         {
    00076         }
    00077
    00083         public GravityEarthGfcModel(string modelFilename)
    00084         {
    00085             this.LoadGfcModel(modelFilename);
    00086         }
    00087
    00091         public int ModelOrder
    00092         {
    00093             get { return this.modelOrder; }
    00094         }
    00095
    00100         public int CalcOrder
    00101         {
    00102             get
    00103             {
    00104                 return this.calcOrder;
    00105             }
    00106
    00107             set
    00108             {
    00109                 if ((value > 0) && (value <= this.modelOrder))
    00110                 {
    00111                     this.calcOrder = value;
    00112                 }
    00113             }
    00114         }
    00115
    00126         public double Calculate(double lat, double lon, double h)
    00127         {
    00128             if (this.calcOrder == 0)
    00129             {
    00130                 throw new ApplicationException(
    00131                                    "Coefficient model not initialized. Load one first.");
    00132             }
    00133
    00134             if ((lat < -90.0) || (lat>90.0))
    00135             {
    00136                 throw new ArgumentOutOfRangeException("lat");
    00137             }
    00138
    00139             if ((lon < 0.0) || (lon > 360.0))
    00140             {
    00141                 throw new ArgumentOutOfRangeException("lon");
    00142             }
    00143
    00144             if ((h < 0.0) || (h > 100000.0))
    00145             {
    00146                 throw new ArgumentOutOfRangeException("h");
    00147             }
    00148
    00149             if (this.calcOrder<1)
    00150             {
    00151                 return 0.0;
    00152             }
    00153
    00154             // Scale input variables
    00155             double r = R + h;
    00156             double phi = Math.PI * lat / 180.0;
    00157             double lambda = Math.PI * lon / 180.0;
    00158
    00159             // Precalculations
    00160             double cos_phi = Math.Cos(phi);
    00161             double sin_phi = Math.Sin(phi);
    00162             double R_r = R / r;
    00163
    00164             // The components of the centrifugal potential
    00165             double cpr = omega * omega * r * cos_phi * cos_phi;
    00166             double cpl = 0.0;
    00167             double cpp = -omega * omega * r * r * cos_phi * sin_phi;
    00168
    00169             // The components of the gradient of the model potential
    00170             double war = 0.0;
    00171             double wal = 0.0;
    00172             double wap = 0.0;
    00173
    00174             // Generated normalized Legendre functions
    00175             double[][] P_lm;
    00176             double[][] dP_lm;
    00177             Legendre.NormalizedAssociatedFunctionAndDerivative(this.calcOrder + 1, sin_phi, out P_lm, out dP_lm);
    00178
    00179             // Sum over L
    00180             for (int l = 0; l <= this.calcOrder; l++)
    00181             {
    00182                 double war_inner = 0.0;
    00183                 double wal_inner = 0.0;
    00184                 double wap_inner = 0.0;
    00185
    00186                 // Sum over M
    00187                 for (int m = 0; m <= l; m++)
    00188                 {
    00189                     double m_lam = (double)m * lambda;
    00190                     double sin_m_lam = Math.Sin(m_lam);
    00191                     double cos_m_lam = Math.Cos(m_lam);
    00192
    00193                     double f1, f2;
    00194                     f1 = (S_lm[l, m] * cos_m_lam - C_lm[l, m] * sin_m_lam);
    00195                     f2 = (C_lm[l, m] * cos_m_lam + S_lm[l, m] * sin_m_lam);
    00196
    00197                     double s1, s2, s3;
    00198                     s1 = P_lm[l][m] * f2;
    00199                     s2 = (double)m * P_lm[l][m] * f1;
    00200                     s3 = dP_lm[l][m] * f2;
    00201
    00202                     war_inner += s1;
    00203                     wal_inner += s2;
    00204                     wap_inner += s3;
    00205                 }
    00206
    00207                 double R_r_l = Math.Pow(R_r, (double)l);
    00208                 war += (R_r_l  * (double)(l + 1) * war_inner);
    00209                 wal += (R_r_l * wal_inner);
    00210                 wap += (R_r_l * wap_inner);
    00211             }
    00212
    00213             double GM_r = GM / r;
    00214             war *= (- GM_r / r);
    00215             wal *= GM_r;
    00216             wap *= GM_r;
    00217
    00218             // Aggregate components
    00219             double c1 = (war + cpr);
    00220             double c2 = (wal + cpl) / (r * cos_phi);
    00221             double c3 = (wap + cpp) / r;
    00222
    00223             double g = Math.Sqrt(c1 * c1 + c2 * c2 + c3 * c3);
    00224
    00225             return g;
    00226         }
    00227
    00234         public void LoadGfcModel(string modelFilename)
    00235         {
    00236             // New list of coefficients
    00237             List<GravityFieldCoefficient> gfcCoefficients =
    00238                 new List<GravityFieldCoefficient>();
    00239
    00240             // Reset model order and coefficient arrays
    00241             this.modelOrder = 0;
    00242             this.C_lm = null;
    00243             this.S_lm = null;
    00244
    00245             // Prep data reader
    00246             Scanf scanf = new Scanf();
    00247             FileStream fileStream = null;
    00248             StreamReader streamReader = null;
    00249             try
    00250             {
    00251                 fileStream = new FileStream(modelFilename, FileMode.Open);
    00252                 streamReader = new StreamReader(fileStream);
    00253                 string line;
    00254                 while ((line = streamReader.ReadLine()) != null)
    00255                 {
    00256                     object[] o = scanf.Scan(line, "gfc    %i    %i  %lf  %lf ");
    00257
    00258                     if ((null != o) && (o.Length == 4))
    00259                     {
    00260                         int l = (int)o[0];
    00261
    00262                         // Track the bigges L as order
    00263                         if (l > this.modelOrder)
    00264                         {
    00265                             this.modelOrder = l;
    00266                         }
    00267
    00268                         GravityFieldCoefficient gcf = new GravityFieldCoefficient(l, (int)o[1], (double)o[2], (double)o[3]);
    00269                         gfcCoefficients.Add(gcf);
    00270                     }
    00271                 }
    00272             }
    00273             finally
    00274             {
    00275                 if (null != streamReader)
    00276                 {
    00277                     streamReader.Close();
    00278                 }
    00279                 if (null != fileStream)
    00280                 {
    00281                     fileStream.Close();
    00282                 }
    00283             }
    00284
    00285             // Repack coefficients into two-dimensional arrays
    00286             // for easy access in iterative calculations
    00287             int imax = modelOrder + 1;
    00288             this.C_lm = new double[imax, imax];
    00289             this.S_lm = new double[imax, imax];
    00290             foreach (GravityFieldCoefficient gfcCoefficient in gfcCoefficients)
    00291             {
    00292                 if ((gfcCoefficient.L < imax) && (gfcCoefficient.M < imax))
    00293                 {
    00294                     C_lm[gfcCoefficient.L, gfcCoefficient.M] = gfcCoefficient.C;
    00295                     S_lm[gfcCoefficient.L, gfcCoefficient.M] = gfcCoefficient.S;
    00296                 }
    00297             }
    00298
    00299             this.calcOrder = this.modelOrder;
    00300         }
    00301
    00305         public void UnloadGfcModel()
    00306         {
    00307             this.C_lm = null;
    00308             this.S_lm = null;
    00309             this.calcOrder = 0;
    00310         }
    00311     }
    00312 }