/* Topocentric Moon Monthly Calendar Keith Burnett (https://k58.uk/tmoon.c) Compiles using... cc tmoon.c -lm -o tmoon Overview -------- This program calculates the approximate values of various numbers of interest to lunar observers for a month at a time. I have written the program using 'standard' C code and libraries. For each day in the month, the program calculates the time of Moon transit, and then calculates for that time the optical libration of the Moon in selenographical longitude and latitude, the position angle of the Moon's rotation axis, the selenographic colongitude and latitude of the sub-Solar point, i.e. the point on the Moon's surface where the Sun is overhead, the position angle of the bright limb of the Moon (just add and subtract 90 degrees from this to get the position angles of the 'horns' of the crescent), and finally the percentage illumination of the Moon. The program also adds a character after the transit time indicating the state of the Sun, blank for night, 'a', 'n', 'c' for astronomical, nautical or civil twilight, and '*' if the Sun is above the horizon at the time of lunar transit. A further character is printed after the colongitude, an 'r' shows that the Sun is rising on the Moon, and an 's' indicates a waning Moon. Topocentric positions are calculated to your position on the Earth. Geocentric positions are calculated to the centre of the Earth. There is an appreciable difference for objects that are near to us, so that your displacement from the centre of the Earth changes your angle to the object by a significant amount. The quantities here are calculated using topocentric positions for the Moon and geocentric positions for the Sun, so that the resulting values should represent the Moon's face that you will see in your eyepiece fairly well. As the Moon is roughly 60 Earth radii away, the topocentric correction has a large effect on the Moon's position in your sky (sometimes up to one degree). The Sun is roughly 140,000 Earth radii away, so I assume that the topocentric correction to the Sun's position is so small as to have no appreciable effect on the illumination of the Moon's face. Your location moves round the Earth's rotation axis in a 'small circle' once a day, and the radius of the small circle depends on your latitude. The geocentric optical libration changes on a scale of the lunar month, as it arises in part from the elliptical nature of the Moon's orbit around the Earth. The size of the difference between the topocentric and geocentric librations will therefore change during a day as you move round your small circle. There will also be a 'constant' or more slowly changing topocentric correction depending on your latitude and the declination of the Moon. I imagine a smooth slow cycle with a small amplitude fast 'wiggle' superimposed on it from the topocentric librations. I have a strong suspicion that the diurnal component of the topocentric correction to the Moon's position is actually zero when the Moon transits, as at that time the centre of the Moon, you and the centre of the Earth are in a plane that includes your meridian. Most published books and articles about Moon observation will record the geocentric librations of the Moon, and these will be different from the figures produced by this program. Accuracy -------- I checked a modified version of this program calculating quantities for 0h UT each day against a topocentric ephemeris generated using the JPL Horizons online ephemeris generator. Quantities agreed to within 0.03 degrees for optical librations (as might be expected as that is the level at which the physical librations become important) and to within better than 0.005 for the colongitude. I chose only one month to check, and a graph of the colongitude errors looked as if there might be a periodic term of amplitude about 0.02 degrees or so. All output has been rounded to one decimal place for the selenographic coordinates, and whole degrees for the position angles. An earlier version of the program used the D46 formulas (0.3 degrees error in position) for the Moon coordinates - and I found the errors in librations were three times as bad. The current moonpos() function is good to 2 arcmin 'most of the time', and 5 arcmin worst case (odd spikes in the error signal over a complete saros). Algorithm --------- The 'logic' of the main() function looks a bit like this: Get inputs Do some error checking print title and column headings for each day in the month find the time of lunar transit using an iteration loop if no transit, print message and skip to next day find the days since J2000.0 corresponding to this time calculate the ecliptic coordinates of the Sun convert to equatorial coordinates find the sun's altitude code sun's altitude using a character ' anc*' find the ecliptic coordinates of the Moon convert to equatorial coordinates find altitude of Moon if Moon below horizon at transit print message and skip to next day apply topocentric correction to Moon equatorial coords convert back to get topo ecliptic coords find librations and pa of pole axis find approximate heliocentric coordinates of Moon use libration procedure to find coords of sub-solar point find if Sun rising 'r' or setting 's' over Moon find pa of bright limb using topocentric Moon coords find % illumination print line of output for that day next day print key to table end I have used 'pass by reference' a lot using pointers in the functions to get round the limitation that C functions can only return a single argument. Perhaps some structures for collections of coordinates would be more elegant. Command line ------------ The program takes a command line with four purely numerical arguments, - the year and month in the format yyyymm, - your time zone as a decimal number of hours before or after Greenwich, before Greenwich taken as negative - your longitude as a signed integer in the format +/-dddmm with West longitudes taken as negative, - your latitude as a signed integer in the format +/-ddmm with North latitudes taken as positive. For example, the command line >tmoon 200009 0 -00155 5230 > sept.txt will calculate a monthly Moon calendar for September 2000, in time zone 0 (UT) for longitude 1 deg 55 minutes West, and 52 deg 30 minutes North. In this example, the symbol '>' tells the operating system to 're-direct' the program output to a file called 'sept.txt'. Redirection of output depends on the compiler and operating system you are using - most support this feature. Error messages are printed to the console using fprintf(stderr, stuff). References ---------- The formulas for the Moon's geocentric coordinates, the siderial time, days since J2000.0, coordinate transformations and all the formulas for calculating the output quantities were modified and simplified from Jean Meeus' book 'Astronomical Algorithms' (1st Edition) published by Willmann-Bell, ISBN 0-943396-35-2. Chapter 51 was especially useful, but I have neglected the physical libration (or 'free space motion') of the Moon and the nutation, resulting in a simplification of the formulas. The formulas for the Sun's geocentric coordinates, and the common sense topocentric correction for the Moon's position (ignoring the polar flattening of Earth), were taken from the Astronomical Almanac pages C24 and D46. The iterative scheme to find the time of Moon transit each day was adapted from the 'Explanatory Supplement to the Astronomical Almanac', 1994, section 9.31 and 9.32. Mealy Mouthed Disclaimer ------------------------ This program was written as a bit of a giggle, and to improve my own understanding of grunge C and selenography. I wouldn't trust it for anything important if I were you. I take no responsibility for anything that might happen as a consequence of your use of this program or any bit of it for anything anywhere at anytime. I would like to hear from you if you find a mistake, or find a nice new application for bits of the code. */ #include #include #include #define RADS 0.0174532925199433 #define DEGS 57.2957795130823 #define TPI 6.28318530717959 #define PI 3.1415927 /* ratio of earth radius to astronomical unit */ #define ER_OVER_AU 0.0000426352325194252 /* all prototypes here */ double getcoord(int coord); void getargs(int argc, char *argv[], int *y, int *m, double *tz, double *glong, double *glat); double range(double y); double rangerad(double y); double days(int y, int m, int dn, double hour); void moonpos(double, double *, double *, double *); void sunpos(double , double *, double *, double *); double moontransit(int y, int m, int d, double timezone, double glat, double glong, int *nt); double atan22(double y, double x); double epsilon(double d); void equatorial(double d, double *lon, double *lat, double *r); void ecliptic(double d, double *lon, double *lat, double *r); double gst(double d); void topo(double lst, double glat, double *alp, double *dec, double *r); double alt(double glat, double ha, double dec); void libration(double day, double lambda, double beta, double alpha, double *l, double *b, double *p); void illumination(double day, double lra, double ldec, double dr, double sra, double sdec, double *pabl, double *ill); int daysinmonth(int y, int m); int isleap(int y); static const char *usage = " Usage: tmoon date[yyyymm] timz[+/-h.hh] long[+/-dddmm] lat[+/-ddmm]\n" "example: tmoon 200009 0 -00155 5230\n"; static const char *keytext = "\n\nKey to column headings\n" "----------------------\n\n" " DY: Date in the month\n" "TRAN: Zone time of Moon transit, followed by a letter indicating\n" " the sun state, * means the Sun is up, blank is full night\n" " and c, n, a stand for civil, nautical and astronomical twilights\n" " ALT: Altitude of the Moon at transit\n" " L: Optical libration in longitude\n" " B: Libration in latitude\n" " PaX: Position angle of Moon's polar axis\n" " Co: Colongitude of the Sun, r means sunrise, s means lunar sunset\n" " Bs: Selenographic latitude of the subsolar point\n" "PaBl: Position angle of the bright limb\n" " %%: Percentage illumination of Moon's earth-facing disc\n"; int main(int argc, char *argv[]) { double mlambda, mbeta, mrv, slambda, sbeta, srv; double tlambda, tbeta, trv; double tz, dt, glong, glat, time; double lst, sunalt, shr, day, Co; double salpha, sdelta, malpha, mdelta, mhr, moonalt, l, b; double ls, bs, hlambda, hbeta, dratio, paxis, dummy, pabl, ill; int nt = 0, date, hour, min, y, m, longitude, latitude; char sunchar, moonchar, eastwest, northsouth; /* get the date, time zone, and observer's position */ getargs(argc, argv, &y, &m, &tz, &glong, &glat); /* print title and column headings */ printf("Topocentric lunar ephemeris\n"); printf("Year: %04d Month: %02d\n", y, m); longitude = atoi(argv[3]); latitude = atoi(argv[4]); if (longitude < 0) eastwest = 'W'; else eastwest = 'E'; if (latitude < 0) northsouth = 'S'; else northsouth = 'N'; printf("Location: %05d%1c, %04d%1c\n\n", abs(longitude), eastwest, abs(latitude), northsouth); printf("DY TRAN ALT L B PaX Co Bs PaBL %%\n"); printf("-------------------------------------------------------\n"); /* this is the main month loop - generates each line in the table */ for(date = 1; date <= daysinmonth(y, m); date++) { /* find the time of transit on this day */ dt = moontransit(y, m, date, tz, glat, glong, &nt); if (nt == 1) { printf("%02d ---- no transit this day\n", date); /* skip to next day */ continue; } day = days(y, m , date, dt * DEGS / 15) - tz/24; lst = gst(day) + glong; /* find Moon topocentric coordinates for libration calculations and check that the Moon is above the horizon at transit */ moonpos(day, &mlambda, &mbeta, &mrv); malpha = mlambda; mdelta = mbeta; equatorial(day, &malpha, &mdelta, &mrv); topo(lst, glat, &malpha, &mdelta, &mrv); mhr = rangerad(lst - malpha); moonalt = alt(glat, mhr, mdelta); if (moonalt < 0 * RADS) { printf(" %02d **** moon below horizon all day\n", date); /* skip to next day */ continue; } /* find sun altitude character */ sunpos(day, &slambda, &sbeta, &srv); salpha = slambda; sdelta = sbeta; equatorial(day, &salpha, &sdelta, &srv); shr = rangerad(lst - salpha); sunalt = alt(glat, shr, sdelta); if (sunalt > 0 * RADS) sunchar = '*'; if (sunalt < 0 * RADS) sunchar = 'c'; if (sunalt < - 6 * RADS) sunchar = 'n'; if (sunalt < -12 * RADS) sunchar = 'a'; if (sunalt < -18 * RADS) sunchar = ' '; /* Optical libration and Position angle of the Pole */ tlambda = malpha; tbeta = mdelta; trv = mrv; ecliptic(day, &tlambda, &tbeta, &trv); libration(day, tlambda, tbeta, malpha, &l, &b, &paxis); /* Selen Colongitude and latitude of sub solar point */ dratio = mrv / srv * ER_OVER_AU; hlambda = slambda + PI + dratio * cos(mbeta) * sin(slambda - mlambda); hbeta = dratio * mbeta; libration(day, hlambda, hbeta, salpha, &ls, &bs, &dummy); ls = rangerad(ls); if(ls < 90 * RADS) Co = 90 * RADS - ls; else Co = 450 * RADS - ls; if(Co < 90 * RADS || Co > 270 * RADS) moonchar = 'r'; else moonchar = 's'; /* PA of bright limb, and percentage illumination */ illumination(day, malpha, mdelta, dratio, salpha, sdelta, &pabl, &ill); /* convert transit time to hhmm */ time = dt * DEGS / 15; hour = floor(time); min = floor((time - hour) * 60 + 0.5); /* Print the line for this day */ printf("%02d %02d%02d", date, hour, min); printf("%1c ", sunchar); printf("%3.0f ", moonalt * DEGS); printf("% 3.1f % 2.1f ", l * DEGS, b * DEGS); printf("%3.0f ", paxis * DEGS); printf("%5.1f%1c % 2.1f ", Co * DEGS, moonchar, bs * DEGS); printf("%4.0f %3.0f\n", pabl * DEGS, ill * 100); } /* end of day in Month loop */ printf(keytext); exit(EXIT_SUCCESS); } /* end of main() */ /* getargs() gets the arguments from the command line, does some basic error checking, and converts arguments into numerical form. Arguments are passed back in pointers. Error messages print to stderr so re-direction of output to file won't leave users blind. Error checking prints list of all errors in a command line before quitting. */ void getargs(int argc, char *argv[], int *y, int *m, double *tz, double *glong, double *glat) { int date, latitude, longitude; int mflag = 0, yflag = 0, longflag = 0, latflag = 0, tzflag = 0; int longminflag = 0, latminflag = 0, dflag = 0; /* if not right number of arguments, then print example command line */ if (argc !=5) { fprintf(stderr, usage); exit(EXIT_FAILURE); } date = atoi(argv[1]); *y = date / 100; *m = date - *y * 100; *tz = (double) atof(argv[2]); longitude = atoi(argv[3]); latitude = atoi(argv[4]); *glong = RADS * getcoord(longitude); *glat = RADS * getcoord(latitude); /* set a flag for each error found */ if (*m > 12 || *m < 1) mflag = 1; if (*y > 2500) yflag = 1; if (date < 150001) dflag = 1; if (fabs((float) *glong) > 180 * RADS) longflag = 1; if (abs(longitude) % 100 > 59) longminflag = 1; if (fabs((float) *glat) > 90 * RADS) latflag = 1; if (abs(latitude) % 100 > 59) latminflag = 1; if (fabs((float) *tz) > 12) tzflag = 1; /* print all the errors found */ if (dflag == 1) { fprintf(stderr, "date: dates must be in form yyyymm, gregorian, and later than 1500 AD\n"); } if (yflag == 1) { fprintf(stderr, "date: too far in future - accurate from 1500 to 2500\n"); } if (mflag == 1) { fprintf(stderr, "date: month must be in range 0 to 12, eg - August 2000 is entered as 200008\n"); } if (tzflag == 1) { fprintf(stderr, "timz: must be in range +/- 12 hours, eg -6 for Chicago\n"); } if (longflag == 1) { fprintf(stderr, "long: must be in range +/- 180 degrees\n"); } if (longminflag == 1) { fprintf(stderr, "long: last two digits are arcmin - max 59\n"); } if (latflag == 1) { fprintf(stderr, " lat: must be in range +/- 90 degrees\n"); } if (latminflag == 1) { fprintf(stderr, " lat: last two digits are arcmin - max 59\n"); } /* quits if one or more flags set */ if (dflag + mflag + yflag + longflag + latflag + tzflag + longminflag + latminflag > 0) { exit(EXIT_FAILURE); } } /* returns coordinates in decimal degrees given the coord as a ddmm value stored in an integer. */ double getcoord(int coord) { int west = 1; double glg, deg; if (coord < 0) west = -1; glg = fabs((double) coord/100); deg = floor(glg); glg = west* (deg + (glg - deg)*100 / 60); return(glg); } /* days() takes the year, month, day in the month and decimal hours in the day and returns the number of days since J2000.0. Assumes Gregorian calendar. */ double days(int y, int m, int d, double h) { int a, b; double day; /* The lines below work from 1900 march to feb 2100 a = 367 * y - 7 * (y + (m + 9) / 12) / 4 + 275 * m / 9 + d; day = (double)a - 730531.5 + hour / 24; */ /* These lines work for any Gregorian date since 0 AD */ if (m ==1 || m==2) { m +=12; y -= 1; } a = y / 100; b = 2 - a + a/4; day = floor(365.25*(y + 4716)) + floor(30.6001*(m + 1)) + d + b - 1524.5 - 2451545 + h/24; return(day); } /* Returns 1 if y a leap year, and 0 otherwise, according to the Gregorian calendar */ int isleap(int y) { int a = 0; if(y % 4 == 0) a = 1; if(y % 100 == 0) a = 0; if(y % 400 == 0) a = 1; return(a); } /* Given the year and the month, function returns the number of days in the month. Valid for Gregorian calendar. */ int daysinmonth(int y, int m) { int b = 31; if(m == 2) { if(isleap(y) == 1) b= 29; else b = 28; } if(m == 4 || m == 6 || m == 9 || m == 11) b = 30; return(b); } /* moonpos() takes days from J2000.0 and returns ecliptic coordinates of moon in the pointers. Note call by reference. This function is within a couple of arcminutes most of the time, and is truncated from the Meeus Ch45 series, themselves truncations of ELP-2000. Returns moon distance in earth radii. Terms have been written out explicitly rather than using the table based method as only a small number of terms is retained. */ void moonpos(double d, double *lambda, double *beta, double *rvec) { double dl, dB, dR, L, D, M, M1, F, e, lm, bm, rm, t; t = d / 36525; L = range(218.3164591 + 481267.88134236 * t) * RADS; D = range(297.8502042 + 445267.1115168 * t) * RADS; M = range(357.5291092 + 35999.0502909 * t) * RADS; M1 = range(134.9634114 + 477198.8676313 * t - .008997 * t * t) * RADS; F = range(93.27209929999999 + 483202.0175273 * t - .0034029 * t * t) * RADS; e = 1 - .002516 * t; dl = 6288774 * sin(M1); dl += 1274027 * sin(2 * D - M1); dl += 658314 * sin(2 * D); dl += 213618 * sin(2 * M1); dl -= e * 185116 * sin(M); dl -= 114332 * sin(2 * F) ; dl += 58793 * sin(2 * D - 2 * M1); dl += e * 57066 * sin(2 * D - M - M1) ; dl += 53322 * sin(2 * D + M1); dl += e * 45758 * sin(2 * D - M); dl -= e * 40923 * sin(M - M1); dl -= 34720 * sin(D) ; dl -= e * 30383 * sin(M + M1) ; dl += 15327 * sin(2 * D - 2 * F) ; dl -= 12528 * sin(M1 + 2 * F); dl += 10980 * sin(M1 - 2 * F); lm = rangerad(L + dl / 1000000 * RADS); dB = 5128122 * sin(F); dB += 280602 * sin(M1 + F); dB += 277693 * sin(M1 - F); dB += 173237 * sin(2 * D - F); dB += 55413 * sin(2 * D - M1 + F); dB += 46271 * sin(2 * D - M1 - F); dB += 32573 * sin(2 * D + F); dB += 17198 * sin(2 * M1 + F); dB += 9266 * sin(2 * D + M1 - F); dB += 8822 * sin(2 * M1 - F); dB += e * 8216 * sin(2 * D - M - F); dB += 4324 * sin(2 * D - 2 * M1 - F); bm = dB / 1000000 * RADS; dR = -20905355 * cos(M1); dR -= 3699111 * cos(2 * D - M1); dR -= 2955968 * cos(2 * D); dR -= 569925 * cos(2 * M1); dR += e * 48888 * cos(M); dR -= 3149 * cos(2 * F); dR += 246158 * cos(2 * D - 2 * M1); dR -= e * 152138 * cos(2 * D - M - M1) ; dR -= 170733 * cos(2 * D + M1); dR -= e * 204586 * cos(2 * D - M); dR -= e * 129620 * cos(M - M1); dR += 108743 * cos(D); dR += e * 104755 * cos(M + M1); dR += 79661 * cos(M1 - 2 * F); rm = 385000.56 + dR / 1000; *lambda = lm; *beta = bm; /* distance to Moon must be in Earth radii */ *rvec = rm / 6378.14; } /* topomoon() takes the local siderial time, the geographical latitude of the observer, and pointers to the geocentric equatorial coordinates. The function overwrites the geocentric coordinates with topocentric coordinates on a simple spherical earth model (no polar flattening). Expects Moon-Earth distance in Earth radii. Formulas scavenged from Astronomical Almanac 'low precision formulae for Moon position' page D46. */ void topo(double lst, double glat, double *alp, double *dec, double *r) { double x, y, z, r1; x = *r * cos(*dec) * cos(*alp) - cos(glat) * cos(lst); y = *r * cos(*dec) * sin(*alp) - cos(glat) * sin(lst); z = *r * sin(*dec) - sin(glat); r1 = sqrt(x*x + y*y + z*z); *alp = atan22(y, x); *dec = asin(z / r1); *r = r1; } /* moontransit() takes date, the time zone and geographic longitude of observer and returns the time (decimal hours) of lunar transit on that day if there is one, and sets the notransit flag if there isn't. See Explanatory Supplement to Astronomical Almanac section 9.32 and 9.31 for the method. */ double moontransit(int y, int m, int d, double tz, double glat, double glong, int *notransit) { double hm, ht, ht1, lon, lat, rv, dnew, lst; int itcount; ht1 = 180 * RADS; ht = 0; itcount = 0; *notransit = 0; do { ht = ht1; itcount++; dnew = days(y, m, d, ht * DEGS/15) - tz/24; lst = gst(dnew) + glong; /* find the topocentric Moon ra (hence hour angle) and dec */ moonpos(dnew, &lon, &lat, &rv); equatorial(dnew, &lon, &lat, &rv); topo(lst, glat, &lon, &lat, &rv); hm = rangerad(lst - lon); ht1 = rangerad(ht - hm); /* if no convergence, then no transit on that day */ if (itcount > 30) { *notransit = 1; break; } } while (fabs(ht - ht1) > 0.04 * RADS); return(ht1); } /* Calculates the selenographic coordinates of either the sub Earth point (optical libration) or the sub-solar point (selen. coords of centre of bright hemisphere). Based on Meeus chapter 51 but neglects physical libration and nutation, with some simplification of the formulas. */ void libration(double day, double lambda, double beta, double alpha, double *l, double *b, double *p) { double i, f, omega, w, y, x, a, t, eps; t = day / 36525; i = 1.54242 * RADS; eps = epsilon(day); f = range(93.2720993 + 483202.0175273 * t - .0034029 * t * t) * RADS; omega = range(125.044555 - 1934.1361849 * t + .0020762 * t * t) * RADS; w = lambda - omega; y = sin(w) * cos(beta) * cos(i) - sin(beta) * sin(i); x = cos(w) * cos(beta); a = atan22(y, x); *l = a - f; /* kludge to catch cases of 'round the back' angles */ if (*l < -90 * RADS) *l += TPI; if (*l > 90 * RADS) *l -= TPI; *b = asin(-sin(w) * cos(beta) * sin(i) - sin(beta) * cos(i)); /* pa pole axis - not used for Sun stuff */ x = sin(i) * sin(omega); y = sin(i) * cos(omega) * cos(eps) - cos(i) * sin(eps); w = atan22(x, y); *p = rangerad(asin(sqrt(x*x + y*y) * cos(alpha - w) / cos(*b))); } /* Takes: days since J2000.0, eq coords Moon, ratio of moon to sun distance, eq coords Sun Returns: position angle of bright limb wrt NCP, percentage illumination of Sun */ void illumination(double day, double lra, double ldec, double dr, double sra, double sdec, double *pabl, double *ill) { double x, y, phi, i; y = cos(sdec) * sin(sra - lra); x = sin(sdec) * cos(ldec) - cos(sdec) * sin(ldec) * cos (sra - lra); *pabl = atan22(y, x); phi = acos(sin(sdec) * sin(ldec) + cos(sdec) * cos(ldec) * cos(sra-lra)); i = atan22(sin(phi) , (dr - cos(phi))); *ill = 0.5*(1 + cos(i)); } /* sunpos() takes days from J2000.0 and returns ecliptic longitude of Sun in the pointers. Latitude is zero at this level of precision, but pointer left in for consistency in number of arguments. This function is within 0.01 degree (1 arcmin) almost all the time for a century either side of J2000.0. This is from the 'low precision fomulas for the Sun' from C24 of Astronomical Alamanac */ void sunpos(double d, double *lambda, double *beta, double *rvec) { double L, g, ls, bs, rs; L = range(280.461 + .9856474 * d) * RADS; g = range(357.528 + .9856003 * d) * RADS; ls = L + (1.915 * sin(g) + .02 * sin(2 * g)) * RADS; bs = 0; rs = 1.00014 - .01671 * cos(g) - .00014 * cos(2 * g); *lambda = ls; *beta = bs; *rvec = rs; } /* this routine returns the altitude given the days since J2000.0 the hour angle and declination of the object and the latitude of the observer. Used to find the Sun's altitude to put a letter code on the transit time, and to find the Moon's altitude at transit just to make sure that the Moon is visible. */ double alt(double glat, double ha, double dec) { return(asin(sin(dec) * sin(glat) + cos(dec) * cos(glat) * cos(ha))); } /* returns an angle in degrees in the range 0 to 360 */ double range(double x) { double a, b; b = x / 360; a = 360 * (b - floor(b)); if (a < 0) a = 360 + a; return(a); } /* returns an angle in rads in the range 0 to two pi */ double rangerad(double x) { double a, b; b = x / TPI; a = TPI * (b - floor(b)); if (a < 0) a = TPI + a; return(a); } /* gets the atan2 function returning angles in the right order and range */ double atan22(double y, double x) { double a; a = atan2(y, x); if (a < 0) a += TPI; return(a); } /* returns mean obliquity of ecliptic in radians given days since J2000.0. */ double epsilon(double d) { double t = d/ 36525; return((23.4392911111111 - (t* (46.8150 + 0.00059*t)/3600)) *RADS); } /* replaces ecliptic coordinates with equatorial coordinates note: call by reference destroys original values R is unchanged. */ void equatorial(double d, double *lon, double *lat, double *r) { double eps, ceps, seps, l, b; l = *lon; b = * lat; eps = epsilon(d); ceps = cos(eps); seps = sin(eps); *lon = atan22(sin(l)*ceps - tan(b)*seps, cos(l)); *lat = asin(sin(b)*ceps + cos(b)*seps*sin(l)); } /* replaces equatorial coordinates with ecliptic ones. Inverse of above, but used to find topocentric ecliptic coords. */ void ecliptic(double d, double *lon, double *lat, double *r) { double eps, ceps, seps, alp, dec; alp = *lon; dec = *lat; eps = epsilon(d); ceps = cos(eps); seps = sin(eps); *lon = atan22(sin(alp)*ceps + tan(dec)*seps, cos(alp)); *lat = asin(sin(dec)*ceps - cos(dec)*seps*sin(alp)); } /* returns the siderial time at greenwich meridian as an angle in radians given the days since J2000.0 */ double gst( double d) { double t = d / 36525; double theta; theta = range(280.46061837 + 360.98564736629 * d + 0.000387933 * t * t); return(theta * RADS); }