home *** CD-ROM | disk | FTP | other *** search
/ HAM Radio 1 / HamRadio.cdr / satel / orbit23 / orbit23.c < prev    next >
C/C++ Source or Header  |  1987-03-24  |  16KB  |  543 lines
  1. /* Copyright (C) 1986,1987 Robert W. Berger
  2.    May be used for non-commercial purposes without prior written permission. */
  3.  
  4. /* Change Log
  5.     3/24/87        v2.3 Adapted for new kepler.dat format.
  6.                  Allow beacon frequencies in mode.dat.
  7.                  Use decay rate for drag compensation.
  8.  
  9.     5/10/86        v2.2 Added single character aliases for satellite
  10.             names.
  11.  
  12.     4/30/86        v2.1 Print blank line if satellite dips below
  13.             horizon and reappears during same orbit and day
  14.  
  15.     4/29/86        v2.0 Changed GetSatelliteParams() to use AMSAT's
  16.             "kepler.dat" file. Moved schedule to "mode.dat" file.
  17.  
  18.         4/22/86         v1.3  Inserted N8FJB's suggestions for variable naming
  19.                         which maintain 8 character uniqueness.
  20.                         Also removed "include" file orbitr.h, which had two
  21.                         definitions of external functions defined in orbit.c
  22.                 -K3MC 
  23.  
  24.         4/1/86          v1.2  Corrected a scanf conversion to %d for an int
  25.                         type.    -K3MC
  26.  
  27.         3/19/86         v1.1  Changed GetSatelliteParams to not pass NULL
  28.                         to sscanf.
  29.  
  30.                                                                         */
  31.  
  32. #define DRAG 1
  33.  
  34. #include <stdio.h>
  35. #include <math.h>
  36.  
  37. extern double Kepler();
  38. extern double GetDay();
  39.  
  40.  
  41. #define PI 3.14159265
  42. #define PI2 (PI*2)
  43. #define MinutesPerDay (24*60.0)
  44. #define SecondsPerDay (60*MinutesPerDay)
  45. #define HalfSecond (0.5/SecondsPerDay)
  46. #define EarthRadius 6378.16             /* Kilometers           */
  47. #define C 2.997925e5                    /* Kilometers/Second    */
  48. #define DegreesPerRadian (180/PI)
  49. #define RadiansPerDegree (PI/180)
  50. #define ABS(x) ((x) < 0 ? (-(x)) : (x))
  51. #define SQR(x) ((x)*(x))
  52.  
  53. #define MaxModes 10
  54. typedef struct {
  55.                 int MinPhase,MaxPhase;
  56.                 char ModeStr[20];
  57.                }  ModeRec;
  58.  
  59. char VersionStr[] = "N3EMO Orbit Simulator  v2.3";
  60.  
  61.     /*  Keplerian Elements and misc. data for the satellite              */
  62.     double  EpochDay;                   /* time of epoch                 */
  63.     double EpochMeanAnomaly;            /* Mean Anomaly at epoch         */
  64.     long EpochOrbitNum;                 /* Integer orbit # of epoch      */
  65.     double EpochRAAN;                   /* RAAN at epoch                 */
  66.     double epochMeanMotion;             /* Revolutions/day               */
  67.     double OrbitalDecay;                /* Revolutions/day^2             */
  68.     double EpochArgPerigee;             /* argument of perigee at epoch  */
  69.     double Eccentricity;
  70.     double Inclination;
  71.     char SatName[100];
  72.     int ElementSet;
  73.     double BeaconFreq;                  /* Mhz, used for doppler calc    */
  74.     double MaxPhase;                    /* Phase units in 1 orbit        */
  75.     double perigeePhase;
  76.     int NumModes;
  77.     ModeRec Modes[MaxModes];
  78.     int PrintApogee;
  79.  
  80.     /* Simulation Parameters */
  81.  
  82.     double StartTime,EndTime, StepTime; /* In Days, 1 = New Year        */
  83.                                         /*      of reference year       */
  84.  
  85.     /* Site Parameters */
  86.     char SiteName[100];
  87.     double SiteLat,SiteLong,SiteAltitude,SiteMinElev;
  88.  
  89.  
  90. /* List the satellites in kepler.dat, and return the number found */
  91. ListSatellites()
  92. {
  93.     char str[100];
  94.     FILE *InFile;
  95.     char satchar;
  96.     int NumSatellites;
  97.  
  98.     printf("Available satellites:\n");
  99.  
  100.     if ((InFile = fopen("kepler.dat","r")) == 0)
  101.         {
  102.     printf("\"kepler.dat\" not found\n");
  103.     exit(-1);
  104.     }
  105.  
  106.     satchar = 'a';
  107.     NumSatellites = 0;
  108.     while (fgets(str,100,InFile))
  109.     if (sscanf(str,"Satellite: %s",str) == 1)
  110.         {
  111.         printf("    %c) %s\n",satchar++,str);
  112.         NumSatellites++;
  113.         }
  114.  
  115.     fclose(InFile);
  116.  
  117.     return NumSatellites;
  118. }
  119.  
  120. /* Match and skip over a string in the input file. Exits on failure. */
  121.  
  122. MatchStr(InFile,FileName,Target)
  123. FILE *InFile;
  124. char *FileName,*Target;
  125. {
  126.     char str[100];
  127.  
  128.     fgets(str,strlen(Target)+1,InFile);
  129.     if (strcmp(Target,str))
  130.        {
  131.        printf("%s: found \"%s\" while expecting \"%s\n\"",FileName,str,Target);
  132.        exit(-1);
  133.        }
  134. }
  135.  
  136. GetSatelliteParams()
  137. {
  138.     FILE *InFile;
  139.     char str[100];
  140.     int EpochYear;
  141.     double EpochHour,EpochMinute,EpochSecond;
  142.     int found;
  143.     int NumSatellites;
  144.     char satchar;
  145.  
  146.     NumSatellites = ListSatellites();
  147.  
  148.     found = 0;
  149.  
  150.     while (!found)
  151.     {
  152.     printf("Letter or satellite name :");
  153.     gets(SatName);
  154.  
  155.     if ((InFile = fopen("kepler.dat","r")) == 0)
  156.         {
  157.         printf("kepler.dat not found\n");
  158.         exit(-1);
  159.         }
  160.  
  161.     if (strlen(SatName) == 1)
  162.         {            /* use single character label */
  163.         satchar = SatName[0];
  164.         if (satchar < 'a' || satchar > 'a'+NumSatellites-1)
  165.             {
  166.             printf("'%c' is out of range\n",satchar);
  167.         fclose(InFile);
  168.         continue;
  169.         }
  170.  
  171.         do    
  172.         do  /* find line beginning with "Satellite: " */
  173.             fgets(str,100,InFile);
  174.         while (sscanf(str,"Satellite: %s",SatName) != 1);
  175.          while (satchar-- > 'a');
  176.         found = 1;
  177.         }
  178.         
  179.      else 
  180.          {
  181.          while (!found)  /* use satellite name */
  182.                 {
  183.             if (! fgets(str,100,InFile))
  184.                 break;    /* EOF */
  185.  
  186.             if (sscanf(str,"Satellite: %s",str) == 1
  187.             && strcmp(SatName,str) == 0)
  188.             found = 1;
  189.             }
  190.  
  191.         if (!found)
  192.         {
  193.         printf("Satellite %s not found\n",SatName);
  194.         fclose(InFile);
  195.         }
  196.         }
  197.     }
  198.  
  199.     BeaconFreq = 146.0;  /* Default value */
  200.  
  201.     fgets(str,100,InFile);    /* Skip line */
  202.  
  203.     MatchStr(InFile,"kepler.dat","Epoch time:");
  204.     fgets(str,100,InFile);
  205.     sscanf(str,"%lf",&EpochDay);
  206.  
  207.     EpochYear = EpochDay / 1000.0;
  208.     EpochDay -= EpochYear*1000.0;
  209.     EpochDay = GetDay(EpochYear,1,EpochDay);
  210.     fgets(str,100,InFile);
  211.  
  212.     if (sscanf(str,"Element set: %ld",&ElementSet) == 0)
  213.        {   /* Old style kepler.dat */
  214.        MatchStr(InFile,"kepler.dat","Element set:");
  215.        fgets(str,100,InFile);
  216.        sscanf(str,"%d",&ElementSet);
  217.        }
  218.  
  219.     MatchStr(InFile,"kepler.dat","Inclination:");
  220.     fgets(str,100,InFile);
  221.     sscanf(str,"%lf",&Inclination);
  222.     Inclination *= RadiansPerDegree;
  223.  
  224.     MatchStr(InFile,"kepler.dat","RA of node:");
  225.     fgets(str,100,InFile);
  226.     sscanf(str,"%lf",&EpochRAAN);
  227.     EpochRAAN *= RadiansPerDegree;
  228.  
  229.     MatchStr(InFile,"kepler.dat","Eccentricity:");
  230.     fgets(str,100,InFile);
  231.     sscanf(str,"%lf",&Eccentricity);
  232.  
  233.     MatchStr(InFile,"kepler.dat","Arg of perigee:");
  234.     fgets(str,100,InFile);
  235.     sscanf(str,"%lf",&EpochArgPerigee);
  236.     EpochArgPerigee *= RadiansPerDegree;
  237.  
  238.     MatchStr(InFile,"kepler.dat","Mean anomaly:");
  239.     fgets(str,100,InFile);
  240.     sscanf(str,"%lf",&EpochMeanAnomaly);
  241.     EpochMeanAnomaly *= RadiansPerDegree;
  242.  
  243.     MatchStr(InFile,"kepler.dat","Mean motion:");
  244.     fgets(str,100,InFile);
  245.     sscanf(str,"%lf",&epochMeanMotion);
  246.  
  247.     MatchStr(InFile,"kepler.dat","Decay rate:");
  248.     fgets(str,100,InFile);
  249.     sscanf(str,"%lf",&OrbitalDecay);
  250.  
  251.     MatchStr(InFile,"kepler.dat","Epoch rev:");
  252.     fgets(str,100,InFile);
  253.     sscanf(str,"%ld",&EpochOrbitNum);
  254.  
  255.     while (1)
  256.         {
  257.         if (! fgets(str,100,InFile))
  258.         break;    /* EOF */
  259.         if (strlen(str) <= 2)
  260.             break;  /* Blank line */
  261.         sscanf(str,"Beacon: %lf",&BeaconFreq);
  262.         }
  263.  
  264.     PrintApogee = (Eccentricity >= 0.3);
  265.  
  266.     perigeePhase = 0; MaxPhase = 360; /* Default values */
  267.     NumModes = 0;
  268.  
  269.     if ((InFile = fopen("mode.dat","r")) == 0)
  270.     return;
  271.  
  272.     found = 0;
  273.     while (!found)
  274.         {
  275.     if (! fgets(str,100,InFile))
  276.         break;    /* EOF */
  277.     if (sscanf(str,"Satellite: %s",str) == 1
  278.         && strcmp(SatName,str) == 0)
  279.         found = 1;
  280.     }
  281.     
  282.     if (found)
  283.     {
  284.     while (1)
  285.         {
  286.         if (! fgets(str,100,InFile))
  287.         break;    /* EOF */
  288.         if (strlen(str) <= 2)
  289.             break;  /* Blank line */
  290.         sscanf(str,"Beacon: %lf",&BeaconFreq);
  291.         sscanf(str,"Perigee phase: %lf",&perigeePhase);
  292.         sscanf(str,"Max phase: %lf",&MaxPhase);
  293.  
  294.         if (sscanf(str,"Mode: %20s from %d to %d",Modes[NumModes].ModeStr,
  295.         &Modes[NumModes].MinPhase,&Modes[NumModes].MaxPhase) == 3
  296.           && NumModes < MaxModes)
  297.           NumModes++;
  298.         }
  299.     fclose(InFile);
  300.     }
  301.  
  302.  
  303.  
  304. }
  305.  
  306.  
  307. GetSiteParams()
  308. {
  309.     FILE *InFile;
  310.     char str[100];
  311.  
  312.     printf("Site name :");
  313.     gets(str);
  314.     strcat(str,".sit");
  315.  
  316.     if ((InFile = fopen(str,"r")) == 0)
  317.         {
  318.         printf("%s not found\n",str);
  319.         exit(-1);
  320.         }
  321.  
  322.     fgets(SiteName,100,InFile);
  323.  
  324.     fgets(str,100,InFile);
  325.     sscanf(str,"%lf",&SiteLat);
  326.     SiteLat *= RadiansPerDegree;
  327.  
  328.     fgets(str,100,InFile);
  329.     sscanf(str,"%lf",&SiteLong);
  330.     SiteLong *= RadiansPerDegree;
  331.  
  332.     fgets(str,100,InFile);
  333.     sscanf(str,"%lf",&SiteAltitude);
  334.     SiteAltitude /= 1000;   /* convert to km */
  335.  
  336.     fgets(str,100,InFile);
  337.     sscanf(str,"%lf",&SiteMinElev);
  338.     SiteMinElev *= RadiansPerDegree;
  339. }
  340.  
  341. GetSimulationParams()
  342. {
  343.     double hour,duration;
  344.     int Month,Day,Year;
  345.  
  346.     printf("Start date (UTC) (Month Day Year) :");
  347.     scanf("%d%d%d",&Month,&Day,&Year);
  348.  
  349.     StartTime = GetDay(Year,Month,(double) Day);
  350.  
  351.     printf("Starting Hour (UTC) :");
  352.     scanf("%lf",&hour);
  353.     StartTime += hour/24;
  354.  
  355.     printf("Duration (Days) :");
  356.     scanf("%lf",&duration);
  357.     EndTime = StartTime + duration;
  358.  
  359.     printf("Time Step (Minutes) :");
  360.     scanf("%lf",&StepTime);
  361.     StepTime /= MinutesPerDay;
  362. }
  363.  
  364. PrintMode(OutFile,Phase)
  365. FILE *OutFile;
  366. {
  367.     int CurMode;
  368.  
  369.     for (CurMode = 0; CurMode < NumModes; CurMode++)
  370.         if ((Phase >= Modes[CurMode].MinPhase
  371.                 && Phase <= Modes[CurMode].MaxPhase)
  372.               || ((Modes[CurMode].MinPhase > Modes[CurMode].MaxPhase)
  373.                   && (Phase >= Modes[CurMode].MinPhase
  374.                         || Phase <= Modes[CurMode].MaxPhase)))
  375.             {
  376.             fprintf(OutFile,"    %s",Modes[CurMode].ModeStr);
  377.             break;
  378.             }
  379. }
  380.  
  381.  
  382. main()
  383. {
  384.     double ReferenceOrbit;      /* Floating point orbit # at epoch */
  385.     double CurrentTime;         /* In Days                         */
  386.     double CurrentOrbit;
  387.     double AverageMotion,       /* Corrected for drag              */
  388.         CurrentMotion;
  389.     double MeanAnomaly,TrueAnomaly;
  390.     double SemiMajorAxis;
  391.     double Radius;              /* From geocenter                  */
  392.     double RAANPrecession,PerigeePrecession;
  393.     double SSPLat,SSPLong;
  394.     long OrbitNum,PrevOrbitNum;
  395.     int Day,PrevDay;
  396.     double Azimuth,Elevation,Range;
  397.     double PrevRange,Velocity,Doppler;
  398.     int Phase;
  399.     char FileName[100];
  400.     FILE *OutFile;
  401.     int DidApogee;
  402.     double TmpTime,PrevTime;
  403.     int PrevVisible;
  404.  
  405.     printf("%s\n",VersionStr);
  406.  
  407.     GetSatelliteParams();
  408.     GetSiteParams();
  409.     GetSimulationParams();
  410.  
  411.     printf("Output file (RETURN for TTY) :");
  412.     gets(FileName);     /* Skip previous RETURN */
  413.     gets(FileName);
  414.  
  415.  
  416.     if (strlen(FileName) > 0)
  417.         {
  418.         if ((OutFile = fopen(FileName,"w")) == 0)
  419.             {
  420.             printf("Can't write to %s\n",FileName);
  421.             exit(-1);
  422.             }
  423.         }
  424.       else OutFile = stdout;
  425.  
  426.     fprintf(OutFile,"%s Element Set %d\n",SatName,ElementSet);
  427.  
  428.     fprintf(OutFile,"%s\n",SiteName);
  429.  
  430.     fprintf(OutFile,"Doppler calculated for freq = %lf MHz\n",BeaconFreq);
  431.  
  432.     SemiMajorAxis = 331.25 * exp(2*log(MinutesPerDay/epochMeanMotion)/3);
  433.     GetPrecession(SemiMajorAxis,Eccentricity,Inclination,&RAANPrecession,
  434.                         &PerigeePrecession);
  435.  
  436.     ReferenceOrbit = EpochMeanAnomaly/PI2 + EpochOrbitNum;
  437.  
  438.     PrevDay = -10000; PrevOrbitNum = -10000;
  439.     PrevTime = StartTime-2*StepTime;
  440.  
  441.     BeaconFreq *= 1E6;          /* Convert to Hz */
  442.  
  443.     DidApogee = 0;
  444.  
  445.     /* Start the loop one step early, to init OldRange */
  446.     for (CurrentTime = StartTime-StepTime; CurrentTime <= EndTime;
  447.                 CurrentTime += StepTime)
  448.         {
  449.  
  450. #if DRAG
  451.         AverageMotion = epochMeanMotion
  452.        + (CurrentTime-EpochDay)*OrbitalDecay/2;
  453.         CurrentMotion = epochMeanMotion
  454.        + (CurrentTime-EpochDay)*OrbitalDecay;
  455. #else
  456.         AverageMotion = epochMeanMotion;
  457.         CurrentMotion = epochMeanMotion;
  458. #endif
  459.         SemiMajorAxis = 331.25 * exp(2*log(MinutesPerDay/CurrentMotion)/3);
  460.  
  461.         CurrentOrbit = ReferenceOrbit +
  462.                         (CurrentTime-EpochDay)*AverageMotion;
  463.         OrbitNum = CurrentOrbit;
  464.  
  465.         MeanAnomaly = (CurrentOrbit-OrbitNum)*PI2;
  466.  
  467.         TmpTime = CurrentTime;
  468.         if (MeanAnomaly < PI)
  469.             DidApogee = 0;
  470.         if (PrintApogee && !DidApogee && MeanAnomaly > PI)
  471.             {                   /* Calculate Apogee */
  472.             TmpTime -= StepTime;   /* So we pick up later where we left off */
  473.             MeanAnomaly = PI;
  474.             CurrentTime=EpochDay+(OrbitNum-ReferenceOrbit+0.5)/AverageMotion;
  475.             }
  476.  
  477.         TrueAnomaly = Kepler(MeanAnomaly,Eccentricity);
  478.         Radius = SemiMajorAxis*(1-SQR(Eccentricity))
  479.                         / (1+Eccentricity*cos(TrueAnomaly));
  480.  
  481.         GetSubSatPoint(EpochDay,EpochRAAN,EpochArgPerigee,
  482.            Inclination,Eccentricity,RAANPrecession,PerigeePrecession,
  483.            CurrentTime,TrueAnomaly,&SSPLat,&SSPLong);
  484.         GetBearings(SiteLat,SiteLong,SiteAltitude,SSPLat,SSPLong,Radius,
  485.             &Azimuth,&Elevation,&Range);
  486.  
  487.         Velocity = (Range-PrevRange)/((CurrentTime-PrevTime)*SecondsPerDay);
  488.         PrevRange = Range;
  489.  
  490.         if (Elevation >= SiteMinElev && CurrentTime >= StartTime)
  491.             {
  492.             Day = CurrentTime + HalfSecond;
  493.             if (((double) Day) > CurrentTime+HalfSecond)
  494.                 Day -= 1;    /* Correct for truncation of negative values */
  495.  
  496.         if (OrbitNum == PrevOrbitNum && Day == PrevDay && !PrevVisible)
  497.             fprintf(OutFile,"\n");    /* Dipped out of sight; print blank */
  498.  
  499.             if (OrbitNum != PrevOrbitNum || Day != PrevDay)
  500.                 {                       /* Print Header */
  501.                 PrintDate(OutFile,Day);
  502.                 fprintf(OutFile,"  ----Orbit # %ld-----\n",OrbitNum);
  503.                 fprintf(OutFile," U.T.C.   Az   El  Doppler  Range");
  504.                 fprintf(OutFile,"  Height  Lat  Long  Phase(%3.0lf)\n",
  505.                                 MaxPhase);
  506.                 }
  507.             PrevOrbitNum = OrbitNum; PrevDay = Day;
  508.             PrintTime(OutFile,CurrentTime + HalfSecond);
  509.  
  510.             Doppler = -BeaconFreq*Velocity/C;
  511.  
  512.             fprintf(OutFile,"  %3.0lf  %3.0lf  %6.0lf  %6.0lf",
  513.             Azimuth*DegreesPerRadian,
  514.                 Elevation*DegreesPerRadian,Doppler,Range);
  515.  
  516.             fprintf(OutFile,"  %6.0lf  %3.0lf  %4.0lf",
  517.                 Radius - EarthRadius,SSPLat*DegreesPerRadian,
  518.                 SSPLong*DegreesPerRadian);
  519.  
  520.             Phase = (MeanAnomaly/PI2*MaxPhase + perigeePhase);
  521.             while (Phase < 0)
  522.                 Phase += MaxPhase;
  523.             while (Phase >= MaxPhase)
  524.                 Phase -= MaxPhase;
  525.  
  526.             fprintf(OutFile,"  %4d", Phase);
  527.             PrintMode(OutFile,Phase);
  528.             if (PrintApogee && (MeanAnomaly == PI))
  529.                 fprintf(OutFile,"    Apogee");
  530.             fprintf(OutFile,"\n");
  531.         PrevVisible = 1;
  532.             }
  533.      else
  534.         PrevVisible = 0;    
  535.         if (PrintApogee && (MeanAnomaly == PI))
  536.             DidApogee = 1;
  537.         PrevTime = CurrentTime;
  538.         CurrentTime = TmpTime;
  539.         }
  540.     fclose(OutFile);
  541.  
  542. }
  543.