Necesito ayuda con el programa de matemáticas en Python para marcar las emisiones de radio jovianas

Hace unos días, decidí asumir el pequeño proyecto de convertir un programa Qbasic en Python (como un proyecto paralelo a mi proyecto Radio Jove), y logré hacerlo funcionar, pero las matemáticas definitivamente no funcionan. Esperaba que un colega astrónomo o astrofísico me ayudara en el aspecto matemático de mi código Python. Obviamente, algo anda mal con el aspecto matemático del código, ya que la mayoría de los valores en el programa Python están muy sesgados. Mis disculpas de antemano por las variables de un solo carácter. Los haría mejores si supiera todo lo que sucede.

Código Qbasic:

'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.405863# * dl + pb
 U2 = 67.81114 + 101.291632# * dl + pb
 z = (2 * (U1 - U2)) MOD 360
 U1 = U1 + .472 * SIN(kr * z)
 U1 = (U1 + 180) MOD 360
RETURN

Outdat:
 dy = INT(th / 24) + 1
 h = th - (dy - 1) * 24
 IF dy > ty THEN
   m = INT((dy - tx) / 30.6) + 3
   da = dy - ty - INT((m - 3) * 30.6 + .5)
 ELSE
   m = INT((dy - 1) / 31) + 1
   da = dy - (m - 1) * 31
 END IF
 mn$ = MID$(m$, (m - 1) * 3 + 1, 3)
     PRINT #1, USING f$; dy; mn$; da; h; U1; L3; dt; s$
RETURN

Salida de muestra:

******************************************************
  JOVIAN IO-DECAMETRIC EMISSION PREDICTIONS FOR 1994 
******************************************************

Day   Date   Hr(UT)  Io_Phase   CML   Dist(AU)  Source

  4  JAN  4   11.0     228      205     5.81     Io-A
  4  JAN  4   11.5     233      224     5.81     Io-A
  4  JAN  4   12.0     237      242     5.81     Io-A
  6  JAN  6    6.0     233      325     5.78     Io-C
  6  JAN  6    6.5     237      343     5.78     Io-C
  7  JAN  7    6.5      81      133     5.76     Io-B
  7  JAN  7    7.0      86      152     5.76     Io-B
  7  JAN  7    7.5      90      170     5.76     Io-B
  9  JAN  9   20.0     241      203     5.73     Io-A
  9  JAN  9   20.5     246      222     5.73     Io-A
 11  JAN 11   12.5     225      232     5.70     Io-A
 11  JAN 11   13.0     229      251     5.70     Io-A
 11  JAN 11   14.5     242      305     5.70     Io-C
 11  JAN 11   15.0     246      323     5.70     Io-C
 12  JAN 12   15.0      90      113     5.68     Io-B
 12  JAN 12   15.5      94      131     5.68     Io-B
 12  JAN 12   16.0      98      150     5.68     Io-B
 14  JAN 14    8.5      82      179     5.66     Io-B
 16  JAN 16   21.0     234      214     5.61     Io-A
 16  JAN 16   21.5     238      232     5.61     Io-A
 16  JAN 16   22.0     242      250     5.61     Io-A
 18  JAN 18   15.5     234      315     5.59     Io-C
 18  JAN 18   16.0     238      333     5.59     Io-C
 19  JAN 19   16.0      82      124     5.57     Io-B
 19  JAN 19   16.5      86      142     5.57     Io-B
 19  JAN 19   17.0      91      160     5.57     Io-B
 19  JAN 19   17.5      95      178     5.57     Io-B

código pitón:

