#### 420 lines 9.6 KiB C Raw Blame History

 ```#include ``` ```#include ``` ```#include "CalcEphem.h" ``` ``` ``` ```int DayofYear(int, int, int); ``` ```int DayofWeek(int, int, int, char[]); ``` ```double kepler(double, double); ``` ```double jd(int, int, int, double); ``` ```double hour24(double); ``` ```double angle2pi(double); ``` ```double angle360(double); ``` ```void Radec_to_Cart(double, double, Vector *); ``` ```int LeapYear(int); ``` ``` ``` ```int DayofYear(year, month, day) ``` ```int year, month, day; ``` ```{ ``` ``` double jd(); ``` ``` return((int)(jd(year, month, day, 0.0) - jd(year, 1, 0, 0.0))); ``` ```} ``` ``` ``` ``` ``` ``` ``` ``` ``` ```int DayofWeek(year, month, day, dowstr) ``` ```int year, month, day; ``` ```char dowstr[]; ``` ```{ ``` ``` double JD, A, Afrac, jd(); ``` ``` int n, iA; ``` ``` ``` ``` JD = jd(year, month, day, 0.0); ``` ``` A = (JD + 1.5)/7.0; ``` ``` iA = (int)A; ``` ``` Afrac = A - (double)iA; ``` ``` n = (int)(Afrac*7.0 + 0.5); ``` ``` switch(n){ ``` ``` case 0: ``` ``` strcpy(dowstr, "Sunday"); ``` ``` break; ``` ``` case 1: ``` ``` strcpy(dowstr, "Monday"); ``` ``` break; ``` ``` case 2: ``` ``` strcpy(dowstr, "Tuesday"); ``` ``` break; ``` ``` case 3: ``` ``` strcpy(dowstr, "Wednesday"); ``` ``` break; ``` ``` case 4: ``` ``` strcpy(dowstr, "Thursday"); ``` ``` break; ``` ``` case 5: ``` ``` strcpy(dowstr, "Friday"); ``` ``` break; ``` ``` case 6: ``` ``` strcpy(dowstr, "Saturday"); ``` ``` break; ``` ``` } ``` ``` return(n); ``` ```} ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ```void ``` ```CalcEphem(date, UT, c) ``` ```long int date; /* integer containing the date (e.g. 960829) */ ``` ```double UT; /* Universal Time */ ``` ```CTrans *c; /* structure containing all the relevent coord trans info */ ``` ```{ ``` ``` int year, month, day; ``` ``` double TU, TU2, TU3, T0, gmst; ``` ``` double varep, varpi; ``` ``` double eccen, epsilon; ``` ``` double days, M, E, nu, lambnew; ``` ``` double r0, earth_sun_distance; ``` ``` double RA, DEC, RA_Moon, DEC_Moon; ``` ``` double TDT; ``` ``` double AGE; ``` ``` double LambdaMoon, BetaMoon, R; ``` ``` double jd(), hour24(), angle2pi(), angle360(), kepler(), Moon(), NewMoon(); ``` ``` double Ta, Tb, Tc, frac(); ``` ``` double SinGlat, CosGlat, SinGlon, CosGlon, Tau, lmst, x, y, z; ``` ``` double SinTau, CosTau, SinDec, CosDec; ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` c->UT = UT; ``` ``` year = (int)(date/10000); ``` ``` month = (int)( (date - year*10000)/100 ); ``` ``` day = (int)( date - year*10000 - month*100 ); ``` ``` c->year = year; ``` ``` c->month = month; ``` ``` c->day = day; ``` ``` ``` ``` c->doy = DayofYear(year, month, day); ``` ``` c->dow = DayofWeek(year, month, day, c->dowstr); ``` ``` ``` ``` ``` ``` ``` ``` /* ``` ``` * Compute Greenwich Mean Sidereal Time (gmst) ``` ``` * The TU here is number of Julian centuries ``` ``` * since 2000 January 1.5 ``` ``` * From the 1996 astronomical almanac ``` ``` */ ``` ``` TU = (jd(year, month, day, 0.0) - 2451545.0)/36525.0; ``` ``` TU2 = TU*TU; ``` ``` TU3 = TU2*TU; ``` ``` T0 = (6.0 + 41.0/60.0 + 50.54841/3600.0) + 8640184.812866/3600.0*TU ``` ``` + 0.093104/3600.0*TU2 - 6.2e-6/3600.0*TU3; ``` ``` T0 = hour24(T0); ``` ``` c->gmst = hour24(T0 + UT*1.002737909); ``` ``` ``` ``` ``` ``` /* convert to radians for ease later on */ ``` ``` gmst = c->gmst*15.0*M_PI/180.0; ``` ``` ``` ``` lmst = 24.0*frac( (c->gmst - c->Glon/15.0) / 24.0 ); ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` /* ``` ``` * ``` ``` * Construct Transformation Matrix from GEI to GSE systems ``` ``` * ``` ``` * ``` ``` * First compute: ``` ``` * mean ecliptic longitude of sun at epoch TU (varep) ``` ``` * elciptic longitude of perigee at epoch TU (varpi) ``` ``` * eccentricity of orbit at epoch TU (eccen) ``` ``` * ``` ``` * The TU here is the number of Julian centuries since ``` ``` * 1900 January 0.0 (= 2415020.0) ``` ``` */ ``` ``` TDT = UT + 59.0/3600.0; ``` ``` TU = (jd(year, month, day, TDT) - 2415020.0)/36525.0; ``` ``` varep = (279.6966778 + 36000.76892*TU + 0.0003025*TU*TU)*RadPerDeg; ``` ``` varpi = (281.2208444 + 1.719175*TU + 0.000452778*TU*TU)*RadPerDeg; ``` ``` eccen = 0.01675104 - 0.0000418*TU - 0.000000126*TU*TU; ``` ``` c->eccentricity = eccen; ``` ``` ``` ``` ``` ``` ``` ``` /* ``` ``` * Compute the Obliquity of the Ecliptic at epoch TU ``` ``` * The TU in this formula is the number of Julian ``` ``` * centuries since epoch 2000 January 1.5 ``` ``` */ ``` ``` TU = (jd(year, month, day, TDT) - jd(2000, 1, 1, 12.0))/36525.0; ``` ``` epsilon = (23.43929167 - 0.013004166*TU - 1.6666667e-7*TU*TU ``` ``` - 5.0277777778e-7*TU*TU*TU)*RadPerDeg; ``` ``` c->epsilon = epsilon; ``` ``` ``` ``` ``` ``` /* ``` ``` * Compute: ``` ``` * Number of Days since epoch 1990.0 (days) ``` ``` * The Mean Anomaly (M) ``` ``` * The True Anomaly (nu) ``` ``` * The Eccentric Anomaly via Keplers equation (E) ``` ``` * ``` ``` * ``` ``` */ ``` ``` days = jd(year, month, day, TDT) - jd(year, month, day, TDT); ``` ``` M = angle2pi(2.0*M_PI/365.242191*days + varep - varpi); ``` ``` E = kepler(M, eccen); ``` ``` nu = 2.0*atan( sqrt((1.0+eccen)/(1.0-eccen))*tan(E/2.0) ); ``` ``` lambnew = angle2pi(nu + varpi); ``` ``` c->lambda_sun = lambnew; ``` ``` ``` ``` ``` ``` /* ``` ``` * Compute distance from earth to the sun ``` ``` */ ``` ``` r0 = 1.495985e8; /* in km */ ``` ``` earth_sun_distance = r0*(1-eccen*eccen)/(1.0 + eccen*cos(nu))/6371.2; ``` ``` c->earth_sun_dist = earth_sun_distance; ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` /* ``` ``` * Compute Right Ascension and Declination of the Sun ``` ``` */ ``` ``` RA = angle360(atan2(sin(lambnew)*cos(epsilon), cos(lambnew))*180.0/M_PI); ``` ``` DEC = asin(sin(epsilon)*sin(lambnew))*180.0/M_PI; ``` ``` c->RA_sun = RA; ``` ``` c->DEC_sun = DEC; ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` /* ``` ``` * Compute Moon Phase and AGE Stuff. The AGE that comes out of Moon() ``` ``` * is actually the Phase converted to days. Since AGE is actually defined ``` ``` * to be time since last NewMoon, we need to figure out what the JD of the ``` ``` * last new moon was. Thats done below.... ``` ``` */ ``` ``` TU = (jd(year, month, day, TDT) - 2451545.0)/36525.0; ``` ``` c->MoonPhase = Moon(TU, &LambdaMoon, &BetaMoon, &R, &AGE); ``` ``` LambdaMoon *= RadPerDeg; ``` ``` BetaMoon *= RadPerDeg; ``` ``` ``` ``` ``` ``` RA_Moon = angle360(atan2(sin(LambdaMoon)*cos(epsilon)-tan(BetaMoon)*sin(epsilon), cos(LambdaMoon))*DegPerRad); ``` ``` DEC_Moon = asin( sin(BetaMoon)*cos(epsilon) + cos(BetaMoon)*sin(epsilon)*sin(LambdaMoon))*DegPerRad; ``` ``` c->RA_moon = RA_Moon; ``` ``` c->DEC_moon = DEC_Moon; ``` ``` ``` ``` ``` ``` /* ``` ``` * Compute Alt/Az coords ``` ``` */ ``` ``` Tau = (15.0*lmst - RA_Moon)*RadPerDeg; ``` ``` CosGlat = cos(c->Glat*RadPerDeg); SinGlat = sin(c->Glat*RadPerDeg); ``` ``` CosGlon = cos(c->Glon*RadPerDeg); SinGlon = sin(c->Glon*RadPerDeg); ``` ``` CosTau = cos(Tau); SinTau = sin(Tau); ``` ``` SinDec = sin(DEC_Moon*RadPerDeg); CosDec = cos(DEC_Moon*RadPerDeg); ``` ``` x = CosDec*CosTau*SinGlat - SinDec*CosGlat; ``` ``` y = CosDec*SinTau; ``` ``` z = CosDec*CosTau*CosGlat + SinDec*SinGlat; ``` ``` c->A_moon = DegPerRad*atan2(y, x); ``` ``` c->h_moon = DegPerRad*asin(z); ``` ``` c->Visible = (c->h_moon < 0.0) ? 0 : 1; ``` ``` ``` ``` ``` ``` ``` ``` /* ``` ``` * Compute accurate AGE of the Moon ``` ``` */ ``` ``` Tb = TU - AGE/36525.0; /* should be very close to minimum */ ``` ``` Ta = Tb - 0.4/36525.0; ``` ``` Tc = Tb + 0.4/36525.0; ``` ``` c->MoonAge = (TU - NewMoon(Ta, Tb, Tc))*36525.0; ``` ``` ``` ``` ``` ``` ``` ``` /* ``` ``` * Compute Earth-Moon distance ``` ``` */ ``` ``` c->EarthMoonDistance = R; ``` ``` ``` ``` ``` ``` ``` ```} ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ```double kepler(M, e) ``` ```double M, e; ``` ```{ ``` ``` int n=0; ``` ``` double E, Eold, eps = 1.0e-8; ``` ``` ``` ``` ``` ``` ``` ``` E = M + e*sin(M); ``` ``` do{ ``` ``` Eold = E; ``` ``` E = Eold + (M-Eold+e*sin(Eold)) ``` ``` /(1.0-e*cos(Eold)); ``` ``` ++n; ``` ``` }while((fabs(E-Eold) > eps) && (n < 100)); ``` ``` return(E); ``` ```} ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ```/* ``` ``` * Compute the Julian Day number for the given date. ``` ``` * Julian Date is the number of days since noon of Jan 1 4713 B.C. ``` ``` */ ``` ```double jd(ny, nm, nd, UT) ``` ```int ny, nm, nd; ``` ```double UT; ``` ```{ ``` ``` double A, B, C, D, JD, day; ``` ``` ``` ``` day = nd + UT/24.0; ``` ``` ``` ``` ``` ``` if ((nm == 1) || (nm == 2)){ ``` ``` ny = ny - 1; ``` ``` nm = nm + 12; ``` ``` } ``` ``` ``` ``` if (((double)ny+nm/12.0+day/365.25)>=(1582.0+10.0/12.0+15.0/365.25)){ ``` ``` A = ((int)(ny / 100.0)); ``` ``` B = 2.0 - A + (int)(A/4.0); ``` ``` } ``` ``` else{ ``` ``` B = 0.0; ``` ``` } ``` ``` ``` ``` if (ny < 0.0){ ``` ``` C = (int)((365.25*(double)ny) - 0.75); ``` ``` } ``` ``` else{ ``` ``` C = (int)(365.25*(double)ny); ``` ``` } ``` ``` ``` ``` D = (int)(30.6001*(double)(nm+1)); ``` ``` ``` ``` ``` ``` JD = B + C + D + day + 1720994.5; ``` ``` return(JD); ``` ``` ``` ```} ``` ``` ``` ```double hour24(hour) ``` ```double hour; ``` ```{ ``` ``` int n; ``` ``` ``` ``` if (hour < 0.0){ ``` ``` n = (int)(hour/24.0) - 1; ``` ``` return(hour-n*24.0); ``` ``` } ``` ``` else if (hour > 24.0){ ``` ``` n = (int)(hour/24.0); ``` ``` return(hour-n*24.0); ``` ``` } ``` ``` else{ ``` ``` return(hour); ``` ``` } ``` ```} ``` ``` ``` ```double angle2pi(angle) ``` ```double angle; ``` ```{ ``` ``` int n; ``` ``` double a; ``` ``` a = 2.0*M_PI; ``` ``` ``` ``` if (angle < 0.0){ ``` ``` n = (int)(angle/a) - 1; ``` ``` return(angle-n*a); ``` ``` } ``` ``` else if (angle > a){ ``` ``` n = (int)(angle/a); ``` ``` return(angle-n*a); ``` ``` } ``` ``` else{ ``` ``` return(angle); ``` ``` } ``` ```} ``` ``` ``` ```double angle360(angle) ``` ```double angle; ``` ```{ ``` ``` int n; ``` ``` ``` ``` if (angle < 0.0){ ``` ``` n = (int)(angle/360.0) - 1; ``` ``` return(angle-n*360.0); ``` ``` } ``` ``` else if (angle > 360.0){ ``` ``` n = (int)(angle/360.0); ``` ``` return(angle-n*360.0); ``` ``` } ``` ``` else{ ``` ``` return(angle); ``` ``` } ``` ```} ``` ``` ``` ```void ``` ```Radec_to_Cart(ra, dec, r) ``` ```double ra, dec; /* RA and DEC */ ``` ```Vector *r; /* returns corresponding cartesian unit vector */ ``` ```{ ``` ``` ``` ``` /* ``` ``` * Convert ra/dec from degrees to radians ``` ``` */ ``` ``` ra *= RadPerDeg; ``` ``` dec *= RadPerDeg; ``` ``` ``` ``` ``` ``` /* ``` ``` * Compute cartesian coordinates (in GEI) ``` ``` */ ``` ``` r->x = cos(dec) * cos(ra); ``` ``` r->y = cos(dec) * sin(ra); ``` ``` r->z = sin(dec); ``` ``` ``` ```} ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ``` ```int LeapYear(year) ``` ```int year; ``` ```{ ``` ``` if ((year%100 == 0)&&(year%400 != 0)) return(0); ``` ``` else if (year%4 == 0) return(1); ``` ``` else return(0); ``` ```} ``` ``` ``` ``` ``` ``` ``` ``` ```