/*----------------------------------------------------------------------------+ | Sunrise/Sunset functions for Heyu. | | | | The solar computation functions below were programmed by Charles W. | | Sullivan utilizing the techniques and astronomical constants published by | | Roger W. Sinnott in the August 1994 issue of Sky & Telescope Magazine, | | Page 84. | +----------------------------------------------------------------------------*/ #include #include #include #include #include #include "sun.h" #include "process.h" #define PI 3.14159265 #define D2R (PI / 180.) /* Ratio, Mean Solar Day / Mean Siderial Day */ #define MSSR 1.0027379 #define KS15 (15. * D2R * MSSR) /* Zenith angles for Sunrise/set and Twilights */ #define ANG_RISESET (90. + 50./60.) #define ANG_CIVILTWI (96.) #define ANG_NAUTITWI (102.) #define ANG_ASTROTWI (108.) static double zenangle[] = { ANG_RISESET, ANG_CIVILTWI, ANG_NAUTITWI, ANG_ASTROTWI }; /* Local functions */ static void sun_position ( double, double *, double *); static double local_sidereal_time (double, long int, double); /*---------------------------------------------------------+ | Calculate local times and Azimuths of Sunrise and | | Sunset on a specified date. | | Input: | | latitude and longitude in degrees are positive for | | North and East of Greenwich respectively. | | JulianDay is the Julian Day number at Greenwich Noon. | | timezone in seconds from Greenwich, positive for | | localities west of Greenwich. | | Return: | | sunrise and sunset times are in minutes after local | | midnight. | | azimuths at rise and set are in degrees (0 - 360). | | return code (defined in sun.h): | | NORMAL_SUN Sun rises and sets on this day. | | DOWN_ALL_DAY Sun below horizon all day. | | UP_ALL_DAY Sun above horizon all day. | | NO_SUNRISE Sun does not rise on this day. | | NO_SUNSET Sun does not set on this day. | | | +---------------------------------------------------------*/ int suntimes ( double latitude, double longitude, long int timezone, long int JulianDay, int sunmode, int *sunrise, int *sunset, double *azrise, double *azset) { double tdays, lst, zendist; double slat, clat,cozend ; double rasc[2], decl[2]; double hj ; int flag_rise, flag_set; int j ; int retcode = NORMAL_SUN; int *srssptr = NULL; double *azptr = NULL; double a, b, a0, a2, d0, d1, d2, dela, deld, p, el0, el2, h0, h1, h2, v0, v1, v2, d, e, t3, h7, n7, d7, az; *sunrise = *sunset = 0 ; /* Elapsed days from 1 Jan 2000 at 00:00 hours UTC0 */ tdays = (double)(JulianDay - 2451545L) - 0.5 ; /* Determine Local Siderial Time. */ lst = local_sidereal_time( tdays, timezone, longitude ); tdays += (double)timezone / (3600. * 24.); /* Get sun's position */ for ( j = 0; j < 2; j++ ) { sun_position ( tdays, &rasc[j], &decl[j] ); tdays += 1.0; } if ( rasc[1] < rasc[0] ) rasc[1] += 2.*PI ; #if 0 /* sunrise and sunset are defined when the sun is 50 minutes below the horizon (at sea level). */ zendist = D2R * (90. + 50./60.) ; #endif zendist = D2R * zenangle[sunmode]; slat = sin(D2R * latitude); clat = cos(D2R * latitude); cozend = cos(zendist); flag_rise = flag_set = 0; a0 = rasc[0]; d0 = decl[0]; v0 = 0.0; dela = rasc[1] - rasc[0]; deld = decl[1] - decl[0]; for ( j = 0; j < 24; j++ ) { hj = (double) j ; p = (1.0 + hj) / 24. ; a2 = rasc[0] + p * dela ; d2 = decl[0] + p * deld ; /* Test an hour for an event */ el0 = lst + hj * KS15 ; el2 = el0 + KS15 ; h0 = el0 - a0; h2 = el2 - a2 ; h1 = 0.5 * (h2 + h0); d1 = 0.5 * (d2 + d0) ; if ( j == 0 ) { v0 = slat * sin(d0) + clat * cos(d0) * cos(h0) - cozend ; } v2 = slat * sin(d2) + clat * cos(d2) * cos(h2) - cozend ; if ( v0 * v2 < 0. ) { v1 = slat * sin(d1) + clat * cos(d1) * cos(h1) - cozend ; a = 2. * v2 - 4. * v1 + 2. * v0 ; b = 4. * v1 - 3. * v0 - v2 ; d = b * b - 4. * a * v0 ; if ( d >= 0. ) { d = sqrt(d); if ( v0 < 0. && v2 > 0. ) { /* Event is Sunrise */ srssptr = sunrise ; azptr = azrise ; flag_rise = 1 ; } if ( v0 > 0. && v2 < 0. ) { /* Event is Sunset */ srssptr = sunset ; azptr = azset ; flag_set = 1; } e = (-b + d) / (2. * a); if ( e > 1. || e < 0. ) e = (-b -d) / (2. * a) ; t3 = hj + e ; *srssptr = (int)(60.* t3 + 0.5) ; /* Round off*/ h7 = h0 + e * (h2 -h0); n7 = - cos(d1) * sin(h7); d7 = clat * sin(d1) - slat * cos(d1) * cos(h7); az = fmod(((atan2(n7, d7) / D2R) + 360.), 360.) ; if ( azptr != NULL ) *azptr = az ; } } if ( flag_rise & flag_set ) break; a0 = a2; d0 = d2; v0 = v2; } if ( !(flag_rise | flag_set) ) { if ( v2 < 0. ) { /* Sun down all day */ retcode = DOWN_ALL_DAY ; } else if ( v2 >= 0. ) { /* Sun up all day */ retcode = UP_ALL_DAY ; } } else if ( !flag_rise ) { /* No sunrise this date */ retcode = NO_SUNRISE ; } else if ( !flag_set ) { /* No sunset this date */ retcode = NO_SUNSET ; } return retcode; } /*---------------------------------------------------------+ | Calculate local sidereal time. | | Input: | | tday - number of days from 2000 Jan 1 at | | 00:00 hours at Greenwich (UTC0). | +---------------------------------------------------------*/ double local_sidereal_time ( double tday, long int timezone, double longitude ) { double s ; s = 24110.5 + 8640184.812999999*tday/36525. ; s += 86636.6*(double)timezone/(3600.*24.) + 86400.*longitude/360. ; s = s/86400. ; s = s - floor(s) ; return (s * 360. * D2R) ; } /*---------------------------------------------------------+ | Calculate sun's Right Ascention and Declination. | +---------------------------------------------------------*/ void sun_position ( double tday, double *rasc, double *decl ) { double l, g, s, u, v, w ; double tcent ; /* Julian centuries from 1900.0 */ tcent = tday / 36525. + 1.0 ; /* Fundamental arguments (Van Flandern & Pulkinnen, 1979) */ l = .779072 + .00273790931 * tday ; g = .993126 + .0027377785 * tday ; l = 2. * PI * (l - floor(l)) ; g = 2. * PI * (g - floor(g)) ; v = .39785 * sin(l) - .01 * sin(l - g) + .00333 * sin(l + g) - .00021 * sin(l) * tcent ; u = 1.0 - .03349 * cos(g) - .00014 * cos(2. * l) + .00008 * cos (l) ; w = - .0001 - .04129 * sin(2. * l) + .03211 * sin(g) + .00104 * sin(2. * l - g) - .00035 * sin( 2. * l + g) - .00008 * sin(g) * tcent ; s = w / sqrt(u - v * v); *rasc = l + atan(s / sqrt(1.0 - s * s)); s = v / sqrt(u); *decl = atan(s / sqrt(1.0 - s * s)); return ; } /*---------------------------------------------------------+ | Gregorian calendar day to Julian Day number. | | Returns the Julian Day number at Greenwich Noon | | (UTC0 12:00) for the year, month, and day arguments. | | | | The Julian Day number is the count of whole days from | | Noon on 1 Jan 4713 B.C.E. in the Julian Proleptic | | Calendar. This calculation is historically valid from | | 15 Oct 1582 onward, however any year, month, and day | | greater than zero are acceptable as arguments and will | | yield the logically correct result. | +---------------------------------------------------------*/ long int greg2jd( int year, int month, int day ) { int count ; if ( month > 12 ) { year += ( month - 1 ) / 12 ; month = ( month - 1 ) % 12 + 1 ; } if ( month > 2 ) count = - ( 4 * month + 23 ) / 10 ; else { count = 365 ; year-- ; } count = count + year / 4 - 3 * ( year / 100 + 1 ) / 4 + 31 * ( month - 1 ) + day ; return ( (long)count + 365L*(long)year + 1721060L ) ; } /*---------------------------------------------------------------------+ | Adjust the times of Dawn and Dusk for abnormal sun conditions, as | | in polar regions, by creating and artificial Dawn and/or Dusk at | | 00:01 or 23:58 as appropriate. | +---------------------------------------------------------------------*/ int abnormal_sun_adjust ( int scode, int *dawnp, int *duskp ) { switch ( scode ) { case NORMAL_SUN : break; case UP_ALL_DAY : *dawnp = 1; *duskp = 23 * 60 + 58; break; case DOWN_ALL_DAY : *dawnp = 23 * 60 + 58; *duskp = 1; break; case NO_SUNRISE : *dawnp = 1; break; case NO_SUNSET : *duskp = 23 * 60 + 58; break; default : break; } return scode; } /*---------------------------------------------------------------------+ | Return the Julian Day corresponding to UTC0 Noon for the UTC0 time | | argument (in seconds from 1/1/1970 at 00:00:00 UTC0). | +---------------------------------------------------------------------*/ long int utc2jd ( long lutc0 ) { extern CONFIG *configp; return (lutc0 - ((lutc0 - configp->tzone) % 86400L)) / 86400L + 2440588L; } /*---------------------------------------------------------------------+ | Compute the UTC0 times of local Dawn and Dusk for the day which | | includes the argument UTC0 time. UTC0 times are all expressed as | | seconds elapsed from 1/1/1970 at 00:00:00 UTC0. | +---------------------------------------------------------------------*/ int local_dawndusk( time_t utc0, time_t *utc0_dawn, time_t *utc0_dusk ) { long int jd; int dawn, dusk, scode; long lutc0, midnight; extern CONFIG *configp; if ( configp->loc_flag != (LATITUDE | LONGITUDE) ) { *utc0_dawn = *utc0_dusk = 0; return -1; } lutc0 = (long)utc0; midnight = lutc0 - ((lutc0 - configp->tzone) % 86400L); jd = utc2jd(lutc0); scode = suntimes(configp->latitude, configp->longitude, configp->tzone, jd, configp->sunmode, &dawn, &dusk, NULL, NULL); /* Adjust for abnormal sun conditions */ abnormal_sun_adjust(scode, &dawn, &dusk); *utc0_dawn = (time_t)(midnight + 60L * (long)dawn); *utc0_dusk = (time_t)(midnight + 60L * (long)dusk); return scode; } /*-------------------------------------------------------------+ | Display a table of Dawn and Dusk for the year in the | | format used by the US Naval Observatory website at the time | | of this writing at: | | http://aa.usno.navy.mil/data/docs/RS_OneYear.html | | (Standard time only is displayed.) | | Argument sunmode indicates the definition of Dawn and Dusk: | | 0 -> Rise and Set | | 1 -> Civil Twilight | | 2 -> Nautical Twilight | | 3 -> Astronomical Twilight | | Printing on letter size or A4 paper requires 8 point fixed | | font and landscape mode. | +-------------------------------------------------------------*/ int display_sun_table_wide ( FILE *fd_sun, int year, long timezone, int sunmode, int timemode, int lat_d, int lat_m, int lon_d, int lon_m ) { static struct tzones { char *name; long seconds; } us_tzones[] = { {"Atlantic", 14400 }, {"Eastern", 18000 }, {"Central", 21600 }, {"Mountain", 25200 }, {"Pacific", 28800 }, {"Alaska", 32400 }, {"Hawaii-Aleutian", 36000 }, {"Samoa", 39600 }, {"Wake Island", -43200 }, {"Guam", -39600 }, }; static int n_tzones = ( sizeof(us_tzones) / sizeof(struct tzones) ); static int mdays[] = {0,31,28,31,30,31,30,31,31,30,31,30,31}; int time_adjust ( int, int, unsigned char ); long julianday, julianday0; int j, month, day, yday ; int rise, set, scode; int retcode; char tmark = ' '; char *timename; double latitude, longitude; get_dst_info(year); timename = ( timemode == TIMEMODE_CIVIL ) ? "Civil" : "Standard"; latitude = (lat_d < 0) ? (double)lat_d - (double)lat_m / 60. : (double)lat_d + (double)lat_m / 60. ; longitude = (lon_d < 0) ? (double)lon_d - (double)lon_m / 60. : (double)lon_d + (double)lon_m / 60. ; julianday0 = greg2jd(year, 1, 1); if ( greg2jd(year + 1, 1, 1) - julianday0 == 366 ) mdays[2] = 29; else mdays[2] = 28; (void) fprintf(fd_sun, "\n o , o , \n"); (void) fprintf(fd_sun, "Location: "); (void) fprintf(fd_sun, "%s%03d %02d,", longitude < 0 ? "W" : "E", abs(lon_d), abs(lon_m)); (void) fprintf(fd_sun, " %s%02d %02d", latitude < 0 ? "S" : "N", abs(lat_d), abs(lat_m)); switch (sunmode) { case RiseSet : (void) fprintf(fd_sun, "%*sSunrise and Sunset for %d", 29, " ", year); break; case CivilTwi : (void) fprintf(fd_sun, "%*sCivil Twilight for %d", 29, " ", year); break; case NautTwi : (void) fprintf(fd_sun, "%*sNautical Twilight for %d", 29, " ", year); break; case AstroTwi : (void) fprintf(fd_sun, "%*sAstronomical Twilight for %d", 26, " ", year); break; }; (void) fprintf(fd_sun, "%*sHEYU ver 2.0 \n\n", 39, " "); for ( j = 0; j < n_tzones; j++ ) { if ( us_tzones[j].seconds == timezone ) break; } if ( j < n_tzones ) (void) fprintf(fd_sun, "%*sUS/%s %s Time\n\n", 54, " ", us_tzones[j].name, timename); else (void) fprintf(fd_sun, "%*sTimezone: %.1fh %s of Greenwich - %s Time\n\n", 51, " ", (double)abs(timezone)/3600., (timezone < 0 ? "East" : "West"), timename ); (void) fprintf(fd_sun, " Jan Feb Mar Apr "); (void) fprintf(fd_sun, " May Jun Jul Aug "); (void) fprintf(fd_sun, " Sep Oct Nov Dec \n"); (void) fprintf(fd_sun, "Day Rise Set"); for ( month = 2; month <= 12; month++ ) (void) fprintf(fd_sun, " Rise Set"); (void) fprintf(fd_sun, "\n "); for ( month = 1; month <= 12; month++ ) (void) fprintf(fd_sun, " h m h m"); scode = NORMAL_SUN; for ( day = 1; day < 32; day++ ) { (void) fprintf(fd_sun, "\n%02d", day); for ( month = 1; month <= 12; month++ ) { if ( day > mdays[month] ) { (void) fprintf(fd_sun, " "); continue; } julianday = greg2jd( year, month, day ); retcode = suntimes(latitude, longitude, timezone, julianday, sunmode, &rise, &set, NULL, NULL ); yday = (int)(julianday - julianday0); if ( timemode == TIMEMODE_CIVIL ) { /* Adjust times for Daylight Time */ if ( retcode == NORMAL_SUN || retcode == NO_SUNSET ) { rise += time_adjust(yday, rise, LGL2STD); if ( rise < 0 || rise > 1439 ) retcode = NO_SUNRISE; } if ( retcode == NORMAL_SUN || retcode == NO_SUNRISE ) { set += time_adjust(yday, set, LGL2STD); if ( set < 0 || set > 1439 ) retcode = NO_SUNSET; } /* Mark day of time change */ tmark = (yday == 0 || time_adjust(yday, 720, LGL2STD) == time_adjust(yday - 1, 720, LGL2STD)) ? ' ' : '*'; } scode |= retcode; switch ( retcode ) { case DOWN_ALL_DAY : (void) fprintf(fd_sun, " %c---- ----", tmark); break ; case UP_ALL_DAY : (void) fprintf(fd_sun, " %c**** ****", tmark); break ; case NO_SUNRISE : (void) fprintf(fd_sun, " %c %02d%02d", tmark, set/60, set % 60); break ; case NO_SUNSET : (void) fprintf(fd_sun, " %c%02d%02d ", tmark, rise/60, rise%60); break ; default : (void) fprintf(fd_sun, " %c%02d%02d %02d%02d", tmark, rise/60, rise%60, set/60,set%60 ); } } } if ( timemode == TIMEMODE_CIVIL ) { fprintf(fd_sun, "\n\n(*) Denotes time change"); } else { fprintf(fd_sun, "\n\n%*sAdd offset for Daylight Time (usually +60 minutes), if and when in effect.", 34, " "); } if ( scode & (DOWN_ALL_DAY | UP_ALL_DAY) ) { fprintf(fd_sun, "\n\n(****) Sun continuously above horizon"); fprintf(fd_sun, "%*s(----) Sun continuously below horizon", 60, " "); } fprintf(fd_sun, "\n"); return 0; } /*-------------------------------------------------------------+ | Display a table of Dawn and Dusk for the year. | | Argument sunmode indicates the definition of Dawn and Dusk: | | 0 -> Rise and Set | | 1 -> Civil Twilight | | 2 -> Nautical Twilight | | 3 -> Astronomical Twilight | +-------------------------------------------------------------*/ int display_sun_table ( FILE *fd_sun, int year, long timezone, int sunmode, int timemode, int lat_d, int lat_m, int lon_d, int lon_m ) { static struct tzones { char *name; long seconds; } us_tzones[] = { {"Atlantic", 14400 }, {"Eastern", 18000 }, {"Central", 21600 }, {"Mountain", 25200 }, {"Pacific", 28800 }, {"Alaska", 32400 }, {"Hawaii-Aleutian", 36000 }, {"Samoa", 39600 }, {"Wake Island", -43200 }, {"Guam", -39600 }, }; static int n_tzones = ( sizeof(us_tzones) / sizeof(struct tzones) ); static int mdays[] = {0,31,28,31,30,31,30,31,31,30,31,30,31}; int time_adjust ( int, int, unsigned char ); long julianday, julianday0; int j, month, day, yday, half ; int rise, set, scode; int retcode; char tmark = ' '; char minibuf[64]; char *timename; double latitude, longitude; get_dst_info(year); timename = ( timemode == TIMEMODE_CIVIL ) ? "Civil" : "Standard"; latitude = (lat_d < 0) ? (double)lat_d - (double)lat_m / 60. : (double)lat_d + (double)lat_m / 60. ; longitude = (lon_d < 0) ? (double)lon_d - (double)lon_m / 60. : (double)lon_d + (double)lon_m / 60. ; julianday0 = greg2jd(year, 1, 1); if ( greg2jd(year + 1, 1, 1) - julianday0 == 366 ) mdays[2] = 29; else mdays[2] = 28; switch (sunmode) { case RiseSet : sprintf(minibuf, "Sunrise and Sunset for %d\n", year); break; case CivilTwi : sprintf(minibuf, "Civil Twilight for %d\n", year); break; case NautTwi : sprintf(minibuf, "Nautical Twilight for %d\n", year); break; case AstroTwi : sprintf(minibuf, "Astronomical Twilight for %d\n", year); break; }; fprintf(fd_sun, "%*s%s", (80 - (int)strlen(minibuf))/2, " ", minibuf); (void) fprintf(fd_sun, "\nLocation: "); (void) fprintf(fd_sun, "%s%03d:%02d,", longitude < 0 ? "W" : "E", abs(lon_d), abs(lon_m)); (void) fprintf(fd_sun, " %s%02d:%02d", latitude < 0 ? "S" : "N", abs(lat_d), abs(lat_m)); for ( j = 0; j < n_tzones; j++ ) { if ( us_tzones[j].seconds == timezone ) break; } if ( j < n_tzones ) (void) fprintf(fd_sun, " US/%s %s Time\n\n", us_tzones[j].name, timename); else (void) fprintf(fd_sun, " Timezone: %.1fh %s of Greenwich - %s Time\n\n", (double)abs(timezone)/3600., (timezone < 0 ? "East" : "West"), timename ); for ( half = 0; half < 2; half++ ) { if ( half == 0 ) (void) fprintf(fd_sun, " Jan Feb Mar Apr May Jun\n"); else (void) fprintf(fd_sun, "\n\n Jul Aug Sep Oct Nov Dec\n"); if ( sunmode == 0 ) { fprintf(fd_sun, "Day Rise Set"); for ( month = 2; month <= 6; month++ ) fprintf(fd_sun, " Rise Set"); fprintf(fd_sun, "\n "); } else { fprintf(fd_sun, "Day Morn Eve"); for ( month = 2; month <= 6; month++ ) fprintf(fd_sun, " Morn Eve"); fprintf(fd_sun, "\n "); } for ( month = 1; month <= 6; month++ ) (void) fprintf(fd_sun, " hh:mm hh:mm"); scode = NORMAL_SUN; for ( day = 1; day < 32; day++ ) { (void) fprintf(fd_sun, "\n%02d", day); for ( month = 6 * half + 1; month <= 6 * half + 6; month++ ) { if ( day > mdays[month] ) { (void) fprintf(fd_sun, " "); continue; } julianday = greg2jd( year, month, day ); yday = (int)(julianday - julianday0); retcode = suntimes(latitude, longitude, timezone, julianday, sunmode, &rise, &set, NULL, NULL ); if ( timemode == TIMEMODE_CIVIL ) { /* Adjust times for Daylight Time */ if ( retcode == NORMAL_SUN || retcode == NO_SUNSET ) { rise += time_adjust(yday, rise, LGL2STD); if ( rise < 0 || rise > 1439 ) retcode = NO_SUNRISE; } if ( retcode == NORMAL_SUN || retcode == NO_SUNRISE ) { set += time_adjust(yday, set, LGL2STD); if ( set < 0 || set > 1439 ) retcode = NO_SUNSET; } /* Mark day of time change */ tmark = (yday == 0 || time_adjust(yday, 720, LGL2STD) == time_adjust(yday - 1, 720, LGL2STD)) ? ' ' : '*'; } scode |= retcode; switch ( retcode ) { case DOWN_ALL_DAY : (void) fprintf(fd_sun, " %c----- -----", tmark); break ; case UP_ALL_DAY : (void) fprintf(fd_sun, " %c***** *****", tmark); break ; case NO_SUNRISE : (void) fprintf(fd_sun, " %c %02d:%02d", tmark, set/60, set % 60); break ; case NO_SUNSET : (void) fprintf(fd_sun, " %c%02d:%02d ", tmark, rise/60, rise%60); break ; default : (void) fprintf(fd_sun, " %c%02d:%02d %02d:%02d", tmark, rise/60, rise%60, set/60,set%60 ); } } } } if ( timemode == TIMEMODE_CIVIL ) { (void) fprintf(fd_sun, "\n\n(*) Denotes time change"); } else { (void) fprintf(fd_sun, "\n\nAdd offset for Daylight Time (usually +60 minutes) if and when in effect."); } if ( scode & (DOWN_ALL_DAY | UP_ALL_DAY) ) { (void) fprintf(fd_sun, "\n\n(*****) Sun continuously above horizon"); (void) fprintf(fd_sun, "%*s(-----) Sun continuously below horizon", 4, " "); } (void) fprintf(fd_sun, "\n"); return 0; }