#!/usr/bin/env python
from __future__ import print_function, division
import math
print('Program to flag Jovian Decametric windows')
month = ['JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG','SEP','OCT','NOV','DEC']
week = 42.46/360
pi = 3.141593
kr = pi / 180
form = "###  \ \ ##   ##.#     ###      ###    ##.##     \  \\"
num1 = open('jovrad.txt', 'w')
yy = int(raw_input(("Year for which predictions are required ")))
e = math.trunc(((yy-1) / 100))
f = 2 - e + math.trunc(e/4)
jd = math.trunc(365.25 * (yy - 1)) + 1721423 + f + .5
d0 = jd - 2435108
incr = 0
dmax = 0
tx = 0
ty = 0
yyly = 0
yylc = 0
if yy / 400 - math.trunc((yy / 400)) == 0:
    incr = 1
    yyly = yy / 4 - math.trunc((yy / 4))
    yylc = yy / 100 - math.trunc((yy / 100))
if yyly == 0 and yylc != 0:
    incr = 1
    ty = 59 + incr
    dmax = 365 + incr
    tx = ty + .5
print("******************************************************", file=num1)
print("  JOVIAN IO-DECAMETRIC EMISSION PREDICTIONS FOR ",yy, file=num1)
print("******************************************************", file=num1)
print("\n", file=num1)
print("Day   Date   Hr(UT)  Io_Phase   CML   Dist(AU)  Source", file=num1)
print("\n", file=num1)
th = 0
def compute(d0, th, dmax):
    global U1, L3, dt, s
    d = d0 + math.trunc(th / 24)
    v = (157.0456 + .0011159 * d) % 360
    m = (357.2148 + .9856003 * d) % 360
    n = (94.3455 + .0830853 * d + .33 * math.sin(kr * v)) % 360
    j = (351.4266 + .9025179 * d - .33 * math.sin(kr * v)) % 360
    a = 1.916 * math.sin(kr * m) + .02 * math.sin(kr * 2 * m)
    b = 5.552 * math.sin(kr * n) + .167 * math.sin(kr * 2 * n)
    k = j + a - b
    r = 1.00014 - .01672 * math.cos(kr * m) - .00014 * math.cos(kr * 2 * m)
    re = 5.20867 - .25192 * math.cos(kr * n) - .0061 * math.cos(kr * 2 * n)
    dt = math.sqrt(re * re + r * r - 2 * re * r * math.cos(kr * k))
    sp = r * math.sin(kr * k) / dt
    ps = sp / .017452
    dl = d - dt / 173
    pb = ps - b
    xi = 150.4529 * math.trunc((dl)) + 870.4529 * (dl - math.trunc((dl)))
    L3 = (274.319 + pb + xi + .01016 * 51) % 360
    U1 = 101.5265 + 203.405863 * dl + pb
    U2 = 67.81114 + 101.291632 * dl + pb
    z = (2 * (U1 - U2)) % 360
    U1 = U1 + .472 * math.sin(kr * z)
    U1 = (U1 + 180) % 360
    s = ""
    while math.trunc(th / 24) + 1 < dmax and math.trunc(th / 24) + 1 == dmax:
        if L3 < 255 and L3 > 200 and U1 < 250 and U1 > 220:
            s = "Io-A"
        if L3 < 180 and L3 > 105 and U1 < 100 and U1 > 80:
            s = "Io-B"
        if L3 < 350 and L3 > 300 and U1 < 250 and U1 > 230:
            s = "Io-C"
        if s != "":
            outdat()
        th = th + .5
    print("Program completed - results in file JOVRAD.TXT", file=num1)
compute(d0, th, dmax)

def outdat(th,tx,ty):
    dy = math.trunc((th / 24)) + 1
    h = th - (dy - 1) * 24
    if(dy > th):
        m = math.trunc((dy - tx) / 30.6) + 3
        da = dy - ty - math.trunc((m - 3) * 30.6 + .5)
    else:
        m = math.trunc((dy - 1) / 31) + 1
        da = dy - (m - 1) * 31
    mn = month[(m-1)*3+1:(m-1)*3-1+3]
    print(dy, mn, da, h, U1, L3, dt, s, file=num1)
outdat(th,tx,ty)

Mi salida actual:

