' sun_riseset.bas March 17, 2017

' rise and set of the sun

' Micromite eXtreme version

'''''''''''''''''''''''''''

option default float

option base 1

' dimension global arrays and variables

dim xsl(50), xsr(50), xsa(50), xsb(50)

dim jdleap(28), leapsec(28)

dim month$(12) as string

dim xnut(11, 13), trr, elev_minima

dim jdsaved, jdprint, obslat, obslong, obsalt

dim jdtdbi, i as integer, j as integer

dim cmonth, cday, cyear, ndays

dim rootflag% as integer

' 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

' astronomical unit (kilometers)

const aunit = 149597870.691

' equatorial radius of the earth (kilometers)

const reqm = 6378.14

' earth flattening factor (nd)

const flat = 1.0 / 298.257

' read solar ephemeris data

for i = 1 to 50

read xsl(i), xsr(i), xsa(i), xsb(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 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 "rise and set of the sun"
print "======================="
print " "

' request initial calendar date (month, day, year)

getdate(cmonth, cday, cyear)

' request observer coordinates

print " "

observer(obslat, obslong, obsalt)

' initial utc julian day

julian(cmonth, cday, cyear, jdutc)

' compute initial tdb julian date

utc2tdb(jdutc, jdtdb)

jdtdbi = jdtdb

' search duration (days)

print " "

print "please input the number of days to simulate"

input ndays

print " "
print "searching for rise and set conditions ..."
print " "

' define search parameters

ti = 0.0

tf = ndays

dt = 0.1

dtsml = 0.01

rootflag% = 0

' find sun rise/set conditions

rs_event(ti, tf, dt, dtsml)

end

'''''''''''''''
'''''''''''''''

sub rsfunc(x, fx)

' rise/set objective function

'''''''''''''''''''''''''''''

local jdtdb, rsun(3), gast, azim, elev

local dis, dref, sdia

' current TDB julian day

jdtdb = jdtdbi + x

' compute eci position vector and distance of the sun

sun(jdtdb, rsun())

dis = vecmag(rsun()) / aunit

' compute topocentric coordinates of the sun

tdb2utc(jdtdb, jdutc)

gast2(jdutc, gast)

eci2topo(gast, rsun(), azim, elev)

if (rootflag% = 1) then

' correct for horizontal refraction

dref = (34.0 / 60.0)

' correct for semidiameter of the sun

tsdia = rtd * asin((0.5 * 696000.0 / aunit) / dis)

else

dref = 0.0

tsdia = 0.0

end if

' evaluate objective function

fx = -(rtd * elev + tsdia + dref)

end sub

''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''

sub rs_event (ti, tf, dt, dtsml)

' predict rise/set 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, iend as integer

local fmin1, tmin1

LOCAL ftemp, df, dflft

local el, er

LOCAL t, ft

local iter1 as integer, iter2 as integer

local iter3 as integer

' initialization

tolm = 1.0e-8

iend = 0

' check the initial time for a minimum

rsfunc(ti, fmin1)

tmin1 = ti

rsfunc(ti + dtsml, ftemp)

df = ftemp - fmin1

dflft = df

el = ti

er = el

t = ti

' if the slope is positive and the minimum is
' negative, calculate event conditions at the initial time

if (df > 0.0 and fmin1 < 0.0) then

events1(ti, tf, tmin1)

end if

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

if (t >= tf) then

' check final time for a minimum

if (iend = 1) then exit sub

rsfunc(tf, fmin1)

rsfunc(tf - dtsml, ftemp)

df = fmin1 - ftemp

' set minimum time to final simulation time

tmin1 = tf

if (df < 0.0) then

' if both the slope and minimum are negative,
' calculate the event conditions at the final
' simulation time

if (fmin1 < 0.0) then

events1(ti, tf, tmin1)

end if

' otherwise, we're finished

exit sub

end if

if (dflft > 0.0) then exit sub

er = tf

iend = 1

end if

rsfunc(t, ft)

