- Main Page
- Classes
- Files

- File List

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

    00001 ﻿// --------------------------------------------------------
    00002 // GravityMarsMroModel.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
    00020     public class GravityMarsMroModel
    00021     {
    00025         private const int LatMin = -90;
    00026
    00030         private const int LatMax = 90;
    00031
    00035         private const int LonMin = 0;
    00036
    00040         private const int LonMax = 360;
    00041
    00045         public const double Gal = 0.01;
    00046
    00050         public const double G = 6.67428E-11;
    00051
    00056         public const double GM = 42828.35796;
    00057
    00062         public const double R = 3397.0;
    00063
    00068         public const double omega = 7.088218081e-5;
    00069
    00074         internal double[,] map;
    00075
    00079         public GravityMarsMroModel()
    00080         {
    00081         }
    00082
    00088         public GravityMarsMroModel(string mapFilename)
    00089         {
    00090             this.LoadAnomalyMap(mapFilename);
    00091         }
    00092
    00102         public double Calculate(double lat, double lon, double h)
    00103         {
    00104             if (null == this.map)
    00105             {
    00106                 throw new ApplicationException(
    00107                     "Anomaly map not initialized. Load one first.");
    00108             }
    00109
    00110             if ((lat < -90.0) || (lat > 90.0))
    00111             {
    00112                 throw new ArgumentOutOfRangeException("lat");
    00113             }
    00114
    00115             if ((lon < 0.0) || (lon > 360.0))
    00116             {
    00117                 throw new ArgumentOutOfRangeException("lon");
    00118             }
    00119
    00120             if ((h < 0.0) || (h > 100000.0))
    00121             {
    00122                 throw new ArgumentOutOfRangeException("h");
    00123             }
    00124
    00125             // Calculate index position
    00126             int intLat = (int)Math.Truncate(lat);
    00127             int intLon = (int)Math.Truncate(lon);
    00128
    00129             // Create bicubic interpolator matrix
    00130             double[,] G = new double[4, 4];
    00131             for (int i = -1; i <= 2; i++)
    00132             {
    00133                 for (int j = -1; j <= 2; j++)
    00134                 {
    00135                     int indexLat = (LatWrap(intLat + i) - LatMin) % 180;
    00136                     int indexLon = (LonWrap(intLon + j) - LonMin) % 360;
    00137                     G[i + 1, j + 1] = map[indexLat, indexLon];
    00138                 }
    00139             }
    00140
    00141             // Interpolate over grid
    00142             Bicubic bicubic = new Bicubic(G);
    00143             double x = Math.Abs(lat - (double)intLat);
    00144             double y = lon - (double)intLon;
    00145             double g = bicubic.Calc(x, y);
    00146
    00147             // Scale from mgal units
    00148             // One milligal equals 0.01 mm/s/s.
    00149             g *= 1e-5;
    00150
    00151             // Add normal gravity
    00152             double g0 = GravityMarsNormalModel.Calculate(lat);
    00153             g += g0;
    00154
    00155             return g;
    00156         }
    00157
    00162         public void LoadAnomalyMap(string mapFilename)
    00163         {
    00164             // Load binary data
    00165             FileStream fileStream = new FileStream(mapFilename, FileMode.Open, FileAccess.Read);
    00166             BinaryReader binaryReader = new BinaryReader(fileStream);
    00167             long totalBytes = new FileInfo(mapFilename).Length;
    00168             byte[] byteArray = binaryReader.ReadBytes((Int32)totalBytes);
    00169             binaryReader.Close();
    00170             fileStream.Close();
    00171
    00172             // GGMRO Image format
    00173             // LINES                      = 180
    00174             // LINE_SAMPLES               = 360
    00175             // SAMPLE_TYPE                = IEEE_REAL
    00176             // SAMPLE_BITS                = 64
    00177             this.map = new double[180, 360];
    00178             int startIndex = 0;
    00179             byte[] tempArray = new byte[8];
    00180             for (int lat = 0; lat < 180; lat++)
    00181             {
    00182                 for (int lon = 0; lon < 360; lon++)
    00183                 {
    00184                     // Get sample
    00185                     // See: http://pds.jpl.nasa.gov/documents/sr/Chapter03.pdf
    00186                     // and http://pds.jpl.nasa.gov/documents/sr/stdref2003/AppendixC.pdf
    00187                     double sample;
    00188                     if (BitConverter.IsLittleEndian)
    00189                     {
    00190                         for (int i = 0; i < 8; i++)
    00191                         {
    00192                             tempArray[7 - i] = byteArray[startIndex + i];
    00193                         }
    00194
    00195                         sample = BitConverter.ToDouble(tempArray, 0);
    00196                     }
    00197                     else
    00198                     {
    00199                         sample = BitConverter.ToDouble(byteArray, startIndex);
    00200                     }
    00201
    00202                     // Data format
    00203                     // OFFSET                     = 0.0E+00
    00204                     // SCALING_FACTOR             = 1.0E+00
    00205                     // The values can be obtained by multiplying the sample in
    00206                     // the map by SCALING_FACTOR and then adding OFFSET.
    00207                     double value = (1.0 * sample + 0.0);
    00208                     this.map[lat, lon] = value ;
    00209
    00210                     // Update source index
    00211                     startIndex += (64 / 8);
    00212                 }
    00213             }
    00214         }
    00215
    00219         public void UnloadAnomalyMap()
    00220         {
    00221             this.map = null;
    00222         }
    00223
    00229         private int LatWrap(int lat)
    00230         {
    00231             while (lat < LatMin)
    00232             {
    00233                 lat =+ (LatMax - LatMin);
    00234             }
    00235
    00236             while (lat > LatMax)
    00237             {
    00238                 lat -= (LatMax - LatMin);
    00239             }
    00240
    00241             return lat;
    00242         }
    00243
    00249         private int LonWrap(int lon)
    00250         {
    00251             while (lon < LonMin)
    00252             {
    00253                 lon += (LonMax - LonMin);
    00254             }
    00255
    00256             while (lon > LonMax)
    00257             {
    00258                 lon -= (LonMax - LonMin);
    00259             }
    00260
    00261             return lon;
    00262         }
    00263     }
    00264 }
    00265