;-------------------------------------------------------------------------------
; This file contains uM-FPU64 user-defined functions for performing date and time
; calculations, calculating the position of celestial objects in the sky, and
; calculating the screen coordinates for displaying the objects on the LCD display.
;
; The compiled code must be loaded to the uM-FPU64 chip over a serial connection
; using the uM-FPU64 IDE software.
;
; author: Geoffrey Gill
;
;---------- register definitions -----------------------------------------------

Year		equ	F10		'FPU registers
Month		equ	F11
Day		equ	F12
Hour		equ	F13
Minute		equ	F14
Second		equ	F15
X		equ	F16		'X screen coordinate
Y		equ	F17		'Y screen coordinate
LX		equ	L18		'int X
LY		equ	L19		'int Y
LMAG		equ	L20		'int Magnitude
Tcube		equ	F134		'T cubed
Tsq		equ	F135		'T squared
Longitude	equ	F136		'Longitude at Las Cruces, NM - radians
Latitude	equ	F137		'Latitude at Las Cruces, NM - radians
A		equ	F138		'Temp register
B		equ	F139		'Temp register
JD		equ	F140
T		equ	F141
SideTimeG	equ	F142
Revolutions	equ	F143
RevRadians	equ	F144
STatGACTUAL	equ	F145		'in radians
STatLCACT	equ	F146
H		equ	F147		'local hour angle of object
RA		equ	F148		'Right Ascencion of object
DEC		equ	F149		'Declination of object 
ALT		equ	F150		'Alt of object in radians
AZ		equ	F151		'Azm of object in radians
C		equ	F152		'Temp reg
D		equ	F153		'Temp reg
MAGNITUDE	equ	F154		'Magnitude of object with sign
SPECTRUM	equ	F155		'Color of object - single letter
PIXPERRAD	equ	F156		'171.8873385 pixels/radian
MeanLSun	equ	F157		'Mean longitutde of the Sun
MeanAnomSun	equ	F158		'Mean Anomaly of the Sun
Ecc		equ	F159		'Eccentricity of earth's orbit
TrueLSun	equ	F160		'True longitutde of Sun - radians
OblEcliptic	equ	F161		'Obliquity of ecliptic - radians
MeanLMoon 	equ	F162		'Mean Longitude of Moon
MeanAnomMoon 	equ	F163		'Mean Anomaly of Moon
MoonMeanElong equ	F164		'Moon Mean Elogation
MeanDisMoonNode equ	F165		'Mean Distance of Moon from Ascending node
MoonGeoCentLong equ	F166		'Geocentric Longitude of Moon
MoonGeoCentLat equ 	F167		'Geocentric Latitude of Moon
CentralAngle	equ	F168		'Central angle of star in radians
EccAnomSun	equ	F181		'Eccentric Anomaly of the Sun
rSun		equ	F182		'Radius vectore of the Sun
L		equ	F229		'Mean Longitude of planet
ecc		equ	F230		'eccentricity of planet orbit
i 		equ	F231		'inclination of planet orbit
O 		equ	F232		'longitude of planet ascending node
MeanAnom	equ	F233		'Mean Anomaly of planet in radians
EccAnom	equ	F234		'Eccentricity Anomaly of planet in radians
TrueAnom	equ	F235		'True Anomaly of planet in radians
ArgLat		equ	F236		'Argument of Latitude of planet
EclLong	equ	F237		'Ecliptical Longitude of planet
EclLat		equ	F238		'Ecliptical Latitude of planet
r		equ	F239		'radius vector of planet in AU
GeoCentLong 	equ 	F240		'GeoCentric Longitude of planet
DistEarth 	equ	F241		'Distance from planet to Earth AU
GeoCentLat 	equ	F242		'GeoCentric Latitude of planet

;---------- function prototypes ------------------------------------------------

JulianDate	func	1		'Function prototypes
CalcT		func	%		'2
STat0hrG	func	%		'3
STatGact	func	%		'4		
STatLoc	func	%		'5
AltAz  	func	%		'6
Sun		func	%		'7
Moon		func 	%		'8
Venus		func	%		'9
Mercury	func	%		'10
Mars		func	%		'11
Jupiter	func	%		'12
Saturn		func	%		'13
Position	func	%		'14

;-------------------------------------------------------------------------------
#FUNCTION	JulianDate		'First function

PIXPERRAD = 171.8873385
Longitude = 1.86562642	'Your logitude in radians should be entered here
Latitude = 0.563741348	'Your latitude in radians should be entered here
Year = RTC(READ_TIME+YEAR) + 2000
Month = RTC(READ_TIME+MONTH)
Day = RTC(READ_TIME+DAY)
Hour = RTC(READ_TIME+HOUR)
Minute = RTC(READ_TIME+MINUTE)