rsfunc(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

if (t >= tf) then

' check final time for a minimum

if (iend = 1) then exit sub

rsfunc(tf, fmin1)

rsfunc(tf - dtsml, ftemp)

df = fmin1 - ftemp

' set minimum time to final simulation time

tmin1 = tf

if (df < 0.0) then

' if both the slope and minimum are negative,
' calculate the event conditions at the final
' simulation time

if (fmin1 < 0.0) then

events1(ti, tf, tmin1)

end if

' otherwise, we're finished

exit sub

end if

if (dflft > 0.0) then exit sub

er = tf

iend = 1

end if

rsfunc(t, ft)

rsfunc(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

dflft = df

' if the minimum is negative,
' calculate event conditions for this minimum

if (fmin1 < 0.0) then

events1(ti, tf, tmin1)

t = trr

end if

next iter1

end sub

''''''''''''''''''''''''
''''''''''''''''''''''''

sub events1 (ti, tf, topt)

' compute and display rise/set events

' input

' ti = initial simulation time
' tf = final simulation time
' topt = extrema time

''''''''''''''''''''''

' define root-bracketing and root-finding control parameters

LOCAL factor = 0.25 ' geometric acceleration factor

LOCAL dxmax = 0.1 ' rectification interval

LOCAL rtol = 1.0e-8 ' root-finding convergence tolerance

LOCAL t1in, t2in

local t1out, t2out

LOCAL troot, froot, jdate

local iter as integer

' compute and display rise conditions

rootflag% = 1

t1in = topt

t2in = t1in - 0.05

broot(t1in, t2in, factor, dxmax, t1out, t2out)

realroot1(t1out, t2out, rtol, troot, froot)

' set to initial time if before ti

if (troot < ti) then

troot = ti

rsfunc(ti, froot)

end if

jdate = jdtdbi + troot

rsprint(1, jdate)

' display maximum elevation conditions

rootflag% = 0

jdate = jdtdbi + topt

rsprint(2, jdate)

' compute and display set conditions

rootflag% = 1

t2in = t1in + 0.05

broot(t1in, t2in, factor, dxmax, t1out, t2out)

realroot1(t1out, t2out, rtol, troot, froot)

' set to final time if after tf

if (troot > tf) then

troot = tf

rsfunc(tf, froot)

end if

trr = troot

jdate = jdtdbi + troot

rsprint(3, jdate)

rootflag% = 0

END sub

'''''''''''''''''''''''
'''''''''''''''''''''''

sub rsprint(iflag, jdtdb)

' print rise and set conditions

'''''''''''''''''''''''''''''''

local rsun(3)

LOCAL jdutc, dms$ as string

local deltat, gast, azimuth, elevation

if (iflag = 1) then

print " "
print "rise conditions"
PRINT "---------------"
print " "

jdprint = jdtdb

end if

if (iflag = 2) then

print " "
print "maximum elevation conditions"
PRINT "----------------------------"
print " "

end if

if (iflag = 3) then

print " "
print "set conditions"
PRINT "--------------"
print " "

end if

' compute and display UTC julian date

tdb2utc(jdtdb, jdutc)

jd2str(jdutc)

PRINT " "

print "UTC julian day ", str$(jdutc, 0, 8)

' compute and display topocentric coordinates of the sun

gast2(jdutc, gast)

sun(jdtdb, rsun())

eci2topo(gast, rsun(), azimuth, elevation)

PRINT " "

print "topocentric coordinates"

PRINT " "

deg2str(rtd * azimuth, dms$)

print "solar azimuth angle ", dms$

PRINT " "

deg2str(rtd * elevation, dms$)

print "solar elevation angle ", dms$

' determine and display event duration

if (iflag = 3) then

deltat = 24.0 * (jdtdb - jdprint)

PRINT " "

hrs2str(deltat)

print " "

end if

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 usr_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

rsfunc(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

rsfunc(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 gsite (angle, rsite())

' ground site position vector

' input

' angle = local sidereal time or east longitude
' (radians; 0 <= angle <= 2*pi)

' input via global

' obslat = geodetic latitude (radians)
' (+north, -south; -pi/2 <= lat <= -pi/2)
' obsalt = geodetic altitude (meters)
' (+ above sea level, - below sea level)

' output

' rsite = ground site position vector (kilometers)

' special notes

' (1) eci coordinates if angle = local sidereal time

' (2) ecf coordinates if angle = east longitude

' (3) geocentric, equatorial coordinates

'''''''''''''''''''''''''''''''''''''''''

LOCAL slat, clat

local sangle, cangle

LOCAL b, c, d

slat = sin(obslat)

clat = cos(obslat)

sangle = sin(angle)

cangle = cos(angle)

' compute geodetic constants

b = sqr(1.0 - (2.0 * flat - flat * flat) * slat * slat)

c = reqm / b + 0.001 * obsalt

d = reqm * (1.0 - flat) * (1.0 - flat) / b + 0.001 * obsalt

' compute x, y and z components of position vector (kilometers)

rsite(1) = c * clat * cangle

rsite(2) = c * clat * sangle

rsite(3) = d * slat

END sub

''''''''''''''''''''
''''''''''''''''''''

sub gast2(jdate, gast)

' greenwich apparent sidereal time

' input

' jdate = julian date

' output

' gast = greenwich apparent sidereal time (radians)

'''''''''''''''''''''''''''''''''''''''''''''''''''

LOCAL t, t2, t3, eqeq, dpsi, deps

LOCAL th, tl, obm, obt, st, x

local tjdh as integer, tjdl

tjdh = int(jdate)

tjdl = jdate - tjdh

th = (tjdh - 2451545.0) / 36525

tl = tjdl / 36525.0

t = th + tl

t2 = t * t

t3 = t2 * t

' obtain equation of the equinoxes

eqeq = 0.0

' obtain nutation parameters in seconds of arc

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

obliq = obm + deps

' compute equation of the equinoxes in seconds of time

eqeq = (dpsi / 15.0) * cos(obliq / seccon)

st = eqeq - 6.2e-6 * t3 + 0.093104 * t2 + 67310.54841 + 8640184.812866 * tl + 3155760000.0 * tl + 8640184.812866 * th + 3155760000.0 * th

' modulo 24 hours

x = st / 3600.0

gast = x - 24.0 * fix(x / 24.0)

if (gast < 0.0) then

gast = gast + 24.0

end if

' convert to radians

gast = pi2 * (gast / 24.0)

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

utc2tdb(jdin, jdwrk)

fx = jdwrk - jdsaved

end sub

'''''''''''''''''
'''''''''''''''''

sub sun(jd, rsun())

' precision ephemeris of the Sun

' input

' jd = julian ephemeris day

' output

' rsun = eci position of the sun

'''''''''''''''''''''''''''''''''

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

obliq = eps

seps = sin(eps)

ceps = cos(eps)

dl = 0.0

dr = 0.0

for i% = 1 to 50

w = xsa(i%) + xsb(i%) * u

dl = dl + xsl(i%) * sin(w)

if (xsr(i%) <> 0.0) then

dr = dr + xsr(i%) * cos(w)

end if

next i%

dl = modulo(dl * 0.0000001 + 4.9353929 + 62833.196168 * u)

dr = aunit * (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)

' geocentric equatorial position vector of the Sun (au)

rsun(1) = dr * cos(rasc) * cos(decl)

rsun(2) = dr * sin(rasc) * cos(decl)

rsun(3) = dr * sin(decl)

end sub

''''''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''''''

sub eci2topo(gast, robj(), azim, elev)

' convert eci position vector to topocentric coordinates

' input

' gast = Greenwich apparent sidereal time (radians)
' robj = eci position vector of object (kilometers)

' output

' azim = topocentric azimuth (radians)
' elev = topocentric elevation (radians)

''''''''''''''''''''''''''''''''''''''''''''''

local rsite(3), rhoijk(3), rhohatijk(3), rhohatsez(3)

local i as integer, tmatrix(3, 3)

LOCAL obslst, srange, sobslat

local cobslat, sobslst, cobslst

' observer local sidereal time

obslst = modulo(gast + obslong)

gsite(obslst, rsite())

' eci vector from observer to moon

for i = 1 to 3

rhoijk(i) = aunit * robj(i) - rsite(i)

next i

' observer-to-object slant range (kilometers)

srange = vecmag(rhoijk())

' compute topocentric unit pointing vector

uvector(rhoijk(), rhohatijk())

' compute eci-to-sez transformation matrix

sobslat = sin(obslat)
cobslat = cos(obslat)

sobslst = sin(obslst)
cobslst = cos(obslst)

tmatrix(1, 1) = sobslat * cobslst
tmatrix(1, 2) = sobslat * sobslst
tmatrix(1, 3) = -cobslat

tmatrix(2, 1) = -sobslst
tmatrix(2, 2) = cobslst
tmatrix(2, 3) = 0.0

tmatrix(3, 1) = cobslat * cobslst
tmatrix(3, 2) = cobslat * sobslst
tmatrix(3, 3) = sobslat

' compute sez coordinates

matxvec(tmatrix(), rhohatijk(), rhohatsez())

' topocentric elevation (radians)

elev = asin(rhohatsez(3))

' topocentric azimuth (radians)

azim = atan3(rhohatsez(2), -rhohatsez(1))

END sub

''''''''''''''''''''''''''''''''''''''''''''''''
''''''''''''''''''''''''''''''''''''''''''''''''

sub broot(x1in, x2in, factor, dxmax, x1out, x2out)

' bracket a single root of a nonlinear equation

' input

' x1in = initial guess for first bracketing x value
' x2in = initial guess for second bracketing x value
' factor = acceleration factor (non-dimensional)
' dxmax = rectification interval

' output

' xout1 = final value for first bracketing x value
' xout2 = final value for second bracketing x value

''''''''''''''''''''''''''''''''''''''''''''''''''''

LOCAL f1, f2

local x3, dx

' evaluate objective function at initial value

rsfunc(x1in, f1)

' save initial value

x3 = x1in

' save initial delta-x

dx = x2in - x1in

' perform bracketing until the product of the
' two function values is negative

do

' geometrically accelerate the second point

x2in = x2in + factor * (x2in - x3)

' evaluate objective function at x2

rsfunc(x2in, f2)

' check to see if rectification is required

if (abs(x2in - x3) > dxmax) then

x3 = x2in - dx

end if

' is the root bracketed?

if ((f1 * f2) < 0.0) then exit do

loop

x1out = x1in

x2out = x2in

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 rsfunc

'''''''''''''''''''''''''''

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

rsfunc(a, fa)

rsfunc(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

rsfunc(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 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 observer(obslat, obslong, obsalt)

' interactive request of latitude, longitude and altitude subroutine

' output

' obslat = latitude (radians)
' obslong = longitude (radians)
' obsalt = altitude (meters)

''''''''''''''''''''''''''''''

do

print "please input the geographic latitude of the observer"
print "(degrees [-90 to +90], minutes [0 - 60], seconds [0 - 60])"
print "(north latitudes are positive, south latitudes are negative)"

input obslat.deg$, obslat.min, obslat.sec

loop until (abs(val(obslat.deg$)) <= 90.0 and (obslat.min >= 0.0 and obslat.min <= 60.0) and (obslat.sec >= 0.0 and obslat.sec <= 60.0))

if (left$(obslat.deg$, 2) = "-0") then

obslat = -dtr * (obslat.min / 60.0 + obslat.sec / 3600.0)

elseif (val(obslat.deg$) = 0.0) then

obslat = dtr * (obslat.min / 60.0 + obslat.sec / 3600.0)

else

obslat = dtr * (sgn(val(obslat.deg$)) * (abs(val(obslat.deg$)) + obslat.min / 60.0 + obslat.sec / 3600.0))

endif

do

print
print "please input the geographic longitude of the observer"
print "(degrees [0 - 360], minutes [0 - 60], seconds [0 - 60])"
print "(east longitude is positive, west longitude is negative)"

input obslong.deg$, obslong.min, obslong.sec

loop until (abs(val(obslong.deg$)) >= 0.0 and abs(val(obslong.deg$)) <= 360.0) and (obslong.min >= 0.0 and obslong.min <= 60.0) and (obslong.sec >= 0.0 and obslong.sec <= 60.0)

if (left$(obslong.deg$, 2) = "-0") then

obslong = -dtr * (obslong.min / 60 + obslong.sec / 3600)

elseif (val(obslong.deg$) = 0.0) then

obslong = dtr * (obslong.min / 60.0 + obslong.sec / 3600.0)

else

obslong = dtr * (sgn(val(obslong.deg$)) * (abs(val(obslong.deg$)) + obslong.min / 60.0 + obslong.sec / 3600.0))

endif

print " "

print "please input the altitude of the observer (meters)"
print "(positive above sea level, negative below sea level)"

input obsalt

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