' seasons.bas March 6, 2017 ' calendar date and UTC time of the seasons ' Micromite eXtreme version ''''''''''''''''''''''''''' ' dimension global arrays dim sl(50) as float, sr(50) as float, sa(50) as float, sb(50) as float dim jdleap(28) as float, leapsec(28) as float, month$(12) as string ' global constants const pi2 = 2.0 * pi, pidiv2 = 0.5 * pi, dtr = pi / 180.0 ' read solar ephemeris data for i% = 1 to 50 read sl(i%), sr(i%), sa(i%), sb(i%) next i% data 403406, 0, 4.721964, 1.621043 data 195207, -97597, 5.937458, 62830.348067 data 119433, -59715, 1.115589, 62830.821524 data 112392, -56188, 5.781616, 62829.634302 data 3891, -1556, 5.5474 , 125660.5691 data 2819, -1126, 1.5120 , 125660.9845 data 1721, -861, 4.1897 , 62832.4766 data 0, 941, 1.163 , 0.813 data 660, -264, 5.415 , 125659.310 data 350, -163, 4.315 , 57533.850 data 334, 0, 4.553 , -33.931 data 314, 309, 5.198 , 777137.715 data 268, -158, 5.989 , 78604.191 data 242, 0, 2.911 , 5.412 data 234, -54, 1.423 , 39302.098 data 158, 0, 0.061 , -34.861 data 132, -93, 2.317 , 115067.698 data 129, -20, 3.193 , 15774.337 data 114, 0, 2.828 , 5296.670 data 99, -47, 0.52 , 58849.27 data 93, 0, 4.65 , 5296.11 data 86, 0, 4.35 , -3980.70 data 78, -33, 2.75 , 52237.69 data 72, -32, 4.50 , 55076.47 data 68, 0, 3.23 , 261.08 data 64, -10, 1.22 , 15773.85 data 46, -16, 0.14 , 188491.03 data 38, 0, 3.44 , -7756.55 data 37, 0, 4.37 , 264.89 data 32, -24, 1.14 , 117906.27 data 29, -13, 2.84 , 55075.75 data 28, 0, 5.96 , -7961.39 data 27, -9, 5.09 , 188489.81 data 27, 0, 1.72 , 2132.19 data 25, -17, 2.56 , 109771.03 data 24, -11, 1.92 , 54868.56 data 21, 0, 0.09 , 25443.93 data 21, 31, 5.98 , -55731.43 data 20, -10, 4.03 , 60697.74 data 18, 0, 4.27 , 2132.79 data 17, -12, 0.79 , 109771.63 data 14, 0, 4.24 , -7752.82 data 13, -5, 2.01 , 188491.91 data 13, 0, 2.65 , 207.81 data 13, 0, 4.98 , 29424.63 data 12, 0, 0.93 , -7.99 data 10, 0, 2.21 , 46941.14 data 10, 0, 3.59 , -68.29 data 10, 0, 1.50 , 21463.25 data 10, -9, 2.55 , 157208.40 ' 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 "time of the seasons" print "===================" print " " ' request calendar year print "please input the calendar year" input year ' process each season for iequsol% = 1 to 4 select case iequsol% case (1) cmonth = 3 day = 15 julian(cmonth, day, year, jdayi) case (2) cmonth = 6 day = 15 julian(cmonth, day, year, jdayi) along2 = 0.5 * pi case (3) cmonth = 9 day = 15 julian(cmonth, day, year, jdayi) case (4) cmonth = 12 day = 15 julian(cmonth, day, year, jdayi) along2 = 1.5 * pi end select ' find event x1 = 0.0 x2 = 10.0 realroot1(x1, x2, 1.0e-8, xroot, froot) ' TDB julian day of event jdtdb = jdayi + xroot ' compute UTC julian day tdb2utc(jdtdb, jdutc) ' print results for this event print " " select case iequsol% case (1) print "spring equinox" print "--------------" case (2) print "summer solstice" print "---------------" case (3) print "fall equinox" print "------------" case (4) print "winter solstice" print "---------------" end select print " " 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) print "UTC time ", str$(thr) + " hours " + str$(tmin) + " minutes " + str$(tsec, 0, 2) + " seconds" next iequsol% print " " end ''''''''''''''' ''''''''''''''' sub esfunc(x, fx) ' equinox/solstice objective function ''''''''''''''''''''''''''''''''''''' local jday, rlsun, rasc, decl jday = jdayi + x solar(jday, rlsun, rasc, decl) if (iequsol% = 1 or iequsol% = 3) then fx = decl else fx = along2 - rlsun end if 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 realroot2(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, tai_utc findleap(jdin, tai_utc) utc2tdb(jdin, tai_utc, jdwrk) fx = jdwrk - jdsaved end sub ''''''''''''''''''''''''''''''' ''''''''''''''''''''''''''''''' sub solar (jd, rlsun, rasc, decl) ' precision ephemeris of the Sun ' input ' jd = julian ephemeris day ' output ' rlsun = ecliptic longitude of the sun ' (0 <= rlsun <= 2 pi) ' rasc = right ascension of the Sun (radians) ' (0 <= rasc <= 2 pi) ' decl = declination of the Sun (radians) ' (-pi/2 <= decl <= pi/2) '''''''''''''''''''''''''''''''''' local u, a1, a2, psi, deps, meps, eps, seps, ceps local dl, dr, w, srl, crl, srb, crb, sra, cra u = (jd - 2451545.0) / 3652500.0 ' compute nutation in longitude a1 = 2.18 + u * (-3375.7 + u * 0.36) a2 = 3.51 + u * (125666.39 + u * 0.1) psi = 0.0000001 * (-834.0 * sin(a1) - 64.0 * sin(a2)) ' compute nutation in obliquity deps = 0.0000001 * u * (-226938 + u * (-75 + u * (96926 + u * (-2491 - u * 12104)))) meps = 0.0000001 * (4090928.0 + 446.0 * cos(a1) + 28.0 * cos(a2)) eps = meps + deps seps = sin(eps) ceps = cos(eps) dl = 0.0 dr = 0.0 for i% = 1 to 50 w = sa(i%) + sb(i%) * u dl = dl + sl(i%) * sin(w) if (sr(i%) <> 0.0) then dr = dr + sr(i%) * cos(w) end if next i% dl = modulo(dl * 0.0000001 + 4.9353929 + 62833.196168 * u) dr = 149597870.691 * (dr * 0.0000001 + 1.0001026) rlsun = modulo(dl + 0.0000001 * (-993.0 + 17.0 * cos(3.1 + 62830.14 * u)) + psi) rb = 0.0 ' compute geocentric declination and right ascension crl = cos(rlsun) srl = sin(rlsun) crb = cos(rb) srb = sin(rb) decl = asin(ceps * srb + seps * crb * srl) sra = -seps * srb + ceps * crb * srl cra = crb * crl rasc = atan3(sra, cra) end sub '''''''''''''''''''''''''''''''''''''' '''''''''''''''''''''''''''''''''''''' sub realroot1(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 esfunc ''''''''''''''''''''''''''' 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 esfunc(a, fa) esfunc(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.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.0) then q = -q p = abs(p) min = abs(e * q) tmp = 3.0 * xm * q - abs(tol1 * q) if (min < tmp) then min = tmp if (2.0 * p < min) 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 esfunc(b, fb) next iter% froot = fb xroot = b end sub '''''''''''''''''''''''''''''''''''''' '''''''''''''''''''''''''''''''''''''' sub realroot2(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 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, tai_utc, jdtdb) ' convert UTC julian date to TDB julian date ' input ' jdutc = UTC julian day ' tai_utc = TAI-UTC (seconds) ' output ' jdtdb = TDB julian day ' Reference Frames in Astronomy and Geophysics ' J. Kovalevsky et al., 1989, pp. 439-442 ''''''''''''''''''''''''''''''''''''''''' local corr, jdtdt, t ' TDT julian date corr = (tai_utc + 32.184) / 86400.0 jdtdt = jdutc + corr ' time argument for correction t = (jdtdt - 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 = jdtdt + 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