IF (Month<2) THEN
	Year = Year-1
	Month = Month+12
ENDIF

A=FLOOR(Year/100)
B=2-A+FLOOR(A/4)

JD=1720994.5 + FLOOR(365.25*Year) + FLOOR(30.6001*(Month+1)) + Day + B

;-------------------------------------------------------------------------------
#FUNCTION	CalcT			'Second function

T=(JD-2415020)/36525
Tsq=T*T
Tcube=Tsq*T

;-------------------------------------------------------------------------------
#FUNCTION	STat0hrG		'3rd function Sidereal Time in radians at 0 HRs GMT

SideTimeG = .27691938 + 100.0021359 * T + .000001075 * Tsq
Revolutions = SideTimeG - FLOOR(SideTimeG)
RevRadians = Revolutions * 2 * pi

;-------------------------------------------------------------------------------
#FUNCTION	STatGact		'4th function

A = 1.002737909350795 * 2 * pi
B = (Hour+(Minute/60))/24 * A
STatGACTUAL = RevRadians + B

;-------------------------------------------------------------------------------
#FUNCTION	STatLoc		'5th function

B = 2 * pi
STatLCACT = STatGACTUAL - Longitude
WHILE (STatLCACT > B)
	STatLCACT = STatLCACT - B
	LOOP
WHILE (STatLCACT < 0)
	STatLCACT = STatLCACT + B
	LOOP

;-------------------------------------------------------------------------------
#FUNCTION	AltAz			'6th function calculates alt-az
					'in screen coordinates
H = STatLCACT - RA			
A = SIN(Latitude) * SIN(DEC)
B = COS(Latitude) * COS(DEC) * COS(H)
ALT = ASIN(A+B)
A = COS(H)*SIN(Latitude)
B = COS(Latitude)*TAN(DEC)
C = SIN(H)
D = A-B
AZ = ATAN2(C,D)
;use stereographic projection from equator 
CentralAngle = ACOS(COS(ALT)*COS(AZ))
Y =CentralAngle*PIXPERRAD*SIN(ALT)/SIN(CentralAngle) 	'get altitude in pixels
Y = 239 - Y				'invert y coordinate
'IF (Y>239.0) THEN
	'Y=999				'999=invisible on screen
'ENDIF
'IF (Y<0.0) THEN
'	Y=999
'ENDIF
X = CentralAngle*PIXPERRAD*COS(ALT)*SIN(AZ)/SIN(CentralAngle)	'get azimuth in pixels
X = X + 159				'center on screen
'IF (X>319) THEN
'	X=999
'ENDIF
'IF (X<0.0) THEN
'	X=999
'ENDIF
LX=FIXR(X)
LY=FIXR(Y)
MAGNITUDE =4.0 - MAGNITUDE
LMAG=FIXR(MAGNITUDE)

;-------------------------------------------------------------------------------
#FUNCTION	Sun			'calculate RA and Dec of Sun

MeanLSun = 279.69668 + 36000.76892 * T + .0003025 * Tsq				'degrees
MeanLSun = MeanLSun/360									'revolutions
MeanLSun = FRAC(MeanLSun) * 2 * PI								'radians
MeanAnomSun = 358.47583 +35999.04975 * T - .00015 * Tsq -.0000033 * Tcube 		'degrees
MeanAnomSun = MeanAnomSun/360								'revolutions
MeanAnomSun = FRAC(MeanAnomSun)*2*PI							'radians
Ecc = .01675104 -.0000418 * T - .000000126 * Tsq
EccAnomSun=ATAN2(SIN(MeanAnomSun),(COS(MeanAnomSun)-Ecc))				'radians
rSun=1.0000002*(1-Ecc*COS(EccAnomSun))							'AU
C = (1.91946 - .004789*T-.000014*Tsq)*SIN(MeanAnomSun) + (.020094 - .0001*T)*SIN(2*MeanAnomSun)+.000293*SIN(3*MeanAnomSun) 'degrees
C = C/360											'revolutions
C = FRAC(C)*2*PI										'radians
TrueLSun = MeanLSun + C									'radians
OblEcliptic = 23.452294 - .0130125*T - .00000164*Tsq + .000000503*Tcube		'degrees
OblEcliptic = (OblEcliptic/180)*PI								'radians
A = COS(OblEcliptic)*SIN(TrueLSun)
B = COS(TrueLSun)
RA = ATAN2(A,B)
DEC = ASIN(SIN(OblEcliptic)*SIN(TrueLSun)	)

