/* * sscalc.c - a sunrise/sunset calculator, adapted from the JavaScript * program on NOAA's site at address: * http://www.srrb.noaa.gov/highlights/sunrise/gen.html * * The original javascript code was written by: * Aaron Horiuchi, Chris Lehman and Chris Cornwall. * This port released with permission of Chris Cornwall. * * This code by: * Keith Walker * 28 November 2000 * $Id: sscalc.c,v 1.3 2000/11/29 18:36:42 kew Exp kew $ */ #include #include #include #include #include static double radToDeg(double angleRad) { return (180.0 * angleRad / M_PI); } static double degToRad(double angleDeg) { return (M_PI * angleDeg / 180.0); } /* Returns the gamma value that is used in the calc for the equation ** of time and the solar declination */ static double calcGamma(double julDate) { return (2.0 * M_PI / 365.0) * (julDate - 1); } /* ** Return the gamma value used to calc equation of time and solar ** declination */ static double calcGamma2(double julDate, double hour) { return (2 * M_PI / 365) * (julDate - 1 + (hour/24)); } /* ** Return the equation of time value for the given date. */ static double calcEqofTime(double gamma) { return (229.18 * (0.000075 + 0.001868 * cos(gamma) - 0.032077 * sin(gamma) - 0.014615 * cos(2 * gamma) - 0.040849 * sin(2 * gamma))); } /* ** Return the solar declination angle (in radians) for the date */ static double calcSolarDec(double gamma) { return (0.006918 - 0.399912 * cos(gamma) + 0.070257 * sin(gamma) - 0.006758 * cos(2 * gamma) + 0.000907 * sin(2 * gamma)); } /* ** The hour angle returned below is only for sunrise/sunset, i.e. when ** the solar zenith angle is 90.8 degrees. The reason why it's not 90 ** degrees is because we need to account for atmospheric refraction. */ static double calcHourAngle(double lat, double solarDec, short calcSunrise) { double latRad = degToRad(lat); if (calcSunrise) return (acos(cos(degToRad(90.833)) / (cos(latRad) * cos(solarDec)) - tan(latRad) * tan(solarDec))); else return -(acos(cos(degToRad(90.833)) / (cos(latRad) * cos(solarDec)) - tan(latRad) * tan(solarDec))); } /* ** Calculate the sunrise at GMT for the lat/long */ double calcSunriseGMT(double julDay, double lat, double lon) { double gamma = calcGamma(julDay); double eqTime = calcEqofTime(gamma); double solarDec = calcSolarDec(gamma); double hourAngle = calcHourAngle(lat, solarDec, 1); double delta = lon - radToDeg(hourAngle); double timeDiff = 4 * delta; double timeGMT = 720 + timeDiff - eqTime; double gammaSunrise = calcGamma2(julDay, timeGMT/60); eqTime = calcEqofTime(gammaSunrise); solarDec = calcSolarDec(gammaSunrise); hourAngle = calcHourAngle(lat, solarDec, 1); delta = lon - radToDeg(hourAngle); timeDiff = 4 * delta; timeGMT = 720 + timeDiff - eqTime; return timeGMT; } double calcSunsetGMT(double julDay, double lat, double lon) { double gamma = calcGamma(julDay + 1); double eqTime = calcEqofTime(gamma); double solarDec = calcSolarDec(gamma); double hourAngle = calcHourAngle(lat, solarDec, 0); double delta = lon - radToDeg(hourAngle); double timeDiff = 4 * delta; double setTimeGMT = 720 + timeDiff - eqTime; double gamma_sunset = calcGamma2(julDay, setTimeGMT/60); eqTime = calcEqofTime(gamma_sunset); solarDec = calcSolarDec(gamma_sunset); hourAngle = calcHourAngle(lat, solarDec, 0); delta = lon - radToDeg(hourAngle); timeDiff = 4 * delta; setTimeGMT = 720 + timeDiff - eqTime; return setTimeGMT; } double findRecentSunrise(double julDay, double lat, double lon) { double jDay = julDay; double t = calcSunriseGMT(jDay, lat, lon); while (isnan(t)) { jDay--; if (jDay < 1) jDay = 365; t = calcSunriseGMT(jDay, lat, lon); } return jDay; } double findRecentSunset(double julDay, double lat, double lon) { double jDay = julDay; double t = calcSunsetGMT(jDay, lat, lon); while (isnan(t)) { jDay--; if (jDay < 1) jDay = 365; t = calcSunsetGMT(jDay, lat, lon); } return jDay; } double findNextSunrise(double julDay, double lat, double lon) { double jDay = julDay; double t = calcSunriseGMT(jDay, lat, lon); while (isnan(t)) { jDay++; if (jDay > 366) jDay = 1; t = calcSunriseGMT(jDay, lat, lon); } return jDay; } double findNextSunset(double julDay, double lat, double lon) { double jDay = julDay; double t = calcSunsetGMT(jDay, lat, lon); while (isnan(t)) { jDay++; if (jDay > 366) jDay = 1; t = calcSunsetGMT(jDay, lat, lon); } return jDay; } /* ** Return a time_t value of the sunrise expressed as GMT. ** It is assumed that you'll use the normal ctime(3) type functions ** to extract the values. */ time_t getSunrise(double lat, double lon, double jDay) { double floatHour; double hour; double floatMinute; double minute; double floatSec; double second; double ssGMT; time_t t = time(0); struct tm tm; struct tm *now = localtime(&t); ssGMT = calcSunriseGMT(jDay, lat, lon); if (isnan(ssGMT)) { /* ** Above the Arctic Circle, during the daylight season, find ** the last date that the sun rose (it's been up ever since ** then). */ if (lat > 66.4 && (jDay > 79 && jDay < 267)) jDay = findRecentSunrise(jDay, lat, lon); /* ** Above the Arctic Circle, during the nighttime season, find ** the date that the sun will come up again and report that. */ else if (lat > 66.4 && (jDay < 83 || jDay > 263)) jDay = findNextSunrise(jDay, lat, lon); /* ** Do the same things as above, only for the southern ** hemisphere */ else if (lat < -66.4 && (jDay < 83 || jDay > 263)) jDay = findRecentSunrise(jDay, lat, lon); else if (lat < -66.4 && (jDay > 79 && jDay < 267)) jDay = findNextSunrise(jDay, lat, lon); ssGMT = calcSunriseGMT(jDay, lat, lon); } floatHour = ssGMT / 60; hour = floor(floatHour); floatMinute = 60 * (floatHour - floor(floatHour)); minute = floor(floatMinute); floatSec = 60 * (floatMinute - floor(floatMinute)); second = floor(floatSec); tm.tm_year = now->tm_year; tm.tm_mon = 0; tm.tm_mday = jDay; tm.tm_hour = (int)hour; tm.tm_min = (int)minute; tm.tm_sec = (int)second; return timegm(&tm); } /* ** Return a time_t value of the sunset expressed as GMT. It is ** assumed that you'll use the normal ctime(3) type functions to ** extract the values. */ time_t getSunset(double lat, double lon, double jDay) { double floatHour; double hour; double floatMinute; double minute; double floatSec; double second; double ssGMT; time_t t = time(0); struct tm tm; struct tm *now = localtime(&t); ssGMT = calcSunsetGMT(jDay, lat, lon); if (isnan(ssGMT)) { if (lat > 66.4 && jDay > 79 && jDay < 267) jDay = findNextSunset(jDay, lat, lon); else if (lat > 66.4 && (jDay < 83 || jDay > 263)) jDay = findRecentSunset(jDay, lat, lon); else if (lat < -66.4 && (jDay < 83 || jDay > 263)) jDay = findNextSunset(jDay, lat, lon); else if (lat < -66.4 && (jDay > 79 && jDay < 267)) jDay = findRecentSunset(jDay, lat, lon); ssGMT = calcSunsetGMT(jDay, lat, lon); } floatHour = ssGMT / 60; hour = floor(floatHour); floatMinute = 60 * (floatHour - floor(floatHour)); minute = floor(floatMinute); floatSec = 60 * (floatMinute - floor(floatMinute)); second = floor(floatSec); tm.tm_year = now->tm_year; tm.tm_mon = 0; tm.tm_mday = jDay; tm.tm_hour = (int)hour; tm.tm_min = (int)minute; tm.tm_sec = (int)second; return timegm(&tm); } int main(int argc, char **argv) { char *p; char *fmt = "%X"; char buf[128]; int c; int i; int printsr = 1; int printss = 1; int reps = 1; time_t t = time(0); struct tm *tm = localtime(&t); double lat = 47.66; /* if nothing else given, use the authors */ double lon = 117.43; /* location of Spokane, Washington, USA */ double jDay = tm->tm_yday + 1; #if defined(LATITUDE) lat = LATITUDE; #endif #if defined(LONGITUDE) lon = LONGITUDE; #endif if ((p = getenv("LATITUDE"))) lat = atof(p); if ((p = getenv("LONGITUDE"))) lon = atof(p); while ((c = getopt(argc, argv, "rsf:m:d:a:o:n:")) != -1) { switch (c) { case 'r': printss = 0; break; case 's': printsr = 0; break; case 'f': fmt = optarg; break; case 'm': tm->tm_mon = atoi(optarg) - 1; break; case 'd': tm->tm_mday = atoi(optarg); break; case 'a': lat = atof(optarg); break; case 'o': lon = atof(optarg); break; case 'n': reps = atoi(optarg); break; default: fprintf(stderr, "usage: sscalc [-n reps][-m mon][-d day]" "[-o lon][-a lat][-rs]\n"); exit(1); break; } } t = mktime(tm); if (!printss && !printsr) printss = printsr = 1; if (lat >= -90 && lat < -89.8) lat = -89.8; if (lat <= 90 && lat > 89.8) lat = 89.8; for (i = 0; i < reps; i++, t += 86400) { tm = localtime(&t); jDay = tm->tm_yday + 1; if (printsr) { t = getSunrise(lat, lon, jDay); strftime(buf, sizeof(buf), fmt, localtime(&t)); printf("%s", buf); } if (printss) { if (printsr) printf(" "); t = getSunset(lat, lon, jDay); strftime(buf, sizeof(buf), fmt, localtime(&t)); printf("%s", buf); } printf("\n"); } return 0; }