58 lines
2.6 KiB
Python
58 lines
2.6 KiB
Python
import astropy.units as unit
|
|
import numpy as np
|
|
from astropy.coordinates import EarthLocation, AltAz
|
|
from astropy.coordinates import get_sun
|
|
|
|
def sun_angles_astropy(lat, lon, h, utc):
|
|
loc = EarthLocation(lat=lat*unit.deg, lon=lon*unit.deg, height=h*unit.m)
|
|
ref = AltAz(obstime=utc, location=loc)
|
|
|
|
sun_pos = get_sun(utc).transform_to(ref)
|
|
|
|
AZ = sun_pos.az.degree
|
|
ELV = sun_pos.alt.degree
|
|
|
|
return AZ, ELV
|
|
|
|
def sun_angles_analytical(lat, lon, utc):
|
|
JD = utc.jd
|
|
JC = (JD - 2451545) / 36525
|
|
GML = (280.46646 + JC * (36000.76983 + JC * 0.0003032)) % 360
|
|
GMA = 357.52911 + JC * (35999.05029 - 0.0001537 * JC)
|
|
EEO = 0.016708634 - JC * (0.000042037 + 0.0000001267 * JC)
|
|
SEC = np.sin(np.deg2rad(GMA)) * (1.914602 - JC * (0.004817 + 0.000014 * JC)) + np.sin(np.deg2rad(2 * GMA)) * (
|
|
0.019993 - 0.000101 * JC) + np.sin(np.deg2rad(3 * GMA)) * 0.000289
|
|
STL = GML + SEC
|
|
# STA = GMA + SEC
|
|
# SRV = (1.000001018 * (1 - EEO ** 2)) / (1 + EEO * np.cos(np.deg2rad(STA)))
|
|
SAL = STL - 0.00569 - 0.00478 * np.sin(np.deg2rad(125.04 - 1934.136 * JC))
|
|
MOE = 23 + (26 + (21.448 - JC * (46.815 + JC * (0.00059 - JC * 0.001813))) / 60) / 60
|
|
OC = MOE + 0.00256 * np.cos(np.deg2rad(125.04 - 1934.136 * JC))
|
|
# SRA = np.rad2deg(np.arctan2(np.cos(np.deg2rad(OC)) * np.sin(np.deg2rad(SAL)), np.cos(np.deg2rad(SAL)))) # radian
|
|
SD = np.rad2deg(np.arcsin(np.sin(np.deg2rad(OC)) * np.sin(np.deg2rad(SAL)))) # radian
|
|
var_y = np.tan(np.deg2rad(OC / 2)) ** 2
|
|
EOT = 4 * np.rad2deg(
|
|
var_y * np.sin(2 * np.deg2rad(GML)) - 2 * EEO * np.sin(np.deg2rad(GMA)) + 4 * EEO * var_y * np.sin(
|
|
np.deg2rad(GMA)) * np.cos(2 * np.deg2rad(GML)) - 0.5 * var_y ** 2 * np.sin(
|
|
4 * np.deg2rad(GML)) - 1.25 * EEO ** 2 * np.sin(2 * np.deg2rad(GMA)))
|
|
TST = (((JD - 0.5) % 1) * 1440 + EOT + 4 * lon) % 1440
|
|
|
|
if TST / 4 < 0:
|
|
HA = TST / 4 + 180
|
|
else:
|
|
HA = TST / 4 - 180
|
|
|
|
SZA = np.rad2deg(np.arccos(
|
|
np.sin(np.deg2rad(lat)) * np.sin(np.deg2rad(SD)) + np.cos(np.deg2rad(lat)) * np.cos(np.deg2rad(SD)) * np.cos(
|
|
np.deg2rad(HA))))
|
|
SEA = 90 - SZA
|
|
|
|
if HA > 0:
|
|
SAA = (np.rad2deg(np.arccos(((np.sin(np.deg2rad(lat)) * np.cos(np.deg2rad(SZA))) - np.sin(np.deg2rad(SD))) / (
|
|
np.cos(np.deg2rad(lat)) * np.sin(np.deg2rad(SZA))))) + 180) % 360
|
|
else:
|
|
SAA = (540 - np.rad2deg(np.arccos(
|
|
((np.sin(np.deg2rad(lat)) * np.cos(np.deg2rad(SZA))) - np.sin(np.deg2rad(SD))) / (
|
|
np.cos(np.deg2rad(lat)) * np.sin(np.deg2rad(SZA)))))) % 360
|
|
|
|
return SAA, SEA |