/* Calculate time of transit * assuming RA and Dec change uniformly with time * * -- S. L. Moshier */ #include "kep.h" extern double glat; /* Earth radii per au */ #define DISFAC 2.3454780e4 /* cosine of 90 degrees 50 minutes: */ #define COSSUN -0.014543897651582657 /* cosine of 90 degrees 34 minutes: */ #define COSZEN -9.8900378587411476e-3 /* Returned transit, rise, and set times in radians (2 pi = 1 day) */ double r_trnsit; double r_rise; double r_set; int f_trnsit; int trnsit(J, lha, dec) double J; /* Julian date */ double lha; /* local hour angle, radians */ double dec; /* declination, radians */ { double x, y, z, N, D, TPI; double lhay, cosdec, sindec, coslat, sinlat; f_trnsit = 0; TPI = 2.0*PI; /* Initialize to no-event flag value. */ r_rise = -10.0; r_set = -10.0; /* observer's geodetic latitude, in radians */ x = glat * DTR; coslat = cos(x); sinlat = sin(x); cosdec = cos(dec); sindec = sin(dec); /* x = (J - floor(J-0.5) - 0.5) * TPI; */ x = floor(J) - 0.5; x = (J - x) * TPI; /* adjust local hour angle */ y = lha; /* printf ("%.7f,", lha); */ while( y < -PI ) y += TPI; while( y > PI ) y -= TPI; lhay = y; y = y/( -dradt/TPI + 1.00273790934); r_trnsit = x - y; /* printf ("rt %.7f ", r_trnsit); */ /* Ordinarily never print here. */ if( prtflg > 1 ) { printf( "approx local meridian transit" ); hms( r_trnsit ); printf( "UT (date + %.5f)\n", r_trnsit/TPI ); } if( (coslat == 0.0) || (cosdec == 0.0) ) goto norise; /* The time at which the upper limb of the body meets the * horizon depends on the body's angular diameter. */ switch( objnum ) { /* Sun */ case 0: N = COSSUN; break; /* Moon, elevation = -34' - semidiameter + parallax * semidiameter = 0.272453 * parallax + 0.0799" */ case 3: N = 1.0/(DISFAC*obpolar[2]); D = asin( N ); /* the parallax */ N = - 9.890199094634534e-3 + (1.0 - 0.2725076)*D - 3.874e-7; N = sin(N); break; /* Other object */ default: N = COSZEN; break; } y = (N - sinlat*sindec)/(coslat*cosdec); if( (y < 1.0) && (y > -1.0) ) { f_trnsit = 1; /* Derivative of y with respect to declination * times rate of change of declination: */ z = -ddecdt*(sinlat + COSZEN*sindec); z /= TPI*coslat*cosdec*cosdec; /* Derivative of acos(y): */ z /= sqrt( 1.0 - y*y); y = acos(y); D = -dradt/TPI + 1.00273790934; r_rise = x - (lhay + y)*(1.0 + z)/D; r_set = x - (lhay - y)*(1.0 - z)/D; /* Ordinarily never print here. */ if( prtflg > 1 ) { printf( "rises approx " ); hms(r_rise); printf( "UT, sets approx " ); hms(r_set); printf( "UT\n" ); } } norise: ; return(0); } extern struct orbit earth; /* Julian dates of rise, transet and set times. */ double t_rise; double t_trnsit; double t_set; /* Compute estimate of lunar rise and set times for iterative solution. */ static void iter_func(t, func) double t; int (* func)(); { int prtsave; prtsave = prtflg; prtflg = 0; JD = t; update(); /* Set UT and TDT */ kepler( TDT, &earth, rearth, eapolar ); (*func)(); prtflg = prtsave; } /* Iterative computation of rise, transit, and set times. */ int iter_trnsit( func ) int (* func)(); { double JDsave, TDTsave, UTsave; double date, date_save, date_trnsit, t0, t1; double rise1, set1, trnsit1, loopctr, retry; double TPI; int prtsave; loopctr = 0; prtsave = prtflg; TPI = PI + PI; JDsave = JD; TDTsave = TDT; UTsave = UT; date = floor(UT - 0.5) + 0.5; retry = 0; date_save = date; t1 = date_save; /* Find transit time. */ do { date = date_save; t0 = t1; iter_func(t0, func); t1 = date + r_trnsit / TPI; if (++loopctr > 10) { printf ("? Transit time did not converge.\n"); goto no_trnsit; } /* Reject transit if it is more than half a day from entered date. */ if (retry == 0) { if ((UTsave - t1) > 0.5) { date_save += 1.0; t1 = t_trnsit + 1.0; retry = 1; /* goto do_retry; */ } if ((UTsave - t1) < -0.5) { date_save -= 1.0; t1 = t_trnsit - 1.0; retry = 1; /* goto do_retry; */ } } } while (fabs(t1 - t0) > .0001); t_trnsit = t1; trnsit1 = r_trnsit; set1 = r_set; if (f_trnsit == 0) goto prtrnsit; /* Set current date to be that of the transit just found. */ date_trnsit = date; t1 = date + r_rise / TPI; /* Choose rising no later than transit. */ if (t1 >= t_trnsit) { date -= 1.0; t1 = date + r_rise / TPI; } loopctr = 0; do { t0 = t1; iter_func(t0, func); /* Skip out if no event found. Don't report rise or set. */ if ((f_trnsit == 0) || (++loopctr > 10)) { if (loopctr > 10) { printf ("? Rise time did not converge.\n"); } f_trnsit = 0; goto prtrnsit; } t1 = date + r_rise / TPI; if (t1 > t_trnsit) { date -= 1; t1 = date + r_rise / TPI; } } while (fabs(t1 - t0) > .0001); rise1 = r_rise; t_rise = t1; /* Set current date to be that of the transit. */ date = date_trnsit; r_set = set1; /* Choose setting no earlier than transit. */ t1 = date + r_set / TPI; if (t1 <= t_trnsit) { date += 1.0; t1 = date + r_set / TPI; } loopctr = 0; do { t0 = t1; iter_func(t0, func); if ((f_trnsit == 0) || (++loopctr > 10)) { if (loopctr > 10) { printf ("? Set time did not converge.\n"); } f_trnsit = 0; goto prtrnsit; } t1 = date + r_set / TPI; if (t1 < t_trnsit) { date += 1.0; t1 = date + r_set / TPI; } } while (fabs(t1 - t0) > .0001); t_set = t1; r_trnsit = trnsit1; r_rise = rise1; prtrnsit: prtflg = 1; printf( "local meridian transit " ); UT = t_trnsit; jtocal (t_trnsit); if (f_trnsit != 0) { printf( "rises " ); UT = t_rise; jtocal (t_rise); printf( "sets " ); UT = t_set; jtocal (t_set); t0 = t_set - t_rise; if ((t0 > 0.0) && (t0 < 1.0)) printf("Visible hours %.4f\n", 24.0 * t0); if ((fabs(JDsave - t_rise) > 0.5) && (fabs(JDsave - t_trnsit) > 0.5) && (fabs(JDsave - t_set) > 0.5)) printf("? wrong event date.\n"); } no_trnsit: JD = JDsave; TDT = TDTsave; UT = UTsave; /* Reset to original input date entry. */ prtflg = 0; update(); prtflg = prtsave; return 0; }