;-------------------------------------------------------------------------------
#FUNCTION 	Moon										'calculate RA and Dec of Moon
;note I changed 270.434164 to 280.434164 to correct moon position
MeanLMoon = 280.434164 + 481267.8831 * T							'degrees						'degrees
MeanLMoon = MeanLMoon/360									'revolutions
MeanLMoon = FRAC(MeanLMoon)*360								'degrees
MeanAnomMoon = 296.104608 + 477198.8491*T							'degrees
MeanAnomMoon = MeanAnomMoon/360								'revolutions
MeanAnomMoon = FRAC(MeanAnomMoon)*2*PI							'radians
MoonMeanElong = 350.737486 + 445267.1142*T						'degrees
MoonMeanElong = MoonMeanElong/360								'revolutions
MoonMeanElong = FRAC(MoonMeanElong)*2*PI							'radians
MeanDisMoonNode = 11.250889 + 483202.0251*T						'degrees
MeanDisMoonNode = MeanDisMoonNode/360							'revolutions
MeanDisMoonNode = FRAC(MeanDisMoonNode)*2*PI						'radians
MoonGeoCentLong = MeanLMoon+6.28875*SIN(MeanAnomMoon)+1.274018*SIN(2*MoonMeanElong-MeanAnomMoon) 'partial
MoonGeoCentLong = MoonGeoCentLong+.658309*SIN(2*MoonMeanElong)+.213616*SIN(2*MeanAnomMoon) 'complete
MoonGeoCentLong=MoonGeoCentLong/360							'revolutions
MoonGeoCentLong=FRAC(MoonGeoCentLong)*2*PI						'radians
MoonGeoCentLat = 5.128189*SIN(MeanDisMoonNode)+.280606*SIN(MeanAnomMoon+MeanDisMoonNode) 'partial
MoonGeoCentLat = MoonGeoCentLat+.277693*SIN(MeanAnomMoon-MeanDisMoonNode)+.173238*SIN(2*MoonMeanElong-MeanDisMoonNode) 'complete
MoonGeoCentLat=MoonGeoCentLat*PI/180								'radians
;convert to RA and DEC
A=SIN(MoonGeoCentLong)*COS(OblEcliptic)-TAN(MoonGeoCentLat)*SIN(OblEcliptic)	'numerator
B=COS(MoonGeoCentLong)									'denominator
RA=ATAN2(A,B)											'radians
A=SIN(MoonGeoCentLat)*COS(OblEcliptic)+COS(MoonGeoCentLat)*SIN(OblEcliptic)*SIN(MoonGeoCentLong)
DEC=ASIN(A)											'radians

;-------------------------------------------------------------------------------
#FUNCTION	Venus			'calculate orbital elements, RA and DEC of Venus

L = 342.767053+58519.21191*T+.0003097*Tsq						'degrees
L = L/360										'revolutions
L = FRAC(L)*2*PI									'radians
ecc = .00682069-.00004774*T+.000000091*Tsq
i = 3.393631+.0010058*T-.000001*Tsq						'degrees
i = i/360										'revolutions
i = FRAC(i)*2*PI									'radians
O = 75.779647+.89985*T+.00041*Tsq							'degrees
O = O/360										'revolutions
O = FRAC(O)*2*PI									'radians
MeanAnom=212.60322+58517.80387*T+.001286*Tsq					'degrees
MeanAnom=MeanAnom/360								'revolutions
MeanAnom=FRAC(MeanAnom)*2*PI							'radians
EccAnom=ATAN2(SIN(MeanAnom),(COS(MeanAnom)-ecc))				'radians
r=.7233316*(1-ecc*(COS(EccAnom)))							'astronomical units
Position

;-------------------------------------------------------------------------------
#FUNCTION	Mercury		'calculate orbital elements, RA and DEC of Mercury

L = 178.179078+149474.07078*T+.0003011*Tsq					'degrees
L = L/360										'revolutions
L = FRAC(L)*2*PI									'radians
ecc = .20561421-.00002046*T-.00000003*Tsq
i = 7.002881+.0018608*T-.0000183*Tsq						'degrees
i = i/360										'revolutions
i = FRAC(i)*2*PI									'radians
O = 47.145944+1.1852083*T+.0001739*Tsq						'degrees
O = O/360										'revolutions
O = FRAC(O)*2*PI									'radians
MeanAnom=102.27938+149472.51529*T+.000007*Tsq					'degrees
MeanAnom=MeanAnom/360								'revolutions
MeanAnom=FRAC(MeanAnom)*2*PI							'radians
EccAnom=ATAN2(SIN(MeanAnom),(COS(MeanAnom)-ecc))				'radians
r=.3870986*(1-ecc*(COS(EccAnom)))							'astronomical units
Position

;-------------------------------------------------------------------------------
#FUNCTION	Mars			'calculate orbital elements, RA and DEC of Mars