******************************************************
  JOVIAN IO-DECAMETRIC EMISSION PREDICTIONS FOR  1998
******************************************************


Day   Date   Hr(UT)  Io_Phase   CML   Dist(AU)  Source


-12.9775 1.0 ['AUG'] -0.5 0.0 135.325782651 10.6054462674 5.7071322535
Es bueno que hayas elegido Python. ¿Has pensado en usar numpy? Entonces tendrás todas las funciones matemáticas, como np.pi, np.sin()y mucho más.
Sí, realmente disfruto de Python. Tengo numpy instalado y scipy, pero no estoy seguro de si necesitaré usarlo. La biblioteca matemática de Python ya tiene una función pi (que debería usar en este programa) y una función de seno y coseno que ya uso.
Ok, pero cual es tu pregunta entonces? ¿Qué está mal con tus matemáticas? me parece bien
Bueno, teniendo en cuenta que obtuve un número negativo para el "Día" y que todos los demás valores son muy diferentes de lo que deberían ser, algo anda mal.
¡Maldad mía, no vi la salida!
@Aabaakawad Sí, el segundo bloque de código lo muestra.
@Macuser lo siento, eliminaré mi comentario.
Tiene una conversión explícita de resultados (que son flotantes) como enteros en el básico pero no en Python. Me preocuparía por las diferencias que resultan de esto. En Python, ¿qué tipos son las variables y cómo se realiza cualquier conversión implícita (si corresponde) a ints?
@Contrad Turner es correcto. Tomemos, por ejemplo, la línea que asigna algo a la variable mi . En el código QBasic, está utilizando el yo norte T ( ) función que trunca la parte decimal de los números y devuelve el número entero. Para hacer eso en python, podrías usar el metro a t h . t r tu norte C ( ) función. Mirando tu código, creo que el yo norte T es lo único que debes cuidar. creo que el METRO O D funciona igual que el módulo de python, si no recuerdo mal.
¡Gran! Actualicé el código de python para tener todos los Ints donde sea necesario. Lo ejecuté de nuevo y parece que sigo obteniendo un número negativo para el valor del día. Pero estoy seguro de que solucionó algunos problemas que estaban presentes antes.
Acabo de actualizar el código aún más. Accidentalmente estaba imprimiendo "f" cuando se suponía que "dy" estaba impreso como lo primero. Ahora solo necesito corregir el formato para poder usar mi variable de "formulario" y luego verificar los valores de salida.
Los números definitivamente todavía no están calculando bien. Solo obtengo una línea de datos para cualquier año que ingrese. Además, la fuente sigue sin aparecer.

Respuestas (1)

Hubo algunos otros pequeños errores en su traducción (como sangrar demasiadas cosas en dos de las ifdeclaraciones), y creo que podría haber estado tratando de ser demasiado inteligente al cambiar SUBROUTINEa funciones definidas cuando era mucho más fácil simplemente eliminar el código en línea ya que solo se usó una vez.

