home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / snip9707.zip / SUNRISET.C < prev    next >
C/C++ Source or Header  |  1997-07-05  |  22KB  |  513 lines
  1. /* +++Date last modified: 05-Jul-1997 */
  2.  
  3. /*
  4.  
  5. SUNRISET.C - computes Sun rise/set times, start/end of twilight, and
  6.              the length of the day at any date and latitude
  7.  
  8. Written as DAYLEN.C, 1989-08-16
  9.  
  10. Modified to SUNRISET.C, 1992-12-01
  11.  
  12. (c) Paul Schlyter, 1989, 1992
  13.  
  14. Released to the public domain by Paul Schlyter, December 1992
  15.  
  16. */
  17.  
  18.  
  19. #include <stdio.h>
  20. #include <math.h>
  21.  
  22.  
  23. /* A macro to compute the number of days elapsed since 2000 Jan 0.0 */
  24. /* (which is equal to 1999 Dec 31, 0h UT)                           */
  25.  
  26. #define days_since_2000_Jan_0(y,m,d) \
  27.     (367L*(y)-((7*((y)+(((m)+9)/12)))/4)+((275*(m))/9)+(d)-730530L)
  28.  
  29. /* Some conversion factors between radians and degrees */
  30.  
  31. #ifndef PI
  32.  #define PI        3.1415926535897932384
  33. #endif
  34.  
  35. #define RADEG     ( 180.0 / PI )
  36. #define DEGRAD    ( PI / 180.0 )
  37.  
  38. /* The trigonometric functions in degrees */
  39.  
  40. #define sind(x)  sin((x)*DEGRAD)
  41. #define cosd(x)  cos((x)*DEGRAD)
  42. #define tand(x)  tan((x)*DEGRAD)
  43.  
  44. #define atand(x)    (RADEG*atan(x))
  45. #define asind(x)    (RADEG*asin(x))
  46. #define acosd(x)    (RADEG*acos(x))
  47. #define atan2d(y,x) (RADEG*atan2(y,x))
  48.  
  49.  
  50. /* Following are some macros around the "workhorse" function __daylen__ */
  51. /* They mainly fill in the desired values for the reference altitude    */
  52. /* below the horizon, and also selects whether this altitude should     */
  53. /* refer to the Sun's center or its upper limb.                         */
  54.  
  55.  
  56. /* This macro computes the length of the day, from sunrise to sunset. */
  57. /* Sunrise/set is considered to occur when the Sun's upper limb is    */
  58. /* 35 arc minutes below the horizon (this accounts for the refraction */
  59. /* of the Earth's atmosphere).                                        */
  60. #define day_length(year,month,day,lon,lat)  \
  61.         __daylen__( year, month, day, lon, lat, -35.0/60.0, 1 )
  62.  
  63. /* This macro computes the length of the day, including civil twilight. */
  64. /* Civil twilight starts/ends when the Sun's center is 6 degrees below  */
  65. /* the horizon.                                                         */
  66. #define day_civil_twilight_length(year,month,day,lon,lat)  \
  67.         __daylen__( year, month, day, lon, lat, -6.0, 0 )
  68.  
  69. /* This macro computes the length of the day, incl. nautical twilight.  */
  70. /* Nautical twilight starts/ends when the Sun's center is 12 degrees    */
  71. /* below the horizon.                                                   */
  72. #define day_nautical_twilight_length(year,month,day,lon,lat)  \
  73.         __daylen__( year, month, day, lon, lat, -12.0, 0 )
  74.  
  75. /* This macro computes the length of the day, incl. astronomical twilight. */
  76. /* Astronomical twilight starts/ends when the Sun's center is 18 degrees   */
  77. /* below the horizon.                                                      */
  78. #define day_astronomical_twilight_length(year,month,day,lon,lat)  \
  79.         __daylen__( year, month, day, lon, lat, -18.0, 0 )
  80.  
  81.  
  82. /* This macro computes times for sunrise/sunset.                      */
  83. /* Sunrise/set is considered to occur when the Sun's upper limb is    */
  84. /* 35 arc minutes below the horizon (this accounts for the refraction */
  85. /* of the Earth's atmosphere).                                        */
  86. #define sun_rise_set(year,month,day,lon,lat,rise,set)  \
  87.         __sunriset__( year, month, day, lon, lat, -35.0/60.0, 1, rise, set )
  88.  
  89. /* This macro computes the start and end times of civil twilight.       */
  90. /* Civil twilight starts/ends when the Sun's center is 6 degrees below  */
  91. /* the horizon.                                                         */
  92. #define civil_twilight(year,month,day,lon,lat,start,end)  \
  93.         __sunriset__( year, month, day, lon, lat, -6.0, 0, start, end )
  94.  
  95. /* This macro computes the start and end times of nautical twilight.    */
  96. /* Nautical twilight starts/ends when the Sun's center is 12 degrees    */
  97. /* below the horizon.                                                   */
  98. #define nautical_twilight(year,month,day,lon,lat,start,end)  \
  99.         __sunriset__( year, month, day, lon, lat, -12.0, 0, start, end )
  100.  
  101. /* This macro computes the start and end times of astronomical twilight.   */
  102. /* Astronomical twilight starts/ends when the Sun's center is 18 degrees   */
  103. /* below the horizon.                                                      */
  104. #define astronomical_twilight(year,month,day,lon,lat,start,end)  \
  105.         __sunriset__( year, month, day, lon, lat, -18.0, 0, start, end )
  106.  
  107.  
  108. /* Function prototypes */
  109.  
  110. double __daylen__( int year, int month, int day, double lon, double lat,
  111.                    double altit, int upper_limb );
  112.  
  113. int __sunriset__( int year, int month, int day, double lon, double lat,
  114.                   double altit, int upper_limb, double *rise, double *set );
  115.  
  116. void sunpos( double d, double *lon, double *r );
  117.  
  118. void sun_RA_dec( double d, double *RA, double *dec, double *r );
  119.  
  120. double revolution( double x );
  121.  
  122. double rev180( double x );
  123.  
  124. double GMST0( double d );
  125.  
  126.  
  127.  
  128. /* A small test program */
  129.  
  130. main()
  131. {
  132.       int year,month,day;
  133.       double lon, lat;
  134.       double daylen, civlen, nautlen, astrlen;
  135.       double rise, set, civ_start, civ_end, naut_start, naut_end,
  136.              astr_start, astr_end;
  137.       int    rs, civ, naut, astr;
  138.       char buf[80];
  139.  
  140.       printf( "Longitude (+ is east) and latitude (+ is north) : " );
  141.       fgets(buf, 80, stdin);
  142.       sscanf(buf, "%lf %lf", &lon, &lat );
  143.  
  144.       for(;;)
  145.       {
  146.             printf( "Input date ( yyyy mm dd ) (ctrl-C exits): " );
  147.             fgets(buf, 80, stdin);
  148.             sscanf(buf, "%d %d %d", &year, &month, &day );
  149.  
  150.             daylen  = day_length(year,month,day,lon,lat);
  151.             civlen  = day_civil_twilight_length(year,month,day,lon,lat);
  152.             nautlen = day_nautical_twilight_length(year,month,day,lon,lat);
  153.             astrlen = day_astronomical_twilight_length(year,month,day,
  154.                   lon,lat);
  155.  
  156.             printf( "Day length:                 %5.2f hours\n", daylen );
  157.             printf( "With civil twilight         %5.2f hours\n", civlen );
  158.             printf( "With nautical twilight      %5.2f hours\n", nautlen );
  159.             printf( "With astronomical twilight  %5.2f hours\n", astrlen );
  160.             printf( "Length of twilight: civil   %5.2f hours\n",
  161.                   (civlen-daylen)/2.0);
  162.             printf( "                  nautical  %5.2f hours\n",
  163.                   (nautlen-daylen)/2.0);
  164.             printf( "              astronomical  %5.2f hours\n",
  165.                   (astrlen-daylen)/2.0);
  166.  
  167.             rs   = sun_rise_set         ( year, month, day, lon, lat,
  168.                                           &rise, &set );
  169.             civ  = civil_twilight       ( year, month, day, lon, lat,
  170.                                           &civ_start, &civ_end );
  171.             naut = nautical_twilight    ( year, month, day, lon, lat,
  172.                                           &naut_start, &naut_end );
  173.             astr = astronomical_twilight( year, month, day, lon, lat,
  174.                                           &astr_start, &astr_end );
  175.  
  176.             printf( "Sun at south %5.2fh UT\n", (rise+set)/2.0 );
  177.  
  178.             switch( rs )
  179.             {
  180.                 case 0:
  181.                     printf( "Sun rises %5.2fh UT, sets %5.2fh UT\n",
  182.                              rise, set );
  183.                     break;
  184.                 case +1:
  185.                     printf( "Sun above horizon\n" );
  186.                     break;
  187.                 case -1:
  188.                     printf( "Sun below horizon\n" );
  189.                     break;
  190.             }
  191.  
  192.             switch( civ )
  193.             {
  194.                 case 0:
  195.                     printf( "Civil twilight starts %5.2fh, "
  196.                             "ends %5.2fh UT\n", civ_start, civ_end );
  197.                     break;
  198.                 case +1:
  199.                     printf( "Never darker than civil twilight\n" );
  200.                     break;
  201.                 case -1:
  202.                     printf( "Never as bright as civil twilight\n" );
  203.                     break;
  204.             }
  205.  
  206.             switch( naut )
  207.             {
  208.                 case 0:
  209.                     printf( "Nautical twilight starts %5.2fh, "
  210.                             "ends %5.2fh UT\n", naut_start, naut_end );
  211.                     break;
  212.                 case +1:
  213.                     printf( "Never darker than nautical twilight\n" );
  214.                     break;
  215.                 case -1:
  216.                     printf( "Never as bright as nautical twilight\n" );
  217.                     break;
  218.             }
  219.  
  220.             switch( astr )
  221.             {
  222.                 case 0:
  223.                     printf( "Astronomical twilight starts %5.2fh, "
  224.                             "ends %5.2fh UT\n", astr_start, astr_end );
  225.                     break;
  226.                 case +1:
  227.                     printf( "Never darker than astronomical twilight\n" );
  228.                     break;
  229.                 case -1:
  230.                     printf( "Never as bright as astronomical twilight\n" );
  231.                     break;
  232.             }
  233.       return 0;
  234.       }
  235. }
  236.  
  237.  
  238. /* The "workhorse" function for sun rise/set times */
  239.  
  240. int __sunriset__( int year, int month, int day, double lon, double lat,
  241.                   double altit, int upper_limb, double *trise, double *tset )
  242. /***************************************************************************/
  243. /* Note: year,month,date = calendar date, 1801-2099 only.             */
  244. /*       Eastern longitude positive, Western longitude negative       */
  245. /*       Northern latitude positive, Southern latitude negative       */
  246. /*       The longitude value IS critical in this function!            */
  247. /*       altit = the altitude which the Sun should cross              */
  248. /*               Set to -35/60 degrees for rise/set, -6 degrees       */
  249. /*               for civil, -12 degrees for nautical and -18          */
  250. /*               degrees for astronomical twilight.                   */
  251. /*         upper_limb: non-zero -> upper limb, zero -> center         */
  252. /*               Set to non-zero (e.g. 1) when computing rise/set     */
  253. /*               times, and to zero when computing start/end of       */
  254. /*               twilight.                                            */
  255. /*        *rise = where to store the rise time                        */
  256. /*        *set  = where to store the set  time                        */
  257. /*                Both times are relative to the specified altitude,  */
  258. /*                and thus this function can be used to compute       */
  259. /*                various twilight times, as well as rise/set times   */
  260. /* Return value:  0 = sun rises/sets this day, times stored at        */
  261. /*                    *trise and *tset.                               */
  262. /*               +1 = sun above the specified "horizon" 24 hours.     */
  263. /*                    *trise set to time when the sun is at south,    */
  264. /*                    minus 12 hours while *tset is set to the south  */
  265. /*                    time plus 12 hours. "Day" length = 24 hours     */
  266. /*               -1 = sun is below the specified "horizon" 24 hours   */
  267. /*                    "Day" length = 0 hours, *trise and *tset are    */
  268. /*                    both set to the time when the sun is at south.  */
  269. /*                                                                    */
  270. /**********************************************************************/
  271. {
  272.       double  d,  /* Days since 2000 Jan 0.0 (negative before) */
  273.       sr,         /* Solar distance, astronomical units */
  274.       sRA,        /* Sun's Right Ascension */
  275.       sdec,       /* Sun's declination */
  276.       sradius,    /* Sun's apparent radius */
  277.       t,          /* Diurnal arc */
  278.       tsouth,     /* Time when Sun is at south */
  279.       sidtime;    /* Local sidereal time */
  280.  
  281.       int rc = 0; /* Return cde from function - usually 0 */
  282.  
  283.       /* Compute d of 12h local mean solar time */
  284.       d = days_since_2000_Jan_0(year,month,day) + 0.5 - lon/360.0;
  285.  
  286.       /* Compute local sidereal time of this moment */
  287.       sidtime = revolution( GMST0(d) + 180.0 + lon );
  288.  
  289.       /* Compute Sun's RA + Decl at this moment */
  290.       sun_RA_dec( d, &sRA, &sdec, &sr );
  291.  
  292.       /* Compute time when Sun is at south - in hours UT */
  293.       tsouth = 12.0 - rev180(sidtime - sRA)/15.0;
  294.  
  295.       /* Compute the Sun's apparent radius, degrees */
  296.       sradius = 0.2666 / sr;
  297.  
  298.       /* Do correction to upper limb, if necessary */
  299.       if ( upper_limb )
  300.             altit -= sradius;
  301.  
  302.       /* Compute the diurnal arc that the Sun traverses to reach */
  303.       /* the specified altitude altit: */
  304.       {
  305.             double cost;
  306.             cost = ( sind(altit) - sind(lat) * sind(sdec) ) /
  307.                   ( cosd(lat) * cosd(sdec) );
  308.             if ( cost >= 1.0 )
  309.                   rc = -1, t = 0.0;       /* Sun always below altit */
  310.             else if ( cost <= -1.0 )
  311.                   rc = +1, t = 12.0;      /* Sun always above altit */
  312.             else
  313.                   t = acosd(cost)/15.0;   /* The diurnal arc, hours */
  314.       }
  315.  
  316.       /* Store rise and set times - in hours UT */
  317.       *trise = tsouth - t;
  318.       *tset  = tsouth + t;
  319.  
  320.       return rc;
  321. }  /* __sunriset__ */
  322.  
  323.  
  324.  
  325. /* The "workhorse" function */
  326.  
  327.  
  328. double __daylen__( int year, int month, int day, double lon, double lat,
  329.                    double altit, int upper_limb )
  330. /**********************************************************************/
  331. /* Note: year,month,date = calendar date, 1801-2099 only.             */
  332. /*       Eastern longitude positive, Western longitude negative       */
  333. /*       Northern latitude positive, Southern latitude negative       */
  334. /*       The longitude value is not critical. Set it to the correct   */
  335. /*       longitude if you're picky, otherwise set to to, say, 0.0     */
  336. /*       The latitude however IS critical - be sure to get it correct */
  337. /*       altit = the altitude which the Sun should cross              */
  338. /*               Set to -35/60 degrees for rise/set, -6 degrees       */
  339. /*               for civil, -12 degrees for nautical and -18          */
  340. /*               degrees for astronomical twilight.                   */
  341. /*         upper_limb: non-zero -> upper limb, zero -> center         */
  342. /*               Set to non-zero (e.g. 1) when computing day length   */
  343. /*               and to zero when computing day+twilight length.      */
  344. /**********************************************************************/
  345. {
  346.       double  d,  /* Days since 2000 Jan 0.0 (negative before) */
  347.       obl_ecl,    /* Obliquity (inclination) of Earth's axis */
  348.       sr,         /* Solar distance, astronomical units */
  349.       slon,       /* True solar longitude */
  350.       sin_sdecl,  /* Sine of Sun's declination */
  351.       cos_sdecl,  /* Cosine of Sun's declination */
  352.       sradius,    /* Sun's apparent radius */
  353.       t;          /* Diurnal arc */
  354.  
  355.       /* Compute d of 12h local mean solar time */
  356.       d = days_since_2000_Jan_0(year,month,day) + 0.5 - lon/360.0;
  357.  
  358.       /* Compute obliquity of ecliptic (inclination of Earth's axis) */
  359.       obl_ecl = 23.4393 - 3.563E-7 * d;
  360.  
  361.       /* Compute Sun's position */
  362.       sunpos( d, &slon, &sr );
  363.  
  364.       /* Compute sine and cosine of Sun's declination */
  365.       sin_sdecl = sind(obl_ecl) * sind(slon);
  366.       cos_sdecl = sqrt( 1.0 - sin_sdecl * sin_sdecl );
  367.  
  368.       /* Compute the Sun's apparent radius, degrees */
  369.       sradius = 0.2666 / sr;
  370.  
  371.       /* Do correction to upper limb, if necessary */
  372.       if ( upper_limb )
  373.             altit -= sradius;
  374.  
  375.       /* Compute the diurnal arc that the Sun traverses to reach */
  376.       /* the specified altitude altit: */
  377.       {
  378.             double cost;
  379.             cost = ( sind(altit) - sind(lat) * sin_sdecl ) /
  380.                   ( cosd(lat) * cos_sdecl );
  381.             if ( cost >= 1.0 )
  382.                   t = 0.0;                      /* Sun always below altit */
  383.             else if ( cost <= -1.0 )
  384.                   t = 24.0;                     /* Sun always above altit */
  385.             else  t = (2.0/15.0) * acosd(cost); /* The diurnal arc, hours */
  386.       }
  387.       return t;
  388. }  /* __daylen__ */
  389.  
  390.  
  391. /* This function computes the Sun's position at any instant */
  392.  
  393. void sunpos( double d, double *lon, double *r )
  394. /******************************************************/
  395. /* Computes the Sun's ecliptic longitude and distance */
  396. /* at an instant given in d, number of days since     */
  397. /* 2000 Jan 0.0.  The Sun's ecliptic latitude is not  */
  398. /* computed, since it's always very near 0.           */
  399. /******************************************************/
  400. {
  401.       double M,         /* Mean anomaly of the Sun */
  402.              w,         /* Mean longitude of perihelion */
  403.                         /* Note: Sun's mean longitude = M + w */
  404.              e,         /* Eccentricity of Earth's orbit */
  405.              E,         /* Eccentric anomaly */
  406.              x, y,      /* x, y coordinates in orbit */
  407.              v;         /* True anomaly */
  408.  
  409.       /* Compute mean elements */
  410.       M = revolution( 356.0470 + 0.9856002585 * d );
  411.       w = 282.9404 + 4.70935E-5 * d;
  412.       e = 0.016709 - 1.151E-9 * d;
  413.  
  414.       /* Compute true longitude and radius vector */
  415.       E = M + e * RADEG * sind(M) * ( 1.0 + e * cosd(M) );
  416.             x = cosd(E) - e;
  417.       y = sqrt( 1.0 - e*e ) * sind(E);
  418.       *r = sqrt( x*x + y*y );              /* Solar distance */
  419.       v = atan2d( y, x );                  /* True anomaly */
  420.       *lon = v + w;                        /* True solar longitude */
  421.       if ( *lon >= 360.0 )
  422.             *lon -= 360.0;                   /* Make it 0..360 degrees */
  423. }
  424.  
  425. void sun_RA_dec( double d, double *RA, double *dec, double *r )
  426. {
  427.       double lon, obl_ecl, x, y, z;
  428.  
  429.       /* Compute Sun's ecliptical coordinates */
  430.       sunpos( d, &lon, r );
  431.  
  432.       /* Compute ecliptic rectangular coordinates (z=0) */
  433.       x = *r * cosd(lon);
  434.       y = *r * sind(lon);
  435.  
  436.       /* Compute obliquity of ecliptic (inclination of Earth's axis) */
  437.       obl_ecl = 23.4393 - 3.563E-7 * d;
  438.  
  439.       /* Convert to equatorial rectangular coordinates - x is unchanged */
  440.       z = y * sind(obl_ecl);
  441.       y = y * cosd(obl_ecl);
  442.  
  443.       /* Convert to spherical coordinates */
  444.       *RA = atan2d( y, x );
  445.       *dec = atan2d( z, sqrt(x*x + y*y) );
  446.  
  447. }  /* sun_RA_dec */
  448.  
  449.  
  450. /******************************************************************/
  451. /* This function reduces any angle to within the first revolution */
  452. /* by subtracting or adding even multiples of 360.0 until the     */
  453. /* result is >= 0.0 and < 360.0                                   */
  454. /******************************************************************/
  455.  
  456. #define INV360    ( 1.0 / 360.0 )
  457.  
  458. double revolution( double x )
  459. /*****************************************/
  460. /* Reduce angle to within 0..360 degrees */
  461. /*****************************************/
  462. {
  463.       return( x - 360.0 * floor( x * INV360 ) );
  464. }  /* revolution */
  465.  
  466. double rev180( double x )
  467. /*********************************************/
  468. /* Reduce angle to within +180..+180 degrees */
  469. /*********************************************/
  470. {
  471.       return( x - 360.0 * floor( x * INV360 + 0.5 ) );
  472. }  /* revolution */
  473.  
  474.  
  475. /*******************************************************************/
  476. /* This function computes GMST0, the Greenwich Mean Sidereal Time  */
  477. /* at 0h UT (i.e. the sidereal time at the Greenwhich meridian at  */
  478. /* 0h UT).  GMST is then the sidereal time at Greenwich at any     */
  479. /* time of the day.  I've generalized GMST0 as well, and define it */
  480. /* as:  GMST0 = GMST - UT  --  this allows GMST0 to be computed at */
  481. /* other times than 0h UT as well.  While this sounds somewhat     */
  482. /* contradictory, it is very practical:  instead of computing      */
  483. /* GMST like:                                                      */
  484. /*                                                                 */
  485. /*  GMST = (GMST0) + UT * (366.2422/365.2422)                      */
  486. /*                                                                 */
  487. /* where (GMST0) is the GMST last time UT was 0 hours, one simply  */
  488. /* computes:                                                       */
  489. /*                                                                 */
  490. /*  GMST = GMST0 + UT                                              */
  491. /*                                                                 */
  492. /* where GMST0 is the GMST "at 0h UT" but at the current moment!   */
  493. /* Defined in this way, GMST0 will increase with about 4 min a     */
  494. /* day.  It also happens that GMST0 (in degrees, 1 hr = 15 degr)   */
  495. /* is equal to the Sun's mean longitude plus/minus 180 degrees!    */
  496. /* (if we neglect aberration, which amounts to 20 seconds of arc   */
  497. /* or 1.33 seconds of time)                                        */
  498. /*                                                                 */
  499. /*******************************************************************/
  500.  
  501. double GMST0( double d )
  502. {
  503.       double sidtim0;
  504.       /* Sidtime at 0h UT = L (Sun's mean longitude) + 180.0 degr  */
  505.       /* L = M + w, as defined in sunpos().  Since I'm too lazy to */
  506.       /* add these numbers, I'll let the C compiler do it for me.  */
  507.       /* Any decent C compiler will add the constants at compile   */
  508.       /* time, imposing no runtime or code overhead.               */
  509.       sidtim0 = revolution( ( 180.0 + 356.0470 + 282.9404 ) +
  510.                           ( 0.9856002585 + 4.70935E-5 ) * d );
  511.       return sidtim0;
  512. }  /* GMST0 */
  513.