GPS Simulator

by
Hans - Henrik Benzon
Atmosphere Ionosphere Remote Sensing Division
Danish Meteorological Institute
e - mail hhb@dmi . dk

<< Geometry`Rotations`
Needs["Graphics`Shapes`"]
<< Miscellaneous`WorldPlot`
<< "Graphics`Graphics3D`"
<< Graphics`MultipleListPlot`
<< Graphics`Legend`
<< Statistics`DataManipulation`
<< Statistics`DescriptiveStatistics`
<< Graphics`

SetDirectory["//priv_1//work//hhb//LocalNotebooksForMathematica"]
Needs["GpsSimulatorV2`"]

ConstellationType = 1;
ViewFromLEOSat = 0;
latitude = 55;
longitude = 10;
CoorYesNo = 1;

PlotGpsAndLeoAndConstellation[ConstellationType, ViewFromLEOSat, 
      latitude, longitude, CoorYesNo]

<< /home/hhb/geomview/Geomview/mathematica/OOGL.m
Geomview[pl3D]

Get["MathGL3d`OpenGLViewer`"]

MathGL3d installed, Version 2.2
for Unix system Silicon Graphics
© by Jens-Peer Kuska 1997,98,99
with MathLink from Wolfram Research.

p3 = 
        Graphics3D[Sphere[6371, 2 10, 2 20], 
          Boxed -> True];
u = 0;
\(\(pl2 = 
        ParametricPlot3D[
\ \ \ \ \ \ \ \ {{0\ Rearth\ Cos[t]\ Cos[u], \ 
              0\ Rearth\ Sin[t]\ Cos[u], \ 0\ Rearth\ Sin[u]}, 
            RotaOcc1\  . \ {roccs\ Cos[t], roccs\ Sin[t], 50\ Abs[u]}, 
            Rota1\  . \ {rs\ Cos[t], rs\ Sin[t], 50\ Abs[u]}, 
            Rota2\  . \ Rota1\  . \ {rs\ Cos[t], rs\ Sin[t], 50\ Abs[u]}, 
            Rota2\  . \ Rota2\  . \ 
              Rota1\  . \ {rs\ Cos[t], rs\ Sin[t], 50\ Abs[u]}, 
            Rota2\  . \ Rota2\  . \ Rota2\  . \ 
              Rota1\  . \ {rs\ Cos[t], rs\ Sin[t], 50\ Abs[u]}, 
            Rota2\  . \ Rota2\  . \ Rota2\  . \ Rota2\  . \ 
              Rota1\  . \ {rs\ Cos[t], rs\ Sin[t], 50\ Abs[u]}, 
            Rota2\  . \ Rota2\  . \ Rota2\  . \ Rota2\  . \ Rota2\  . \ 
              Rota1\  . \ {rs\ Cos[t], rs\ Sin[t], 
                50\ Abs[u]}}, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {t, \ 0, \ 
            2  Pi}, PlotRange -> All, 
          AxesLabel -> {"km", "km", "km"}, Boxed -> False, 
          DisplayFunction -> Identity];\)\)
shadefunc[country_] := RGBColor[0, 1, 0]
f = 
        WorldPlot[{World, shadefunc}, WorldFrame -> None, 
          WorldGrid -> None, 
          WorldRange -> {{-90, 90}, {-180, 180}}, 
          WorldBackground -> RGBColor[0, 0, 1], 
          DisplayFunction -> Identity];
gl = 
        Show[Graphics[f], 
          PlotRange -> {{-10007.5, 10007.5}, Automatic}, 
          ImageSize -> 600, 
          DisplayFunction -> Identity];