Hubo ajustes menores y, por lo que sé, esto parece funcionar (querrá jugar con la impresión formateada para limpiarla, y esto se ejecuta en Python 3; si necesita que se ejecute en Python 2, usted ' Probablemente necesite cambiar las declaraciones printy num1.write):

#!/usr/bin/env python

import math

print('Program to flag Jovian Decametric windows')
month = ['JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG','SEP','OCT','NOV','DEC']
week = 42.46/360
pi = 3.141593
kr = pi / 180
num1 = open('jovrad.txt', 'w')
yy = int(input(("Year for which predictions are required ")))
e = math.trunc(((yy-1) / 100))
print(e)
f = 2 - e + math.trunc(e/4)
print(f)
jd = math.trunc(365.25 * (yy - 1)) + 1721423 + f + .5
print(jd)
d0 = jd - 2435108
print(d0)
incr = 0
dmax = 0
tx = 0
ty = 0
yyly = 0
yylc = 0
if yy / 400 - math.trunc((yy / 400)) == 0:
    incr = 1
    print("in if-1")
yyly = yy / 4 - math.trunc((yy / 4))
yylc = yy / 100 - math.trunc((yy / 100))
if yyly == 0 and yylc != 0:
    print("in if-2")
    incr = 1
ty = 59 + incr
dmax = 365 + incr
tx = ty + .5
num1.write("******************************************************\n")
pout = "  JOVIAN IO-DECAMETRIC EMISSION PREDICTIONS FOR " + str(yy) + "\n"
num1.write(pout)
num1.write("******************************************************\n")
num1.write("\n")
num1.write("Day   Date   Hr(UT)  Io_Phase   CML   Dist(AU)  Source")
num1.write("\n")
th = 0

while int(th / 24) + 1 <= dmax:
    d = d0 + (th / 24)
    v = (157.0456 + .0011159 * d) % 360
    m = (357.2148 + .9856003 * d) % 360
    n = (94.3455 + .0830853 * d + .33 * math.sin(kr * v)) % 360
    j = (351.4266 + .9025179 * d - .33 * math.sin(kr * v)) % 360
    a = 1.916 * math.sin(kr * m) + .02 * math.sin(kr * 2 * m)
    b = 5.552 * math.sin(kr * n) + .167 * math.sin(kr * 2 * n)
    k = j + a - b
    r = 1.00014 - .01672 * math.cos(kr * m) - .00014 * math.cos(kr * 2 * m)
    re = 5.20867 - .25192 * math.cos(kr * n) - .0061 * math.cos(kr * 2 * n)
    dt = math.sqrt(re * re + r * r - 2 * re * r * math.cos(kr * k))
    sp = r * math.sin(kr * k) / dt
    ps = sp / .017452
    dl = d - dt / 173
    pb = ps - b
    xi = 150.4529 * math.trunc((dl)) + 870.4529 * (dl - math.trunc((dl)))
    L3 = (274.319 + pb + xi + .01016 * 51) % 360
    U1 = 101.5265 + 203.405863 * dl + pb
    U2 = 67.81114 + 101.291632 * dl + pb
    z = (2 * (U1 - U2)) % 360
    U1 = U1 + .472 * math.sin(kr * z)
    U1 = (U1 + 180) % 360
    s = ""
    if L3 < 255 and L3 > 200 and U1 < 250 and U1 > 220:
        s = "Io-A"
    if L3 < 180 and L3 > 105 and U1 < 100 and U1 > 80:
        s = "Io-B"
    if L3 < 350 and L3 > 300 and U1 < 250 and U1 > 230:
        s = "Io-C"
    if s != "":
        dy = math.trunc((th / 24)) + 1
        h = th - (dy - 1) * 24
        if(dy > th):
            m = math.trunc((dy - tx) / 30.6) + 3
            da = dy - ty - math.trunc((m - 3) * 30.6 + .5)
        else:
            m = math.trunc((dy - 1) / 31) + 1
            da = dy - (m - 1) * 31
        mn = month[(m-1)]
#        mn = month[(m-1)*3+1:(m-1)*3-1+3]
        outstring = "%i  %s  %i  %5.3f  %5.3f  %5.3f  %5.3f  %s\n" % (dy, mn, da, h, U1, L3, dt, s)
        num1.write(outstring)
    th = th + .5

num1.close()
print("Program Complete  - results in file JOVRAD.TXT")
¡Muchas gracias! Realmente lo aprecio. Por cierto, esto funciona al 100% en Python 2. Y sí, tendré que trabajar un poco en el formato de impresión.
Además, eliminé la primera variable en la "cadena exterior" que escribiste. Simplemente comienza con el mes y luego el día (variable "da").
Necesito ayuda con el programa de matemáticas en Python para marcar las emisiones de radio jovianas
RespuestasAquíPolítica de privacidad1y, 0v