1

There were 2 questions showing different solutions: "Pyephem calculate current solar time" and "Local solar time function".

I tried both solutions using a timestamp of a transit. And I thought, that should result in 12:00:00 local solar time. But the results differed:

sun transit                 2021-03-15 12:16:42.157316
                             radians             hours
sun.az,sun.alt to tau,delta   3.1416   1.0000  12.0000
observer_ra-sun.a_ra          3.1416   1.0000  12.0000
observer_ra-sun.g_ra          3.1370   0.9985  11.9825
observer_ra-sun.ra            3.1370   0.9985  11.9825
sidereal_time-sun.ra          3.1416   1.0000  12.0000

The question is:

The code testing those calculation methods and producing the above output is:

sun=ephem.Sun()
hier = ephem.Observer()
hier.lat,hier.lon='51.123','13.040'
hier.elevation=171
ts=1615762800
dt=datetime.datetime.fromtimestamp(ts,datetime.timezone(datetime.timedelta(),'UTC'))
hier.date=dt.astimezone(datetime.timezone(datetime.timedelta(),'UTC'))
sun.compute(hier)
hier.date = hier.next_transit(sun)
sun.compute(hier)
hier_ra,hier_dec = hier.radec_of('0','-90')
print("sun transit                 %s" % ephem.localtime(hier.date))
print("                             radians             hours")
tau,delta=horiz_to_equ(float(hier.lat),float(sun.az),float(sun.alt))
print("sun.az,sun.alt to tau,delta %8.4f %8.4f %8.4f" % (tau,tau/math.pi,tau/math.pi*12))
solartime = ephem.hours((hier_ra-sun.a_ra)%(2*math.pi))
print("observer_ra-sun.a_ra        %8.4f %8.4f %8.4f" % (solartime,solartime/math.pi,solartime/math.pi*12))
solartime = ephem.hours((hier_ra-sun.g_ra)%(2*math.pi))
print("observer_ra-sun.g_ra        %8.4f %8.4f %8.4f" % (solartime,solartime/math.pi,solartime/math.pi*12))
solartime = ephem.hours((hier_ra-sun.ra)%(2*math.pi))
print("observer_ra-sun.ra          %8.4f %8.4f %8.4f" % (solartime,solartime/math.pi,solartime/math.pi*12))
solartime = ephem.hours((hier.sidereal_time()-sun.ra)+ephem.hours('12:00'))
print("sidereal_time-sun.ra        %8.4f %8.4f %8.4f" % (solartime,solartime/math.pi,solartime/math.pi*12))

The function horiz_to_equ does a coordinate transformation from azimuth and altitude to hour angle and declination. It is not the question here and given for information, only.

def horiz_to_equ(lat,a,h):
    a-=math.pi
    # Hilfswerte
    p34=math.pi*3.0/4.0
    p14=math.pi/4.0
    p38=math.pi*3.0/8.0
    p58=math.pi*5.0/8.0
    # Deklination
    delta=math.asin(math.sin(lat)*math.sin(h)-math.cos(lat)*math.cos(h)*math.cos(a))
    # Stundenwinkel
    uu=math.sin(a)
    vv=math.sin(lat)*math.cos(a)+math.cos(lat)*math.tan(h)
    if a>-p58 and a<-p38:
        uu*=2
        vv*=2
    tau=math.atan2(uu,vv)
    return tau+math.pi,delta

Thanks in advance.

roe
  • 11
  • 2

0 Answers0