cdeagle Senior Member
Joined: 22/06/2014 Location: United StatesPosts: 261 |
Posted: 10:02am 05 Apr 2017 |
Copy link to clipboard |
Print this post |
|
This post describes a MicroMite eXtreme computer program that can be used to compute the circumstances of apogee and perigee of the Moon. Apogee corresponds to greatest distance from the Earth and perigee is the instant of closest approach to the Earth.
Here's a typical interaction with apmoon.bas and the computed results.
apogee and perigee of the moon ============================== please input the calendar date (month [1 - 12], day [1 - 31], year [yyyy]) < for example, october 21, 1986 is input as 10,21,1986 > < b.c. dates are negative, a.d. dates are positive > < the day of the month may also include a decimal part > ? 1,1,2017 please input the number of days to search ? 60 searching for apogee and perigee conditions ... time and conditions at perigee ============================== calendar date January 10 2017 UTC time 6 hours 0 minutes 53.82 seconds UTC julian day 2457763.75062294 geocentric distance 363238.6560 kilometers time and conditions at apogee ============================= calendar date January 22 2017 UTC time 0 hours 13 minutes 47.68 seconds UTC julian day 2457775.50957959 geocentric distance 404914.5482 kilometers time and conditions at perigee ============================== calendar date February 6 2017 UTC time 14 hours 2 minutes 8.54 seconds UTC julian day 2457791.08482106 geocentric distance 368816.0306 kilometers time and conditions at apogee ============================= calendar date February 18 2017 UTC time 21 hours 13 minutes 47.53 seconds UTC julian day 2457803.38457784 geocentric distance 404376.3116 kilometers time and conditions at perigee ============================== calendar date March 3 2017 UTC time 7 hours 33 minutes 5.70 seconds UTC julian day 2457815.81464930 geocentric distance 369062.0120 kilometers Here's the MMBASIC source code for this application.
' apmoon.bas April 5, 2017 ' apogee and perigee of the moon ' Micromite eXtreme version ''''''''''''''''''''''''''' option default float option base 1 ' dimension global arrays and variables dim jdleap(28), leapsec(28), xnut(11, 13) dim month$(12) as string dim trr, iap_flg%, jdtdbi dim cmonth, cday, cyear, ndays ' global constants const pi2 = 2.0 * pi, pidiv2 = 0.5 * pi, dtr = pi / 180.0, rtd = 180.0 / pi const atr = pi / 648000.0, seccon = 206264.8062470964 ' read subset of IAU2000 nutation data for j% = 1 to 13 for i% = 1 to 11 read xnut(i%, j%) next i% next j% data 0, 0, 0, 0, 1, -172064161, -174666, 92052331, 9086, 33386, 15377 data 0, 0, 2, -2, 2, -13170906, -1675, 5730336, -3015, -13696, -4587 data 0, 0, 2, 0, 2, -2276413, -234, 978459, -485, 2796, 1374 data 0, 0, 0, 0, 2, 2074554, 207, -897492, 470, -698, -291 data 0, 1, 0, 0, 0, 1475877, -3633, 73871, -184, 11817, -1924 data 0, 1, 2, -2, 2, -516821, 1226, 224386, -677, -524, -174 data 1, 0, 0, 0, 0, 711159, 73, -6750, 0, -872, 358 data 0, 0, 2, 0, 1, -387298, -367, 200728, 18, 380, 318 data 1, 0, 2, 0, 2, -301461, -36, 129025, -63, 816, 367 data 0, -1, 2, -2, 2, 215829, -494, -95929, 299, 111, 132 data 0, 0, 2, -2, 1, 128227, 137, -68982, -9, 181, 39 data -1, 0, 2, 0, 2, 123457, 11, -53311, 32, 19, -4 data -1, 0, 0, 2, 0, 156994, 10, -1235, 0, -168, 82 ' read leap second data for i% = 1 to 28 read jdleap(i%), leapsec(i%) next i% data 2441317.5, 10.0 data 2441499.5, 11.0 data 2441683.5, 12.0 data 2442048.5, 13.0 data 2442413.5, 14.0 data 2442778.5, 15.0 data 2443144.5, 16.0 data 2443509.5, 17.0 data 2443874.5, 18.0 data 2444239.5, 19.0 data 2444786.5, 20.0 data 2445151.5, 21.0 data 2445516.5, 22.0 data 2446247.5, 23.0 data 2447161.5, 24.0 data 2447892.5, 25.0 data 2448257.5, 26.0 data 2448804.5, 27.0 data 2449169.5, 28.0 data 2449534.5, 29.0 data 2450083.5, 30.0 data 2450630.5, 31.0 data 2451179.5, 32.0 data 2453736.5, 33.0 data 2454832.5, 34.0 DATA 2456109.5, 35.0 data 2457204.5, 36.0 data 2457754.5, 37.0 ' calendar months month$(1) = "January" month$(2) = "February" month$(3) = "March" month$(4) = "April" month$(5) = "May" month$(6) = "June" month$(7) = "July" month$(8) = "August" month$(9) = "September" month$(10) = "October" month$(11) = "November" month$(12) = "December" '''''''''''''''''' ' begin simulation '''''''''''''''''' print " " print "apogee and perigee of the moon" print "==============================" ' request initial calendar date (month, day, year) getdate(cmonth, cday, cyear) ' initial utc julian day julian(cmonth, cday, cyear, jdutc) ' compute initial tdb julian date utc2tdb(jdutc, jdtdb) jdtdbi = jdtdb ' request search duration (days) print " " print "please input the number of days to search" input ndays print " " print "searching for apogee and perigee conditions ..." print " " ' initialize apogee/perigee "flag" iap_flg% = 1 ' define search parameters ti = 0.0 tf = ndays tisaved = ti dt = 1.0 dtsml = 0.1 ' find moon rise/set conditions ap_event(ti, tf, dt, dtsml) end '''''''''''''''''''''''''''''' '''''''''''''''''''''''''''''' sub ap_event (ti, tf, dt, dtsml) ' predict apogee/perigee events ' input ' ti = initial simulation time ' tf = final simulation time ' dt = step size used for bounding minima ' dtsml = small step size used to determine whether ' the function is increasing or decreasing ''''''''''''''''''''''''''''''''''''''''''''''''''' LOCAL tolm local fmin1, tmin1 LOCAL ftemp, df, dflft local el, er LOCAL t, ft local iter1%, iter2%, iter3% ' initialization tolm = 1.0e-10 df = 1.0 for iter1% = 1 to 1000 ' find where function first starts decreasing for iter2% = 1 to 1000 if (df <= 0.0) then exit for end if t = t + dt ap_func(t, ft) ap_func(t- dtsml, ftemp) df = ft - ftemp next iter2% ' function decreasing - find where function ' first starts increasing for iter3% = 1 to 1000 el = t dflft = df t = t + dt ap_func(t, ft) ap_func(t - dtsml, ftemp) df = ft - ftemp if (df > 0.0) then exit for next iter3% er = t ' calculate minimum using Brent's method minima(el, er, tolm, tmin1, fmin1) el = er ' print current conditions ap_print(tmin1) if (t >= tf) then exit for next iter1% end sub '''''''''''''''' '''''''''''''''' sub ap_func(x, fx) ' apogee/perigee objective function ''''''''''''''''''''''''''''''''''' local jdtdb, rmoon(3), srange ' current TDB julian day jdtdb = jdtdbi + x ' compute eci position vector and distance of the moon moon(jdtdb, rmoon()) srange = iap_flg% * vecmag(rmoon()) ' current objective function fx = srange end sub '''''''''''''''' '''''''''''''''' sub ap_print(topt) ' print apogee and perigee conditions ''''''''''''''''''''''''''''''''''''' LOCAL jdutc, rmoon(3) if (iap_flg% = 1) then print " " print "time and conditions at perigee" print "==============================" print " " else print " " print "time and conditions at apogee" print "=============================" print " " end if ' TDB julian day jdtdb = jdtdbi + topt ' compute and display UTC julian date tdb2utc(jdtdb, jdutc) jd2str(jdutc) PRINT " " print "UTC julian day ", str$(jdutc, 0, 8) PRINT " " moon(jdtdb, rmoon()) print "geocentric distance ", str$(vecmag(rmoon()), 0, 4), " kilometers" print " " ' toggle type of event flag iap_flg% = -iap_flg% END sub '''''''''''''''''''''''''''''''' '''''''''''''''''''''''''''''''' sub minima(a, b, tolm, xmin, fmin) ' one-dimensional minimization ' Brent's method ' input ' a = initial x search value ' b = final x search value ' tolm = convergence criterion ' output ' xmin = minimum x value ' note ' user-defined objective subroutine ' coded as ap_func(x, fx) ' remember: a maximum is simply a minimum ' with a negative attitude! ''''''''''''''''''''''''''''''''''''' ' machine epsilon LOCAL epsm = 2.23e-16 ' golden number LOCAL c = 0.38196601125 LOCAL iter as integer, d, e LOCAL t2, p, q local r, u, fu LOCAL x, xm, w local v, fx, fw LOCAL tol1, fv x = a + c * (b - a) w = x v = w e = 0.0 p = 0.0 q = 0.0 r = 0.0 ap_func(x, fx) fw = fx fv = fw for iter = 1 to 100 if (iter > 50) then print ("error in function minima!") print ("(more than 50 iterations)") end if xm = 0.5 * (a + b) tol1 = tolm * abs(x) + epsm t2 = 2.0 * tol1 if (abs(x - xm) <= (t2 - 0.5 * (b - a))) then xmin = x fmin = fx exit sub end if if (abs(e) > tol1) then r = (x - w) * (fx - fv) q = (x - v) * (fx - fw) p = (x - v) * q - (x - w) * r q = 2.0 * (q - r) if (q > 0.0) then p = -p end if q = abs(q) r = e e = d end if if ((abs(p) >= abs(0.5 * q * r)) or (p <= q * (a - x)) or (p >= q * (b - x))) then if (x >= xm) then e = a - x else e = b - x end if d = c * e else d = p / q u = x + d if ((u - a) < t2) or ((b - u) < t2) then d = sgn(xm - x) * tol1 end if end if if (abs(d) >= tol1) then u = x + d else u = x + sgn(d) * tol1 end if ap_func(u, fu) if (fu <= fx) then if (u >= x) then a = x else b = x end if v = w fv = fw w = x fw = fx x = u fx = fu else if (u < x) then a = u else b = u end if if ((fu <= fw) Or (w = x)) then v = w fv = fw w = u fw = fu elseif ((fu <= fv) Or (v = x) Or (v = w)) then v = u fv = fu end if end if next iter end sub '''''''''''''''''''''''' '''''''''''''''''''''''' sub tdb2utc (jdtdb, jdutc) ' convert TDB julian day to UTC julian day subroutine ' input ' jdtdb = TDB julian day ' output ' jdutc = UTC julian day ''''''''''''''''''''''''' local x1, x2, xroot, froot jdsaved = jdtdb ' set lower and upper bounds x1 = jdsaved - 0.1 x2 = jdsaved + 0.1 ' solve for UTC julian day using Brent's method realroot(x1, x2, 1.0e-8, xroot, froot) jdutc = xroot end sub ''''''''''''''''''' ''''''''''''''''''' sub jdfunc (jdin, fx) ' objective function for tdb2utc ' input ' jdin = current value for UTC julian day ' output ' fx = delta julian day '''''''''''''''''''''''' local jdwrk utc2tdb(jdin, jdwrk) fx = jdwrk - jdsaved end sub '''''''''''''''''''''' '''''''''''''''''''''' sub moon(jdate, rmoon()) ' geocentric position of the moon subroutine ' input ' jdate = tdb julian day ' output ' rmoon = eci position vector of the moon (kilometers) ''''''''''''''''''''''''''''''''''''''''''''''''''''''' local t1, t2, t3 local t4, dpsi, deps LOCAL ll, d, lp local l, f, t LOCAL ve, ma, ju local rm, dv, pl local plat, plon LOCAL a, b local rasc, decl, obliq ' get nutations and obliquity obliq_lp(jdate, dpsi, deps, obliq) ' evaluate lunar ephemeris t1 = (jdate - 2451545.0) / 36525.0 t2 = t1 * t1 t3 = t1 * t1 * t1 t4 = t1 * t1 * t1 * t1 ' calculate fundamental arguments (radians) ll = dtr * (218 + (18 * 60 + 59.95571) / 3600) ll = modulo(ll + atr * (1732564372.83264 * t1 - 4.7763 * t2 + .006681 * t3 - 0.00005522 * t4)) d = dtr * (297 + (51 * 60 + .73512) / 3600) d = modulo(d + atr * (1602961601.4603 * t1 - 5.8681 * t2 + .006595 * t3 - 0.00003184 * t4)) lp = dtr * (357 + (31 * 60 + 44.79306) / 3600) lp = modulo(lp + atr * (129596581.0474 * t1 - .5529 * t2 + 0.000147 * t3)) l = dtr * (134 + (57 * 60 + 48.28096) / 3600) l = modulo(l + atr * (1717915923.4728 * t1 + 32.3893 * t2 + .051651 * t3 - 0.0002447 * t4)) f = dtr * (93 + (16 * 60 + 19.55755) / 3600) f = modulo(f + atr * (1739527263.0983 * t1 - 12.2505 * t2 - .001021 * t3 + 0.00000417 * t4)) t = dtr * (100 + (27 * 60 + 59.22059) / 3600) t = modulo(t + atr * (129597742.2758 * t1 - .0202 * t2 + .000009 * t3 + 0.00000015 * t4)) ve = dtr * (181 + (58 * 60 + 47.28305) / 3600) ve = modulo(ve + atr * 210664136.43355 * t1) ma = dtr * (355 + (25 * 60 + 59.78866) / 3600) ma = modulo(ma + atr * 68905077.59284 * t1) ju = dtr * (34 + (21 * 60 + 5.34212) / 3600) ju = modulo(ju + atr * 10925660.42861 * t1) ' compute geocentric distance (kilometers) ' a(c,0,r) series rm = 385000.52899 rm = rm - 20905.35504 * sin(l + pidiv2) rm = rm - 3699.11092 * sin(2 * d - l + pidiv2) rm = rm - 2955.96756 * sin(2 * d + pidiv2) rm = rm - 569.92512 * sin(2 * l + pidiv2) rm = rm + 246.15848 * sin(2 * d - 2 * l + pidiv2) rm = rm - 204.58598 * sin(2 * d - lp + pidiv2) rm = rm - 170.73308 * sin(2 * d + l + pidiv2) rm = rm - 152.13771 * sin(2 * d - lp - l + pidiv2) rm = rm - 129.62014 * sin(lp - l + pidiv2) rm = rm + 108.7427 * sin(d + pidiv2) rm = rm + 104.75523 * sin(lp + l + pidiv2) rm = rm + 79.66056 * sin(l - 2 * f + pidiv2) rm = rm + 48.8883 * sin(lp + pidiv2) rm = rm - 34.78252 * sin(4 * d - l + pidiv2) rm = rm + 30.82384 * sin(2 * d + lp + pidiv2) rm = rm + 24.20848 * sin(2 * d + lp - l + pidiv2) rm = rm - 23.21043 * sin(3 * l + pidiv2) rm = rm - 21.63634 * sin(4 * d - 2 * l + pidiv2) rm = rm - 16.67471 * sin(d + lp + pidiv2) rm = rm + 14.40269 * sin(2 * d - 3 * l + pidiv2) rm = rm - 12.8314 * sin(2 * d - lp + l + pidiv2) rm = rm - 11.64995 * sin(4 * d + pidiv2) rm = rm - 10.44476 * sin(2 * d + 2 * l + pidiv2) rm = rm + 10.32111 * sin(2 * d - 2 * f + pidiv2) rm = rm + 10.0562 * sin(2 * d - lp - 2 * l + pidiv2) rm = rm - 9.88445 * sin(2 * d - 2 * lp + pidiv2) rm = rm + 8.75156 * sin(2 * d - l - 2 * f + pidiv2) rm = rm - 8.37911 * sin(d - l + pidiv2) rm = rm - 7.00269 * sin(lp - 2 * l + pidiv2) rm = rm + 6.322 * sin(d + l + pidiv2) rm = rm + 5.75085 * sin(lp + 2 * l + pidiv2) rm = rm - 4.95013 * sin(2 * d - 2 * lp - l + pidiv2) rm = rm - 4.42118 * sin(2 * l - 2 * f + pidiv2) rm = rm + 4.13111 * sin(2 * d + l - 2 * f + pidiv2) rm = rm - 3.95798 * sin(4 * d - lp - l + pidiv2) rm = rm + 3.25824 * sin(3 * d - l + pidiv2) rm = rm - 3.1483 * sin(2 * f + pidiv2) rm = rm + 2.61641 * sin(2 * d + lp + l + pidiv2) rm = rm + 2.35363 * sin(2 * d + 2 * lp - l + pidiv2) rm = rm - 2.11713 * sin(2 * lp - l + pidiv2) rm = rm - 1.89704 * sin(4 * d - lp - 2 * l + pidiv2) rm = rm - 1.73853 * sin(d - 2 * l + pidiv2) rm = rm - 1.57139 * sin(4 * d - lp + pidiv2) rm = rm - 1.42255 * sin(4 * d + l + pidiv2) rm = rm - 1.41893 * sin(3 * d + pidiv2) rm = rm + 1.16553 * sin(2 * lp + l + pidiv2) rm = rm - 1.11694 * sin(4 * l + pidiv2) rm = rm + 1.06567 * sin(2 * lp + pidiv2) rm = rm - .93332 * sin(d + lp + l + pidiv2) rm = rm + .86243 * sin(3 * d - 2 * l + pidiv2) rm = rm + .85124 * sin(d + lp - l + pidiv2) rm = rm - .8488 * sin(2 * d - lp + 2 * l + pidiv2) rm = rm - .79563 * sin(d - 2 * f + pidiv2) rm = rm + .77854 * sin(2 * d - 4 * l + pidiv2) rm = rm + .77404 * sin(2 * d - 2 * l + 2 * f + pidiv2) rm = rm - .66968 * sin(2 * d + 3 * l + pidiv2) rm = rm - .65753 * sin(2 * d - 2 * lp + l + pidiv2) rm = rm + .65706 * sin(2 * d - lp - 2 * f + pidiv2) rm = rm + .59632 * sin(2 * d - l + 2 * f + pidiv2) rm = rm + .57879 * sin(4 * d + lp - l + pidiv2) rm = rm - .51423 * sin(4 * d - 3 * l + pidiv2) rm = rm - .50792 * sin(4 * d - 2 * f + pidiv2) rm = rm + .49755 * sin(d - lp + pidiv2) rm = rm + .49504 * sin(2 * d - lp - 3 * l + pidiv2) rm = rm + .47262 * sin(2 * d - 2 * l - 2 * f + pidiv2) rm = rm - .4225 * sin(6 * d - 2 * l + pidiv2) rm = rm - .42241 * sin(lp - 3 * l + pidiv2) rm = rm - .41071 * sin(2 * d - 3 * lp + pidiv2) rm = rm + .37852 * sin(d + 2 * l + pidiv2) rm = rm + .35508 * sin(lp + 3 * l + pidiv2) rm = rm + .34302 * sin(2 * d - 2 * lp - 2 * l + pidiv2) rm = rm + .33463 * sin(lp - l + 2 * f + pidiv2) rm = rm + .33225 * sin(d + lp - 2 * l + pidiv2) rm = rm + .32334 * sin(2 * d - lp - l - 2 * f + pidiv2) rm = rm - .32176 * sin(4 * d - l - 2 * f + pidiv2) rm = rm - .28663 * sin(6 * d - l + pidiv2) rm = rm + .28399 * sin(2 * d + 2 * l - 2 * f + pidiv2) rm = rm - .27904 * sin(4 * d - 2 * lp - l + pidiv2) rm = rm + .2556 * sin(3 * d - lp - l + pidiv2) rm = rm - .2481 * sin(lp + l - 2 * f + pidiv2) rm = rm + .24452 * sin(4 * d + lp + pidiv2) rm = rm + .23695 * sin(4 * d + lp - 2 * l + pidiv2) rm = rm - .21258 * sin(3 * d + lp - l + pidiv2) rm = rm + .21251 * sin(2 * d + lp + 2 * l + pidiv2) rm = rm + .20941 * sin(2 * d - lp + l - 2 * f + pidiv2) rm = rm - .20285 * sin(4 * d - lp + l + pidiv2) rm = rm + .20099 * sin(3 * d - 2 * f + pidiv2) rm = rm - .18567 * sin(lp - 2 * f + pidiv2) rm = rm - .18316 * sin(6 * d - 3 * l + pidiv2) rm = rm + .16857 * sin(2 * d + lp - 3 * l + pidiv2) rm = rm - .15802 * sin(lp + 2 * f + pidiv2) rm = rm - .15707 * sin(3 * d - lp + pidiv2) rm = rm - .14806 * sin(2 * d - 3 * lp - l + pidiv2) rm = rm + .14763 * sin(2 * d + 2 * lp + pidiv2) rm = rm + .14368 * sin(2 * d + lp - 2 * l + pidiv2) rm = rm - .13922 * sin(4 * d + 2 * l + pidiv2) rm = rm - .13617 * sin(2 * lp - 2 * l + pidiv2) rm = rm - .13571 * sin(2 * d + lp - 2 * f + pidiv2) rm = rm - .12805 * sin(4 * d - 2 * lp + pidiv2) rm = rm + .11411 * sin(d - lp - l + pidiv2) rm = rm + .10998 * sin(d - lp + l + pidiv2) rm = rm - .10887 * sin(2 * d + 2 * lp - 2 * l + pidiv2) rm = rm - .10833 * sin(4 * d - 2 * lp - 2 * l + pidiv2) rm = rm - .10766 * sin(3 * d + lp + pidiv2) rm = rm - .10326 * sin(l + 2 * f + pidiv2) rm = rm - .09938 * sin(d - 3 * l + pidiv2) rm = rm - .08587 * sin(6 * d + pidiv2) rm = rm - .07982 * sin(4 * d - 2 * l - 2 * f + pidiv2) rm = rm - 6.678e-02 * sin(6 * d - lp - 2 * l + pidiv2) rm = rm - 6.545e-02 * sin(3 * d + l + pidiv2) rm = rm + .06055 * sin(d + l - 2 * f + pidiv2) rm = rm - .05904 * sin(d + lp + 2 * l + pidiv2) rm = rm - .05888 * sin(5 * l + pidiv2) rm = rm - .0585 * sin(2 * d - lp + 3 * l + pidiv2) rm = rm - .05789 * sin(4 * d - lp - 2 * f + pidiv2) rm = rm - .05527 * sin(2 * d + lp + l - 2 * f + pidiv2) rm = rm + .05293 * sin(3 * d - lp - 2 * l + pidiv2) rm = rm - .05191 * sin(6 * d - lp - l + pidiv2) rm = rm + .05072 * sin(2 * lp + 2 * l + pidiv2) rm = rm - .0502 * sin(lp - 2 * l + 2 * f + pidiv2) rm = rm - .04843 * sin(3 * d - 3 * l + pidiv2) rm = rm + .0474 * sin(2 * d - 5 * l + pidiv2) rm = rm - .04736 * sin(2 * d + lp - l - 2 * f + pidiv2) rm = rm - .04608 * sin(2 * d - 2 * lp + 2 * l + pidiv2) rm = rm + .04591 * sin(5 * d - 2 * l + pidiv2) rm = rm - .04422 * sin(2 * d + 4 * l + pidiv2) rm = rm - .04316 * sin(4 * d - lp - 3 * l + pidiv2) rm = rm - .04232 * sin(d - l - 2 * f + pidiv2) rm = rm - .03894 * sin(3 * lp - l + pidiv2) rm = rm + .0381 * sin(3 * d + lp - 2 * l + pidiv2) rm = rm + .03734 * sin(2 * d - lp - l + 2 * f + pidiv2) rm = rm + .03729 * sin(d + 2 * lp + pidiv2) rm = rm + .03682 * sin(4 * d + lp + l + pidiv2) rm = rm + .03379 * sin(d + lp - 2 * f + pidiv2) rm = rm + .03265 * sin(lp + 2 * l - 2 * f + pidiv2) rm = rm + .03143 * sin(2 * d + 2 * f + pidiv2) rm = rm + .03024 * sin(2 * d - lp - 2 * l + 2 * f + pidiv2) rm = rm - .02948 * sin(d - 2 * lp + pidiv2) rm = rm - .02939 * sin(4 * d - 4 * l + pidiv2) rm = rm + .0291 * sin(2 * d - 3 * l - 2 * f + pidiv2) rm = rm - .02855 * sin(2 * d - 3 * lp + l + pidiv2) rm = rm + .02839 * sin(2 * d - 2 * lp - 2 * f + pidiv2) rm = rm - .02698 * sin(4 * d - lp - l - 2 * f + pidiv2) rm = rm - .02674 * sin(lp - 4 * l + pidiv2) rm = rm + .02658 * sin(4 * d + 2 * lp - 2 * l + pidiv2) rm = rm - .02471 * sin(d - l + 2 * f + pidiv2) rm = rm - .02436 * sin(6 * d - lp - 3 * l + pidiv2) rm = rm - .02399 * sin(4 * d + lp - 3 * l + pidiv2) rm = rm + .02368 * sin(d + 3 * l + pidiv2) rm = rm + .02334 * sin(2 * d - lp - 4 * l + pidiv2) rm = rm + .02304 * sin(lp + 4 * l + pidiv2) rm = rm + .02127 * sin(3 * lp + pidiv2) rm = rm - .02079 * sin(4 * d - lp + 2 * l + pidiv2) rm = rm - .02008 * sin(2 * d - 3 * l + 2 * f + pidiv2) ' a(p,0,r) series rm = rm + 1.0587 * sin(2 * t - 2 * ju + 2 * d - l + 90.11969000000001 * dtr) rm = rm + .72783 * sin(18 * ve - 16 * t - 2 * l + 116.54311 * dtr) rm = rm + .68256 * sin(18 * ve - 16 * t + 296.54574 * dtr) rm = rm + .59827 * sin(3 * ve - 3 * t + 2 * d - l + 89.98187 * dtr) rm = rm + .45648 * sin(ll + l - f + 270.00126 * dtr) rm = rm + .45276 * sin(ll - l - f + 90.00128 * dtr) rm = rm + .41011 * sin(2 * t - 3 * ju + 2 * d - l + 280.06924 * dtr) rm = rm + .20497 * sin(t - ju - 2 * d + 91.79862 * dtr) rm = rm + .20473 * sin(18 * ve - 16 * t - 2 * d - l + 116.54222 * dtr) rm = rm + .20367 * sin(18 * ve - 16 * t + 2 * d - l + 296.54299 * dtr) rm = rm + .16644 * sin(2 * ve - 2 * t - 2 * d + 90.36386 * dtr) rm = rm + .1578 * sin(4 * t - 8 * ma + 3 * ju + l + 194.98833 * dtr) rm = rm + .1578 * sin(4 * t - 8 * ma + 3 * ju - l + 14.98841 * dtr) rm = rm + .15751 * sin(t - ju - 2 * d + l + 91.74578 * dtr) rm = rm + .1445 * sin(2 * t - 2 * ju - l + 89.97863 * dtr) rm = rm + .13811 * sin(ve - t - l + 270.00993 * dtr) rm = rm + .13477 * sin(18 * ve - 16 * t - 2 * d + 116.53978 * dtr) rm = rm + .12671 * sin(18 * ve - 16 * t + 2 * d - 2 * l + 296.54238 * dtr) rm = rm + .12666 * sin(t - ju - l + 91.22751 * dtr) rm = rm + .12362 * sin(ve - t - 2 * d + 269.98523 * dtr) rm = rm + .12047 * sin(2 * ve - 2 * t + 2 * d - l + 269.99692 * dtr) rm = rm + .11998 * sin(ve - t + l + 90.01606 * dtr) rm = rm + .11617 * sin(2 * ve - 2 * t - 2 * d + l + 90.31081 * dtr) rm = rm + .11256 * sin(4 * t - 8 * ma + 3 * ju + 2 * d - l + 197.11421 * dtr) rm = rm + .11251 * sin(4 * t - 8 * ma + 3 * ju - 2 * d + l + 17.11263 * dtr) rm = rm + .11226 * sin(4 * t - 8 * ma + 3 * ju + 2 * d + 196.69224 * dtr) rm = rm + .11216 * sin(4 * t - 8 * ma + 3 * ju - 2 * d + 16.68897 * dtr) rm = rm + .10689 * sin(ll + 2 * d - f + 270.00092 * dtr) rm = rm + .10504 * sin(t - ju + l + 271.06726 * dtr) rm = rm + .1006 * sin(ve - t - 2 * d + l + 269.98452 * dtr) rm = rm + 9.932e-02 * sin(3 * ve - 3 * t - l + 90.1054 * dtr) rm = rm + .09554 * sin(ll - 2 * d - f + 90.00096 * dtr) rm = rm + .08508 * sin(ll - f + 270.00061 * dtr) rm = rm + 7.9450e-02 * sin(4 * ve - 4 * t + 2 * d - l + 89.99224 * dtr) rm = rm + .07725 * sin(2 * t - 3 * ju - l + 280.16516 * dtr) rm = rm + 7.054e-02 * sin(6 * ve - 8 * t + 2 * d - l + 77.22087 * dtr) rm = rm + 6.313e-02 * sin(2 * ju - 5 * l + l + 256.56163 * dtr) ' a(p,1,r) series rm = rm + .51395 * t1 * sin(2 * d - lp + pidiv2) rm = rm + .38245 * t1 * sin(2 * d - lp - l + pidiv2) rm = rm + .32654 * t1 * sin(lp - l + pidiv2) rm = rm + .26396 * t1 * sin(lp + l + 270 * dtr) rm = rm + .12302 * t1 * sin(lp + 270 * dtr) rm = rm + .07754 * t1 * sin(2 * d + lp + 270 * dtr) rm = rm + .06068 * t1 * sin(2 * d + lp - l + 270 * dtr) rm = rm + .0497 * t1 * sin(2 * d - 2 * lp + pidiv2) rm = rm + .04194 * t1 * sin(d + lp + pidiv2) rm = rm + .03222 * t1 * sin(2 * d - lp + l + pidiv2) rm = rm + .02529 * t1 * sin(2 * d - lp - 2 * l + 270 * dtr) rm = rm + .0249 * t1 * sin(2 * d - 2 * lp - l + pidiv2) rm = rm + .00149 * t2 * sin(2 * d - lp + pidiv2) rm = rm + .00111 * t2 * sin(2 * d - lp - l + pidiv2) rmm = rm ' compute geocentric ecliptic longitude (arc seconds) ' a(c,0,v) series dv = 22639.58578 * sin(l) dv = dv + 4586.4383 * sin(2 * d - l) dv = dv + 2369.91394 * sin(2 * d) dv = dv + 769.02571 * sin(2 * l) dv = dv - 666.4171 * sin(lp) dv = dv - 411.59567 * sin(2 * f) dv = dv + 211.65555 * sin(2 * d - 2 * l) dv = dv + 205.43582 * sin(2 * d - lp - l) dv = dv + 191.9562 * sin(2 * d + l) dv = dv + 164.72851 * sin(2 * d - lp) dv = dv - 147.32129 * sin(lp - l) dv = dv - 124.98812 * sin(d) dv = dv - 109.38029 * sin(lp + l) dv = dv + 55.17705 * sin(2 * d - 2 * f) dv = dv - 45.0996 * sin(l + 2 * f) dv = dv + 39.53329 * sin(l - 2 * f) dv = dv + 38.42983 * sin(4 * d - l) dv = dv + 36.12381 * sin(3 * l) dv = dv + 30.77257 * sin(4 * d - 2 * l) dv = dv - 28.39708 * sin(2 * d + lp - l) dv = dv - 24.35821 * sin(2 * d + lp) dv = dv - 18.58471 * sin(d - l) dv = dv + 17.95446 * sin(d + lp) dv = dv + 14.53027 * sin(2 * d - lp + l) dv = dv + 14.3797 * sin(2 * d + 2 * l) dv = dv + 13.89906 * sin(4 * d) dv = dv + 13.19406 * sin(2 * d - 3 * l) dv = dv - 9.67905 * sin(lp - 2 * l) dv = dv - 9.36586 * sin(2 * d - l + 2 * f) dv = dv + 8.60553 * sin(2 * d - lp - 2 * l) dv = dv - 8.45310 * sin(d + l) dv = dv + 8.05016 * sin(2 * d - 2 * lp) dv = dv - 7.63015 * sin(lp + 2 * l) dv = dv - 7.44749 * sin(2 * lp) dv = dv + 7.37119 * sin(2 * d - 2 * lp - l) dv = dv - 6.38315 * sin(2 * d + l - 2 * f) dv = dv - 5.74161 * sin(2 * d + 2 * f) dv = dv + 4.37401 * sin(4 * d - lp - l) dv = dv - 3.99761 * sin(2 * l + 2 * f) dv = dv - 3.20969 * sin(3 * d - l) dv = dv - 2.91454 * sin(2 * d + lp + l) dv = dv + 2.73189 * sin(4 * d - lp - 2 * l) dv = dv - 2.56794 * sin(2 * lp - l) dv = dv - 2.5212 * sin(2 * d + 2 * lp - l) dv = dv + 2.48889 * sin(2 * d + lp - 2 * l) dv = dv + 2.14607 * sin(2 * d - lp - 2 * f) dv = dv + 1.97773 * sin(4 * d + l) dv = dv + 1.93368 * sin(4 * l) dv = dv + 1.87076 * sin(4 * d - lp) dv = dv - 1.75297 * sin(d - 2 * l) dv = dv - 1.43716 * sin(2 * d + lp - 2 * f) dv = dv - 1.37257 * sin(2 * l - 2 * f) dv = dv + 1.26182 * sin(d + lp + l) dv = dv - 1.22412 * sin(3 * d - 2 * l) dv = dv + 1.18683 * sin(4 * d - 3 * l) dv = dv + 1.177 * sin(2 * d - lp + 2 * l) dv = dv - 1.16169 * sin(2 * lp + l) dv = dv + 1.07769 * sin(d + lp - l) dv = dv + 1.0595 * sin(2 * d + 3 * l) dv = dv - .99022 * sin(2 * d + l + 2 * f) dv = dv + .94828 * sin(2 * d - 4 * l) dv = dv + .75168 * sin(2 * d - 2 * lp + l) dv = dv - .66938 * sin(lp - 3 * l) dv = dv - .63521 * sin(4 * d + lp - l) dv = dv - .58399 * sin(d + 2 * l) dv = dv - .58331 * sin(d - 2 * f) dv = dv + .57156 * sin(6 * d - 2 * l) dv = dv - .56064 * sin(2 * d - 2 * l - 2 * f) dv = dv - .55692 * sin(d - lp) dv = dv - .54592 * sin(lp + 3 * l) dv = dv - .53571 * sin(2 * d - 2 * l + 2 * f) dv = dv + .4784 * sin(2 * d - lp - 3 * f) dv = dv - .45379 * sin(2 * d + 2 * l - 2 * f) dv = dv - .42622 * sin(2 * d - lp - l + 2 * f) dv = dv + .42033 * sin(4 * f) dv = dv + .4134 * sin(lp + 2 * f) dv = dv + .40423 * sin(3 * d) dv = dv + .39451 * sin(6 * d - l) dv = dv - .38213 * sin(2 * d - lp + 2 * f) dv = dv - .37451 * sin(2 * d - lp + l - 2 * f) dv = dv - .35758 * sin(4 * d + lp - 2 * l) dv = dv + .34965 * sin(d + lp - 2 * l) dv = dv + .33979 * sin(2 * d - 3 * lp) dv = dv - .32866 * sin(3 * l + 2 * f) dv = dv + .30872 * sin(4 * d - 2 * lp - l) dv = dv + .30155 * sin(lp - l - 2 * f) dv = dv + .30086 * sin(4 * d - l - 2 * f) dv = dv + .2942 * sin(2 * d - 2 * lp - 2 * l) dv = dv + .29255 * sin(6 * d - 3 * l) dv = dv - .29022 * sin(2 * d + lp + 2 * l) ' a(p,2,v) series dv = dv + .00487 * t2 * sin(lp) dv = dv - .0015 * t2 * sin(2 * d - lp - l + pi) dv = dv - .0012 * t2 * sin(2 * d - lp + pi) dv = dv + .00108 * t2 * sin(lp - l) dv = dv + .0008 * t2 * sin(lp + l) ' a(p,0,v) series dv = dv + 14.24883 * sin(18 * ve - 16 * t - l + dtr * 26.54261) dv = dv + 7.06304 * sin(ll - f + dtr * .00094) dv = dv + 1.14307 * sin(2 * t - 2 * ju + 2 * d - l + dtr * 180.11977) dv = dv + .901140 * sin(4 * t - 8 * ma + 3 * ju + dtr * 285.98707) dv = dv + .82155 * sin(ve - t + dtr * 180.00988) dv = dv + .78811 * sin(18 * ve - 16 * t - 2 * l + dtr * 26.54324) dv = dv + .7393 * sin(18 * ve - 16 * t + dtr * 26.5456) dv = dv + .64371 * sin(3 * ve - 3 * t + 2 * d - l + dtr * 179.98144) dv = dv + .6388 * sin(t - ju + dtr * 1.2289) dv = dv + .56341 * sin(10 * ve - 3 * t - l + dtr * 333.30551) dv = dv + .49331 * sin(ll + l - f + .00127 * dtr) dv = dv + .49141 * sin(ll - l - f + .00127 * dtr) dv = dv + .44532 * sin(2 * t - 3 * ju + 2 * d - l + 10.07001 * dtr) dv = dv + .36061 * sin(ll + f + .00071 * dtr) dv = dv + .34355 * sin(2 * ve - 3 * t + 269.95393 * dtr) dv = dv + .32455 * sin(t - 2 * ma + 318.13776 * dtr) dv = dv + .30155 * sin(2 * ve - 2 * t + .20448 * dtr) dv = dv + .28938 * sin(t + d - f + 95.13523 * dtr) dv = dv + .28281 * sin(2 * t - 3 * ju + 2 * d - 2 * l + 10.03835 * dtr) dv = dv + .24515 * sin(2 * t - 2 * ju + 2 * d - 2 * l + .08642 * dtr) ' a(p,1,v) series dv = dv + 1.6768 * t1 * sin(lp) dv = dv + .51642 * t1 * sin(2 * d - lp - l + pi) dv = dv + .41383 * t1 * sin(2 * d - lp + pi) dv = dv + .37115 * t1 * sin(lp - l) dv = dv + .2756 * t1 * sin(lp + l) dv = dv + .25425 * t1 * sin(18 * ve - 16 * t - l + 114.5655 * dtr) dv = dv + 7.1178e-02 * t1 * sin(2 * d + lp - l) dv = dv + .06128 * t1 * sin(2 * d + lp) dv = dv + .04516 * t1 * sin(d + lp + pi) dv = dv + .04048 * t1 * sin(2 * d - 2 * lp + pi) dv = dv + .03747 * t1 * sin(2 * lp) dv = dv + .03707 * t1 * sin(2 * d - 2 * lp - l + pi) dv = dv + .03649 * t1 * sin(2 * d - lp + l + pi) dv = dv + .02438 * t1 * sin(lp - 2 * l) dv = dv + .02165 * t1 * sin(2 * d - lp - 2 * l + pi) dv = dv + .01923 * t1 * sin(lp + 2 * l) plon = modulo(ll + atr * dv + dpsi) ' compute geocentric ecliptic latitude (arc seconds) ' a(c,0,u) series pl = 18461.23868 * sin(f) pl = pl + 1010.16707 * sin(l + f) pl = pl + 999.69358 * sin(l - f) pl = pl + 623.65243 * sin(2 * d - f) pl = pl + 199.48374 * sin(2 * d - l + f) pl = pl + 166.5741 * sin(2 * d - l - f) pl = pl + 117.26069 * sin(2 * d + f) pl = pl + 61.91195 * sin(2 * l + f) pl = pl + 33.3572 * sin(2 * d + l - f) pl = pl + 31.75967 * sin(2 * l - f) pl = pl + 29.57658 * sin(2 * d - lp - f) pl = pl + 15.56626 * sin(2 * d - 2 * l - f) pl = pl + 15.12155 * sin(2 * d + l + f) pl = pl - 12.09414 * sin(2 * d + lp - f) pl = pl + 8.86814 * sin(2 * d - lp - l + f) pl = pl + 7.95855 * sin(2 * d - lp + f) pl = pl + 7.43455 * sin(2 * d - lp - l - f) pl = pl - 6.73143 * sin(lp - l - f) pl = pl + 6.57957 * sin(4 * d - l - f) pl = pl - 6.46007 * sin(lp + f) pl = pl - 6.29648 * sin(3 * f) pl = pl - 5.63235 * sin(lp - l + f) pl = pl - 5.3684 * sin(d + f) pl = pl - 5.31127 * sin(lp + l + f) pl = pl - 5.07591 * sin(lp + l - f) pl = pl - 4.83961 * sin(lp - f) pl = pl - 4.80574 * sin(d - f) pl = pl + 3.98405 * sin(3 * l + f) pl = pl + 3.67446 + sin(4 * d - f) pl = pl + 2.99848 * sin(4 * d - l + f) pl = pl + 2.79864 * sin(l - 3 * f) pl = pl + 2.41388 * sin(4 * d - 2 * l + f) pl = pl + 2.18631 * sin(2 * d - 3 * f) pl = pl + 2.14617 * sin(2 * d + 2 * l - f) pl = pl + 1.76598 * sin(2 * d - lp + l - f) pl = pl - 1.62442 * sin(2 * d - 2 * l + f) pl = pl + 1.5813 * sin(3 * l - f) pl = pl + 1.51975 * sin(2 * d + 2 * l + f) pl = pl - 1.51563 * sin(2 * d - 3 * l - f) pl = pl - 1.31782 * sin(2 * d + lp - l + f) pl = pl - 1.26427 * sin(2 * d + lp + f) pl = pl + 1.19187 * sin(4 * d + f) pl = pl + 1.13461 * sin(2 * d - lp + l + f) pl = pl + 1.08578 * sin(2 * d - 2 * lp - f) pl = pl - 1.01938 * sin(l + 3 * f) pl = pl - .822710 * sin(2 * d + lp + l - f) pl = pl + .80422 * sin(d + lp - f) pl = pl + .80259 * sin(d + lp + f) pl = pl - .79319 * sin(lp - 2 * l - f) pl = pl - .79101 * sin(2 * d + lp - l - f) pl = pl - .66741 * sin(d + l + f) pl = pl + .65022 * sin(2 * d - lp - 2 * l - f) pl = pl - .63881 * sin(lp + 2 * l + f) pl = pl + .63371 * sin(4 * d - 2 * l - f) pl = pl + .59577 * sin(4 * d - lp - l - f) pl = pl - .58893 * sin(d + l - f) pl = pl + .47338 * sin(4 * d + l - f) pl = pl - .42989 * sin(d - l - f) pl = pl + .41494 * sin(4 * d - lp - f) pl = pl + .3835 * sin(2 * d - 2 * lp + f) pl = pl - .35183 * sin(3 * d - f) pl = pl + .33881 * sin(4 * d - lp - l + f) pl = pl + .32906 * sin(2 * d - l - 3 * f) ' a(p,0,u) series pl = pl + 8.04508 * sin(ll + 180.00071 * dtr) pl = pl + 1.51021 * sin(t + d + 276.68007 * dtr) pl = pl + .63037 * sin(18 * ve - 16 * t - l + f + 26.54287 * dtr) pl = pl + .63014 * sin(18 * ve - 16 * t - l - f + 26.54272 * dtr) pl = pl + .45586 * sin(ll - l + .00075 * dtr) pl = pl + .41571 * sin(ll + l + 180.00069 * dtr) pl = pl + .32622 * sin(ll - 2 * f + .00086 * dtr) pl = pl + .29854 * sin(ll - 2 * d + .00072 * dtr) ' a(p,1,u) series pl = pl + .0743 * t1 * sin(2 * d - lp - f + pi) pl = pl + .03043 * t1 * sin(2 * d + lp - f) pl = pl + .02229 * t1 * sin(2 * d - lp - l + f + pi) pl = pl + .01999 * t1 * sin(2 * d - lp + f + pi) pl = pl + .01869 * t1 * sin(2 * d - lp - l - f + pi) pl = pl + .01696 * t1 * sin(lp - l - f) pl = pl + .01623 * t1 * sin(lp + f) plat = atr * pl ' compute geocentric right ascension and declination a = sin(plon) * cos(obliq) - tan(plat) * sin(obliq) b = cos(plon) rasc = atan3(a, b) decl = asin(sin(plat) * cos(obliq) + cos(plat) * sin(obliq) * sin(plon)) ' compute the geocentric position vector of the moon (kilometers) rmoon(1) = rmm * cos(rasc) * cos(decl) rmoon(2) = rmm * sin(rasc) * cos(decl) rmoon(3) = rmm * sin(decl) end sub ''''''''''''''''''''''''''''''''''''' ''''''''''''''''''''''''''''''''''''' sub realroot(x1, x2, tol, xroot, froot) ' real root of a single non-linear function subroutine ' input ' x1 = lower bound of search interval ' x2 = upper bound of search interval ' tol = convergence criter%ia ' output ' xroot = real root of f(x) = 0 ' froot = function value ' note: requires sub jdfunc ''''''''''''''''''''''''''' local eps, a, b, c, d, e, fa, fb, fcc, tol1 local xm, p, q, r, s, xmin, tmp eps = 2.23e-16 e = 0.0 a = x1 b = x2 jdfunc(a, fa) jdfunc(b, fb) fcc = fb for iter% = 1 to 50 if (fb * fcc > 0.0) then c = a fcc = fa d = b - a e = d end if if (abs(fcc) < abs(fb)) then a = b b = c c = a fa = fb fb = fcc fcc = fa end if tol1 = 2.0 * eps * abs(b) + 0.5 * tol xm = 0.5 * (c - b) if (abs(xm) <= tol1 or fb = 0) then exit for if (abs(e) >= tol1 and abs(fa) > abs(fb)) then s = fb / fa if (a = c) then p = 2.0 * xm * s q = 1.0 - s else q = fa / fcc r = fb / fcc p = s * (2.0 * xm * q * (q - r) - (b - a) * (r - 1.0)) q = (q - 1.0) * (r - 1.0) * (s - 1.0) end if if (p > 0) then q = -q p = abs(p) xmin = abs(e * q) tmp = 3.0 * xm * q - abs(tol1 * q) if (xmin < tmp) then xmin = tmp if (2.0 * p < xmin) then e = d d = p / q else d = xm e = d end if else d = xm e = d end if a = b fa = fb if (abs(d) > tol1) then b = b + d else b = b + sgn(xm) * tol1 end if jdfunc(b, fb) next iter% froot = fb xroot = b end sub '''''''''''''''''''''''''''''''''''' '''''''''''''''''''''''''''''''''''' sub obliq_lp(jdate, dpsi, deps, obliq) ' nutations and true obliquity ' input ' jdate = julian date ' output ' dpsi = nutation in longitude in radians ' deps = nutation in obliquity in radians ' obliq = true obliquity of the ecliptic in radians '''''''''''''''''''''''''''''''''''''''''''''''''''' LOCAL t, t2, t3, eqeq LOCAL th, tl, obm, obt, st LOCAL tjdh, tjdl ' split the julian date tjdh = int(jdate) tjdl = jdate - tjdh ' fundamental time units th = (tjdh - 2451545.0) / 36525.0 tl = tjdl / 36525.0 t = th + tl t2 = t * t t3 = t2 * t ' obtain equation of the equinoxes eqeq = 0.0 ' obtain nutations nut2000_lp(jdate, dpsi, deps) ' compute mean obliquity of the ecliptic in seconds of arc obm = 84381.4480 - 46.8150 * t - 0.00059 * t2 + 0.001813 * t3 ' compute true obliquity of the ecliptic in seconds of arc obt = obm + deps ' return elements in radians deps = atr * deps dpsi = atr * dpsi obliq = atr * obt END sub ''''''''''''''''''''''''''''''' ''''''''''''''''''''''''''''''' sub nut2000_lp(jdate, dpsi, deps) ' low precison nutation based on iau 2000a ' this function evaluates a short nutation series and returns approximate ' values for nutation in longitude and nutation in obliquity for a given ' tdb julian date. in this mode, only the largest 13 terms of the iau 2000a ' nutation series are evaluated. ' input ' jdate = tdb julian date ' output ' dpsi = nutation in longitude in arcseconds ' deps = nutation in obliquity in arcseconds ''''''''''''''''''''''''''''''''''''''''''''' LOCAL rev = 360.0 * 3600.0 LOCAL el, elp, f, d, omega LOCAL i%, arg LOCAL t = (jdate - 2451545.0) / 36525.0 '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' ' computation of fundamental (delaunay) arguments from simon et al. (1994) '''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''' ' mean anomaly of the moon el = (485868.249036 + t * (1717915923.2178 + t * (31.8792 + t * (0.051635 + t * (-0.00024470)))) mod rev) / seccon ' mean anomaly of the sun elp = (1287104.79305 + t * (129596581.0481 + t * (-0.5532 + t * (0.000136 + t * (- 0.00001149)))) mod rev) / seccon ' mean argument of the latitude of the moon f = (335779.526232 + t * (1739527262.8478 + t * (-12.7512 + t * (-0.001037 + t * (0.00000417)))) mod rev) / seccon ' mean elongation of the moon from the sun d = (1072260.70369 + t * (1602961601.2090 + t * (- 6.3706 + t * (0.006593 + t * (- 0.00003169)))) mod rev) / seccon ' mean longitude of the ascending node of the moon (from simon section 3.4(b.3), precession = 5028.8200 arcsec/cy) omega = (450160.398036 + t * (- 6962890.5431 + t * (7.4722 + t * (0.007702 + t * (- 0.00005939)))) mod rev) / seccon dpsi = 0.0 deps = 0.0 ' sum nutation series terms for i% = 13 to 1 step -1 arg = xnut(1, i%) * el + xnut(2, i%) * elp + xnut(3, i%) * f + xnut(4, i%) * d + xnut(5, i%) * omega dpsi = (xnut(6, i%) + xnut(7, i%) * t) * sin(arg) + xnut(10, i%) * cos(arg) + dpsi deps = (xnut(8, i%) + xnut(9, i%) * t) * cos(arg) + xnut(11, i%) * sin(arg) + deps next i% dpsi = 1.0e-7 * dpsi deps = 1.0e-7 * deps ' add in out-of-phase component of principal (18.6-year) term ' (to avoid small but long-term bias in results) dpsi = dpsi + 0.0033 * cos(omega) deps = deps + 0.0015 * sin(omega) END sub '''''''''''''''''''''''''''' '''''''''''''''''''''''''''' sub getdate (month, day, year) ' request calendar date subroutine do print " " print "please input the calendar date" print " " print "(month [1 - 12], day [1 - 31], year [yyyy])" print "< for example, october 21, 1986 is input as 10,21,1986 >" print "< b.c. dates are negative, a.d. dates are positive >" print "< the day of the month may also include a decimal part >" print " " input month, day, year loop until (month >= 1 and month <= 12) and (day >= 1 and day <= 31) end sub '''''''''''''''''''''''''''''''' '''''''''''''''''''''''''''''''' sub julian(month, day, year, jday) ' Gregorian date to julian day subroutine ' input ' month = calendar month ' day = calendar day ' year = calendar year (all four digits) ' output ' jday = julian day ' special notes ' (1) calendar year must include all digits ' (2) will report October 5, 1582 to October 14, 1582 ' as invalid calendar dates and exit ''''''''''''''''''''''''''''''''''''''''' local a, b, c, m, y y = year m = month b = 0.0 c = 0.0 if (m <= 2.0) then y = y - 1.0 m = m + 12.0 end if if (y < 0.0) then c = -0.75 if (year < 1582.0) then ' null elseif (year > 1582.0) then a = fix(y / 100.0) b = 2.0 - a + fix(a / 4.0) elseif (month < 10.0) then ' null elseif (month > 10.0) then a = fix(y / 100.0) b = 2.0 - a + fix(a / 4.0) elseif (day <= 4.0) then ' null elseif (day > 14.0) then a = fix(y / 100.0) b = 2.0 - a + fix(a / 4.0) else print "this date does not exist!!" exit end if jday = fix(365.25 * y + c) + fix(30.6001 * (m + 1.0)) + day + b + 1720994.5 end sub '''''''''''''''''''''''''''''''' '''''''''''''''''''''''''''''''' sub gdate (jday, month, day, year) ' Julian day to Gregorian date subroutine ' input ' jday = julian day ' output ' month = calendar month ' day = calendar day ' year = calendar year '''''''''''''''''''''''' local a, b, c, d, e, f, z, alpha z = fix(jday + 0.5) f = jday + 0.5 - z if (z < 2299161) then a = z else alpha = fix((z - 1867216.25) / 36524.25) a = z + 1.0 + alpha - fix(alpha / 4.0) end if b = a + 1524.0 c = fix((b - 122.1) / 365.25) d = fix(365.25 * c) e = fix((b - d) / 30.6001) day = b - d - fix(30.6001 * e) + f if (e < 13.5) then month = e - 1.0 else month = e - 13.0 end if if (month > 2.5) then year = c - 4716.0 else year = c - 4715.0 end if end sub '''''''''''''''''''''''' '''''''''''''''''''''''' sub utc2tdb (jdutc, jdtdb) ' convert UTC julian date to TDB julian date ' input ' jdutc = UTC julian day ' output ' jdtdb = TDB julian day ' Reference Frames in Astronomy and Geophysics ' J. Kovalevsky et al., 1989, pp. 439-442 ''''''''''''''''''''''''''''''''''''''''' local corr, jdtt, t, leapsecond ' find current number of leap seconds findleap(jdutc, leapsecond) ' compute TDT julian date corr = (leapsecond + 32.184) / 86400.0 jdtt = jdutc + corr ' time argument for correction t = (jdtt - 2451545.0) / 36525.0 ' compute correction in microseconds corr = 1656.675 * sin(dtr * (35999.3729 * t + 357.5287)) corr = corr + 22.418 * sin(dtr * (32964.467 * t + 246.199)) corr = corr + 13.84 * sin(dtr * (71998.746 * t + 355.057)) corr = corr + 4.77 * sin(dtr * ( 3034.906 * t + 25.463)) corr = corr + 4.677 * sin(dtr * (34777.259 * t + 230.394)) corr = corr + 10.216 * t * sin(dtr * (35999.373 * t + 243.451)) corr = corr + 0.171 * t * sin(dtr * (71998.746 * t + 240.98 )) corr = corr + 0.027 * t * sin(dtr * ( 1222.114 * t + 194.661)) corr = corr + 0.027 * t * sin(dtr * ( 3034.906 * t + 336.061)) corr = corr + 0.026 * t * sin(dtr * ( -20.186 * t + 9.382)) corr = corr + 0.007 * t * sin(dtr * (29929.562 * t + 264.911)) corr = corr + 0.006 * t * sin(dtr * ( 150.678 * t + 59.775)) corr = corr + 0.005 * t * sin(dtr * ( 9037.513 * t + 256.025)) corr = corr + 0.043 * t * sin(dtr * (35999.373 * t + 151.121)) ' convert corrections to days corr = 0.000001 * corr / 86400.0 ' TDB julian date jdtdb = jdtt + corr end sub '''''''''''''''''''''''''''' '''''''''''''''''''''''''''' sub findleap(jday, leapsecond) ' find number of leap seconds for utc julian day ' input ' jday = utc julian day ' input via global ' jdleap = array of utc julian dates ' leapsec = array of leap seconds ' output ' leapsecond = number of leap seconds '''''''''''''''''''''''''''''''''''''' if (jday <= jdleap(1)) then ' date is <= 1972; set to first data element leapsecond = leapsec(1) exit sub end if if (jday >= jdleap(28)) then ' date is >= end of current data ' set to last data element leapsecond = leapsec(28) exit sub end if ' find data within table for i% = 1 to 27 if (jday >= jdleap(i%) and jday < jdleap(i% + 1)) then leapsecond = leapsec(i%) exit sub end if next i% end sub ''''''''''''''''''''''''' ''''''''''''''''''''''''' function modulo(x) as float ' modulo 2 pi function '''''''''''''''''''''' local a a = x - pi2 * fix(x / pi2) if (a < 0.0) then a = a + pi2 end if modulo = a end function ''''''''''''''''''''''''''' ''''''''''''''''''''''''''' function atan3(a, b) as float ' four quadrant inverse tangent function ' input ' a = sine of angle ' b = cosine of angle ' output ' atan3 = angle (0 =< atan3 <= 2 * pi; radians) '''''''''''''''''''''''''''''''''''''''''''''''' local c if (abs(a) < 1.0e-10) then atan3 = (1.0 - sgn(b)) * pidiv2 exit function else c = (2.0 - sgn(a)) * pidiv2 endif if (abs(b) < 1.0e-10) then atan3 = c exit function else atan3 = c + sgn(a) * sgn(b) * (abs(atn(a / b)) - pidiv2) endif end function '''''''''''''''''' '''''''''''''''''' function vecmag(a()) ' vector magnitude function ' input ' { a } = column vector ( 3 rows by 1 column ) ' output ' vecmag = scalar magnitude of vector { a } vecmag = sqr(a(1) * a(1) + a(2) * a(2) + a(3) * a(3)) end function '''''''''''''''''''' '''''''''''''''''''' sub uvector (a(), b()) ' unit vector subroutine ' input ' a = column vector (3 rows by 1 column) ' output ' b = unit vector (3 rows by 1 column) ''''''''''''''''''''''''''''''''''''''' local i as integer, amag amag = vecmag(a()) for i = 1 to 3 if (amag <> 0.0) then b(i) = a(i) / amag else b(i) = 0.0 end if next i end sub ''''''''''''''''''''' ''''''''''''''''''''' function vdot(a(), b()) ' vector dot product function ' c = { a } dot { b } ' input ' n% = number of rows ' { a } = column vector with n rows ' { b } = column vector with n rows ' output ' vdot = dot product of { a } and { b } '''''''''''''''''''''''''''''''''''''''' local c = 0.0 for i% = 1 to 3 c = c + a(i%) * b(i%) next i% vdot = c end function ''''''''''''''''''''''' ''''''''''''''''''''''' sub vcross(a(), b(), c()) ' vector cross product subroutine ' { c } = { a } x { b } ' input ' { a } = vector a ( 3 rows by 1 column ) ' { b } = vector b ( 3 rows by 1 column ) ' output ' { c } = { a } x { b } ( 3 rows by 1 column ) c(1) = a(2) * b(3) - a(3) * b(2) c(2) = a(3) * b(1) - a(1) * b(3) c(3) = a(1) * b(2) - a(2) * b(2) end sub '''''''''''''''''''''''' '''''''''''''''''''''''' sub matxvec(a(), b(), c()) ' matrix/vector multiplication subroutine ' { c } = [ a ] * { b } ' input ' a = matrix a ( 3 rows by 3 columns ) ' b = vector b ( 3 rows ) ' output ' c = vector c ( 3 rows ) '''''''''''''''''''''''''' local s, i%, j% for i% = 1 to 3 s = 0.0 for j% = 1 to 3 s = s + a(i%, j%) * b(j%) next j% c(i%) = s next i% end sub '''''''''''''''''''''' '''''''''''''''''''''' sub transpose (a(), b()) ' matrix traspose subroutine ' input ' m = number of rows in matrix [ a ] ' n = number of columns in matrix [ a ] ' a = matrix a ( 3 rows by 3 columns ) ' output ' b = matrix transpose ( 3 rows by 3 columns ) ''''''''''''''''''''''''''''''''''''''''''''''' local i%, j% for i% = 1 to 3 for j% = 1 to 3 b(i%, j%) = a(j%, i%) next j% next i% end sub ''''''''''''''' ''''''''''''''' sub jd2str(jdutc) ' convert julian day to calendar date and UTC time '''''''''''''''''''''''''''''''''''''''''''''''''' gdate (jdutc, cmonth, day, year) print "calendar date ", month$(cmonth), " ", STR$(int(day)), " ", str$(year) print " " thr0 = 24.0 * (day - int(day)) thr = int(thr0) tmin0 = 60.0 * (thr0 - thr) tmin = int(tmin0) tsec = 60.0 * (tmin0 - tmin) ' fix seconds and minutes for rollover if (tsec >= 60.0) then tsec = 0.0 tmin = tmin + 1.0 end if ' fix minutes for rollover if (tmin >= 60.0) then tmin = 0.0 thr = thr + 1.0 end if print "UTC time ", str$(thr) + " hours " + str$(tmin) + " minutes " + str$(tsec, 0, 2) + " seconds" end sub ''''''''''''''''''' ''''''''''''''''''' sub deg2str(dd, dms$) ' convert decimal degrees to degrees, ' minutes, seconds string ' input ' dd = angle in decimal degrees ' output ' dms$ = string equivalent ''''''''''''''''''''''''''' local d1, d, m, s d1 = abs(dd) d = fix(d1) d1 = (d1 - d) * 60.0 m = fix(d1) s = (d1 - m) * 60.0 if (dd < 0.0) then if (d <> 0.0) then d = -d elseif (m <> 0.0) then m = -m else s = -s end if end if dms$ = str$(d) + " deg " + str$(m) + " min " + str$(s, 0, 2) + " sec" end sub '''''''''''''''' '''''''''''''''' sub hrs2str(hours) ' convert hours to equivalent string '''''''''''''''''''''''''''''''''''' local thr, tmin0, tmin, tsec thr = fix(hours) tmin0 = 60.0 * (hours - thr) tmin = fix(tmin0) tsec = 60.0 * (tmin0 - tmin) ' fix seconds and minutes for rollover if (tsec >= 60.0) then tsec = 0.0 tmin = tmin + 1.0 end if ' fix minutes for rollover if (tmin >= 60.0) then tmin = 0.0 thr = thr + 1 end if print "event duration ", str$(thr) + " hours " + str$(tmin) + " minutes " + str$(tsec, 0, 2) + " seconds" end sub
|