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