asolardistancelo.py

Created on December 07, 2023
2.08 KB
Calculates solar distance (centre Sun-centre Earth) given a JDE. Based on low-precision Jean Meeus algorithm (precision maybe 1000’s of km)
"""Solar distance.

A low-accuracy algorithm.
Calculate solar distance given JDE.
Results in AU km or light seconds.
"""
__version__="0.1"
crid="Astronomy:Solar distance V"+__version__+"""
© 2023 SteveG1CMZ"""

#Jean Meeus low accuracy algorithm

from math import *
#from aeclipse import * #just for example dates

showme=True
#select edition1 or edition2
edition1=None #not implemented
edition2=1

#24=25 Solar Coordinates
#low accuracy (5sf au, 1000's km, 1/10s)

rads=pi/180
aum = 149597870700 #au  to m (exact)
kms = 3.33564e-6 #km to light seconds

def solar_distau(JDE):

  #=25.1 T centuries
  T=(JDE-2451545.0)/36525
  #=25.2 mean long of Sun (degrees)
  L0=(edition2*280.46646) + 36000.76983*T + 0.0003032*T**2
  #=25.3 mean anomaly of Sun (degrees)
  M=edition2*(357.52911 + 35999.05029*T - 0.0001537*T**2)
  #=25.4 eccentricity (Earth orbit)
  ecc=edition2*(0.016708634 - 0.000042037*T -0.0000001267*T**2)
  Mr=M*rads
  
  #Suns equation of centre (degrees)
  C= edition2*( (1.914602 - 0.004817*T - 0.000014*T**2)*sin(Mr) +(0.019993 - 0.000101*T)*sin(2*Mr) +0.000289*sin(3*Mr))
  #Sun true
  Suntruelong=L0+C #degrees ???
  Suntrueanomaly=M+C #degrees
  vr=Suntrueanomaly*rads
  #print("T L0 M ecc C long anom")
  #print(T,L0,M,ecc)
  #print(C,Suntruelong,Suntrueanomaly)

  #25.5 Sun-Earth distance (centre-centre)
  den=(1+ecc*cos(vr))
  R=(1.000001018*(1-ecc**2))/den #au
  return round(R,5) #au
  
def solar_distkm(JDE):
  Rau=solar_distau(JDE)
  return round(Rau*aum/1000) #km
  
def solar_dists(JDE):
  Rkm=solar_distkm(JDE)
  return round(Rkm*kms,1) #light seconds
  
def examples():
  JDE=2448908.5
  print(JDE)
  au=solar_distau(JDE)
  print(au,"au")
  km=solar_distkm(JDE)
  print(km,"km")
  print(solar_dists(JDE),"s")
  #today
  today=2460246 #when I wrote this line
  au=solar_distau(today)
  print(au,"au")
  km=solar_distkm(today)
  print(km,"km")
  print(solar_dists(today),"s")
  
if showme or __name__=="__main__":
  print(crid)
  print("distance Earth-Sun (centre-centre)")
  #print("au: given a JDE")
  #print("see aeclipse to change Date to JDE")
  #25.a
  examples()