Code:
## -------------------------- ORIGINAL LICENSE -------------------
##
##This program is free software: you can redistribute it and/or modify
##it under the terms of the GNU General Public License as published by the Free Software
##Foundation, either version 3 of the license or any later version.<br><br>
##This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
##without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
##See the GNU General Public License
##
##Check this web site for updated versions:
## https://www.celnav.de/index.htm</a></p>
# by rmax 2018
# Ported into python from Javascript
# Original is from from https://www.celnav.de/
import math
import sys
dtr = math.pi / 180;
# set-up global variables. Kludgy but will work for the moment
#inputs
year = 0;
month = 0;
day =0;
hour =0;
minute =0;
second =0;
#time
JD=0;
JD0h = 0;
TE=0;
TE2=0;
TE3=0;
TE4=0;
TE5=0;
Tau=0;
Tau2 = 0;
Tau3 = 0;
Tau4 = 0;
Tau5 = 0;
T2=0;
T3=0;
T4=0;
T5=0;
# Aries
delta_psi=0;
eps =0;
GHAAtrue=0;
# Sun
Dsun = 0;
RAsun = 0;
lmbda = 0;
Lsun_true = 0;
# moon
lambdaMapp = 0;
k = 0;
DECmoon = 0;
RAmoon =0;
# dataOutputList
dataOutput = [];
# Main function
def main():
ReadData();
if (ValidateData()):
print(str(day) + "/" + str(month)+ "/" +str(year) + " " + str(hour) +":"+str(minute) + ":" + str(second) + ",")
TimeMeasures();
Nutation();
Aries();
Sun();
Moon();
MoonPhase();
Weekday();
Polaris();
# star loop
x=0;
while (x<58):
Star(x);
x+=1;
# Auxiliary functions
# Sine of angles in degrees
def sind(x):
return math.sin(dtr*x);
# Cosine of angles in degrees
def cosd(x):
return math.cos(dtr*x);
#Tangent of angles in degrees
def tand(x):
return math.tan(dtr*x);
# Normalize large angles
def norm_360_deg(x):
while (x<0):
x=x+360;
return x%360;
# Output Hour Angle
def OutHA(x):
GHAdeg = math.floor(x);
GHAmin = math.floor(60*(x-GHAdeg));
GHAsec = round(3600*(x-GHAdeg-GHAmin/60));
if (GHAsec==60):
GHAsec=0;
GHAmin+=1;
if (GHAmin==60):
GHAmin=0;
GHAdeg+=1;
if (GHAdeg==0):
GHAdeg="000";
elif (GHAdeg <10):
GHAdeg="00" + str(GHAdeg);
elif (GHAdeg <100):
GHAdeg="0" + str(GHAdeg);
if (GHAmin==0):
GHAmin="00";
elif (GHAmin<10):
GHAmin="0" + str(GHAmin);
if (GHAsec<10):
GHAsec="0" + str(GHAsec);
printtext = " " + str(GHAdeg) + u'\xB0' + " " + str(GHAmin) + "'" + " " + str(GHAsec) + "''";
return printtext.encode('utf-8');
#Output Mean Obliqity of Ecliptic
def OutECL(x):
ECLdeg = math.floor(x);
ECLmin = math.floor(60*(x-ECLdeg));
ECLsec = round(3600000*(x-ECLdeg-ECLmin/60))/1000;
if (ECLsec==60):
ECLsec=0;
ECLmin+=1;
if (ECLmin==60):
ECLmin=0;
ECLdeg+=1;
if (ECLmin==0):
ECLmin="00";
elif (ECLmin<10):
ECLmin="0" + str(ECLmin);
if (ECLsec<10):
ECLsec="0" + str(ECLsec);
#kludge didn't want to return unicode
printtext = " " + str(ECLdeg) + u'\xB0' + " " + str(ECLmin) + "'" + " " + str(ECLsec) + "''";
return printtext.encode('utf-8');
# Output Sidereal Time
def OutSidTime(x):
GMSTdecimal = x/15;
GMSTh = math.floor(GMSTdecimal);
GMSTmdecimal = 60*(GMSTdecimal-GMSTh);
GMSTm = math.floor(GMSTmdecimal);
GMSTsdecimal = 60*(GMSTmdecimal-GMSTm);
GMSTs = round(1000*GMSTsdecimal)/1000;
GMSTs_string = str(GMSTs)
if (GMSTs-math.floor(GMSTs)==0):
GMSTs_string = GMSTs_string + ".000";
elif (10*GMSTs-math.floor(10*GMSTs)==0):
GMSTs_string = GMSTs_string + "00";
elif (100*GMSTs-math.floor(100*GMSTs)==0):
GMSTs_string = GMSTs_string + "0";
return " " + str(GMSTh) + "h" + " "+ str(GMSTm) + "m" + " " + GMSTs_string + "s";
# Output Declination
def OutDec(x):
if (x<0):
signDEC=-1;
name="S";
else:
signDEC=1;
name="N";
DEC = abs(x);
DECdeg = math.floor(DEC);
DECmin = math.floor(60*(DEC-DECdeg));
DECsec = round(3600*(DEC-DECdeg-DECmin/60));
if (DECsec==60):
DECsec=0;
DECmin+=1;
if (DECmin==60):
DECmin=0;
DECdeg+=1;
if (DECdeg==0):
DECdeg="00";
elif (DECdeg<10):
DECdeg="0" + str(DECdeg);
if (DECmin==0):
DECmin="00";
elif (DECmin<10):
DECmin="0" + str(DECmin);
if (DECsec<10):
DECsec="0"+ str(DECsec);
printtext = " " + name+" "+ str(DECdeg) + u'\xB0' + " "+ str(DECmin) + "'" + " " + str(DECsec) + "''";
return printtext.encode('utf-8');
# Output SD and HP
def OutSDHP(x):
x = round(10*x)/10;
printx = str(x);
if (x-math.floor(x)== 0):
printx = printx + ".0";
return " " + printx + "''";
# Input data conversion
def ReadData():
global year;
global month;
global day;
global hour;
global minute;
global second;
startDate = "";
startTime = "";
startDate = input("Enter Date for GHA (dd/mm/yyyy > ");
startTime = input("Enter time for GHA (hh:mm:ss> ");
if (startDate == ""):
startDate = "31/07/1973";
# divide the dateformat
# startDate = startDate.strip();
print ("Entered " + startDate);
delimChar="";
if (startDate.find('/')>0):
delimChar='/';
if (startDate.find('-')>0):
delimChar='-';
if (startDate.find('?')>0):
delimChar='?';
if (startDate.find(' ')>0):
delimChar=' ';
startDateSplit=startDate.split(delimChar);
day = eval(startDateSplit[0]);
month = eval(startDateSplit[1]);
year = eval(startDateSplit[2]);
if (startTime == ""):
startTime="00:00:00";
delimChar="";
if (startTime.find(':')>0):
delimChar=':';
if (startTime.find('-')>0):
delimChar='-';
if (startTime.find(' ')>0):
delimChar=' ';
startTimeSplit=startTime.split(delimChar);
hour = eval(startTimeSplit[0]);
minute = eval(startTimeSplit[1]);
second = eval(startTimeSplit[2]);
# Extract data conversion
def getData(date, inputTime):
global year;
global month;
global day;
global hour;
global minute;
global second;
delimChar="/";
startDateSplit=date.split(delimChar);
if(startDateSplit[0][:1]=="0"):
day = eval(startDateSplit[0][1:2]);
else:
day = eval(startDateSplit[0]);
if(startDateSplit[1][:1]=="0"):
month= eval(startDateSplit[1][1:2]);
else:
month = eval(startDateSplit[1]);
year = eval(startDateSplit[2]);
delimChar=":";
startTimeSplit = inputTime.split(delimChar);
print (startTimeSplit[0])
if(startTimeSplit[0][:1]=="0"):
hour= eval(startTimeSplit[0][1:2]);
else:
hour = eval(startTimeSplit[0]);
if(startTimeSplit[1][:1]=="0"):
minute = eval(startTimeSplit[1][1:2]);
else:
minute = eval(startTimeSplit[1]);
if(startTimeSplit[2][:1]=="0"):
second = eval(startTimeSplit[2][1:2]);
else:
second = eval(startTimeSplit[2]);
def ValidateData():
## if(InForm.year.value == "") alert("Missing year! Restart calculation.");
## else year = eval(InForm.year.value);
## if(InForm.month.value == "") alert("Missing month! Restart calculation.");
## else month = eval(InForm.month.value);
## if(InForm.day.value == "") alert("Missing day! Restart calculation.");
## else day = eval(InForm.day.value);
validData = True;
if(month < 1 | month > 12):
print("Month out of range! Restart calculation.");
validData = False;
if(day < 1 | day > 31):
print("Day out of range! Restart calculation.");
validData = False;
if(year/4- math.floor(year/4) == 0):
schj=1;
if(year/100 - math.floor(year/100) == 0):
schj=0;
if(year/400 - math.floor(year/400) == 0):
schj=1;
if(month == 2 & day > 28 & schj == 0):
print("February has only 28 days! Restart calculation.");
validData = False;
if(month == 2 & day > 29 & schj == 1):
print("February has only 29 days in a leap year! Restart calculation.");
validData = False;
if(month == 4 & day > 30):
validData = False;
print("April has only 30 days! Restart calculation.");
validData = False;
if(month == 6 & day > 30):
print("June has only 30 days! Restart calculation.");
validData = False;
if(month == 9 & day > 30):
print("September has only 30 days! Restart calculation.");
validData = False;
if(month == 11 & day > 30):
print("November has only 30 days! Restart calculation.");
validData = False;
## if(InForm.hour.value == "") var hour = 0;
## else var hour = eval(InForm.hour.value);
## if(InForm.minute.value == "") var minute = 0;
## else var minute = eval(InForm.minute.value);
## if(InForm.second.value == "") var second = 0;
## else var second = eval(InForm.second.value);
dayfraction = (hour + minute/60 + second/3600)/24;
if(dayfraction < 0 or dayfraction > 1):
print("Time out of range! Restart calculation.");
validData = False;
## if(InForm.delta.value == "") deltaT = 0;
## else deltaT = eval(InForm.delta.value);
return validData;
# Calculating Julian date, century, and millennium
def TimeMeasures():
global year;
global month;
global day;
global hour;
global minute;
global second;
##### note set DeltaT as 0 for the moment #####
# Julian day (UT1)
if ( month <= 2):
year -=1;
month += 12;
dayfraction =0 ;
dayfraction = (float(hour) + minute/60 + second/3600)/24;
A = math.floor(year/100);
B = 2 - A + math.floor(A/4);
global JD0h;
JD0h = math.floor(365.25*(year+4716)) + math.floor(30.6001*(month+1))+day+B-1524.5;
global JD;
JD = JD0h + dayfraction;
print ("JD " + str(JD) +",")
# Julian centuries (GMT) since 2000 January 0.5
global T2;
global T3;
global T4;
global T5;
T = (JD - 2451545)/36525;
T2 = T*T;
T3 = T*T2;
T4 = T*T3;
T5 = T*T4;
#### Fixed delta T for the moment #####
deltaT = 0;
# Julian ephemeris day (TDT)
JDE = JD + deltaT / 86400;
# Julian centuries (TDT) from 2000 January 0.5
global TE;
global TE2;
global TE3;
global TE4;
global TE5;
TE = (JDE - 2451545) / 36525;
TE2 = TE*TE;
TE3 = TE*TE2;
TE4 = TE*TE3;
TE5 = TE*TE4;
# Julian millenniums (TDT) from 2000 January 0.5
global Tau;
global Tau2;
global Tau3;
global Tau4;
global Tau5;
Tau = 0.1 * TE;
Tau2 = Tau * Tau;
Tau3 = Tau * Tau2;
Tau4 = Tau * Tau3;
Tau5 = Tau * Tau4;
# Nutation, obliquity of the ecliptic
def Nutation():
# IAU 1980 nutation theory:
# Mean anomaly of the moon
Mm = 134.962981389+198.867398056 * TE + norm_360_deg(477000 * TE) + 0.008697222222*TE2 + TE3/56250;
# Mean anomaly of the sun
M = 357.527723333+359.05034*TE+norm_360_deg(35640*TE)-0.0001602777778*TE2-TE3/300000;
# Mean distance of the moon from the ascending node
F = 93.271910277+82.017538055*TE+norm_360_deg(483120*TE)-0.0036825*TE2+TE3/327272.7273;
#Mean elongation of the moon
D = 297.850363055+307.11148*TE+norm_360_deg(444960*TE)-0.001914166667*TE2+TE3/189473.6842;
# Longitude of the ascending node of the moon
omega = 125.044522222 - 134.136260833 * TE - norm_360_deg(1800*TE)+0.002070833333*TE2+TE3/450000;
# Periodic terms for nutation
nut = [];
nut.append(" 0 0 0 0 1-171996-174.2 92025 8.9 ");
nut.append(" 0 0 2-2 2 -13187 -1.6 5736-3.1 ")
nut.append(" 0 0 2 0 2 -2274 -0.2 977-0.5 ")
nut.append(" 0 0 0 0 2 2062 0.2 -895 0.5 ")
nut.append(" 0-1 0 0 0 -1426 3.4 54-0.1 ")
nut.append(" 1 0 0 0 0 712 0.1 -7 0.0 ")
nut.append(" 0 1 2-2 2 -517 1.2 224-0.6 ")
nut.append(" 0 0 2 0 1 -386 -0.4 200 0.0 ")
nut.append(" 1 0 2 0 2 -301 0.0 129-0.1 ")
nut.append(" 0-1 2-2 2 217 -0.5 -95 0.3 ")
nut.append("-1 0 0 2 0 158 0.0 -1 0.0 ")
nut.append(" 0 0 2-2 1 129 0.1 -70 0.0 ")
nut.append("-1 0 2 0 2 123 0.0 -53 0.0 ")
nut.append(" 1 0 0 0 1 63 0.1 -33 0.0 ")
nut.append(" 0 0 0 2 0 63 0.0 -2 0.0 ")
nut.append("-1 0 2 2 2 -59 0.0 26 0.0 ")
nut.append("-1 0 0 0 1 -58 -0.1 32 0.0 ")
nut.append(" 1 0 2 0 1 -51 0.0 27 0.0 ")
nut.append("-2 0 0 2 0 -48 0.0 1 0.0 ")
nut.append("-2 0 2 0 1 46 0.0 -24 0.0 ")
nut.append(" 0 0 2 2 2 -38 0.0 16 0.0 ")
nut.append(" 2 0 2 0 2 -31 0.0 13 0.0 ")
nut.append(" 2 0 0 0 0 29 0.0 -1 0.0 ")
nut.append(" 1 0 2-2 2 29 0.0 -12 0.0 ")
nut.append(" 0 0 2 0 0 26 0.0 -1 0.0 ")
nut.append(" 0 0 2-2 0 -22 0.0 0 0.0 ")
nut.append("-1 0 2 0 1 21 0.0 -10 0.0 ")
nut.append(" 0 2 0 0 0 17 -0.1 0 0.0 ")
nut.append(" 0 2 2-2 2 -16 0.1 7 0.0 ")
nut.append("-1 0 0 2 1 16 0.0 -8 0.0 ")
nut.append(" 0 1 0 0 1 -15 0.0 9 0.0 ")
nut.append(" 1 0 0-2 1 -13 0.0 7 0.0 ")
nut.append(" 0-1 0 0 1 -12 0.0 6 0.0 ")
nut.append(" 2 0-2 0 0 11 0.0 0 0.0 ")
nut.append("-1 0 2 2 1 -10 0.0 5 0.0 ")
nut.append(" 1 0 2 2 2 -8 0.0 3 0.0 ")
nut.append(" 0-1 2 0 2 -7 0.0 3 0.0 ")
nut.append(" 0 0 2 2 1 -7 0.0 3 0.0 ")
nut.append(" 1 1 0-2 0 -7 0.0 0 0.0 ")
nut.append(" 0 1 2 0 2 7 0.0 -3 0.0 ")
nut.append("-2 0 0 2 1 -6 0.0 3 0.0 ")
nut.append(" 0 0 0 2 1 -6 0.0 3 0.0 ")
nut.append(" 2 0 2-2 2 6 0.0 -3 0.0 ")
nut.append(" 1 0 0 2 0 6 0.0 0 0.0 ")
nut.append(" 1 0 2-2 1 6 0.0 -3 0.0 ")
nut.append(" 0 0 0-2 1 -5 0.0 3 0.0 ")
nut.append(" 0-1 2-2 1 -5 0.0 3 0.0 ")
nut.append(" 2 0 2 0 1 -5 0.0 3 0.0 ")
nut.append(" 1-1 0 0 0 5 0.0 0 0.0 ")
nut.append(" 1 0 0-1 0 -4 0.0 0 0.0 ")
nut.append(" 0 0 0 1 0 -4 0.0 0 0.0 ")
nut.append(" 0 1 0-2 0 -4 0.0 0 0.0 ")
nut.append(" 1 0-2 0 0 4 0.0 0 0.0 ")
nut.append(" 2 0 0-2 1 4 0.0 -2 0.0 ")
nut.append(" 0 1 2-2 1 4 0.0 -2 0.0 ")
nut.append(" 1 1 0 0 0 -3 0.0 0 0.0 ")
nut.append(" 1-1 0-1 0 -3 0.0 0 0.0 ")
nut.append("-1-1 2 2 2 -3 0.0 1 0.0 ")
nut.append(" 0-1 2 2 2 -3 0.0 1 0.0 ")
nut.append(" 1-1 2 0 2 -3 0.0 1 0.0 ")
nut.append(" 3 0 2 0 2 -3 0.0 1 0.0 ")
nut.append("-2 0 2 0 2 -3 0.0 1 0.0 ")
nut.append(" 1 0 2 0 0 3 0.0 0 0.0 ")
nut.append("-1 0 2 4 2 -2 0.0 1 0.0 ")
nut.append(" 1 0 0 0 2 -2 0.0 1 0.0 ")
nut.append("-1 0 2-2 1 -2 0.0 1 0.0 ")
nut.append(" 0-2 2-2 1 -2 0.0 1 0.0 ")
nut.append("-2 0 0 0 1 -2 0.0 1 0.0 ")
nut.append(" 2 0 0 0 1 2 0.0 -1 0.0 ")
nut.append(" 3 0 0 0 0 2 0.0 0 0.0 ")
nut.append(" 1 1 2 0 2 2 0.0 -1 0.0 ")
nut.append(" 0 0 2 1 2 2 0.0 -1 0.0 ")
nut.append(" 1 0 0 2 1 -1 0.0 0 0.0 ")
nut.append(" 1 0 2 2 1 -1 0.0 1 0.0 ")
nut.append(" 1 1 0-2 1 -1 0.0 0 0.0 ")
nut.append(" 0 1 0 2 0 -1 0.0 0 0.0 ")
nut.append(" 0 1 2-2 0 -1 0.0 0 0.0 ")
nut.append(" 0 1-2 2 0 -1 0.0 0 0.0 ")
nut.append(" 1 0-2 2 0 -1 0.0 0 0.0 ")
nut.append(" 1 0-2-2 0 -1 0.0 0 0.0 ")
nut.append(" 1 0 2-2 0 -1 0.0 0 0.0 ")
nut.append(" 1 0 0-4 0 -1 0.0 0 0.0 ")
nut.append(" 2 0 0-4 0 -1 0.0 0 0.0 ")
nut.append(" 0 0 2 4 2 -1 0.0 0 0.0 ")
nut.append(" 0 0 2-1 2 -1 0.0 0 0.0 ")
nut.append("-2 0 2 4 2 -1 0.0 1 0.0 ")
nut.append(" 2 0 2 2 2 -1 0.0 0 0.0 ")
nut.append(" 0-1 2 0 1 -1 0.0 0 0.0 ")
nut.append(" 0 0-2 0 1 -1 0.0 0 0.0 ")
nut.append(" 0 0 4-2 2 1 0.0 0 0.0 ")
nut.append(" 0 1 0 0 2 1 0.0 0 0.0 ")
nut.append(" 1 1 2-2 2 1 0.0 -1 0.0 ")
nut.append(" 3 0 2-2 2 1 0.0 0 0.0 ")
nut.append("-2 0 2 2 2 1 0.0 -1 0.0 ")
nut.append("-1 0 0 0 2 1 0.0 -1 0.0 ")
nut.append(" 0 0-2 2 1 1 0.0 0 0.0 ")
nut.append(" 0 1 2 0 1 1 0.0 0 0.0 ")
nut.append("-1 0 4 0 2 1 0.0 0 0.0 ")
nut.append(" 2 1 0-2 0 1 0.0 0 0.0 ")
nut.append(" 2 0 0 2 0 1 0.0 0 0.0 ")
nut.append(" 2 0 2-2 1 1 0.0 -1 0.0 ")
nut.append(" 2 0-2 0 1 1 0.0 0 0.0 ")
nut.append(" 1-1 0-2 0 1 0.0 0 0.0 ")
nut.append("-1 0 0 1 1 1 0.0 0 0.0 ")
nut.append("-1-1 0 2 1 1 0.0 0 0.0 ")
nut.append(" 0 1 0 1 0 1 0.0 0 0.0 ")
# Reading periodic terms
fMm=0;
fM=0;
fF=0;
fD=0;
f_omega=0;
dp=0.0;
de=0.0;
x=0;
while(x<=105):
tempstring = nut[x];
fMm = eval(tempstring[0:2]);
fM = eval(tempstring[2:4]);
fF = eval(tempstring[4:6]);
fD = eval(tempstring[6:8]);
f_omega = eval(tempstring[8:10]);
dp = dp + (eval(tempstring[10:17])+ TE * eval(tempstring[17:23]))* sind(fD*D + fM*M + fMm*Mm + fF*F + f_omega*omega);
de = de + (eval(tempstring[23:29])+ TE * eval(tempstring[29:33]))* cosd(fD*D + fM*M + fMm*Mm + fF*F + f_omega*omega);
x=x+1;
# Corrections (Herring, 1987) - not implemented
## /*
## var corr = new Array(4);
## corr[0] = " 0 0 0 0 1-725 417 213 224 ";
## corr[1] = " 0 1 0 0 0 523 61 208 -24 ";
## corr[2] = " 0 0 2-2 2 102-118 -41 -47 ";
## corr[3] = " 0 0 2 0 2 -81 0 32 0 ";
##
## for (x=0; x<4; x++)
## {
## fMm = eval(corr[x].substring(0,2));
## fM = eval(corr[x].substring(2,4));
## fF = eval(corr[x].substring(4,6));
## fD = eval(corr[x].substring(6,8));
## f_omega = eval(corr[x].substring(8,10));
## dp += 0.1*(eval(corr[x].substring(10,14))*sind(fD*D+fM*M+fMm*Mm+fF*F+f_omega*omega)+eval(corr[x].substring(14,18))*cosd(fD*D+fM*M+fMm*Mm+fF*F+f_omega*omega));
## de += 0.1*(eval(corr[x].substring(18,22))*cosd(fD*D+fM*M+fMm*Mm+fF*F+f_omega*omega)+eval(corr[x].substring(22,26))*sind(fD*D+fM*M+fMm*Mm+fF*F+f_omega*omega));
## }
## */
global delta_psi;
#Nutation in longitude
delta_psi = dp/36000000;
print ("Delta psi " + str(round(3600000*delta_psi)/1000)+"''"+",");
#Nutation in obliquity
delta_eps = de/36000000;
print ("Delta eps " + str(round(3600000*delta_eps)/1000)+"''"+",");
# Mean obliquity of the ecliptic
eps0 = (84381.448 - 46.815 * TE - 0.00059 * TE2 + 0.001813 * TE3)/3600;
global eps;
# True obliquity of the ecliptic
eps = eps0 + delta_eps;
# GHA Aries, GAST, GMST, equation of the equinoxes
def Aries():
# Mean GHA Aries
GHAAmean = norm_360_deg(280.46061837 + 360.98564736629 * (JD - 2451545) + 0.000387933 * T2 - T3 / 38710000);
print ("GHAAmean " + str(GHAAmean));
# GMST
SidTm = OutSidTime(GHAAmean);
print("GMST " + str(SidTm));
global GHAAtrue;
# True GHA Aries
GHAAtrue = norm_360_deg(GHAAmean + delta_psi * cosd(eps));
# GAST
SidTa = OutSidTime(GHAAtrue);
print ("GAST " + SidTa);
print ("True obl of Ecl " +str(OutECL(eps)));
# Equation of the equinoxes
EoE = 240*delta_psi*cosd(eps);
EoEout = round(1000*EoE)/1000;
EoEout = " " + str(EoEout) + "s"+",";
print ("- Aries -");
print ("GHA true " + str(OutHA(GHAAtrue)));
# Calculations for the sun
def Sun():
# Periodic terms for the sun
# Longitude
L0= 175347046;
L0+=3341656*math.cos(4.6692568+6283.0758500*Tau);
L0+=34894*math.cos(4.62610+12566.15170*Tau);
L0+=3497*math.cos(2.7441+5753.3849*Tau);
L0+=3418*math.cos(2.8289+3.5231*Tau);
L0+=3136*math.cos(3.6277+77713.7715*Tau);
L0+=2676*math.cos(4.4181+7860.4194*Tau);
L0+=2343*math.cos(6.1352+3930.2097*Tau);
L0+=1324*math.cos(0.7425+11506.7698*Tau);
L0+=1273*math.cos(2.0371+529.6910*Tau);
L0+=1199*math.cos(1.1096+1577.3435*Tau);
L0+=990*math.cos(5.233+5884.927*Tau);
L0+=902*math.cos(2.045+26.298*Tau);
L0+=857*math.cos(3.508+398.149*Tau);
L0+=780*math.cos(1.179+5223.694*Tau);
L0+=753*math.cos(2.533+5507.553*Tau);
L0+=505*math.cos(4.583+18849.228*Tau);
L0+=492*math.cos(4.205+775.523*Tau);
L0+=357*math.cos(2.920+0.067*Tau);
L0+=317*math.cos(5.849+11790.629*Tau);
L0+=284*math.cos(1.899+796.298*Tau);
L0+=271*math.cos(0.315+10977.079*Tau);
L0+=243*math.cos(0.345+5486.778*Tau);
L0+=206*math.cos(4.806+2544.314*Tau);
L0+=205*math.cos(1.869+5573.143*Tau);
L0+=202*math.cos(2.458+6069.777*Tau);
L0+=156*math.cos(0.833+213.299*Tau);
L0+=132*math.cos(3.411+2942.463*Tau);
L0+=126*math.cos(1.083+20.775*Tau);
L0+=115*math.cos(0.645+0.980*Tau);
L0+=103*math.cos(0.636+4694.003*Tau);
L0+=102*math.cos(0.976+15720.839*Tau);
L0+=102*math.cos(4.267+7.114*Tau);
L0+=99*math.cos(6.21+2146.17*Tau);
L0+=98*math.cos(0.68+155.42*Tau);
L0+=86*math.cos(5.98+161000.69*Tau);
L0+=85*math.cos(1.30+6275.96*Tau);
L0+=85*math.cos(3.67+71430.70*Tau);
L0+=80*math.cos(1.81+17260.15*Tau);
L0+=79*math.cos(3.04+12036.46*Tau);
L0+=75*math.cos(1.76+5088.63*Tau);
L0+=74*math.cos(3.50+3154.69*Tau);
L0+=74*math.cos(4.68+801.82*Tau);
L0+=70*math.cos(0.83+9437.76*Tau);
L0+=62*math.cos(3.98+8827.39*Tau);
L0+=61*math.cos(1.82+7084.90*Tau);
L0+=57*math.cos(2.78+6286.60*Tau);
L0+=56*math.cos(4.39+14143.50*Tau);
L0+=56*math.cos(3.47+6279.55*Tau);
L0+=52*math.cos(0.19+12139.55*Tau);
L0+=52*math.cos(1.33+1748.02*Tau);
L0+=51*math.cos(0.28+5856.48*Tau);
L0+=49*math.cos(0.49+1194.45*Tau);
L0+=41*math.cos(5.37+8429.24*Tau);
L0+=41*math.cos(2.40+19651.05*Tau);
L0+=39*math.cos(6.17+10447.39*Tau);
L0+=37*math.cos(6.04+10213.29*Tau);
L0+=37*math.cos(2.57+1059.38*Tau);
L0+=36*math.cos(1.71+2352.87*Tau);
L0+=36*math.cos(1.78+6812.77*Tau);
L0+=33*math.cos(0.59+17789.85*Tau);
L0+=30*math.cos(0.44+83996.85*Tau);
L0+=30*math.cos(2.74+1349.87*Tau);
L0+=25*math.cos(3.16+4690.48*Tau);
#L1
L1=628331966747;
L1+=206059*math.cos(2.678235+6283.075850*Tau);
L1+=4303*math.cos(2.6351+12566.1517*Tau);
L1+=425*math.cos(1.590+3.523*Tau);
L1+=119*math.cos(5.796+26.298*Tau);
L1+=109*math.cos(2.966+1577.344*Tau);
L1+=93*math.cos(2.59+18849.23*Tau);
L1+=72*math.cos(1.14+529.69*Tau);
L1+=68*math.cos(1.87+398.15*Tau);
L1+=67*math.cos(4.41+5507.55*Tau);
L1+=59*math.cos(2.89+5223.69*Tau);
L1+=56*math.cos(2.17+155.42*Tau);
L1+=45*math.cos(0.40+796.30*Tau);
L1+=36*math.cos(0.47+775.52*Tau);
L1+=29*math.cos(2.65+7.11*Tau);
L1+=21*math.cos(5.34+0.98*Tau);
L1+=19*math.cos(1.85+5486.78*Tau);
L1+=19*math.cos(4.97+213.30*Tau);
L1+=17*math.cos(2.99+6275.96*Tau);
L1+=16*math.cos(0.03+2544.31*Tau);
L1+=16*math.cos(1.43+2146.17*Tau);
L1+=15*math.cos(1.21+10977.08*Tau);
L1+=12*math.cos(2.83+1748.02*Tau);
L1+=12*math.cos(3.26+5088.63*Tau);
L1+=12*math.cos(5.27+1194.45*Tau);
L1+=12*math.cos(2.08+4694.00*Tau);
L1+=11*math.cos(0.77+553.57*Tau);
L1+=10*math.cos(1.30+6286.60*Tau);
L1+=10*math.cos(4.24+1349.87*Tau);
L1+=9*math.cos(2.70+242.73*Tau);
L1+=9*math.cos(5.64+951.72*Tau);
L1+=8*math.cos(5.30+2352.87*Tau);
L1+=6*math.cos(2.65+9437.76*Tau);
L1+=6*math.cos(4.67+4690.48*Tau);
# L2
L2=52919;
L2+=8720*math.cos(1.0721+6283.0758*Tau);
L2+=309*math.cos(0.867+12566.152*Tau);
L2+=27*math.cos(0.05+3.52*Tau);
L2+=16*math.cos(5.19+26.30*Tau);
L2+=16*math.cos(3.68+155.42*Tau);
L2+=10*math.cos(0.76+18849.23*Tau);
L2+=9*math.cos(2.06+77713.77*Tau);
L2+=7*math.cos(0.83+775.52*Tau);
L2+=5*math.cos(4.66+1577.34*Tau);
L2+=4*math.cos(1.03+7.11*Tau);
L2+=4*math.cos(3.44+5573.14*Tau);
L2+=3*math.cos(5.14+796.30*Tau);
L2+=3*math.cos(6.05+5507.55*Tau);
L2+=3*math.cos(1.19+242.73*Tau);
L2+=3*math.cos(6.12+529.69*Tau);
L2+=3*math.cos(0.31+398.15*Tau);
L2+=3*math.cos(2.28+553.57*Tau);
L2+=2*math.cos(4.38+5223.69*Tau);
L2+=2*math.cos(3.75+0.98*Tau);
#L3
L3=289*math.cos(5.844+6283.076*Tau);
L3+=35;
L3+=17*math.cos(5.49+12566.15*Tau);
L3+=3*math.cos(5.20+155.42*Tau);
L3+=1*math.cos(4.72+3.52*Tau);
L3+=1*math.cos(5.30+18849.23*Tau);
L3+=1*math.cos(5.97+242.73*Tau);
# L4
L4=114*math.cos(3.142);
L4+=8*math.cos(4.13+6283.08*Tau);
L4+=1*math.cos(3.84+12566.15*Tau);
#L5
L5 = 1*math.cos(3.14);
# Mean longitude of the sun
Lsun_mean = norm_360_deg(280.4664567+360007.6982779*Tau+0.03032028*Tau2+Tau3/49931-Tau4/15299-Tau5/1988000);
# Heliocentric longitude
Lhelioc = norm_360_deg((L0+L1*Tau+L2*Tau2+L3*Tau3+L4*Tau4+L5*Tau5)/1e8/dtr);
# Geocentric longitude
global Lsun_true;
Lsun_true = norm_360_deg(Lhelioc+180-0.000025);
# Latitude
#B0
B0=280*math.cos(3.199+84334.662*Tau);
B0+=102*math.cos(5.422+5507.553*Tau);
B0+=80*math.cos(3.88+5223.69*Tau);
B0+=44*math.cos(3.70+2352.87*Tau);
B0+=32*math.cos(4.00+1577.34*Tau);
#B1
B1=9*math.cos(3.90+5507.55*Tau);
B1+=6*math.cos(1.73+5223.69*Tau);
# Heliocentric latitude
B = (B0+B1*Tau)/1e8/dtr;
# Geocentric latitude
beta = norm_360_deg(-B);
# Corrections
Lsun_prime = norm_360_deg(Lhelioc+180-1.397*TE-0.00031*TE2);
beta = beta+0.000011*(cosd(Lsun_prime)-sind(Lsun_prime));
# Distance earth-sun
R0=100013989;
R0+=1670700*math.cos(3.0984635+6283.0758500*Tau);
R0+=13956*math.cos(3.05525+12566.15170*Tau);
R0+=3084*math.cos(5.1985+77713.7715*Tau);
R0+=1628*math.cos(1.1739+5753.3849*Tau);
R0+=1576*math.cos(2.8469+7860.4194*Tau);
R0+=925*math.cos(5.453+11506.770*Tau);
R0+=542*math.cos(4.564+3930.210*Tau);
R0+=472*math.cos(3.661+5884.927*Tau);
R0+=346*math.cos(0.964+5507.553*Tau);
R0+=329*math.cos(5.900+5223.694*Tau);
R0+=307*math.cos(0.299+5573.143*Tau);
R0+=243*math.cos(4.273+11790.629*Tau);
R0+=212*math.cos(5.847+1577.344*Tau);
R0+=186*math.cos(5.022+10977.079*Tau);
R0+=175*math.cos(3.012+18849.228*Tau);
R0+=110*math.cos(5.055+5486.778*Tau);
R0+=98*math.cos(0.89+6069.78*Tau);
R0+=86*math.cos(5.69+15720.84*Tau);
R0+=86*math.cos(1.27+161000.69*Tau);
R0+=65*math.cos(0.27+17260.15*Tau);
R0+=63*math.cos(0.92+529.69*Tau);
R0+=57*math.cos(2.01+83996.85*Tau);
R0+=56*math.cos(5.24+71430.70*Tau);
R0+=49*math.cos(3.25+2544.31*Tau);
R0+=47*math.cos(2.58+775.52*Tau);
R0+=45*math.cos(5.54+9437.76*Tau);
R0+=43*math.cos(6.01+6275.96*Tau);
R0+=39*math.cos(5.36+4694.00*Tau);
R0+=38*math.cos(2.39+8827.39*Tau);
R0+=37*math.cos(0.83+19651.05*Tau);
R0+=37*math.cos(4.90+12139.55*Tau);
R0+=36*math.cos(1.67+12036.46*Tau);
R0+=35*math.cos(1.84+2942.46*Tau);
R0+=33*math.cos(0.24+7084.90*Tau);
R0+=32*math.cos(0.18+5088.63*Tau);
R0+=32*math.cos(1.78+398.15*Tau);
R0+=28*math.cos(1.21+6286.60*Tau);
R0+=28*math.cos(1.90+6279.55*Tau);
R0+=26*math.cos(4.59+10447.39*Tau);
# R1
R1=103019*math.cos(1.107490+6283.075850*Tau);
R1+=1721*math.cos(1.0644+12566.1517*Tau);
R1+=702*math.cos(3.142);
R1+=32*math.cos(1.02+18849.23*Tau);
R1+=31*math.cos(2.84+5507.55*Tau);
R1+=25*math.cos(1.32+5223.69*Tau);
R1+=18*math.cos(1.42+1577.34*Tau);
R1+=10*math.cos(5.91+10977.08*Tau);
R1+=9*math.cos(1.42+6275.96*Tau);
R1+=9*math.cos(0.27+5486.78*Tau);
#R2
R2=4359*math.cos(5.7846+6283.0758*Tau);
R2+=124*math.cos(5.579+12566.152*Tau);
R2+=12*math.cos(3.14);
R2+=9*math.cos(3.63+77713.77*Tau);
R2+=6*math.cos(1.87+5573.14*Tau);
R2+=3*math.cos(5.47+18849.23*Tau);
#R3
R3=145*math.cos(4.273+6283.076*Tau);
R3+=7*math.cos(3.92+12566.15*Tau);
#R4
R4 = 4*math.cos(2.56+6283.08*Tau);
R = (R0+R1*Tau+R2*Tau2+R3*Tau3+R4*Tau4)/1e8;
# Apparent longitude of the sun
global lmbda;
lmbda = norm_360_deg(Lsun_true+delta_psi-0.005691611/R);
# Right ascension of the sun, apparent
global RAsun;
RAsun = norm_360_deg(math.atan2((sind(lmbda)*cosd(eps)-tand(beta)*sind(eps)),cosd(lmbda))/dtr);
# Sidereal hour angle of the sun, apparent
SHAsun = 360-RAsun;
# Declination of the sun, apparent
DECsun = math.asin(sind(beta)*cosd(eps)+cosd(beta)*sind(eps)*sind(lmbda))/dtr;
global Dsun;
Dsun = DECsun;
# GHA of the sun
GHAsun = norm_360_deg(GHAAtrue-RAsun);
# Semidiameter of the sun
SDsun = 959.63/R;
# Horizontal parallax of the sun
HPsun = 8.794/R;
dayfraction = (hour + minute/60 + second/3600)/24;
# Equation of time
# EOT = 4*(Lsun_mean-0.0057183-0.0008-RAsun+delta_psi*cosd(eps));
EOT = 4*GHAsun+720-1440*dayfraction;
if (EOT>20):
EOT-=1440;
if (EOT<-20):
EOT+=1440;
print ("- Sun -");
print ("GHA " + str(OutHA(GHAsun)));
print ("SHA " + str(OutHA(SHAsun)));
print ("Dec " + str(OutDec(Dsun)));
print ("SD " + str(OutSDHP(SDsun)));
print ("HP " + str(OutSDHP(HPsun)));
# Calculation of ephemerides for the moon
def Moon():
# Mean longitude of the moon
Lmoon_mean = norm_360_deg(218.3164591 + 481267.88134236*TE - 0.0013268*TE2 + TE3/538841 - TE4/65194000);
# Mean elongation of the moon
D = norm_360_deg(297.8502042+445267.1115168*TE-0.00163*TE2+TE3/545868-TE4/113065000);
# Mean anomaly of the sun
Msun_mean = norm_360_deg(357.5291092+35999.0502909*TE-0.0001536*TE2+TE3/24490000);
# Mean anomaly of the moon
Mmoon_mean = norm_360_deg(134.9634114+477198.8676313*TE+0.008997*TE2+TE3/69699-TE4/14712000);
# Mean distance of the moon from her ascending node
F = norm_360_deg(93.2720993+483202.0175273*TE-0.0034029*TE2-TE3/3526000+TE4/863310000);
# Corrections
A1 = 119.75+131.849*TE;
A1 = 360*(A1/360 - math.floor(A1/360));
A2 = 53.09 + 479264.29 * TE;
A2 = 360 * (A2/360 - math.floor(A2/360));
A3 = 313.45 + 481266.484 * TE;
A3 = 360*(A3/360 - math.floor(A3/360));
fE = 1-0.002516*TE-0.0000074*TE2;
fE2 = fE*fE;
# Periodic terms for the moon
# Longitude and distance
fD = [];
fMms = [];
fMmm = [];
fF = [];
coeffs = [];
coeffc = [];
fD.append(0);
fMms.append(0);
fMmm.append(1);
fF.append(0);
coeffs.append(6288774);
coeffc.append(-20905355);
fD.append(2);
fMms.append(0);
fMmm.append(-1);
fF.append(0);
coeffs.append(1274027);
coeffc.append(-3699111);
fD.append(2);
fMms.append(0);
fMmm.append(0);
fF.append(0);
coeffs.append(658314);
coeffc.append(-2955968);
fD.append(0);
fMms.append(0);
fMmm.append(2);
fF.append(0);
coeffs.append(213618);
coeffc.append(-569925);
fD.append(0);
fMms.append(1);
fMmm.append(0);
fF.append(0);
coeffs.append(-185116);
coeffc.append(48888);
fD.append(0);
fMms.append(0);
fMmm.append(0);
fF.append(2);
coeffs.append(-114332);
coeffc.append(-3149);
fD.append(2);
fMms.append(0);
fMmm.append(-2);
fF.append(0);
coeffs.append(58793);
coeffc.append(246158);
fD.append(2);
fMms.append(-1);
fMmm.append(-1);
fF.append(0);
coeffs.append(57066);
coeffc.append(-152138);
fD.append(2);
fMms.append(0);
fMmm.append(1);
fF.append(0);
coeffs.append(53322);
coeffc.append(-170733);
fD.append(2);
fMms.append(-1);
fMmm.append(0);
fF.append(0);
coeffs.append(45758);
coeffc.append(-204586);
fD.append(0);
fMms.append(1);
fMmm.append(-1);
fF.append(0);
coeffs.append(-40923);
coeffc.append(-129620);
fD.append(1);
fMms.append(0);
fMmm.append(0);
fF.append(0);
coeffs.append(-34720);
coeffc.append(108743);
fD.append(0);
fMms.append(1);
fMmm.append(1);
fF.append(0);
coeffs.append(-30383);
coeffc.append(104755);
fD.append(2);
fMms.append(0);
fMmm.append(0);
fF.append(-2);
coeffs.append(15327);
coeffc.append(10321);
fD.append(0);
fMms.append(0);
fMmm.append(1);
fF.append(2);
coeffs.append(-12528);
coeffc.append(0);
fD.append(0);
fMms.append(0);
fMmm.append(1);
fF.append(-2);
coeffs.append(10980);
coeffc.append(79661);
fD.append(4);
fMms.append(0);
fMmm.append(-1);
fF.append(0);
coeffs.append(10675);
coeffc.append(-34782);
fD.append(0);
fMms.append(0);
fMmm.append(3);
fF.append(0);
coeffs.append(10034);
coeffc.append(-23210);
fD.append(4);
fMms.append(0);
fMmm.append(-2);
fF.append(0);
coeffs.append(8548);
coeffc.append(-21636);
fD.append(2);
fMms.append(1);
fMmm.append(-1);
fF.append(0);
coeffs.append(-7888);
coeffc.append(24208);
fD.append(2);
fMms.append(1);
fMmm.append(0);
fF.append(0);
coeffs.append(-6766);
coeffc.append(30824);
fD.append(1);
fMms.append(0);
fMmm.append(-1);
fF.append(0);
coeffs.append(-5163);
coeffc.append(-8379);
fD.append(1);
fMms.append(1);
fMmm.append(0);
fF.append(0);
coeffs.append(4987);
coeffc.append(-16675);
fD.append(2);
fMms.append(-1);
fMmm.append(1);
fF.append(0);
coeffs.append(4036);
coeffc.append(-12831);
fD.append(2);
fMms.append(0);
fMmm.append(2);
fF.append(0);
coeffs.append(3994);
coeffc.append(-10445);
fD.append(4);
fMms.append(0);
fMmm.append(0);
fF.append(0);
coeffs.append(3861);
coeffc.append(-11650);
fD.append(2);
fMms.append(0);
fMmm.append(-3);
fF.append(0);
coeffs.append(3665);
coeffc.append(14403);
fD.append(0);
fMms.append(1);
fMmm.append(-2);
fF.append(0);
coeffs.append(-2689);
coeffc.append(-7003);
fD.append(2);
fMms.append(0);
fMmm.append(-1);
fF.append(2);
coeffs.append(-2602);
coeffc.append(0);
fD.append(2);
fMms.append(-1);
fMmm.append(-2);
fF.append(0);
coeffs.append(2390);
coeffc.append(10056);
fD.append(1);
fMms.append(0);
fMmm.append(1);
fF.append(0);
coeffs.append(-2348);
coeffc.append(6322);
fD.append(2);
fMms.append(-2);
fMmm.append(0);
fF.append(0);
coeffs.append(2236);
coeffc.append(-9884);
fD.append(0);
fMms.append(1);
fMmm.append(2);
fF.append(0);
coeffs.append(-2120);
coeffc.append(5751);
fD.append(0);
fMms.append(2);
fMmm.append(0);
fF.append(0);
coeffs.append(-2069);
coeffc.append(0);
fD.append(2);
fMms.append(-2);
fMmm.append(-1);
fF.append(0);
coeffs.append(2048);
coeffc.append(-4950);
fD.append(2);
fMms.append(0);
fMmm.append(1);
fF.append(-2);
coeffs.append(-1773);
coeffc.append(4130);
fD.append(2);
fMms.append(0);
fMmm.append(0);
fF.append(2)
coeffs.append(-1595);
coeffc.append(0);
fMms.append(0);
fD.append(4);
fMms.append(-1);
fMmm.append(-1);
fF.append(0);
coeffs.append(1215);
coeffc.append(-3958);
fD.append(0);
fMms.append(0);
fMmm.append(2);
fF.append(2);
coeffs.append(-1110);
coeffc.append(0);
fD.append(3);
fMms.append(0);
fMmm.append(-1);
fF.append(0);
coeffs.append(-892);
coeffc.append(3258);
fD.append(2);
fMms.append(1);
fMmm.append(1);
fF.append(0);
coeffs.append(-810);
coeffc.append(2616);
fD.append(4);
fMms.append(-1);
fMmm.append(-2);
fF.append(0);
coeffs.append(759);
coeffc.append(-1897);
fD.append(0);
fMms.append(2);
fMmm.append(-1);
fF.append(0);
coeffs.append(-713);
coeffc.append(-2117);
fD.append(2);
fMms.append(2);
fMmm.append(-1);
fF.append(0);
coeffs.append(-700);
coeffc.append(2354);
fD.append(2);
fMms.append(1);
fMmm.append(-2);
fF.append(0);
coeffs.append(691);
coeffc.append(0);
fD.append(2);
fMms.append(-1);
fMmm.append(0);
fF.append(-2);
coeffs.append(596);
coeffc.append(0);
fD.append(4);
fMms.append(0);
fMmm.append(1);
fF.append(0);
coeffs.append(549);
coeffc.append(-1423);
fD.append(0);
fMms.append(0);
fMmm.append(4);
fF.append(0);
coeffs.append(537);
coeffc.append(-1117);
fD.append(4);
fMms.append(-1);
fMmm.append(0);
fF.append(0);
coeffs.append(520);
coeffc.append(-1571);
fD.append(1);
fMms.append(0);
fMmm.append(-2);
fF.append(0);
coeffs.append(-487);
coeffc.append(-1739);
fD.append(2);
fMms.append(1);
fMmm.append(0);
fF.append(-2);
coeffs.append(-399);
coeffc.append(0);
fD.append(0);
fMms.append(0);
fMmm.append(2);
fF.append(-2);
coeffs.append(-381);
coeffc.append(-4421);
fD.append(1);
fMms.append(1);
fMmm.append(1);
fF.append(0);
coeffs.append(351);
coeffc.append(0);
fD.append(3);
fMms.append(0);
fMmm.append(-2);
fF.append(0);
coeffs.append(-340);
coeffc.append(0);
fD.append(4);
fMms.append(0);
fMmm.append(-3);
fF.append(0);
coeffs.append(330);
coeffc.append(0);
fD.append(2);
fMms.append(-1);
fMmm.append(2);
fF.append(0);
coeffs.append(327);
coeffc.append(0);
fD.append(0);
fMms.append(2);
fMmm.append(1);
fF.append(0);
coeffs.append(-323);
coeffc.append(1165);
fD.append(1);
fMms.append(1);
fMmm.append(-1);
fF.append(0);
coeffs.append(299);
coeffc.append(0);
fD.append(2);
fMms.append(0);
fMmm.append(3);
fF.append(0);
coeffs.append(294);
coeffc.append(0);
fD.append(2);
fMms.append(0);
fMmm.append(-1);
fF.append(-2);
coeffs.append(0);
coeffc.append(8752);
# Latitude
fD2 = [];
fMms2 = [];
fMmm2 = [];
fF2 = [];
coeffs2 = [];
fD2.append(0);
fMms2.append(0);
fMmm2.append(0);
fF2.append(1);
coeffs2.append(5128122);
fD2.append(0);
fMms2.append(0);
fMmm2.append(1);
fF2.append(1);
coeffs2.append(280602);
fD2.append(0);
fMms2.append(0);
fMmm2.append(1);
fF2.append(-1);
coeffs2.append(277693);
fD2.append(2);
fMms2.append(0);
fMmm2.append(0);
fF2.append(-1);
coeffs2.append(173237);
fD2.append(2);
fMms2.append(0);
fMmm2.append(-1);
fF2.append(1);
coeffs2.append(55413);
fD2.append(2);
fMms2.append(0);
fMmm2.append(-1);
fF2.append(-1);
coeffs2.append(46271);
fD2.append(2);
fMms2.append(0);
fMmm2.append(0);
fF2.append(1);
coeffs2.append(32573);
fD2.append(0);
fMms2.append(0);
fMmm2.append(2);
fF2.append(1);
coeffs2.append(17198);
fD2.append(2);
fMms2.append(0);
fMmm2.append(1);
fF2.append(-1);
coeffs2.append(9266);
fD2.append(0);
fMms2.append(0);
fMmm2.append(2);
fF2.append(-1);
coeffs2.append(8822);
fD2.append(2);
fMms2.append(-1);
fMmm2.append(0);
fF2.append(-1);
coeffs2.append(8216);
fD2.append(2);
fMms2.append(0);
fMmm2.append(-2);
fF2.append(-1);
coeffs2.append(4324);
fD2.append(2);
fMms2.append(0);
fMmm2.append(1);
fF2.append(1);
coeffs2.append(4200);
fD2.append(2);
fMms2.append(1);
fMmm2.append(0);
fF2.append(-1);
coeffs2.append(-3359);
fD2.append(2);
fMms2.append(-1);
fMmm2.append(-1);
fF2.append(1);
coeffs2.append(2463);
fD2.append(2);
fMms2.append(-1);
fMmm2.append(0);
fF2.append(1);
coeffs2.append(2211);
fD2.append(2);
fMms2.append(-1);
fMmm2.append(-1);
fF2.append(-1);
coeffs2.append(2065);
fD2.append(0);
fMms2.append(1);
fMmm2.append(-1);
fF2.append(-1);
coeffs2.append(-1870);
fD2.append(4);
fMms2.append(0);
fMmm2.append(-1);
fF2.append(-1);
coeffs2.append(1828);
fD2.append(0);
fMms2.append(1);
fMmm2.append(0);
fF2.append(1);
coeffs2.append(-1794);
fD2.append(0);
fMms2.append(0);
fMmm2.append(0);
fF2.append(3);
coeffs2.append(-1749);
fD2.append(0);
fMms2.append(1);
fMmm2.append(-1);
fF2.append(1);
coeffs2.append(-1565);
fD2.append(1);
fMms2.append(0);
fMmm2.append(0);
fF2.append(1);
coeffs2.append(-1491);
fD2.append(0);
fMms2.append(1);
fMmm2.append(1);
fF2.append(1);
coeffs2.append(-1475);
fD2.append(0);
fMms2.append(1);
fMmm2.append(1);
fF2.append(-1);
coeffs2.append(-1410);
fD2.append(0);
fMms2.append(1);
fMmm2.append(0);
fF2.append(-1);
coeffs2.append(-1344);
fD2.append(1);
fMms2.append(0);
fMmm2.append(0);
fF2.append(-1);
coeffs2.append(-1335);
fD2.append(0);
fMms2.append(0);
fMmm2.append(3);
fF2.append(1);
coeffs2.append(1107);
fD2.append(4);
fMms2.append(0);
fMmm2.append(0);
fF2.append(-1);
coeffs2.append(1021);
fD2.append(4);
fMms2.append(0);
fMmm2.append(-1);
fF2.append(1);
coeffs2.append(833);
fD2.append(0);
fMms2.append(0);
fMmm2.append(1);
fF2.append(-3);
coeffs2.append(777);
fD2.append(4);
fMms2.append(0);
fMmm2.append(-2);
fF2.append(1);
coeffs2.append(671);
fD2.append(2);
fMms2.append(0);
fMmm2.append(0);
fF2.append(-3);
coeffs2.append(607);
fD2.append(2);
fMms2.append(0);
fMmm2.append(2);
fF2.append(-1);
coeffs2.append(596);
fD2.append(2);
fMms2.append(-1);
fMmm2.append(1);
fF2.append(-1);
coeffs2.append(491);
fD2.append(2);
fMms2.append(0);
fMmm2.append(-2);
fF2.append(1);
coeffs2.append(-451);
fD2.append(0);
fMms2.append(0);
fMmm2.append(3);
fF2.append(-1);
coeffs2.append(439);
fD2.append(2);
fMms2.append(0);
fMmm2.append(2);
fF2.append(1);
coeffs2.append(422);
fD2.append(2);
fMms2.append(0);
fMmm2.append(-3);
fF2.append(-1);
coeffs2.append(421);
fD2.append(2);
fMms2.append(1);
fMmm2.append(-1);
fF2.append(1);
coeffs2.append(-366);
fD2.append(2);
fMms2.append(1);
fMmm2.append(0);
fF2.append(1);
coeffs2.append(-351);
fD2.append(4);
fMms2.append(0);
fMmm2.append(0);
fF2.append(1);
coeffs2.append(331);
fD2.append(2);
fMms2.append(-1);
fMmm2.append(1);
fF2.append(1);
coeffs2.append(315);
fD2.append(2);
fMms2.append(-2);
fMmm2.append(0);
fF2.append(-1);
coeffs2.append(302);
fD2.append(0);
fMms2.append(0);
fMmm2.append(1);
fF2.append(3);
coeffs2.append(-283);
fD2.append(2);
fMms2.append(1);
fMmm2.append(1);
fF2.append(-1);
coeffs2.append(-229);
fD2.append(1);
fMms2.append(1);
fMmm2.append(0);
fF2.append(-1);
coeffs2.append(223);
fD2.append(1);
fMms2.append(1);
fMmm2.append(0);
fF2.append(1);
coeffs2.append(223);
fD2.append(0);
fMms2.append(1);
fMmm2.append(-2);
fF2.append(-1);
coeffs2.append(-220);
fD2.append(2);
fMms2.append(1);
fMmm2.append(-1);
fF2.append(-1);
coeffs2.append(-220);
fD2.append(1);
fMms2.append(0);
fMmm2.append(1);
fF2.append(1);
coeffs2.append(-185);
fD2.append(2);
fMms2.append(-1);
fMmm2.append(-2);
fF2.append(-1);
coeffs2.append(181);
fD2.append(0);
fMms2.append(1);
fMmm2.append(2);
fF2.append(1);
coeffs2.append(-177);
fD2.append(4);
fMms2.append(0);
fMmm2.append(-2);
fF2.append(-1);
coeffs2.append(176);
fD2.append(4);
fMms2.append(-1);
fMmm2.append(-1);
fF2.append(-1);
coeffs2.append(166);
fD2.append(1);
fMms2.append(0);
fMmm2.append(1);
fF2.append(-1);
coeffs2.append(-164);
fD2.append(4);
fMms2.append(0);
fMmm2.append(1);
fF2.append(-1);
coeffs2.append(132);
fD2.append(1);
fMms2.append(0);
fMmm2.append(-1);
fF2.append(-1);
coeffs2.append(-119);
fD2.append(4);
fMms2.append(-1);
fMmm2.append(0);
fF2.append(-1);
coeffs2.append(115);
fD2.append(2);
fMms2.append(-2);
fMmm2.append(0);
fF2.append(1);
coeffs2.append(107);
sumL = 0;
sumr = 0;
sumB = 0;
x=0;
while(x<60):
f = 1;
if(abs(fMms[x])==1):
f=fE;
if(abs(fMms[x])==2):
f=fE2;
# print str(x) +"," + str(fD[x])+"," + str(fMms[x])+"," + str(fMmm[x])+"," + str(fF[x])+"," + str(coeffs[x])+"," + str(coeffc[x])+"," + str(fD2[x])+"," +str(fMms2[x])+"," + str(fMmm2[x])+"," + str(coeffs2[x])+"," + str(fF2[x]);
sumL += f*(coeffs[x]*sind(fD[x]*D+fMms[x]*Msun_mean+fMmm[x]*Mmoon_mean+fF[x]*F));
sumr += f*(coeffc[x]*cosd(fD[x]*D+fMms[x]*Msun_mean+fMmm[x]*Mmoon_mean+fF[x]*F));
f = 1;
if(abs(fMms2[x])==1):
f=fE;
if(abs(fMms2[x])==2):
f=fE2;
sumB += f*(coeffs2[x]*sind(fD2[x]*D+fMms2[x]*Msun_mean+fMmm2[x]*Mmoon_mean+fF2[x]*F));
x=x+1
# Corrections
sumL = sumL+3958*sind(A1)+1962*sind(Lmoon_mean-F)+318*sind(A2);
sumB = sumB - 2235 * sind(Lmoon_mean) + 382 * sind(A3) + 175 * sind(A1-F) + 175 * sind(A1+F) + 127 * sind(Lmoon_mean-Mmoon_mean)- 115 * sind(Lmoon_mean + Mmoon_mean);
# Longitude of the moon
lambdaMm = norm_360_deg(Lmoon_mean + sumL/1000000);
# Latitude of the moon
betaM = sumB/1000000;
# Distance earth-moon
dEM = 385000.56 + sumr / 1000;
# Apparent longitude of the moon
global lambdaMapp;
lambdaMapp = lambdaMm + delta_psi;
global RAmoon;
# Right ascension of the moon, apparent
RAmoon = norm_360_deg(math.atan2((sind(lambdaMapp)*cosd(eps)-tand(betaM)*sind(eps)),cosd(lambdaMapp))/dtr);
# Sidereal hour angle of the moon, apparent
SHAmoon = 360-RAmoon;
global DECMoon;
# Declination of the moon
DECmoon = math.asin(sind(betaM)*cosd(eps)+cosd(betaM)*sind(eps)*sind(lambdaMapp))/dtr;
Dmoon = DECmoon;
# GHA of the moon
GHAmoon = norm_360_deg(GHAAtrue - RAmoon);
print ("- Moon -");
print ("GHA " + str(OutHA(GHAmoon)));
print ("SHA " + str(OutHA(SHAmoon)));
print ("Dec " + str(OutDec(Dmoon)));
# Horizontal parallax of the moon
HPmoon = 3600*math.asin(6378.14/dEM)/dtr;
print ("HP Moon " + str(OutSDHP(HPmoon)));
# Semidiameter of the moon
SDmoon = 3600*math.asin(1738/dEM)/dtr;
print ("SD Moon " + str(OutSDHP(SDmoon)));
# Geocentric angular distance between moon and sun
LDist = math.acos(sind(Dmoon)*sind(Dsun)+cosd(Dmoon)*cosd(Dsun)*cosd(RAmoon-RAsun))/dtr;
# Phase of the moon
i = lambdaMapp-lmbda;
global k;
k = 100*(1-cosd(i))/2;
k = round(10*k)/10;
# Ephemerides of Polaris
def Polaris():
# Equatorial coordinates of Polaris at 2000.0 (mean equinox and equator 2000.0)
RApol0 = 37.95293333;
DECpol0 = 89.26408889;
# Proper motion per year
dRApol = 2.98155/3600;
dDECpol = -0.0152/3600;
# Equatorial coordinates at Julian Date T (mean equinox and equator 2000.0)
RApol1 = RApol0+100*TE*dRApol;
DECpol1 = DECpol0+100*TE*dDECpol;
# Mean obliquity of ecliptic at 2000.0 in degrees
eps0_2000 = 23.439291111;
# Transformation to ecliptic coordinates in radians (mean equinox and equator 2000.0)
lambdapol1 = math.atan2((sind(RApol1)*cosd(eps0_2000)+tand(DECpol1)*sind(eps0_2000)),cosd(RApol1));
betapol1 = math.asin(sind(DECpol1)*cosd(eps0_2000)-cosd(DECpol1)*sind(eps0_2000)*sind(RApol1));
# Precession
eta = (47.0029*TE-0.03302*TE2+0.00006*TE3)*dtr/3600;
PI0 = (174.876384-(869.8089*TE+0.03536*TE2)/3600)*dtr;
p0 = (5029.0966*TE+1.11113*TE2-0.0000006*TE3)*dtr/3600;
A1 = math.cos(eta)*math.cos(betapol1)*math.sin(PI0-lambdapol1)-math.sin(eta)*math.sin(betapol1);
B1 = math.cos(betapol1)*math.cos(PI0-lambdapol1);
C1 = math.cos(eta)*math.sin(betapol1)+math.sin(eta)*math.cos(betapol1)*math.sin(PI0-lambdapol1);
lambdapol2 = p0+PI0-math.atan2(A1,B1);
betapol2 = math.asin(C1);
# Nutation in longitude
lambdapol2 += dtr*delta_psi;
# Aberration
kappa = dtr*20.49552/3600;
pi0 = dtr*(102.93735+1.71953*TE+0.00046*TE2);
e = 0.016708617-0.000042037*TE-0.0000001236*TE2;
dlambdapol = (e*kappa*math.cos(pi0-lambdapol2)-kappa*math.cos(dtr*Lsun_true-lambdapol2))/math.cos(betapol2);
dbetapol = -kappa*math.sin(betapol2)*(math.sin(dtr*Lsun_true-lambdapol2)-e*math.sin(pi0-lambdapol2));
lambdapol2 += dlambdapol;
betapol2 += dbetapol;
# Transformation back to equatorial coordinates in radians
RApol2 = math.atan2((math.sin(lambdapol2)*cosd(eps)-math.tan(betapol2)*sind(eps)),math.cos(lambdapol2));
DECpol2 = math.asin(math.sin(betapol2)*cosd(eps)+math.cos(betapol2)*sind(eps)*math.sin(lambdapol2));
# Finals
GHApol = GHAAtrue-RApol2/dtr;
GHApol = norm_360_deg(GHApol);
SHApol = 360-RApol2/dtr;
SHApol = norm_360_deg(SHApol);
DECpol = DECpol2/dtr;
print ("- Polaris -");
print ("GHA " + str(OutHA(GHApol)));
print ("SHA " + str(OutHA(SHApol)));
print ("Dec " + str(OutDec(DECpol)));
# Calculation of the phase of the moon
def MoonPhase():
x = lambdaMapp-lmbda;
x = norm_360_deg(x);
if(x>0 and x<180):
quarter = " (+)";
if(x>180 and x<360):
quarter = " (-)";
if (k==0):
quarter = ", new";
if (k==100):
quarter = ", full";
print ("Moon Phase " + str(k) + "%" + str(quarter));
# Day of the week
def Weekday():
global JD0h;
JD0h = JD0h + 1.5;
res = JD0h-7 * math.floor(JD0h/7);
if (res == 0): DoW = " SUN";
if (res == 1): DoW = " MON";
if (res == 2): DoW = " TUE";
if (res == 3): DoW = " WED";
if (res == 4): DoW = " THU";
if (res == 5): DoW = " FRI";
if (res == 6): DoW = " SAT";
print ("Day of week = " + DoW);
def Star(starname):
# Star catalog
navstar =[];
starcat=[];
# Acamar
navstar.append( " 2.9710266670 -40.3047138890 -0.3910 1.9400 0.0280")
# Achernar
navstar.append( " 1.6285700000 -57.2367166670 1.1730 -3.4700 0.0230")
#Acrux
navstar.append( "12.4432975000 -63.0990500000 -0.5240 -1.2100 0.0000")
# Adhara
navstar.append( " 6.9770966670 -28.9720833330 0.0310 0.2800 0.0000")
# Aldebaran
navstar.append( " 4.5986769440 16.5092750000 0.4390 -18.9700 0.0480")
# Alioth
navstar.append( "12.9004855560 55.9598527780 1.3280 -0.5800 0.0090")
# Alkaid
navstar.append( "13.7923427780 49.3133194440 -1.2490 -1.0900 0.0350")
# Al Na'ir
navstar.append( "22.1372222220 -46.9609972220 1.2590 -15.1000 0.0510")
# Alnilam
navstar.append( " 5.6035580560 -1.2019500000 0.0060 -0.2400 0.0000")
# Alphard
navstar.append( " 9.4597908330 -8.6586527780 -0.0930 3.2800 0.0170")
# Alphecca
navstar.append("15.5781322220 26.7147055560 0.9060 -8.8600 0.0430")
# Alpheratz
navstar.append(" 0.1397958330 29.0904388890 1.0390 -16.3300 0.0240")
# Altair
navstar.append("19.8463894440 8.8683416670 3.6290 38.6300 0.1981")
# Ankaa
navstar.append(" 0.4380638890 -42.3060583330 1.8330 -39.5700 0.0350")
# Antares
navstar.append("16.4901219440 -26.4319861110 -0.0710 -2.0300 0.0190")
# Arcturus
navstar.append("14.2610213890 19.1824194440 -7.7140 -199.8400 0.0900")
# Atria
navstar.append("16.8110747220 -69.0277277780 0.2600 -3.4000 0.0240")
# Avior
navstar.append(" 8.3752313890 -59.5095861110 -0.3460 1.4400 0.0000")
# Bellatrix
navstar.append(" 5.4188491670 6.3496500000 -0.0590 -1.3900 0.0260")
# Betelgeuse
navstar.append(" 5.9195297220 7.4070416670 0.1730 0.8700 0.0050")
# Canopus
navstar.append(" 6.3991997220 -52.6956944440 0.2450 2.0700 0.0180")
# Capella
navstar.append(" 5.2781536110 45.9980277780 0.7280 -42.4700 0.0730")
# Deneb
navstar.append("20.6905325000 45.2803638890 0.0270 0.2300 0.0000")
# Denebola
navstar.append("11.8176611110 14.5720416670 -3.4220 -11.4100 0.0760")
# Diphda
navstar.append(" 0.7264922220 -17.9866166670 1.6370 3.2500 0.0570")
# Dubhe
navstar.append("11.0621294440 61.7508944440 -1.6750 -6.6500 0.0310")
# Elnath
navstar.append(" 5.4381975000 28.6074083330 0.1690 -17.5100 0.0180")
# Eltanin
navstar.append("17.9434352780 51.4889472220 -0.0810 -1.9400 0.0170")
# Enif
navstar.append("21.7364344440 9.8749777780 0.2070 -0.0600 0.0060")
# Fomalhaut
navstar.append("22.9608486110 -29.6222500000 2.5510 -16.4700 0.1440")
# Gacrux
navstar.append("12.5194247220 -57.1131944440 0.2850 -26.2300 0.0000")
# Gienah
navstar.append("12.2634350000 -17.5419361110 -1.1240 2.3300 0.0000")
# Hadar
navstar.append("14.0637244440 -60.3729972220 -0.4260 -1.9300 0.0160")
# Hamal
navstar.append(" 2.1195563890 23.4624055560 1.3830 -14.8300 0.0430")
# Kaus Aust.
navstar.append("18.4028686110 -34.3846472220 -0.3090 -12.4100 0.0150")
# Kochab
navstar.append("14.8450961110 74.1554944440 -0.7630 1.2200 0.0310")
# Markab
navstar.append("23.0793494440 15.2052500000 0.4360 -4.2500 0.0300")
# Menkar
navstar.append(" 3.0379925000 4.0897027780 -0.0630 -7.8000 0.0090")
# Menkent
navstar.append("14.1113752780 -36.3700083330 -4.2930 -51.9000 0.0590")
# Miaplacidus
navstar.append(" 9.2199880560 -69.7172083330 -3.1080 10.7800 0.0380")
# Mirfak
navstar.append(" 3.4053791670 49.8612055560 0.2460 -2.4600 0.0290")
# Nunki
navstar.append("18.9210900000 -26.2967305560 0.0990 -5.4200 0.0000")
# Peacock
navstar.append("20.4274588890 -56.7351055560 0.0820 -8.9100 0.0000")
# Polaris
navstar.append(" 2.5301955560 89.2640888890 19.8770 -1.520 0.0070")
# Pollux
navstar.append(" 7.7552627780 28.0261833330 -4.7400 -4.5900 0.0930")
# Procyon
navstar.append(" 7.6550313890 5.2250166670 -4.7550 -102.2900 0.2880")
# Rasalhague
navstar.append("17.5822433330 12.5600388890 0.8220 -22.6400 0.0560")
# Regulus
navstar.append("10.1395319440 11.9671916670 -1.6930 0.6400 0.0390")
# Rigel
navstar.append(" 5.2422966670 -8.2016611110 0.0030 -0.1300 0.0130")
# Rigil Kent.
navstar.append("14.6599680560 -60.8354000000 -49.8260 69.9300 0.7516")
# Sabik
navstar.append("17.1729669440 -15.7249194400 0.2600 9.5000 0.0520")
# Schedar
navstar.append(" 0.6751250000 56.5373500000 0.6360 -3.1900 0.0160")
# Shaula
navstar.append("17.5601483330 -37.1038111110 -0.0110 -2.9200 0.0000")
# Sirius
navstar.append(" 6.7524641670 -16.7161083330 -3.8470 -120.5300 0.3751")
# Spica
navstar.append("13.4198852780 -11.1613083330 -0.2780 -2.8300 0.0210")
# Suhail
navstar.append(" 9.1332711110 -43.4326055560 -0.1720 1.2700 0.0150")
# Vega
navstar.append("18.6156477780 38.7836583330 1.7260 28.6100 0.1230")
# Zubenelgenubi
navstar.append("14.8479758330 -16.0417833330 -0.7340 -6.6800 0.0490")
starcat.append("Acamar");
starcat.append("Achernar");
starcat.append("crux");
starcat.append("Adhara");
starcat.append("Aldebaran");
starcat.append("Alioth");
starcat.append("Alkaid");
starcat.append("Al Na'ir");
starcat.append("Alnilam");
starcat.append("Alphard");
starcat.append("Alphecca");
starcat.append("Alpheratz");
starcat.append("Altair");
starcat.append("Ankaa");
starcat.append("Antares");
starcat.append("Arcturus");
starcat.append("Atria");
starcat.append("Avior");
starcat.append("Bellatrix");
starcat.append("Betelgeuse");
starcat.append("Canopus");
starcat.append("Capella");
starcat.append("Deneb");
starcat.append("Denebola");
starcat.append("Diphda");
starcat.append("Dubhe");
starcat.append("Elnath");
starcat.append("Eltanin");
starcat.append("Enif");
starcat.append("Fomalhaut");
starcat.append("Gacrux");
starcat.append("Gienah");
starcat.append("Hadar");
starcat.append("Hamal");
starcat.append("Kaus Aust.");
starcat.append("Kochab");
starcat.append("Markab");
starcat.append("Menkar");
starcat.append("Menkent");
starcat.append("Miaplacidus");
starcat.append("Mirfak");
starcat.append("Nunki");
starcat.append("Peacock");
starcat.append("Polaris");
starcat.append("Pollux");
starcat.append("Procyon");
starcat.append("Rasalhague");
starcat.append("Regulus");
starcat.append("Rigel");
starcat.append("Rigil Kent.");
starcat.append("Sabik");
starcat.append("Schedar");
starcat.append("Shaula");
starcat.append("Sirius");
starcat.append("Spica");
starcat.append("Suhail");
starcat.append("Vega");
starcat.append("Zubenelgenubi");
# starname=1;
tempstring = navstar[starname]
# Read catalog
RAstar0 = 15*eval(tempstring[0:13]);
DECstar0 = eval(tempstring[14:28]);
dRAstar = 15*eval(tempstring[29:37])/3600;
dDECstar = eval(tempstring[38:47])/3600;
par = eval( tempstring[49:55])/3600;
# Equatorial coordinates at Julian Date T (mean equinox and equator 2000.0)
RAstar1 = RAstar0+TE*dRAstar;
DECstar1 = DECstar0+TE*dDECstar;
# Mean obliquity of ecliptic at 2000.0 in degrees
eps0_2000 = 23.439291111;
# Transformation to ecliptic coordinates in radians (mean equinox and equator 2000.0)
lambdastar1 = math.atan2((sind(RAstar1)*cosd(eps0_2000)+tand(DECstar1)*sind(eps0_2000)),cosd(RAstar1));
betastar1 = math.asin(sind(DECstar1)*cosd(eps0_2000)-cosd(DECstar1)*sind(eps0_2000)*sind(RAstar1));
# Precession
eta = (47.0029*TE-0.03302*TE2+0.00006*TE3)*dtr/3600;
PI0 = (174.876384-(869.8089*TE+0.03536*TE2)/3600)*dtr;
p0 = (5029.0966*TE+1.11113*TE2-0.0000006*TE3)*dtr/3600;
A1 = math.cos(eta)*math.cos(betastar1)*math.sin(PI0-lambdastar1)-math.sin(eta)*math.sin(betastar1);
B1 = math.cos(betastar1)*math.cos(PI0-lambdastar1);
C1 = math.cos(eta)*math.sin(betastar1)+math.sin(eta)*math.cos(betastar1)*math.sin(PI0-lambdastar1);
lambdastar2 = p0+PI0-math.atan2(A1,B1);
betastar2 = math.asin(C1);
# Annual parallax
par_lambda = dtr*par*math.sin(dtr*Lsun_true-lambdastar2)/math.cos(betastar2);
par_beta = -dtr*par*math.sin(betastar2)*math.cos(dtr*Lsun_true-lambdastar2);
lambdastar2 += par_lambda;
betastar2 += par_beta;
# Nutation in longitude
lambdastar2 += dtr*delta_psi;
# Aberration
kappa = dtr*20.49552/3600;
pi0 = dtr*(102.93735+1.71953*TE+0.00046*TE2);
e = 0.016708617-0.000042037*TE-0.0000001236*TE2;
dlambdastar = (e*kappa*math.cos(pi0-lambdastar2)-kappa*math.cos(dtr*Lsun_true-lambdastar2))/math.cos(betastar2);
dbetastar = -kappa*math.sin(betastar2)*(math.sin(dtr*Lsun_true-lambdastar2)-e*math.sin(pi0-lambdastar2));
lambdastar2 += dlambdastar;
betastar2 += dbetastar;
# Transformation back to equatorial coordinates in radians
RAstar2 = math.atan2((math.sin(lambdastar2)*cosd(eps)-math.tan(betastar2)*sind(eps)),math.cos(lambdastar2));
DECstar2 = math.asin(math.sin(betastar2)*cosd(eps)+math.cos(betastar2)*sind(eps)*math.sin(lambdastar2));
global DECmoon;
global RAmoon;
# Lunar distance of star
LDist = math.acos(sind(DECmoon)*math.sin(DECstar2)+cosd(DECmoon)*math.cos(DECstar2)*cosd(RAmoon-RAstar2/dtr))/dtr;
# Finals
GHAstar = norm_360_deg(GHAAtrue-RAstar2/dtr);
SHAstar = norm_360_deg(360-RAstar2/dtr);
DECstar = DECstar2/dtr;
print ("- " + starcat[starname] + " -");
print ("GHA = " + str(OutHA(GHAstar)));
print ("SHA = " + str(OutHA(SHAstar)));
print ("Dec = " + str(OutDec(DECstar)));
# Data output
def Output():
# Sun
GHAsun = OutHA(GHAsun);
SHAsun = OutHA(SHAsun);
DECsun = OutDec(DECsun);
SDsun = OutSDHP(SDsun);
HPsun = OutSDHP(HPsun);
# Moon
GHAmoon = OutHA(GHAmoon);
SHAmoon = OutHA(SHAmoon);
DECmoon = OutDec(DECmoon);
SDmoon = OutSDHP(SDmoon);
HPmoon = OutSDHP(HPmoon);
# Aries
AR = OutHA(GHAAtrue);
#Polaris
GHApolaris = OutHA(GHApol);
SHApolaris = OutHA(SHApol);
DECpolaris = OutDec(DECpol);
# Obliquity of Ecliptic
OoE = OutECL(eps0);
tOoE = OutECL(eps);
# Equation of time
if (EOT<0):
sign="-";
else: sign="+";
EOT = abs(EOT);
EOTmin = math.floor(EOT);
EOTsec = round(600*(EOT-EOTmin))/10;
if (EOTsec-math.floor(EOTsec)==0):
EOTsec+=".0";
if (EOTmin==0):
EOT = " " + sign + " " + EOTsec + "s";
else:
EOT = " "+sign+" "+EOTmin+"m "+EOTsec+"s";
# Illumination
illum=" "+k+"%"+quarter;
# Lunar Distance of Sun
LDist = OutHA(LDist);
print (GHAsun);
if __name__ =="__main__":
main()