home *** CD-ROM | disk | FTP | other *** search
/ PsL Monthly 1993 December / PSL Monthly Shareware CD-ROM (December 1993).iso / prgmming / dos / c / snippet.exe / DAYLEN.C < prev    next >
Encoding:
C/C++ Source or Header  |  1991-09-24  |  9.3 KB  |  241 lines
  1. /*
  2.  
  3. DAYLEN.C - computes the length of the day at any date and latitude
  4.  
  5. Paul Schlyter, 1989-08-16
  6.  
  7. (c) Paul Schlyter, 1989
  8.  
  9. This program may be used by anyone for any purpose, iff:
  10.    1. it is not being sold for profit
  11.    2. this notice is not removed
  12.  
  13. */
  14.  
  15.  
  16. #include <stdio.h>
  17. #include <math.h>
  18.  
  19.  
  20. /* A macro to compute the number of days elapsed since 2000 Jan 0.0 */
  21. /* (which is equal to 1999 Dec 31, 0h UT)                           */
  22.  
  23. #define days_since_2000_Jan_0(y,m,d) \
  24.     (367L*(y)-((7*((y)+(((m)+9)/12)))/4)+((275*(m))/9)+(d)-730530L)
  25.  
  26. /* Some conversion factors between radians and degrees */
  27.  
  28. #ifndef PI
  29.  #define PI        3.1415926535897932384
  30. #endif
  31.  
  32. #define RADEG     ( 180.0 / PI )
  33. #define DEGRAD    ( PI / 180.0 )
  34.  
  35. /* The trigonometric functions in degrees */
  36.  
  37. #define sind(x)  sin((x)*DEGRAD)
  38. #define cosd(x)  cos((x)*DEGRAD)
  39. #define tand(x)  tan((x)*DEGRAD)
  40.  
  41. #define atand(x)    (RADEG*atan(x))
  42. #define asind(x)    (RADEG*asin(x))
  43. #define acosd(x)    (RADEG*acos(x))
  44. #define atan2d(y,x) (RADEG*atan2(y,x))
  45.  
  46.  
  47. /* Following are some macros around the "workhorse" function __daylen__ */
  48. /* They mainly fill in the desired values for the reference altitude    */
  49. /* below the horizon, and also selects whether this altitude should     */
  50. /* refer to the Sun's center or its upper limb.                         */
  51.  
  52.  
  53. /* This macro computes the length of the day, from sunrise to sunset. */
  54. /* Sunrise/set is considered to occur when the Sun's upper limb is    */
  55. /* 35 arc minutes below the horizon (this accounts for the refraction */
  56. /* of the Earth's atmosphere).                                        */
  57. #define day_length(year,month,day,lon,lat)  \
  58.         __daylen__( year, month, day, lon, lat, -35.0/60.0, 1 )
  59.  
  60. /* This macro computes the length of the day, including civil twilight. */
  61. /* Civil twilight starts/ends when the Sun's center is 6 degrees below  */
  62. /* the horizon.                                                         */
  63. #define day_civil_twilight_length(year,month,day,lon,lat)  \
  64.         __daylen__( year, month, day, lon, lat, -6.0, 0 )
  65.  
  66. /* This macro computes the length of the day, incl. nautical twilight.  */
  67. /* Nautical twilight starts/ends when the Sun's center is 12 degrees    */
  68. /* below the horizon.                                                   */
  69. #define day_nautical_twilight_length(year,month,day,lon,lat)  \
  70.         __daylen__( year, month, day, lon, lat, -12.0, 0 )
  71.  
  72. /* This macro computes the length of the day, incl. astronomical twilight. */
  73. /* Astronomical twilight starts/ends when the Sun's center is 18 degrees   */
  74. /* below the horizon.                                                      */
  75. #define day_astronomical_twilight_length(year,month,day,lon,lat)  \
  76.         __daylen__( year, month, day, lon, lat, -18.0, 0 )
  77.  
  78.  
  79. /* Function prototypes */
  80.  
  81. double __daylen__( int year, int month, int day, double lon, double lat,
  82.                    double altit, int upper_limb );
  83.  
  84. void sunpos( double d, double *lon, double *r );
  85.  
  86. double revolution( double x );
  87.  
  88.  
  89.  
  90. /* A small test program */
  91.  
  92. void main(void)
  93. {
  94.       int year,month,day;
  95.       double lon, lat;
  96.       double daylen, civlen, nautlen, astrlen;
  97.  
  98.       printf( "Longitude (+ is east) and latitude (+ is north) : " );
  99.       scanf( "%lf %lf", &lon, &lat );
  100.  
  101.       for(;;)
  102.       {
  103.             printf( "Input date ( yyyy mm dd ) (ctrl-C exits): " );
  104.             scanf( "%d %d %d", &year, &month, &day );
  105.  
  106.             daylen  = day_length(year,month,day,lon,lat);
  107.             civlen  = day_civil_twilight_length(year,month,day,lon,lat);
  108.             nautlen = day_nautical_twilight_length(year,month,day,lon,lat);
  109.             astrlen = day_astronomical_twilight_length(year,month,day,
  110.                   lon,lat);
  111.  
  112.             printf( "Day length:                 %5.2f hours\n", daylen );
  113.             printf( "With civil twilight         %5.2f hours\n", civlen );
  114.             printf( "With nautical twilight      %5.2f hours\n", nautlen );
  115.             printf( "With astronomical twilight  %5.2f hours\n", astrlen );
  116.             printf( "Length of twilight: civil   %5.2f hours\n",
  117.                   (civlen-daylen)/2.0);
  118.             printf( "                  nautical  %5.2f hours\n",
  119.                   (nautlen-daylen)/2.0);
  120.             printf( "              astronomical  %5.2f hours\n",
  121.                   (astrlen-daylen)/2.0);
  122.       }
  123. }
  124.  
  125.  
  126. /* The "workhorse" function */
  127.  
  128.  
  129. double __daylen__( int year, int month, int day, double lon, double lat,
  130.                    double altit, int upper_limb )
  131. /**********************************************************************/
  132. /* Note: year,month,date = calendar date, 1801-2099 only.             */
  133. /*       Eastern longitude positive, Western longitude negative       */
  134. /*       Northern latitude positive, Southern latitude negative       */
  135. /*       The longitude value is not critical. Set it to the correct   */
  136. /*       longitude if you're picky, otherwise set to to, say, 0.0     */
  137. /*       The latitude however IS critical - be sure to get it correct */
  138. /*       altit = the latitude where the Sun should cross              */
  139. /*               Set to -35/60 degrees for rise/set, -6 degrees       */
  140. /*               for civil, -12 degrees for nautical and -18          */
  141. /*               degrees for astronomical twilight.                   */
  142. /*         upper_limb: non-zero -> upper limb, zero -> center         */
  143. /*               Set to non-zero (e.g. 1) when computing day length   */
  144. /*               and to zero when computing day+twilight length.      */
  145. /**********************************************************************/
  146. {
  147.       double  d,  /* Days since 2000 Jan 0.0 (negative before) */
  148.       obl_ecl,    /* Obliquity (inclination) of Earth's axis */
  149.       sr,         /* Solar distance, astronomical units */
  150.       slon,       /* True solar longitude */
  151.       sin_sdecl,  /* Sine of Sun's declination */
  152.       cos_sdecl,  /* Cosine of Sun's declination */
  153.       sradius,    /* Sun's apparent radius */
  154.       t;          /* Diurnal arc */
  155.  
  156.       /* Compute d of 12h local mean solar time */
  157.       d = days_since_2000_Jan_0(year,month,day) - lon/360.0;
  158.  
  159.       /* Compute obliquity of ecliptic (inclination of Earth's axis */
  160.       obl_ecl = 23.4393 - 3.563E-7 * d;
  161.  
  162.       /* Compute Sun's position */
  163.       sunpos( d, &slon, &sr );
  164.  
  165.       /* Compute sine and cosine of Sun's declination */
  166.       sin_sdecl = sind(obl_ecl) * sind(slon);
  167.       cos_sdecl = sqrt( 1.0 - sin_sdecl * sin_sdecl );
  168.  
  169.       /* Compute the Sun's apparent radius, degrees */
  170.       sradius = 0.2666 / sr;
  171.  
  172.       /* Do correction to upper limb, if necessary */
  173.       if ( upper_limb )
  174.             altit -= sradius;
  175.  
  176.       /* Compute the diurnal arc that the Sun traverses to reach */
  177.       /* the specified altitide altit: */
  178.       {
  179.             double cost;
  180.             cost = ( sind(altit) - sind(lat) * sin_sdecl ) /
  181.                   ( cosd(lat) * cos_sdecl );
  182.             if ( cost >= 1.0 )
  183.                   t = 0.0;                      /* Sun always below altit */
  184.             else if ( cost <= -1.0 )
  185.                   t = 24.0;                     /* Sun always above altit */
  186.             else  t = (2.0/15.0) * acosd(cost); /* The diurnal arc, hours */
  187.       }
  188.       return t;
  189. }
  190.  
  191.  
  192. /* This function computes the Sun's position at any instant */
  193.  
  194. void sunpos( double d, double *lon, double *r )
  195. /******************************************************/
  196. /* Computes the Sun's ecliptic longitude and distance */
  197. /* at an instant given in d, number of days since     */
  198. /* 2000 Jan 0.0.  The Sun's ecliptic latitude is not  */
  199. /* computed, since it's always very near 0.           */
  200. /******************************************************/
  201. {
  202.       double M,         /* Mean anomaly of the Sun */
  203.              w,         /* Mean longitude of perihelion */
  204.                         /* Note: Sun's mean longitude = M + w */
  205.              e,         /* Eccentricity of Earth's orbit */
  206.              E,         /* Eccentric anomaly */
  207.              x, y,      /* x, y coordinates in orbit */
  208.              v;         /* True anomaly */
  209.  
  210.       /* Compute mean elements */
  211.       M = revolution( 356.0470 + 0.9856002585 * d );
  212.       w = 282.9404 + 4.70935E-5 * d;
  213.       e = 0.016709 - 1.151E-9 * d;
  214.  
  215.       /* Compute true longitude and radius vector */
  216.       E = M + e * RADEG * sind(M) * ( 1.0 + e * cosd(M) );
  217.             x = cosd(E) - e;
  218.       y = sqrt( 1.0 - e*e ) * sind(E);
  219.       *r = sqrt( x*x + y*y );              /* Solar distance */
  220.       v = atan2d( y, x );                  /* True anomaly */
  221.       *lon = v + w;                        /* True solar longitude */
  222.       if ( *lon >= 360.0 )
  223.             *lon -= 360.0;                   /* Make it 0..360 degrees */
  224. }
  225.  
  226. /******************************************************************/
  227. /* This function reduces any angle to within the first revolution */
  228. /* by subtracting or adding even multiples of 360.0 until the     */
  229. /* result is >= 0.0 and < 360.0                                   */
  230. /******************************************************************/
  231.  
  232. #define INV360    ( 1.0 / 360.0 )
  233.  
  234. double revolution( double x )
  235. /*****************************************/
  236. /* Reduce angle to within 0..360 degrees */
  237. /*****************************************/
  238. {
  239.       return( x - 360.0 * floor( x * INV360 ) );
  240. }
  241.