home *** CD-ROM | disk | FTP | other *** search
/ rtsi.com / 2014.01.www.rtsi.com.tar / www.rtsi.com / OS9 / OSK / TELECOM / OSKBox.lzh / MAILBOX / SAT / sun.c < prev   
C/C++ Source or Header  |  1990-10-28  |  8KB  |  281 lines
  1.  
  2. #define DRAG 1
  3.  
  4. #include <stdio.h>
  5. #include <math.h>
  6. #include <ctype.h>
  7.  
  8. extern double Kepler();
  9. extern long GetDayNum();
  10.  
  11. #define LC(c) (isupper(c) ? tolower(c) : (c))
  12.  
  13. #ifndef PI
  14. #define PI 3.14159265
  15. #endif
  16.  
  17. #ifdef PI2
  18. #undef PI2
  19. #endif
  20.  
  21. #define PI2 (PI*2)
  22.  
  23. #define MinutesPerDay (24*60.0)
  24. #define SecondsPerDay (60*MinutesPerDay)
  25. #define HalfSecond (0.5/SecondsPerDay)
  26. #define EarthRadius 6378.16             /* Kilometers           */
  27. #define C 2.997925e5                    /* Kilometers/Second    */
  28. #define TropicalYear 365.24199        /* Mean solar days    */
  29. #define EarthEccentricity 0.016713
  30. #define DegreesPerRadian (180/PI)
  31. #define RadiansPerDegree (PI/180)
  32. #define ABS(x) ((x) < 0 ? (-(x)) : (x))
  33. #define SQR(x) ((x)*(x))
  34.  
  35. #define SunRadius 695000        
  36. #define SunSemiMajorAxis  149598845.0          /* Kilometers            */
  37.  
  38. #define MaxModes 10
  39. typedef struct {
  40.                 int MinPhase,MaxPhase;
  41.                 char ModeStr[20];
  42.                }  ModeRec;
  43.  
  44. char VersionStr[] = "Sun ephemeris 1.0";
  45.  
  46.     /*  Keplerian Elements and misc. data for the satellite              */
  47.     double  EpochDay;                   /* time of epoch                 */
  48.     double EpochMeanAnomaly;            /* Mean Anomaly at epoch         */
  49.     long EpochOrbitNum;                 /* Integer orbit # of epoch      */
  50.     double EpochRAAN;                   /* RAAN at epoch                 */
  51.     double epochMeanMotion;             /* Revolutions/day               */
  52.     double OrbitalDecay;                /* Revolutions/day^2             */
  53.     double EpochArgPerigee;             /* argument of perigee at epoch  */
  54.     double Eccentricity;
  55.     double Inclination;
  56.     char SatName[100];
  57.     int ElementSet;
  58.     double BeaconFreq;                  /* Mhz, used for doppler calc    */
  59.     double MaxPhase;                    /* Phase units in 1 orbit        */
  60.     double perigeePhase;
  61.     int NumModes;
  62.     ModeRec Modes[MaxModes];
  63.     int PrintApogee;
  64.     int PrintEclipses;
  65.     int Flip;
  66.  
  67.     /* Simulation Parameters */
  68.  
  69.     double StartTime,EndTime, StepTime; /* In Days, 1 = New Year        */
  70.                                         /*      of reference year       */
  71.  
  72.     /* Site Parameters */
  73.     char SiteName[100];
  74.     double SiteLat,SiteLong,SiteAltitude,SiteMinElev;
  75.  
  76.  
  77.  
  78. #define SITELAT 37.943
  79. #define SITELONG 122.318
  80. #define SITEALT 73
  81.  
  82. GetSiteParams()
  83. {
  84.  
  85.     SiteLat = SITELAT*RadiansPerDegree;
  86.     SiteLong = SITELONG*RadiansPerDegree;
  87.     SiteAltitude = SITEALT/1000.0;   /* convert to km */
  88.     SiteMinElev = 0;
  89.     Flip = PrintEclipses = 0;
  90.     }
  91.  
  92. #define ONEDAY 86400    /* in seconds */
  93. #define FONEDAY 86400.0
  94. #define UTC 8
  95.  
  96. double LocalTime (Time)
  97. double Time;
  98. {
  99.     return Time-UTC/24.0;
  100.     }
  101.     
  102. SPrintTime2(Str,Time)
  103. char *Str;
  104. double Time;
  105. {
  106.     int day,hours,minutes,seconds;
  107.     char ampm;
  108.  
  109.     day = Time;
  110.     Time -= day;
  111.     if (Time < 0)
  112.         Time += 1.0;   /* Correct for truncation problems with negatives */
  113.  
  114.     hours = Time*24;
  115.     Time -=  hours/24.0;
  116.  
  117.     minutes = Time*MinutesPerDay;
  118.     Time -= minutes/MinutesPerDay;
  119.  
  120.     seconds = Time*SecondsPerDay;
  121.     seconds -= seconds/SecondsPerDay;
  122.  
  123.     ampm = 'A';
  124.     if (hours > 12) {
  125.         hours -= 12;
  126.         ampm = 'P';
  127.         }
  128.  
  129.     sprintf(Str,"%d:%02d%cM",hours,minutes,ampm);
  130. }
  131.  
  132. PrintTime2(OutFile,Time)
  133. FILE *OutFile;
  134. double Time;
  135. {
  136.     char str[100];
  137.  
  138.     SPrintTime2(str,Time);
  139.     fprintf(OutFile,"%s",str);
  140. }
  141.  
  142. GetSimulationParams()
  143. {
  144.     int time, date, day, tick;
  145.  
  146.     _sysdate (1, &time, &date, &day, &tick);
  147.     if ((time += UTC*60*60) > ONEDAY) {time -= ONEDAY; date++;}
  148.     StartTime = date - 2415019.0 + time / FONEDAY;
  149. }
  150.  
  151. SunPos (CurrentTime, Azimuth, Elevation, SSPLat, SSPLong, Height)
  152. double CurrentTime;
  153. double *Azimuth, *Elevation;
  154. double *SSPLat, *SSPLong, *Height;
  155. {
  156.     extern double SunEpochTime,SunInclination,SunRAAN,SunEccentricity,
  157.         SunArgPerigee,SunMeanAnomaly,SunMeanMotion;
  158.     double MeanAnomaly, TrueAnomaly;
  159.     double SatX,SatY,SatZ;    /* In Right Ascension based system */
  160.     double SatVX,SatVY,SatVZ;   /* Kilometers/second           */
  161.     double SiteX,SiteY,SiteZ;
  162.     double SiteVX,SiteVY;
  163.     double SiteMatrix[3][3];
  164.     double Radius;              /* From geocenter                  */
  165.  
  166.     MeanAnomaly = SunMeanAnomaly+ (CurrentTime-SunEpochTime)*SunMeanMotion*PI2;
  167.     TrueAnomaly = Kepler(MeanAnomaly,SunEccentricity);
  168.     GetSatPosition(SunEpochTime,SunRAAN,SunArgPerigee,SunSemiMajorAxis,
  169.         SunInclination,SunEccentricity,0.0,0.0,CurrentTime,
  170.         TrueAnomaly,&SatX,&SatY,&SatZ,&Radius,&SatVX,&SatVY,&SatVZ);
  171.     GetSitPosition(SiteLat,SiteLong,SiteAltitude,CurrentTime,
  172.         &SiteX,&SiteY,&SiteZ,&SiteVX,&SiteVY,SiteMatrix);
  173.     GetBearings(SatX,SatY,SatZ,SiteX,SiteY,SiteZ,SiteMatrix,
  174.         Azimuth,Elevation);
  175.     if (SSPLat && SSPLong && Height)
  176.         GetSubSatPoint(SatX,SatY,SatZ,CurrentTime,SSPLat,SSPLong,Height);
  177.     } 
  178.  
  179. double SunSet (CurrentTime)
  180. double CurrentTime;
  181. {
  182.     double Azimuth, Elevation;
  183.     double StartTime, EndTime;
  184.     
  185.     StartTime = CurrentTime;
  186.     EndTime = StartTime + .125;
  187.     SunPos (EndTime, &Azimuth, &Elevation, 0, 0, 0);
  188.     while (Elevation > 0.0) {
  189.         EndTime += .125;
  190.         SunPos (EndTime, &Azimuth, &Elevation, 0, 0, 0);
  191.         }
  192.     do {
  193.         CurrentTime = StartTime + ((EndTime - StartTime) / 2);
  194.         SunPos (CurrentTime, &Azimuth, &Elevation, 0, 0, 0);
  195.         if (Elevation > 0.0)
  196.             StartTime = CurrentTime;
  197.         else
  198.             EndTime = CurrentTime;
  199.         }
  200.     while (EndTime - StartTime > 1.0 / ONEDAY);
  201.     return CurrentTime;
  202.     }
  203.  
  204. double SunRise (CurrentTime)
  205. double CurrentTime;
  206. {
  207.     double Azimuth, Elevation;
  208.     double StartTime, EndTime;
  209.     
  210.     StartTime = CurrentTime;
  211.     EndTime = StartTime + .125;
  212.     SunPos (EndTime, &Azimuth, &Elevation, 0, 0, 0);
  213.     while (Elevation < 0.0) {
  214.         EndTime += .125;
  215.         SunPos (EndTime, &Azimuth, &Elevation, 0, 0, 0);
  216.         }
  217.     do {
  218.         CurrentTime = StartTime + ((EndTime - StartTime) / 2);
  219.         SunPos (CurrentTime, &Azimuth, &Elevation, 0, 0, 0);
  220.         if (Elevation < 0.0)
  221.             StartTime = CurrentTime;
  222.         else
  223.             EndTime = CurrentTime;
  224.         }
  225.     while (EndTime - StartTime > 1.0 / ONEDAY);
  226.     return CurrentTime;
  227.     }
  228.  
  229. main()
  230. {
  231.     double ReferenceOrbit;      /* Floating point orbit # at epoch */
  232.     double CurrentTime;         /* In Days                         */
  233.     double CurrentOrbit;
  234.     double AverageMotion,       /* Corrected for drag              */
  235.     CurrentMotion;
  236.     double MeanAnomaly,TrueAnomaly;
  237.     double SemiMajorAxis;
  238.     double Height;
  239.     double RAANPrecession,PerigeePrecession;
  240.     double SSPLat,SSPLong;
  241.     long OrbitNum,PrevOrbitNum;
  242.     long Day,PrevDay;
  243.     double Azimuth,Elevation,Range;
  244.     double RangeRate,Doppler;
  245.     int Phase;
  246.     char FileName[100];
  247.     FILE *OutFile;
  248.     int DidApogee;
  249.     double TmpTime,PrevTime;
  250.     int PrevVisible;
  251.  
  252.  
  253.     GetSiteParams();
  254.     GetSimulationParams();
  255.     InitOrbitRoutines((StartTime+EndTime)/2);
  256.     CurrentTime = StartTime;
  257.     SunPos (CurrentTime, &Azimuth, &Elevation, &SSPLat, &SSPLong, &Height);
  258.     printf("The sun is currently above %0.2lf%c degrees lattitude, %0.0lf%c degrees longitude.\n",
  259.         fabs (SSPLat*DegreesPerRadian),
  260.         ((SSPLat*DegreesPerRadian > 0) ? 'N' : 'S'),
  261.         ((SSPLong > PI) ? PI2-SSPLong : SSPLong)*DegreesPerRadian,
  262.         (SSPLong > PI) ? 'E' : 'W');
  263.     if (Elevation > 0.0) {
  264.         printf("It's azimuth is %0.0lf degrees and its elevation is %0.0lf degrees.\n",
  265.             Azimuth*DegreesPerRadian,Elevation*DegreesPerRadian);
  266.         printf ("The sun will set tonight at ");
  267.         CurrentTime = SunSet (CurrentTime);
  268.         PrintTime2 (stdout, LocalTime (CurrentTime));
  269.         printf (" and will rise tomorrow morning at ");
  270.         PrintTime2 (stdout, LocalTime (SunRise (CurrentTime)));
  271.         }
  272.     else {
  273.         printf ("The sun will rise at ");
  274.         CurrentTime = SunRise (CurrentTime);
  275.         PrintTime2 (stdout, LocalTime (CurrentTime));
  276.         printf (" and will set tomorrow evening at ");
  277.         PrintTime2 (stdout, LocalTime (SunSet (CurrentTime)));
  278.         }
  279.     printf (".\n");
  280.     }
  281.