Need help with the math in Python program to flag Jovian radio emissions
A few days ago I decided to take on the little project of converting a Qbasic program into Python (as a side project to my Radio Jove project), and I've managed to get it to run, but the math is definitely off. I was hoping for a fellow astronomer or astrophysicist to assist in the mathematical side of my Python code. 394
Obviously something is wrong with the math aspect of the code since most of the values in the python program are very skewed. My apologies in advance for the single character variables. I would make them better if I knew everything going on.
Qbasic code:
'Program to flag Jovian Decametric windows
m$ = "JANFEBMARAPRMAYJUNJULAUGSEPOCTNOVDEC"
wi = 42.46 / 360
pi = 3.141593
kr = pi / 180
f$ = "### \ \ ## ##.# ### ### ##.## \ \"
OPEN "jovrad.txt" FOR OUTPUT AS #1
INPUT "Year for which predictions are required"; yy
e = INT((yy - 1) / 100)
f = 2 - e + INT(e / 4)
jd = INT(365.25 * (yy - 1)) + 1721423 + f + .5
d0 = jd - 2435108
incr = 0
IF yy / 400 - INT(yy / 400) = 0 THEN incr = 1
yyly = yy / 4 - INT(yy / 4)
yylc = yy / 100 - INT(yy / 100)
IF yyly = 0 AND yylc <> 0 THEN incr = 1
ty = 59 + incr
dmax = 365 + incr
tx = ty + .5
PRINT #1, "******************************************************"
PRINT #1, " JOVIAN IO-DECAMETRIC EMISSION PREDICTIONS FOR"; yy
PRINT #1, "******************************************************"
PRINT #1,
PRINT #1, "Day Date Hr(UT) Io_Phase CML Dist(AU) Source"
PRINT #1,
th = 0
DO
GOSUB Compute
s$ = ""
IF L3 < 255 AND L3 > 200 AND U1 < 250 AND U1 > 220 THEN s$ = "Io-A"
IF L3 < 180 AND L3 > 105 AND U1 < 100 AND U1 > 80 THEN s$ = "Io-B"
IF L3 < 350 AND L3 > 300 AND U1 < 250 AND U1 > 230 THEN s$ = "Io-C"
IF s$ <> "" THEN GOSUB Outdat
th = th + .5
LOOP UNTIL INT(th / 24) + 1 > dmax
PRINT "Program completed - results in file JOVRAD.TXT"
END
Compute:
d = d0 + th / 24
v = (157.0456 + .0011159# * d) MOD 360
m = (357.2148 + .9856003# * d) MOD 360
n = (94.3455 + .0830853# * d + .33 * SIN(kr * v)) MOD 360
j = (351.4266 + .9025179# * d - .33 * SIN(kr * v)) MOD 360
a = 1.916 * SIN(kr * m) + .02 * SIN(kr * 2 * m)
b = 5.552 * SIN(kr * n) + .167 * SIN(kr * 2 * n)
k = j + a - b
r = 1.00014 - .01672 * COS(kr * m) - .00014 * COS(kr * 2 * m)
re = 5.20867 - .25192 * COS(kr * n) - .0061 * COS(kr * 2 * n)
dt = SQR(re * re + r * r - 2 * re * r * COS(kr * k))
sp = r * SIN(kr * k) / dt
ps = sp / .017452
dl = d - dt / 173
pb = ps - b
xi = 150.4529 * INT(dl) + 870.4529 * (dl - INT(dl))
L3 = (274.319 + pb + xi + .01016 * 51) MOD 360
U1 = 101.5265 + 203. 613