Source code for hetutils
'''The hetutils.py package provides functions to do calculations specific
to the Hobby Eberly Telescope at McDonald Observatory. These functions
provide information about telescope pointing as well as tracker information.
'''
from __future__ import print_function
from math import *
#
# Some astronomical routines
#
Latitude_deg = 30.681436
Latitude_rad = radians(Latitude_deg)
LongitudeWest_deg = 104.014742
LongitudeWest_rad = radians(LongitudeWest_deg)
LongitudeEast_deg = 360.0 - LongitudeWest_deg
LongitudeEast_rad = radians(LongitudeEast_deg)
Altitude_deg = 55.055223
Altitude_rad = radians(Altitude_deg)
Height_meters = 2003.0
Fs_meters = 13.37884180
Fs_millimeters = Fs_meters * 1000.0
het_P = cos(Altitude_rad) * cos(Latitude_rad) # radians
het_Q = sin(Altitude_rad) * sin(Latitude_rad) # radians
# limit on the total angular travel of the tracker
BetaLimit_deg = 8.5
BetaLimit_rad = radians(BetaLimit_deg)
X_offset_mm = -3.030
Y_offset_mm = 25.690
Z_offset_mm = 0.025
T_offset_mm = 0.050
P_offset_mm = -0.100
R_offset_mm = 0.000
# Some basic math functions
[docs]def degrees_to_hours(deg):
'''Convert decimal degrees to decimal hours.'''
return (deg / 15.0)
[docs]def hours_to_degrees(hr):
'''Convert decimal hours to decimal degrees.'''
return (hr * 15.0)
[docs]def radians_to_hours(rad):
'''Convert decimal radians to decimal hours.'''
return (degrees_to_hours(degrees(deg)))
[docs]def hours_to_radians(hr):
'''Convert decimal hours to decimal radians.'''
return radians(hours_to_degrees(hr))
[docs]def modulo(x, y):
'''A general purpose modulo (x mod y) function.
Definition of 'mod'
x mod y = x - y * floor( x/y) for y != 0
x mod 0 = x
if 0 < y then 0 <= (x mod y) < y
if 0 > y then 0 >= (x mod y) > y
'''
if 0.0 == y:
return x
return x - y * floor(x/y)
[docs]def clipToCircle_deg(x):
''' Return the value (dX mod 360). This has the result of clipping the
value between 0 <= dR < 360.0.
Inputs dX - the value you want clipped
Outputs none
Return (dX mod 360)
'''
return modulo(x, 360.0)
[docs]def clipToCircle_rad(x):
'''Return the value (dX mod 2PI). This has the result of clipping the
value between 0 <= dR < 2PI.
Inputs dX - the value you want clipped
Outputs none
Return (dX mod 2PI)
'''
return modulo(x, 2*pi)
# azimuth with the longest track time as a function of declination
[docs]def BestAzimuth(dec_deg):
'''This function calculates the best azimuth for a given declination.
The function returns an azimuth only on the east side. If you want
the west azimuth, simply add 180 degrees to the return value.
The input declination is in degrees, The return value is the azimuth
in degrees or None is the declinination is not visible at HET.'''
t = type(dec_deg)
if (t is not float and t is not int and t is not long):
# should throw an exception here
#print('error BestAzimuth dec_deg is %s' % str(dec_deg))
return None
dec_rad = radians(dec_deg)
tde_n = asin(het_Q + het_P)
tde_s = asin(het_Q - het_P)
if dec_rad > tde_s and dec_rad < tde_n:
return degrees(acos((sin(dec_rad) - het_Q) / het_P))
elif dec_rad <= tde_s and dec_rad > (tde_s - BetaLimit_rad):
# valid south trajectory
return 180.0
elif dec_rad >= tde_n and dec_rad < (tde_n + BetaLimit_rad):
# valid north trajectory
return 0.0
else:
return None
# telecentric declination as a function of azimuth
[docs]def Tde(azimuth_deg):
'''Calculate the telecentric declination as a function of azimuth.
The input azimuth is in degrees. Returns the declination in degrees
or None if the input is not between 0 <= az < 360.0'''
t = type(azimuth_deg)
if (t is not float and t is not int and t is not long) or azimuth_deg < 0.0 or azimuth_deg >= 360.0:
# should throw an exception here
#print('error Tde az is %s' % str(azimuth_deg))
return None
az_rad = radians(azimuth_deg)
return degrees(asin(het_P * cos(az_rad) + het_Q))
# telecentric hour angle as a function of azimuth
[docs]def Ho(azimuth_deg): #returns degrees
'''Calculate the telecentric hour angle for a given azimuth.
Negative values are east hour angles. The input azimuth is in
degrees. Returns the hour angle in degrees.'''
t = type(azimuth_deg)
if (t is not float and t is not int and t is not long) or azimuth_deg < 0.0 or azimuth_deg >= 360.0:
# should throw an exception here
#print('error Ho az is %s' % str(azimuth_deg))
return None
tde_rad = radians(Tde(azimuth_deg))
az_rad = radians(azimuth_deg)
return degrees(asin(-cos(Altitude_rad) * sin(az_rad) / cos(tde_rad)))
# paralatic angle as a function of azimuth
[docs]def Pa(azimuth_deg): # returns degrees
'''Calculate the zero-crossing paralactic angle as a function of azimuth.
The input azimuth is in degrees. Returns the paralactic angle in degrees.'''
t = type(azimuth_deg)
if (t is not float and t is not int and t is not long) or azimuth_deg < 0.0 or azimuth_deg >= 360.0:
# should throw an exception here
#print('error Pa az is %s' % str(azimuth_deg))
return None
Ho_rad = radians(Ho(azimuth_deg))
az_rad = radians(azimuth_deg)
Pa_rad = acos( cos(Ho_rad) * cos(az_rad)
+ sin(Ho_rad) * sin(az_rad) * sin(Latitude_rad))
if azimuth_deg > 180.0:
Pa_rad += pi
return degrees(Pa_rad)
# available time for track as a function of azimuth and declination
[docs]def Att(az_deg, dec_deg):
'''Calculate the total possible track time for a give azimuth and declination.
Input arguments are in degrees. The function returns the track time in minutes or
0 if no trajectory is available.'''
t = type(az_deg)
if (t is not float and t is not int and t is not long) or az_deg < 0.0 or az_deg >= 360.0:
# should throw an exception here
#print('error Att az_deg is %s' % str(azimuth_deg))
return 0.0
t = type(dec_deg)
if (t is not float and t is not int and t is not long):
# should throw an exception here
#print('error Att dec_deg is %s' % str(azimuth_deg))
return 0.0
az_rad = radians(az_deg)
dec_rad = radians(dec_deg)
tde_rad = Tde(az_deg)
Ho_rad = Ho(az_deg)
# C1-C5, a, b, c, t1, t2 are intermidiate variables
C1 = sin(dec_rad - tde_rad)
C2 = sin(tde_rad)
C3 = cos(dec_rad)
C4 = C2 * C3
C5 = C1 + C4
C6 = cos(tde_rad)
a = C3 * Fs_millimeters
b = C1 * C3 * Fs_millimeters
c = Fs_millimeters * (C5 - C4)
Rmax = Fs_millimeters * sin(BetaLimit_rad)
t1 = a**2 * (Rmax**2 - c**2 - (2 * b * c))
t2 = sqrt(a**4 - t1 + (b* Rmax)**2)
t1 = a**2 * (b + c) - b * t2
Dlim = t2 / (a**2 - b**2)
if abs(Dlim) > abs(Rmax):
# avoids a negative square root in the calculation of Rlim
#print('Error Att: no trajectory available for dec %.4f and azimuth %.4d' % (dec_deg, az_deg))
return 0.0
Rlim = sqrt(Rmax**2 - Dlim**2)
Hmax = asin( Rlim / (Fs_millimeters *C3))
# att in minutes = distance / velocity
att = (2 * Hmax) * (86164.091 / (2 * pi * 60))
return att
#
# Test functions in this module
#
if __name__ == '__main__':
def space():
print('')
print('test modulo() functions')
print('clipToCircle_deg(361.0) is %f'%(clipToCircle_deg(361.0)))
print('clipToCircle_deg(360.0) is %f'%(clipToCircle_deg(360.0)))
print('clipToCircle_deg(359.0) is %f'%(clipToCircle_deg(359.0)))
space()
print('clipToCircle_deg(1.0) is %f'%(clipToCircle_deg(1.0)))
print('clipToCircle_deg(0.0) is %f'%(clipToCircle_deg(0.0)))
print('clipToCircle_deg(-1.0) is %f'%(clipToCircle_deg(-1.0)))
space(); space()
TwoPi = 2.0 * pi
print('clipToCircle_rad(TwoPi + 0.001) is %f'%(clipToCircle_rad(TwoPi + 0.001)))
print('clipToCircle_rad(TwoPi) is %f'%(clipToCircle_rad(TwoPi)))
print('clipToCircle_rad(TwoPi - 0.001 is %f'%(clipToCircle_rad(TwoPi - 0.001)))
space()
print('clipToCircle_rad(1.0) is %f'%(clipToCircle_rad(0.0001)))
print('clipToCircle_rad(0.0) is %f'%(clipToCircle_rad(0.0)))
print('clipToCircle_rad(-1.0) is %f'%(clipToCircle_rad(-0.001)))
space(); space()
print(' az Tde Ho Pa')
for az in [0.0, 45, 64, 90, 180, 270, 's']:
if type(az) is not float and type(az) is not int:
az_str = str(az)
else:
az_str = '%8.4f' % az
t = Tde(az)
if type(t) is float:
t_str = '%8.4f' % t
else:
t_str = str(t)
h = Ho(az)
if type(h) is float:
h_str = '%8.4f' % h
else:
h_str = str(h)
p = Pa(az)
if type(p) is float:
p_str = '%8.4f' % p
else:
p_str = str(p)
print( '%8s %8s %8s %8s' % (az_str, t_str, h_str, p_str ))
print('\n Dec Az Att')
for dec in [-20, -13, -12, -10, -8, -5, -4, 0, 63, 64, 65, 66, 70, 71, 80, 's']:
b = BestAzimuth(dec)
att = Att( b, dec )
if b is not None:
b_str = '%8.4f' % b
else:
b_str = str(b)
if type(dec) is not float and type(dec) is not int:
dec_str = str(dec)
else:
dec_str = '%8.4f' % dec
att_str = '%8.4f' % att
print( '%8s %8s %8s' % (dec_str, b_str, att_str) )