home *** CD-ROM | disk | FTP | other *** search
/ ftp.parl.clemson.edu / 2015-02-07.ftp.parl.clemson.edu.tar / ftp.parl.clemson.edu / pub / portedOneB.tar / OneB / SolarTools.c < prev    next >
C/C++ Source or Header  |  1999-07-13  |  10KB  |  370 lines
  1. #include "RAConeb.h"
  2.  
  3. const LDOUBLE
  4. e_Re = 6378149.196,
  5. e_ecc_sq = 6.694378665e-3;
  6.  
  7. //------------------------------------------------------------------------------
  8. // 
  9. // Solar Zenith Angle Computation Program 
  10. //
  11. //      coded by C. Hoisington  SSAI
  12. //      19 Feb 98   version     1.0
  13. //
  14. //      modified by Scott Mitchell
  15. //      20 Feb 98   version     1.1
  16. //
  17. //    This contains the procedures and a test main program to compute
  18. //    Solar Zenith Angles
  19. //
  20. //------------------------------------------------------------------------------
  21.  
  22. /*
  23.  
  24.  CC -g SolarAngleProgram.c -I../mystuff2 ../mystuff2/mylib.a myorbit.a \
  25.    /usr/lib/libm.a -lm -o SolarAngleProgram
  26.  
  27. */
  28.  
  29. LDOUBLE PrincipleAngleDeg  (LDOUBLE angle);
  30. LDOUBLE JulianDate         (time_t unixTimeSec);
  31. LDOUBLE chGreenwichHourAngle (LDOUBLE t1970);
  32.  
  33. //------------------------------------------------------------------------------
  34. //
  35. // PrincipleAngleDeg
  36. //
  37. //    Convert angle to principle positive angle. (degrees)
  38. //
  39. //------------------------------------------------------------------------------
  40.  
  41. LDOUBLE PrincipleAngleDeg (LDOUBLE angle) 
  42. {
  43.     LDOUBLE base;
  44.  
  45.     base = fmod (angle, 360L);
  46.     if (base < 0.L) base += 360L;
  47.     return base;
  48. }
  49.  
  50.  
  51. //------------------------------------------------------------------------------
  52. //
  53. // AngleBetweenVectors
  54. //
  55. //    Compute angle between to vectors in rectangular coordinates. 
  56. //    (returns angle in radians)
  57. //
  58. //------------------------------------------------------------------------------
  59.  
  60.  
  61. LDOUBLE AngleBetweenVectors 
  62.     (
  63.     LDOUBLE x1, 
  64.     LDOUBLE y1, 
  65.     LDOUBLE z1,
  66.  
  67.     LDOUBLE x2, 
  68.     LDOUBLE y2, 
  69.     LDOUBLE z2
  70.     )
  71. {
  72.  
  73.     LDOUBLE dotProduct, absV1, absV2, aTemp;
  74.  
  75.     dotProduct = x1*x2 + y1*y2 + z1*z2;
  76.     absV1      = sqrt (x1*x1 + y1*y1 + z1*z1);
  77.     absV2      = sqrt (x2*x2 + y2*y2 + z2*z2);
  78.  
  79.     /* Perform an FINIT if defined in platform.h as needed */
  80.         FINIT
  81.  
  82.     aTemp = acos ( dotProduct / (absV1*absV2) );
  83.  
  84.     return aTemp;
  85. }
  86.      
  87.  
  88. //------------------------------------------------------------------------------
  89. //
  90. // chGreenwichHourAngle 
  91. //
  92. //    - compute Greenwich Hour Angle from time in seconds from 0 Jan 1970.
  93. //
  94. //
  95. //------------------------------------------------------------------------------
  96.  
  97. LDOUBLE chGreenwichHourAngle (LDOUBLE t1970) 
  98.  
  99. {
  100.     const LDOUBLE
  101.         DEL= +.22737064e0;           // correct GHA to 91/10/02 
  102.                                      //                01:13:18.030
  103.     const LDOUBLE 
  104.         G0 = 0.1746647027e+1,        // GHA on 01 Jan 1950  (Rad)
  105.         G1 = 0.1720279100e-1,        // GHA revolutions/Day (Rad/Day)
  106.         RE = 0.7292115850e-4,        // Earth rotation rate (Rad/Sec)
  107.         CA = 0.7671432057e-4,        // Nutation amplitude  (Rad)
  108.         A1 = 0.2114043493e+0,        // Nutation component 01 Jan 50
  109.         A2 = 0.9242186639e-3,        // Nutation rate       (Rad/Day)
  110.         CB = 0.5662286171e-5,        // Earth orbit rev comp.   (Rad)
  111.         B1 = 0.3493493163e+1,        // Earth rev comp. 01 Jan  (Rad)
  112.         B2 = 0.3440558258e-1;        // Earth rev rate      (Rad/day)
  113.     
  114.     LDOUBLE SEC,SD,RAD;
  115.     long DY;
  116.  
  117.     SEC = t1970 + 631152000;
  118.     DY  = SEC/86400e0;        // days since 1950
  119.     SD  = SEC-DY*86400 + DEL; // second of day
  120.  
  121.     RAD = G0 + G1 * DY + RE * SD         // Earth Rotation
  122.         - CA * sin (A1 - DY*A2)          // Nutation component
  123.         - CB * sin (B1 + DY*B2);         // revolution component
  124.  
  125.     RAD = fmod (RAD,m_twoPi);                   // convert to principle angle
  126.  
  127.     return RAD;
  128. }
  129.  
  130. //------------------------------------------------------------------------------
  131. //
  132. // RectInertialToFixed
  133. //
  134. //    Convert Geocentric Rectangular Inertial coordinates to Earth Fixed 
  135. //    Coordinates
  136. //
  137. //------------------------------------------------------------------------------
  138.  
  139. void RectInertialToFixed 
  140.     (
  141.     LDOUBLE greenwichHourAngle,
  142.     LDOUBLE inertxp,            // inertial coordinates input
  143.     LDOUBLE inertyp,
  144.     LDOUBLE inertzp,
  145.  
  146.     LDOUBLE *fixedxp,        // earth fixed coordinates output
  147.     LDOUBLE *fixedyp,
  148.     LDOUBLE *fixedzp
  149.     )
  150. {
  151.  
  152.     LDOUBLE sinGha,cosGha;
  153.  
  154.         sinGha = sin (greenwichHourAngle);
  155.         cosGha = cos (greenwichHourAngle);
  156.  
  157.         *fixedxp =  inertxp * cosGha + inertyp * sinGha;      // position
  158.         *fixedyp = -inertxp * sinGha + inertyp * cosGha;
  159.         *fixedzp =  inertzp;
  160. }
  161.  
  162.  
  163. //------------------------------------------------------------------------------
  164. //
  165. // GeodeticToFixedRect
  166. //
  167. //    Convert Geodetic (latitude and longitude) coordinates to earth fixed
  168. //    rectangular coordinates
  169. //
  170. //------------------------------------------------------------------------------
  171.  
  172.  
  173. void GeodeticToFixedRect 
  174.     (
  175.     LDOUBLE phi, 
  176.     LDOUBLE lambda,
  177.  
  178.     LDOUBLE *xp,
  179.     LDOUBLE *yp,
  180.     LDOUBLE *zp
  181.     ) 
  182. {
  183.  
  184.     LDOUBLE sinPhi,cosPhi,sinLambda,cosLambda,Rs;
  185.  
  186.     sinPhi = sin (phi);              // Sine phi
  187.     cosPhi = cos (phi);              // Cos  phi
  188.     sinLambda = sin (lambda);        // Sine   lambda
  189.     cosLambda = cos (lambda);        // Cosine lambda
  190.     Rs = e_Re / sqrt (1. - e_ecc_sq * sinPhi*sinPhi);
  191.  
  192.                                            // rectangular coordinates
  193.     *xp = (Rs) * cosPhi *  cosLambda;
  194.     *yp = (Rs) * cosPhi *  sinLambda;
  195.     *zp = (Rs * (1.-e_ecc_sq)) * sinPhi;
  196. }
  197.  
  198.  
  199. //------------------------------------------------------------------------------
  200. //
  201. // GeoInertialSunLocation
  202. //
  203. //    Compute sun location in geocentric inertial rectangular coordinates
  204. //
  205. //    Note: See "The Astronomical Almanac for the Year 1996", U.S. Government
  206. //          Printing Office, page C24
  207. //
  208. //------------------------------------------------------------------------------
  209.  
  210. void GeoInertialSunLocation 
  211.     (
  212.       LDOUBLE JulianDate,
  213.       LDOUBLE *xp,
  214.       LDOUBLE *yp,
  215.       LDOUBLE *zp
  216.     )
  217. {
  218.  
  219.     LDOUBLE au = 149598845.0L;         // astronomical unit in kilometers
  220.     LDOUBLE n, L, g, gr, eta, etar, lambda, lambdar, R, x, y, z;
  221.  
  222.     n = JulianDate - 2451545.0L;         // days from J2000.0
  223.     L =  280.461L + 0.9856474L * n;      // mean longitude of sun
  224.     L =  PrincipleAngleDeg(L);
  225.     g =  357.528L + 0.9856003L * n;      // mean anomaly (deg)
  226.     g =  PrincipleAngleDeg(g);
  227.     gr=  g*m_dToR;
  228.     eta = 23.439L - 0.0000004L * n;      // Obliquity of ecliptic
  229.     etar= eta * m_dToR;
  230.     lambda = L + 1.915L * sin(gr) + 0.020 * sin (2*gr); // Ecliptic longitude
  231.     lambdar= lambda *m_dToR;
  232.     R = 1.00014L - 0.01671L * cos(gr) - 0.00014L * cos(2*gr); // distance of sun from earth in a.u
  233.     x = R * cos(lambdar);
  234.     y = R * cos(etar   ) * sin(lambdar);
  235.     z = R * sin(etar   ) * sin(lambdar);
  236.  
  237.     *xp = x * au;
  238.     *yp = y * au;
  239.     *zp = z * au;
  240. }
  241.  
  242.  
  243. //------------------------------------------------------------------------------
  244. //
  245. // GeoFixedSunLocation
  246. //
  247. //    Compute sun location in geocentric fixed rectangular coordinates
  248. //
  249. //------------------------------------------------------------------------------
  250.  
  251. void GeoFixedSunLocation 
  252.     (
  253.       LDOUBLE unixSec,
  254.       LDOUBLE *sunXfixed,
  255.       LDOUBLE *sunYfixed,
  256.       LDOUBLE *sunZfixed
  257.     )
  258. {
  259.     LDOUBLE sunXinertial, sunYinertial, sunZinertial;
  260.     LDOUBLE julianDays, greenwichHourAngle;
  261.  
  262.     julianDays =  JulianDate (unixSec);
  263.  
  264.     GeoInertialSunLocation 
  265.     (
  266.         julianDays, 
  267.        &sunXinertial, &sunYinertial, &sunZinertial
  268.     );
  269.  
  270.     greenwichHourAngle = chGreenwichHourAngle ( (LDOUBLE)unixSec );
  271.  
  272.     RectInertialToFixed 
  273.     (
  274.         greenwichHourAngle, 
  275.         sunXinertial, sunYinertial, sunZinertial,
  276.         sunXfixed,    sunYfixed,    sunZfixed
  277.     );
  278. }
  279.  
  280. //------------------------------------------------------------------------------
  281. //
  282. // JulianDate
  283. //
  284. //------------------------------------------------------------------------------
  285.  
  286. LDOUBLE JulianDate (time_t unixTimeSec)
  287. {
  288.     return 2440587.5L + (LDOUBLE)unixTimeSec / 86400.L;
  289. }
  290.  
  291. //------------------------------------------------------------------------------
  292. //
  293. // SolarZenithAngle
  294. //
  295. //    Compute the solar zenith angle at a given time for the specified 
  296. //    latitude and longitude. 
  297. //
  298. //    (angle is returned in degrees.)
  299. //
  300. //------------------------------------------------------------------------------
  301.  
  302.  
  303. LDOUBLE SolarZenithAngle
  304.     (
  305.     char *buffer, int fourDigitYear,  int dayOfYear, long milliSecondOfDay, time_t unixTime,
  306.  
  307.     LDOUBLE latitudeDeg, LDOUBLE longitudeDeg
  308.     )
  309.  
  310.  
  311. //------------------------------------------------------------------------
  312. //  procedure body
  313. //
  314. //------------------------------------------------------------------------
  315. {
  316.     LDOUBLE sunXfixed, sunYfixed, sunZfixed;
  317.     LDOUBLE pixelXfixed, pixelYfixed, pixelZfixed;    
  318.     time_t unixSec; 
  319.     LDOUBLE aTemp;
  320.  
  321.     unixSec = unixTime;
  322.  
  323.         //--------------------------------------------------------------------
  324.     //  compute Sun location
  325.     //--------------------------------------------------------------------
  326.     
  327.     GeoFixedSunLocation 
  328.     (
  329.       (LDOUBLE)unixSec, 
  330.       &sunXfixed, &sunYfixed, &sunZfixed
  331.     );
  332.  
  333.     //--------------------------------------------------------------------
  334.     //  compute pixel location
  335.     //--------------------------------------------------------------------
  336.    
  337.     GeodeticToFixedRect 
  338.     (
  339.         latitudeDeg*m_dToR, longitudeDeg*m_dToR,  
  340.  
  341.         &pixelXfixed, &pixelYfixed, &pixelZfixed 
  342.     );
  343.  
  344.     //--------------------------------------------------------------------
  345.     //  compute angle between pixel and sun
  346.     //--------------------------------------------------------------------
  347.  
  348. // printf("sunXfixed = %f, sunYfixed = %f, sunZfixed = %f\n", sunXfixed, sunYfixed, sunZfixed);
  349. // printf("pixelXfixed = %f, pixelYfixed = %f, pixelZfixed = %f\n", pixelXfixed, pixelYfixed, pixelZfixed);
  350.  
  351.     aTemp = (AngleBetweenVectors
  352.                 (
  353.                         sunXfixed,  sunYfixed,   sunZfixed,
  354.                         pixelXfixed, pixelYfixed, pixelZfixed
  355.  
  356.                 ) * m_rToD);
  357.  
  358.  
  359.     aTemp *= 2.0;
  360. //        printf("(uchar)aTemp = %c\n", (uchar)aTemp);
  361.     buffer[0] = (uchar)aTemp;
  362.         // printf("buffer[0] = %c\n", buffer[0]);
  363.  
  364.         // printf("aTemp = %Lf\n", aTemp);
  365.  
  366.     return aTemp;
  367. }
  368.  
  369.  
  370.