From 39902fe0d6545eb16f3bad0d5e008a71e943578a Mon Sep 17 00:00:00 2001 From: curt Date: Fri, 1 Aug 1997 15:27:56 +0000 Subject: [PATCH] Initial revision. --- Time/sunpos.c | 260 ++++++++++++++++++++++++++++++++++++++++++++++++++ Time/sunpos.h | 53 ++++++++++ 2 files changed, 313 insertions(+) create mode 100644 Time/sunpos.c create mode 100644 Time/sunpos.h diff --git a/Time/sunpos.c b/Time/sunpos.c new file mode 100644 index 000000000..c022b4b73 --- /dev/null +++ b/Time/sunpos.c @@ -0,0 +1,260 @@ +/* + * sunpos.c + * kirk johnson + * july 1993 + * + * code for calculating the position on the earth's surface for which + * the sun is directly overhead (adapted from _practical astronomy + * with your calculator, third edition_, peter duffett-smith, + * cambridge university press, 1988.) + * + * RCS $Id$ + * + * Copyright (C) 1989, 1990, 1993, 1994, 1995 Kirk Lauritz Johnson + * + * Parts of the source code (as marked) are: + * Copyright (C) 1989, 1990, 1991 by Jim Frost + * Copyright (C) 1992 by Jamie Zawinski + * + * Permission to use, copy, modify and freely distribute xearth for + * non-commercial and not-for-profit purposes is hereby granted + * without fee, provided that both the above copyright notice and this + * permission notice appear in all copies and in supporting + * documentation. + * + * The author makes no representations about the suitability of this + * software for any purpose. It is provided "as is" without express or + * implied warranty. + * + * THE AUTHOR DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, + * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, + * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, INDIRECT + * OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM + * LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, + * NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN + * CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + * + * $Id$ + * (Log is kept at end of this file) + */ + + +#include +#include +#include + +#include "sunpos.h" +#include "../constants.h" +#undef E + + +/* + * the epoch upon which these astronomical calculations are based is + * 1990 january 0.0, 631065600 seconds since the beginning of the + * "unix epoch" (00:00:00 GMT, Jan. 1, 1970) + * + * given a number of seconds since the start of the unix epoch, + * DaysSinceEpoch() computes the number of days since the start of the + * astronomical epoch (1990 january 0.0) + */ + +#define EpochStart (631065600) +#define DaysSinceEpoch(secs) (((secs)-EpochStart)*(1.0/(24*3600))) + +/* + * assuming the apparent orbit of the sun about the earth is circular, + * the rate at which the orbit progresses is given by RadsPerDay -- + * FG_2PI radians per orbit divided by 365.242191 days per year: + */ + +#define RadsPerDay (FG_2PI/365.242191) + +/* + * details of sun's apparent orbit at epoch 1990.0 (after + * duffett-smith, table 6, section 46) + * + * Epsilon_g (ecliptic longitude at epoch 1990.0) 279.403303 degrees + * OmegaBar_g (ecliptic longitude of perigee) 282.768422 degrees + * Eccentricity (eccentricity of orbit) 0.016713 + */ + +#define Epsilon_g (279.403303*(FG_2PI/360)) +#define OmegaBar_g (282.768422*(FG_2PI/360)) +#define Eccentricity (0.016713) + +/* + * MeanObliquity gives the mean obliquity of the earth's axis at epoch + * 1990.0 (computed as 23.440592 degrees according to the method given + * in duffett-smith, section 27) + */ +#define MeanObliquity (23.440592*(FG_2PI/360)) + +static double solve_keplers_equation(double); +static double sun_ecliptic_longitude(time_t); +static void ecliptic_to_equatorial(double, double, double *, double *); +static double julian_date(int, int, int); +static double GST(time_t); + +/* + * solve Kepler's equation via Newton's method + * (after duffett-smith, section 47) + */ +static double solve_keplers_equation(double M) { + double E; + double delta; + + E = M; + while (1) { + delta = E - Eccentricity*sin(E) - M; + if (fabs(delta) <= 1e-10) break; + E -= delta / (1 - Eccentricity*cos(E)); + } + + return E; +} + + +/* compute ecliptic longitude of sun (in radians) (after + * duffett-smith, section 47) */ + +static double sun_ecliptic_longitude(time_t ssue) { + /* time_t ssue; seconds since unix epoch */ + double D, N; + double M_sun, E; + double v; + + D = DaysSinceEpoch(ssue); + + N = RadsPerDay * D; + N = fmod(N, FG_2PI); + if (N < 0) N += FG_2PI; + + M_sun = N + Epsilon_g - OmegaBar_g; + if (M_sun < 0) M_sun += FG_2PI; + + E = solve_keplers_equation(M_sun); + v = 2 * atan(sqrt((1+Eccentricity)/(1-Eccentricity)) * tan(E/2)); + + return (v + OmegaBar_g); +} + + +/* convert from ecliptic to equatorial coordinates (after + * duffett-smith, section 27) */ + +static void ecliptic_to_equatorial(double lambda, double beta, + double *alpha, double *delta) { + /* double lambda; ecliptic longitude */ + /* double beta; ecliptic latitude */ + /* double *alpha; (return) right ascension */ + /* double *delta; (return) declination */ + + double sin_e, cos_e; + + sin_e = sin(MeanObliquity); + cos_e = cos(MeanObliquity); + + *alpha = atan2(sin(lambda)*cos_e - tan(beta)*sin_e, cos(lambda)); + *delta = asin(sin(beta)*cos_e + cos(beta)*sin_e*sin(lambda)); +} + + +/* computing julian dates (assuming gregorian calendar, thus this is + * only valid for dates of 1582 oct 15 or later) (after duffett-smith, + * section 4) */ + +static double julian_date(int y, int m, int d) { + /* int y; year (e.g. 19xx) */ + /* int m; month (jan=1, feb=2, ...) */ + /* int d; day of month */ + + int A, B, C, D; + double JD; + + /* lazy test to ensure gregorian calendar */ + if (y < 1583) { + printf("WHOOPS! Julian dates only valid for 1582 oct 15 or later\n"); + } + + if ((m == 1) || (m == 2)) { + y -= 1; + m += 12; + } + + A = y / 100; + B = 2 - A + (A / 4); + C = 365.25 * y; + D = 30.6001 * (m + 1); + + JD = B + C + D + d + 1720994.5; + + return JD; +} + + +/* compute greenwich mean sidereal time (GST) corresponding to a given + * number of seconds since the unix epoch (after duffett-smith, + * section 12) */ +static double GST(time_t ssue) { + /* time_t ssue; seconds since unix epoch */ + + double JD; + double T, T0; + double UT; + struct tm *tm; + + tm = gmtime(&ssue); + + JD = julian_date(tm->tm_year+1900, tm->tm_mon+1, tm->tm_mday); + T = (JD - 2451545) / 36525; + + T0 = ((T + 2.5862e-5) * T + 2400.051336) * T + 6.697374558; + + T0 = fmod(T0, 24.0); + if (T0 < 0) T0 += 24; + + UT = tm->tm_hour + (tm->tm_min + tm->tm_sec / 60.0) / 60.0; + + T0 += UT * 1.002737909; + T0 = fmod(T0, 24.0); + if (T0 < 0) T0 += 24; + + return T0; +} + + +/* given a particular time (expressed in seconds since the unix + * epoch), compute position on the earth (lat, lon) such that sun is + * directly overhead. (lat, lon are reported in radians */ + +void fgSunPosition(time_t ssue, double *lon, double *lat) { + /* time_t ssue; seconds since unix epoch */ + /* double *lat; (return) latitude */ + /* double *lon; (return) longitude */ + + double lambda; + double alpha, delta; + double tmp; + + lambda = sun_ecliptic_longitude(ssue); + ecliptic_to_equatorial(lambda, 0.0, &alpha, &delta); + + tmp = alpha - (FG_2PI/24)*GST(ssue); + if (tmp < -FG_PI) { + do tmp += FG_2PI; + while (tmp < -FG_PI); + } else if (tmp > FG_PI) { + do tmp -= FG_2PI; + while (tmp < -FG_PI); + } + + *lon = tmp; + *lat = delta; +} + + +/* $Log$ +/* Revision 1.1 1997/08/01 15:27:56 curt +/* Initial revision. +/* + */ diff --git a/Time/sunpos.h b/Time/sunpos.h new file mode 100644 index 000000000..d153042de --- /dev/null +++ b/Time/sunpos.h @@ -0,0 +1,53 @@ +/* + * sunpos.h + * kirk johnson + * july 1993 + * + * code for calculating the position on the earth's surface for which + * the sun is directly overhead (adapted from _practical astronomy + * with your calculator, third edition_, peter duffett-smith, + * cambridge university press, 1988.) + * + * RCS $Id$ + * + * Copyright (C) 1989, 1990, 1993, 1994, 1995 Kirk Lauritz Johnson + * + * Parts of the source code (as marked) are: + * Copyright (C) 1989, 1990, 1991 by Jim Frost + * Copyright (C) 1992 by Jamie Zawinski + * + * Permission to use, copy, modify and freely distribute xearth for + * non-commercial and not-for-profit purposes is hereby granted + * without fee, provided that both the above copyright notice and this + * permission notice appear in all copies and in supporting + * documentation. + * + * The author makes no representations about the suitability of this + * software for any purpose. It is provided "as is" without express or + * implied warranty. + * + * THE AUTHOR DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, + * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS, + * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, INDIRECT + * OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM + * LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, + * NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN + * CONNECTION WITH THE USE OR PERFORMANCE OF THIS SOFTWARE. + */ + + +#ifndef SUNPOS_H +#define SUNPOS_H + + +#include + + +/* given a particular time (expressed in seconds since the unix + * epoch), compute position on the earth (lat, lon) such that sun is + * directly overhead. (lat, lon are reported in radians */ + +void fgSunPosition(time_t ssue, double *lon, double *lat); + + +#endif /* SUNPOS_H */ -- 2.39.2