/* $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 <assert.h>
#include <iostream.h>
#include <math.h>
#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;
}
syntax highlighted by Code2HTML, v. 0.9.1