|
|
SUBSIM: The Web's #1 resource for all submarine & naval simulations since 1997
|
SUBSIM: The Web's #1 resource for all submarine & naval simulations since 1997 |
03-27-19, 03:56 AM
|
#46 | |
|
Count Dracula
Join Date: Aug 2012
Location: 46°13′N, 24°47′E
Posts: 1,816
Downloads: 405
Uploads: 0
|
Quote:
__________________
|
|
|
|
|
03-27-19, 12:43 PM
|
#47 | |
|
Stowaway
Posts: n/a
Downloads:
Uploads:
|
Quote:
I hope your research is more productive than mine. I searched the contents of all textual files (*.ini, *.cfg, *.py, etc.) for "clock" (as much as I dared) and found no switch I could identify. I really hope this Local Time issue can be resolved. Without it, identifying accurate longitude will be impossible, AFAICT. It may just turn out I am asking the simulator to do more than it was ever designed to do. Unfortunately, this would be a show-stopper for me as my interest in Silent Hunter is to sail the Seven Seas like the mariners did in the past. I know there could be some utility in using latitude (via sextant) with DR as if daytime was solid overcast and nights were clear. I'll give it a run and see how it goes. Another possibility I thought about was to use a real-world watch/clock for Local Time, but success would require the game clock and the RW clock to run at the same speed. Another thing to test out.  If anyone has any comment on the above, please share your thoughts. EDIT... I don't give up easy and SH5 is gorgeous! I don't want to hive that up. I also like to focus on what can be done as opposed to complaining about what is lacking. Soooo, I am going to test out using DR with Latitude via sextant for a TC1 voyage from Memal Harbor to Keil. I will also be double-plotting in SH5 and Google Earth Pro. Last edited by Capt Codetrucker; 03-27-19 at 06:46 PM. |
|
|
|
|
03-28-19, 01:01 PM
|
#48 |
|
Stowaway
Posts: n/a
Downloads:
Uploads:
|
Hmmm... One thing which keeps going around in my head is how has anyone been able to calculate longitude given the clock minutes lock between GMT and Local/Naut in SH5? I can't be the only one since the beginning of the SH5 franchise that has sailed without any help from the automatic Navigator?
|
|
|
|
03-28-19, 09:25 PM
|
#49 |
|
Seaman
 Join Date: Mar 2016
Posts: 33
Downloads: 12
Uploads: 0
|
My 2 pfennigs
I have done sextant sights sometime ago in SH5. I seem to remember being in the first patrol around Gdansk, and my position was in Berlin or some such.
AFAIR: You only need a clock on GMT to sight any celestial body. You then use the SHA and the site reduction tables to generate your first position. You then run, plot your new position using DR, and then take a sight again (after say an hour). This gives you your longitude and latitude (but it is fiddly - and as you can tell from the story I was not that accurate). I wrote a python programme that would generate a perpetual calendar for sun, moon and major stars (although there were some slight errors with the star SHA I think). If anyone wants it I can add here somewhere (possibly). Let me know. Finally: I thought that as the earth was flat in SH5, you could not do CelNav. hence I stopped. Is this not the case? Can someone confirm. Danke. Rmax. |
|
|
|
|
03-29-19, 02:43 AM
|
#50 | ||
|
Stowaway
Posts: n/a
Downloads:
Uploads:
|
Quote:
Yes, I would like to try your program. I have been doing many experiments trying to derive longitude manually, but without success. Quote:
The biggest hurdle in SH5 for CelNav is the SH5 Local/Naut Time's minutes do not reflect real world local time. The SH5 time's minutes are always synced to the GMT's minutes. I can't wrap my brain around this since it makes sense, but when I try to use this for calculating longitude I always am stuck with a multiple of 15 with zero minutes. I hope this makes sense? |
||
|
|
|
03-29-19, 06:20 AM
|
#51 |
|
Seaman
Join Date: Mar 2016
Posts: 33
Downloads: 12
Uploads: 0
|
Perpetual Calendar in Python
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()
Last edited by rmax; 04-01-19 at 04:36 AM. Reason: Updated to Python 3 |
|
|
|
|
03-29-19, 06:25 AM
|
#52 |
|
Seaman
Join Date: Mar 2016
Posts: 33
Downloads: 12
Uploads: 0
|
Perpetual Calendar in Python install
Instructions:
Download python: https://www.python.org/downloads/ Open a new file in IDLE Copy and paste the code into the window Hit f5 It will ask you to save it. Give it a filename It will give a perpetual calendar so you can always get GHA and SHA. This is the first step. It was one of my first python programmes, so the code is pretty hairy. It works for the Sun well, but for some reasons there seem to be some errors when I compare to professional tables. I don't think it will give much error, but it is a start at least. |
|
|
|
|
03-31-19, 04:02 PM
|
#53 |
|
Stowaway
Posts: n/a
Downloads:
Uploads:
|
Thank you very much. I will look forward to making use of your contribution.
|
|
|
|
04-01-19, 02:10 AM
|
#54 |
|
Seaman
Join Date: Mar 2016
Posts: 33
Downloads: 12
Uploads: 0
|
Cel Nav PDF
You are welcome. I converted to Python as I was going to generate a SH5 Almanac (very few people are concerned with Navigation from 1939 - 1945...can't think why).
Anyhow. I found this gives a good description of Cel Nav, and how you use SHA and LHA. http://www.dacust.com/navigation/pdf...lNavigator.pdf It is a time since I tried this in SH5 (and the results were not good - but that could have been my lousy sun shots). At present I am stalking a Merchant ship off the coast of the UK, so not going to be getting into this in the short-term. If I do, I'll try and help as much as possible. |
|
|
|
|
04-01-19, 02:26 AM
|
#55 | |
|
Stowaway
Posts: n/a
Downloads:
Uploads:
|
Quote:
Thank you. FYI - I created a thread to log my CelNav experiences. You can find it at "Navigating Using Only Pilotage, DR, Compass, Speed, Chronometer, and In-Game Sextant." I'm also going to be keeping a sharp eye out for MagDec, current influences, etc. if any. It's a long-term project, but it might be my way to "give back" to this fantastic community. |
|
|
|
|
04-01-19, 03:29 AM
|
#56 |
|
Seaman
Join Date: Mar 2016
Posts: 33
Downloads: 12
Uploads: 0
|
Your CelNav adventure
I saw the thread, but didn't want to comment there and pollute it!
Just looking through my notes, I made a lot of use of these http://www.navigation-spreadsheets.com/ Some of the calculations end up being a little fiddly to do so this helped me a lot. |
|
|
|
|
04-01-19, 04:37 AM
|
#57 |
|
Seaman
Join Date: Mar 2016
Posts: 33
Downloads: 12
Uploads: 0
|
Updated Python Code
For various reasons I originally wrote the code in Python 2.x. I have updated it to Python 3, so now should work!
|
|
|
|
 |
| Thread Tools | |
| Display Modes | |
|
|
Linear Mode