MVClear[];
\(MVShow3D[
	Show[p3], 
	MVTexture -> Graphics[gl], 
	
      MVTextureMapType\  -> \ 
        MVSphereMapping, 
\ \ \ \ \ \ MVScaleTextureTarget\  -> \ {1, 1, 
          1}, 
\ \ \ \ \ \ MVRotateTextureTarget\  -> \ {{0, 0, 1}, Pi}, 
MVScaleTexture -> {1, 1}
]\)
\(MVShow3D[
	
      Show[{pl2, plcoorNEW1, plcoorNEW2, plcoorNEW3, leosat, gpssat1, 
          gpssat2, gpssat3, gpssat4, gpssat5, gpssat6, gpssat7, gpssat8, 
          gpssat9, gpssat10, gpssat11, gpssat12, gpssat13, gpssat14, 
          gpssat15, gpssat16, gpssat17, gpssat18, gpssat19, gpssat20, 
          gpssat21, gpssat22, gpssat23, gpssat24}], 
	
      MVLineTubeSize -> 2\ 50.0
]\)

MVPasteGraphics[]

MVClose[];

PlotAntennaPattern[5, 175];

ConstellationType = 1;
ViewFromLEOSat = 0;
latitude = 55;
longitude = 10;
ViewAngle = 175;
NrOfFrames = 50;
EarthRotateYesNo = 0;

AnimateConstellation[ConstellationType, ViewFromLEOSat, latitude, 
      longitude, ViewAngle, NrOfFrames, EarthRotateYesNo]

ConstellationType = 1;
NrOfLat = 181/(5 3);
NrOfLon = 361/(5 3);
iiLat = 4;
jjLon = 4;
TimeResolution = 70/7;
ViewAngle = 175;
EarthRotateYesNo = 1;

ConstellationCover[ConstellationType, NrOfLat, NrOfLon, iiLat, jjLon, 
      TimeResolution, ViewAngle, EarthRotateYesNo]

(* The average number (averaged over the time for one orbit) of GPS satellites as a function of the position on the earth is presented in the following 3D plot .  The ExtendGraphics`Geometry3D` package must be loaded before the plot can be generated. *)

SetDirectory[
      ToFileName[ {$TopDirectory,  "AddOns",  "Applications"}]]
Needs[ "ExtendGraphics`SurfaceGraphics3D`"]

plpos = 
        ListSurfacePlot3D[ListOfSats1, 
          AxesLabel  ->  {"X",  "Y", "Z"}, 
          ColorFunction -> Hue, 
          PlotLabel -> "Satellite coverage distribution", 
          Mesh -> True, Boxed -> False, 
          DisplayFunction -> Identity];
plneg = 
        ListSurfacePlot3D[ListOfSats2, 
          AxesLabel  ->  {"X",  "Y", "Z"}, Mesh -> True, 
          Boxed -> False, ColorFunction -> Hue, 
          DisplayFunction -> Identity];
pldis = 
        Show[{Graphics3D[plpos], Graphics3D[plneg]}, 
          PlotRange -> {{-7000, 7000}, {-7000, 
                7000}, {-7000, 7000}}, AspectRatio -> Automatic, 
          BoxRatios -> {1, 1, 1}];
ShowLegend[
      pldis, {Hue[If[(#) >  .9,  .9, (#)], 1, 1] &, 10, 
        ToString[Max[Column[ListOfSats1, 4]]], 
        ToString[Min[Column[ListOfSats1, 4]]], 
        LegendPosition -> {0.9, - .4}}, 
      ImageSize -> 600]

Max[Column[ListOfSats1, 4]]
Max[Column[ListOfSats2, 4]]
Min[Column[ListOfSats1, 4]]
Min[Column[ListOfSats2, 4]]

Mean[ListCoverAtOBS] // N
Variance[ListCoverAtOBS] // N

TextBarC = {};
AppendTo[
        TextBarC, "Nr of satellites in 12 intervals at (lat,lon)="];
AppendTo[TextBarC, ToString[OBSlatitude]];
AppendTo[TextBarC, ToString[OBSlongitude]];
freq  =  
        BinCounts[ListCoverAtOBS,  {0,  12.0,  1.00}];
midpoints = Table[1.00/2 + 1.00*(i - 1), {i, 12}];
\(BarChart[Transpose[{freq, \ midpoints}], PlotLabel\  -> 
TextBarC, 
      PlotRange -> All\ , ImageSize -> 450]\)