/* $Id: astro.cc,v 1.6 2000/01/13 20:11:22 mac Exp $ */ /* * glbiff -- A Mesa/OpenGL-based `xbiff' substitute * Copyright (C) 2000 Maciej Kalisiak * * This program is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License * as published by the Free Software Foundation; either version 2 * of the License, or (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. * */ //////////////////////////////////////////////////////////// // // HUGE WARNING!!!!!! // // The following code has been conceived somewhere in the // whee hours of a morning, certainly past the 4am mark. // It is terribly ugly, and brain damage is evident. // // Read the following at your own risk. // // You've been warned. // // P.S. The "magic numbers" are not my fault. That's what // my source (the book mentioned below) had... // //-------------- // // the routines and formulae present here are based on: // Practical Astronomy with your Calculator // Third Edition // Peter Duffett-Smith // Cambridge University Press // ISBN 0 521 35699 7 (paperback) // #include #include #include #include "astro.h" const int month_length[12] = {31,29,31,30,31,30,31,31,30,31,30,31}; char* month_name[12] = { "Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec", }; char* weekday_name[7] = { "Sunday", "Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", }; ///////////// prototypes /////////////// bool is_DST( const Date & ); // this is the reference point in time for all calculations, // just like in the book (note: 0.0 jan, 1990 really refers to // the midnight between 30 and 31 dec, 1989 Date epoch(1990,JAN,0.0); ostream& operator<<( ostream & sout, Date & dt ) { sout << dt.day << " " << month_name[dt.month-1] << ", " << dt.year; return sout; } //////////////////////////////////////////////////////////// // class Date stuff bool Date::operator>( const Date & dt ) const { if( year > dt.year ) return true; if( year == dt.year && month > dt.month ) return true; if( year == dt.year && month == dt.month && day > dt.day ) return true; return false; } // bool operator>( const Date & dt1, const Date & dt2 ) { // if( dt1.year > dt2.year ) // return true; // if( dt1.year == dt2.year && dt1.month > dt2.month ) // return true; // if( dt1.year == dt2.year && dt1.month == dt2.month && dt1.day > dt2.day ) // return true; // return false; // } bool Date::operator<( const Date & dt ) const { if( year < dt.year ) return true; if( year == dt.year && month < dt.month ) return true; if( year == dt.year && month == dt.month && day < dt.day ) return true; return false; } //////////////////////////////////////////////////////////// // Date::valid() // // checks whether the Date is a valid one //////////////////////////////////////////////////////////// bool Date::valid() { if( month==0 || month > 12 ) return false; if( day<0 ) return false; if( day > (double)month_length[ month-1 ] + 1 ) return false; // perhaps we should be more specific here for February, // and actually figure out how many days it had on given // year, and use that number // if made it here, everything is ok return true; } double norm_angle( double d ) { while( d<0 ) d += 360; while( d>=360 ) d -= 360; return d; } double norm_hour( double h ) { while( h<0 ) h += 24; while( h>=24 ) h -= 24; return h; } //////////////////////////////////////////////////////////// // returns answer in radians double inv_tan( double y, double x ) { double res = atan( y/x ); if( x < 0 ) res += M_PI; return res; } //////////////////////////////////////////////////////////// // Date_to_JD() // // given a `Date', returns the Julian Day //////////////////////////////////////////////////////////// julian_day Date_to_JD( Date dt ) { if( !dt.valid() ) { cerr << "Date_to_JD(): invalid Date passed in (" << dt << ")" << endl; return INVALID_JD; } if( dt.month == 1 || dt.month == 2 ) { dt.month += 12; dt.year--; } int B, C, D; if( dt >= Date(1582, OCT, 15) ) B = 2 - dt.year/100 + dt.year/400; else B = 0; if( dt.year < 0 ) C = (int)(365.25 * dt.year - 0.75); else C = (int)(365.25 * dt.year); D = (int)(30.6001 * (dt.month+1)); return (julian_day)(B + C + D + dt.day + 1720994.5); } //////////////////////////////////////////////////////////// // result is in degrees double mean_anomaly( const Date & dt ) { double D = Date_to_JD( dt ) - Date_to_JD( epoch ); return norm_angle((360.0/SOLAR_D_PER_Y)*D + EPSILON_G - OMEGA_G); } //////////////////////////////////////////////////////////// // result is in degrees double true_anomaly( const Date & dt ) { double Ecc; // in text it is `E' double M = mean_anomaly( dt )/180 * M_PI; // first, Kepler's equation (using what looks like Netwon's Method) Ecc = M; const double max_error = 1e-12; double delta; for(;;) { delta = Ecc - E * sin(Ecc) - M; if( fabs(delta) <= max_error ) break; Ecc -= delta/(1 - E*cos(Ecc)); } double tmp = sqrt( (1+E)/(1-E) ) * tan(Ecc/2); double v = atan(tmp) * 2 * (180 / M_PI); return v; } //////////////////////////////////////////////////////////// // returns answer in degrees double obliquity( const Date & dt ) { double T = Date_to_JD( dt ) - Date_to_JD( Date(2000,JAN,1.5) ); T /= 36525.0; double dEps = 46.815*T + 0.0006*T*T - 0.00181*T*T*T; dEps /= 3600; return 23.439292 - dEps; } coord ecl_to_equ( const coord & ecl, const Date & dt ) { coord res; double epsilon = obliquity( dt )/180*M_PI; // declination res.b = asin( sin(ecl.b)*cos(epsilon) + cos(ecl.b)*sin(epsilon)*sin(ecl.a) ); // RA double y = sin(ecl.a)*cos(epsilon)-tan(ecl.b)*sin(epsilon); double x = cos(ecl.a); res.a = inv_tan( y, x ); return res; } coord sun_pos( const Date & dt ) { double true_ano = true_anomaly( dt ); double lambda = (OMEGA_G + true_ano)/180*M_PI; return ecl_to_equ( coord(lambda,0), dt ); } double rise_time( const coord & equ, const coord & geo ) { double H = ((1.0/15.0) * acos( -tan(geo.b)*tan(equ.b)))*(180.0/M_PI); double res=24 - H + equ.a*12/M_PI; if( res >= 24 ) res -= 24; return res; } double set_time( const coord & equ, const coord & geo ) { double H = ((1.0/15.0) * acos( -tan(geo.b)*tan(equ.b)))*(180.0/M_PI); double res=H + equ.a*12/M_PI; if( res >= 24 ) res -= 24; return res; } double lst_to_gst( double lst, const coord & geo ) { double res = lst + geo.a*12.0/M_PI; if( res < 0 ) res += 24; if( res >= 24 ) res -= 24; return res; } double ut_to_gst( double ut, const Date & dt ) { julian_day jd = Date_to_JD( dt ); julian_day S = jd - 2451545.0; double T = S/36525.0; double T0 = 6.697374558 + (2400.051336*T) + (0.000025862*T*T); T0 = norm_hour( T0 ); ut *= 1.002737909; return norm_hour( ut + T0 ); } double gst_to_ut( double gst, const Date & dt ) { julian_day jd = Date_to_JD( dt ); double S = jd - 2451545.0; double T = S/36525.0; double T0 = 6.697374558 + (2400.051336*T) + (0.000025862*T*T); T0 = norm_hour( T0 ); return norm_hour(gst - T0) * 0.9972695663; } double ut_to_civil( double ut ) { return norm_hour( ut + TimeZone ); } // vert_displacement is to be in rads // // returns dt in seconds double corr_for_rise_set( double vert_displacement, coord geo, coord equ ) { double dt; double phi = acos( sin(geo.b)/cos(equ.b) ); //rads double y = asin( sin(vert_displacement)/sin(phi) ); //rads dt = 240 * (y/M_PI*180) / cos(equ.b); return dt; } coord sun_rise_set( Date dt, coord geo ) { coord pos_1 = sun_pos( dt ); double lambda = (OMEGA_G + true_anomaly(dt))/180*M_PI; coord pos_2 = ecl_to_equ( coord(lambda+0.985647/180*M_PI,0), dt ); double lst1_r = rise_time( pos_1, geo ); double lst1_s = set_time( pos_1, geo ); double lst2_r = rise_time( pos_2, geo ); double lst2_s = set_time( pos_2, geo ); double gst1_r = lst_to_gst( lst1_r, geo ); double gst1_s = lst_to_gst( lst1_s, geo ); double gst2_r = lst_to_gst( lst2_r, geo ); double gst2_s = lst_to_gst( lst2_s, geo ); if( gst1_r > gst2_r ) gst2_r += 24; if( gst1_s > gst2_s ) gst2_s += 24; double T00 = ut_to_gst( 0, dt ); double T00_ = T00 - (-geo.a/M_PI*12) * 1.002738; if( T00_ < 0 ) T00_ += 24; if( gst1_r < T00_ ) { gst1_r += 24; gst2_r += 24; } if( gst1_s < T00_ ) { gst1_s += 24; gst2_s += 24; } double gst_r = (24.07*gst1_r - T00*(gst2_r-gst1_r))/(24.07+gst1_r-gst2_r); double gst_s = (24.07*gst1_s - T00*(gst2_s-gst1_s))/(24.07+gst1_s-gst2_s); ////// do some corrections now // for vertical displacement: using an arbitrary, relatively ok value for // refraction, and totally ignoring parallax (it's small) double x = (THETA_O/2)/180.0*M_PI + (34.0/60)/180.0*M_PI; coord equ( 0, (pos_1.b+pos_2.b)/2 ); // RA doesn't matter double delta_t = corr_for_rise_set( x, geo, equ ) / 3600.0; gst_r -= delta_t; gst_s += delta_t; if( is_DST( dt ) ) { gst_r += 1; gst_s += 1; } return coord(ut_to_civil(gst_to_ut(gst_r,dt)), ut_to_civil(gst_to_ut(gst_s,dt)) ); } // declination of the sun should be at midday (this is rough though) double twilight_dur( coord equ_sun, coord geo ) { double H=acos( -tan(geo.b)*tan(equ_sun.b) ); double H_=acos((cos(108.0/180.0*M_PI)-sin(geo.b)*sin(equ_sun.b))/(cos(geo.b)*cos(equ_sun.b)) ); double dt= (H_ - H)/(15.0/180.0*M_PI); // sidereal time dt *= 0.9973; // UT return dt; } double twilight_dur( Date dt, coord geo ) { coord pos_1 = sun_pos( dt ); double lambda = (OMEGA_G + true_anomaly(dt))/180*M_PI; coord pos_2 = ecl_to_equ( coord(lambda+0.985647/180*M_PI,0), dt ); coord pos_mid = coord( (pos_1.a+pos_2.a)/2, (pos_1.b+pos_2.b)/2 ); return twilight_dur( pos_mid, geo ); } //////////////////////////////////////////////////////////// // day_of_week() // // returns the day of the week // 0==Sunday, 1==Monday, etc... //////////////////////////////////////////////////////////// int day_of_week( const Date & dt ) { julian_day jd = Date_to_JD( dt ); double A = (jd+1.5)/7.0; A -= ((int)A); return (int)(A*7+0.5); } //////////////////////////////////////////////////////////// // is_DST() // // returns true if on a given Date the DST is in effect // // Note: DST // * begins at 2 a.m. on the first Sunday of April // * ends at 2 a.m. on the last Sunday of October //////////////////////////////////////////////////////////// bool is_DST( const Date & dt ) { Date start(dt.year,APR,1), end(dt.year,OCT,31); int tmp = day_of_week( start ); if( tmp ) start.day += 7-tmp; tmp = day_of_week( end ); if( tmp ) end.day -= tmp; if( dt < start || dt > end ) return false; return true; } time_hms time_convert( double decimal_hour ) { time_hms result; assert( decimal_hour >= 0 ); result.h = (unsigned)decimal_hour; result.m = (unsigned)((decimal_hour - result.h) * 60); result.s = (unsigned)rint((((decimal_hour - result.h)*60) - result.m)*60); return result; } ostream& operator<<( ostream& os, const time_hms& tm ) { os << tm.h << ":" << tm.m << ":" << tm.s; return os; } void astro_test(void) { Date today(1998,SEP,11); cout.precision(24); cout << "Julian day for 17.25 Feb, 1985 is " << Date_to_JD( Date(1985,FEB,17.25) ) << endl; cout << "mean anomaly for 27 July, 1980 is " << mean_anomaly( Date(1980,JUL,27) ) << endl; cout << "real anomaly for 27 July, 1988 is " << true_anomaly( Date(1988,JUL,27) ) << endl; cout << "obliquity on 1980.0 is " << obliquity( Date(1980,JAN,0) ) << endl; cout << "obliquity on 1950.0 is " << obliquity( Date(1950,JAN,0) ) << endl; coord res = ecl_to_equ( coord(139.686111/180*M_PI,4.875278/180*M_PI), Date(1980,JAN,0) ); cout << "long.=139o41'10\" and lat.=4o52'31\": RA is " << res.a*180/M_PI << ", decl is " << res.b*180/M_PI << endl; res = sun_pos( Date(1988,JUL,27) ); cout << "sun position on 27 Jul, 1988 is " << res.a*180/M_PI << " RA, " << res.b*180/M_PI << " decl." << endl; cout << "X-rise of ... " << rise_time( coord(23.655556/12.0*M_PI,21.7/180*M_PI), coord(64.0/180.0*M_PI,30/180.0*M_PI) ) << endl; cout << "X-set of ... " << set_time( coord(23.655556/12.0*M_PI,21.7/180*M_PI), coord(64.0/180.0*M_PI,30/180.0*M_PI) ) << endl; cout << "GST at 64oW if LST= 0h24m5.23s is " << lst_to_gst( 0.401453, coord( 64/180.0*M_PI, 0 ) ) << endl; cout << "when UT = 14.614353 on Apr 22, 1980, GST is " << ut_to_gst( 14.614353, Date(1980,APR,22) ) << endl; cout << "when gst = 4.668119 on Apr 22, 1980, ut is " << gst_to_ut( 4.668119, Date(1980,APR,22) ) << endl; cout << "Big Sunrise/Sunset calc from book" << endl; res = sun_rise_set( Date(1986,MAR,10), coord(71.05/180*M_PI,42.37/180*M_PI) ); cout << "sunrise is " << time_convert( res.a ) << endl; cout << "sunset is " << time_convert( res.b ) << endl; res = sun_rise_set( today, coord(79.5/180*M_PI, 43.75/180*M_PI) ); cout << "TODAY, HERE" << endl; cout << "sunrise is " << time_convert( res.a ) << endl; cout << "sunset is " << time_convert( res.b ) << endl; int day = day_of_week( Date(1998,MAY,10) ); cout << "Today is " << weekday_name[day] << endl; cout << "Today we are " << ((is_DST( Date(1998,MAY,10) ))?"":"not") << " under DST" << endl; Date test_date(1980,AUG,24); res = sun_rise_set( test_date, coord(30/180*M_PI,-64/180*M_PI) ); cout << "rise = " << time_convert( res.a ) << " set = " << time_convert( res.b ) << endl; test_date = Date(1979,SEP,7); cout << "S50: twilight duration: " << twilight_dur( test_date, coord(0,52.0/180.0*M_PI) ) << endl; }