L = 293.737334+19141.69551*T+.0003107*Tsq						'degrees
L = L/360										'revolutions
L = FRAC(L)*2*PI									'radians
ecc = .0933129+.000092064*T-.000000077*Tsq
i = 1.85033-.000675*T+.0000126*Tsq							'degrees
i = i/360										'revolutions
i = FRAC(i)*2*PI									'radians
O = 48.786442+.7709917*T-.0000014*Tsq	-.00000533*Tcube			'degrees
O = O/360										'revolutions
O = FRAC(O)*2*PI									'radians
MeanAnom=319.51913+19139.85475*T+.000181*Tsq        				'degrees
MeanAnom=MeanAnom/360								'revolutions
MeanAnom=FRAC(MeanAnom)*2*PI							'radians
EccAnom=ATAN2(SIN(MeanAnom),(COS(MeanAnom)-ecc))				'radians
r=1.5236883*(1-ecc*(COS(EccAnom)))							'astronomical units
Position

;-------------------------------------------------------------------------------										'radians
#FUNCTION	Jupiter			'calculate orbital elements, RA and DEC of Jupiter

L = 238.049257+3036.301986*T+.0003347*Tsq-.00000165*Tcube			'degrees
L = L/360										'revolutions
L = FRAC(L)*2*PI									'radians
ecc = .04833475+.00016418*T-.0000004676*Tsq-.0000000017*Tcube
i = 1.308736-.0056961*T+.0000039*Tsq						'degrees
i = i/360										'revolutions
i = FRAC(i)*2*PI									'radians
O = 99.443414+1.01053*T+.00035222*Tsq-.00000851*Tcube				'degrees
O = O/360										'revolutions
O = FRAC(O)*2*PI									'radians
MeanAnom=225.32833+3034.69202*T-.000722*Tsq					'degrees
MeanAnom=MeanAnom/360								'revolutions
MeanAnom=FRAC(MeanAnom)*2*PI							'radians
EccAnom=ATAN2(SIN(MeanAnom),(COS(MeanAnom)-ecc))				'radians
r=5.202561*(1-ecc*(COS(EccAnom)))							'astronomical units
Position

;-------------------------------------------------------------------------------
#FUNCTION	Saturn			'calculate orbital elements, RA and DEC of Saturn

L = 266.564377+1223.509884*T+.0003245*Tsq-.0000058*Tcube			'degrees
L = L/360										'revolutions
L = FRAC(L)*2*PI									'radians
ecc = .05589232-.0003455*T-.000000728*Tsq+.00000000074*Tcube
i = 2.492519-.0039189*T-.00001549*Tsq						'degrees
i = i/360										'revolutions
i = FRAC(i)*2*PI									'radians
O = 112.790414+.8731951*T-.00015218*Tsq-.00000531*Tcube			'degrees
O = O/360										'revolutions
O = FRAC(O)*2*PI									'radians
MeanAnom=175.46622+1221.55147*T-.000502*Tsq					'degrees
MeanAnom=MeanAnom/360								'revolutions
MeanAnom=FRAC(MeanAnom)*2*PI							'radians
EccAnom=ATAN2(SIN(MeanAnom),(COS(MeanAnom)-ecc))				'radians
r=9.554747*(1-ecc*(COS(EccAnom)))							'astronomical units
Position

;-------------------------------------------------------------------------------
#FUNCTION	Position		'calculate geocentric long, lat, RA, DEC	
		
TrueAnom=2*ATAN(SQRT((1+ecc)/(1-ecc))*TAN(EccAnom/2))				'radians
ArgLat=L+TrueAnom-MeanAnom-O							'radians
EclLong=O+ATAN2((COS(i)*SIN(ArgLat)),COS(ArgLat))				'radians
EclLat=ASIN(SIN(ArgLat)*SIN(i))							'radians
A=r*COS(EclLat)*SIN(EclLong-TrueLSun)
B=r*COS(EclLat)*COS(EclLong-TrueLSun)+rSUN
GeoCentLong=ATAN2(A,B)+ TrueLSun							'radians
DistEarth=SQRT(A*A+B*B+r*r*SIN(EclLat)*SIN(EclLat))				'AU
GeoCentLat=ASIN(r*SIN(EclLat)/DistEarth)						'radians
;convert to RA and DEC
A=SIN(GeoCentLong)*COS(OblEcliptic)-TAN(GeoCentLat)*SIN(OblEcliptic)		'numerator
B=COS(GeoCentLong)									'denominator
RA=ATAN2(A,B)										'radians
A=SIN(GeoCentLat)*COS(OblEcliptic)+COS(GeoCentLat)*SIN(OblEcliptic)*SIN(GeoCentLong)
DEC=ASIN(A)								
#END




