1+ #!/usr/bin/env python3
2+ import math
3+ from astropy .coordinates import SkyCoord , EarthLocation , AltAz , Angle
4+ from astropy import units as u
5+ from astropy .time import Time
6+ from astropy .utils import iers
7+ import sys , subprocess
8+
9+ def polarcalc (mylat , mylong , myelev , time , p1RA , p1DEC , p2RA , p2DEC , p3RA , p3DEC ):
10+ #iers.conf.auto_download = False
11+ #iers.conf.auto_max_age = None
12+
13+ #Create time object based on given time
14+ observing_time = Time ('2020-08-10 3:00:05' )
15+
16+ #Create location object based on lat/long/elev
17+ observing_location = EarthLocation (lat = mylat * u .deg , lon = mylong * u .deg , height = 177 * u .m )
18+
19+ #Create coordinate objects for each point
20+ p1 = SkyCoord (p1RA , p1DEC , unit = 'deg' )
21+ p2 = SkyCoord (p2RA , p2DEC , unit = 'deg' )
22+ p3 = SkyCoord (p3RA , p3DEC , unit = 'deg' )
23+ p1X = (90 - p1 .dec .degree ) * math .cos (p1 .ra .radian )
24+ p1Y = (90 - p1 .dec .degree ) * math .sin (p1 .ra .radian )
25+ p2X = (90 - p2 .dec .degree ) * math .cos (p2 .ra .radian )
26+ p2Y = (90 - p2 .dec .degree ) * math .sin (p2 .ra .radian )
27+ p3X = (90 - p3 .dec .degree ) * math .cos (p3 .ra .radian )
28+ p3Y = (90 - p3 .dec .degree ) * math .sin (p3 .ra .radian )
29+
30+ #Calculate center of circle using three points in the complex plane. DEC is treated as unitless for the purposes of the calculation.
31+ x , y , z = complex (p1X ,p1Y ), complex (p2X ,p2Y ), complex (p3X ,p3Y )
32+ w = z - x
33+ w /= y - x
34+ c = (x - y )* (w - abs (w )** 2 )/ 2j / w .imag - x
35+ resultX = - c .real
36+ resultY = c .imag
37+
38+ #Convert X/Y values of circle into RA/DEC
39+ resultDEC = (90 - math .sqrt (resultX ** 2 + resultY ** 2 ))
40+ resultRA = math .atan2 (resultY , resultX )* 360 / (2 * math .pi )
41+ if resultRA < 0 :
42+ resultRA = (180 - abs (resultRA ))+ 180
43+
44+ #Create coordinate object for current alignment offset
45+ offset = SkyCoord (resultRA , resultDEC , frame = 'icrs' , unit = 'deg' )
46+ print (f"Current alignment in RA/DEC: { offset .ra .to_string (u .hour , precision = 2 )} /{ offset .dec .to_string (u .deg , precision = 2 )} ." )
47+
48+ #Create coordinate object for pole
49+ pole = SkyCoord (0 , 90 , frame = 'icrs' , unit = 'deg' )
50+
51+ #Create coordinate object for pole
52+ poleAzAlt = pole .transform_to (AltAz (obstime = observing_time ,location = observing_location ))
53+ print (f"True polar alignment in Az./Alt.: 0h00m00s/{ Angle (mylat * u .degree ).to_string (u .degree , sep = ('d' , 'm' , 's' ), precision = 2 )} ." )
54+
55+ #Transform current alignment to Alt/Az coordinate system
56+ #offsetAzAlt = offset.transform_to(AltAz(obstime=observing_time,location=observing_location, temperature=25*u.deg_C, pressure=101325*u.Pa, relative_humidity=0.5))
57+ offsetAzAlt = offset .transform_to (AltAz (obstime = observing_time ,location = observing_location ))
58+ print (f"Current alignment in Az./Alt.: { offsetAzAlt .az .to_string (u .hour , precision = 2 )} /{ offsetAzAlt .alt .to_string (u .deg , precision = 2 )} ." )
59+
60+ #Calculate offset deltas from pole
61+ if offsetAzAlt .az .deg < 180 :
62+ errorAz = - offsetAzAlt .az .deg * 60
63+ else :
64+ errorAz = (360 - offsetAzAlt .az .deg )* 60
65+ print (f"Azimuth error correction is: { errorAz :.4f} arcminutes." )
66+ errorAlt = (mylat - offsetAzAlt .alt .deg )* 60
67+ print (f"Altitude error correction is: { errorAlt :.4f} arcminutes." )
68+
69+ return errorAz , errorAlt
70+
71+ #Latitude in degrees
72+ mylat = float (sys .argv [1 ])
73+
74+ #Longitude in degrees
75+ mylong = float (sys .argv [2 ])
76+
77+ #Elevation in meters
78+ myelev = sys .argv [3 ]
79+
80+ #YYYY-MM-DD HH:MM:SS format
81+ time = sys .argv [4 ]
82+
83+ #All RA/DEC values must be in compatible format to Astropy.coordinates library.
84+ #Preferrably degrees, but 00h00m00.0s and 00d00m00.0s should also work
85+ p1RA = sys .argv [5 ]
86+ p1DEC = sys .argv [6 ]
87+ p2RA = sys .argv [7 ]
88+ p2DEC = sys .argv [8 ]
89+ p3RA = sys .argv [9 ]
90+ p3DEC = sys .argv [10 ]
91+
92+ serialport = '/dev/ttyACM0'
93+
94+ result = polarcalc (mylat , mylong , myelev , time , p1RA , p1DEC , p2RA , p2DEC , p3RA , p3DEC )
95+
96+ #Verify error correction can be handled by AutoPA hardware (assuming it is in home/centered position)
97+ moveAz = "N"
98+ if abs (result [0 ]) > 192 :
99+ moveAz = input ("Azimuth error may be out of bounds of hardware capabilities if not in home position. Continue? (Y/N): " )
100+ else :
101+ moveAz = "Y"
102+ if moveAz .upper () == "Y" :
103+ #Call process to move azimuth using elevated privileges to override any existing serial connection
104+ subprocess .call (['sudo' , './altaz.py' , "az" , str (result [0 ]), serialport ])
105+
106+ moveAlt = "N"
107+ if result [1 ] > 168 :
108+ moveAz = input ("Altitude error may be out of bounds of hardware capabilities if not in home position. Continue? (Y/N): " )
109+ elif result [1 ] > 432 :
110+ moveAz = input ("Altitude error may be out of bounds of hardware capabilities if not in home position. Continue? (Y/N): " )
111+ else :
112+ moveAlt = "Y"
113+ if moveAlt .upper () == "Y" :
114+ #Call process to move altitude using elevated privileges to override any existing serial connection
115+ subprocess .call (['sudo' , './altaz.py' , "alt" , str (result [1 ]), serialport ])
0 commit comments