mmbasic:mmx_time_of_the_seasons
MMX - time of the seasons
' 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
mmbasic/mmx_time_of_the_seasons.txt · Last modified: 2024/01/19 09:30 by 127.0.0.1