Документ взят из кэша поисковой машины. Адрес оригинального документа : http://www.naic.edu/~phil/hardware/byuPhasedAr/rfischercode/array_pointingV0.py Дата изменения: Wed Jan 27 16:46:07 2010 Дата индексирования: Sat Jun 26 01:51:18 2010 Кодировка: Поисковые слова: arp 220

from slalib import *
import numpy as np
import string as st
import time as tm

# HA, Dec -> Az, El
def hadec_to_azel ( ha_hrs, dec_deg ) :
ha_rad = np.pi * ha_hrs / 12.0
dec_rad = np.pi * dec_deg / 180.0
vec = sla_dcs2c(ha_rad, dec_rad)
desig, name, of_long, of_lat, elev = sla_obs(0, 'GBVA140')
#print of_lat * 180 / np.pi
rmat = sla_deuler('yxx', np.pi / 2.0 - of_lat, 0.0, 0.0)
new_vec = sla_dmxv(rmat, vec)
az, el = sla_dcc2s(new_vec)
return {'az':180.0 * (1.0 + az / np.pi), 'el':180.0 * el / np.pi}

#print hadec_to_azel(0.0, 0.0)
#print hadec_to_azel(-6.0, 0.0)

# 'xx:yy:zz' to xx.xxxx
def str2frac ( dt_str ) :
parsed = st.split(dt_str, ':')
sign = 1.0
if parsed[0][0] == '-' :
sign = -1.0
return sign * (abs(st.atof(parsed[0])) + (st.atof(parsed[1]) +
st.atof(parsed[2]) / 60.0) / 60.0)

#print str2frac('11:22:33')
#print str2frac('-11:22:33')

# J2000 to epoch
def j2000_to_epoch ( ra_str, dec_str, mjd ) :
rmat = sla_prenut(2000, mjd)
ra = np.pi * str2frac(ra_str) / 12.0
dec = np.pi * str2frac(dec_str) / 180.0
vec = sla_dcs2c(ra, dec)
new_vec = sla_dmxv(rmat, vec)
pra, pdec = sla_dcc2s(new_vec)
return {'ra_hrs':12.0 * pra / np.pi, 'dec_deg':180.0 * pdec / np.pi}

#print 'Precessed:', j2000_to_epoch('1:00:00', '1:00:00', 54630)

def get_current_mjd ( ) :
"""
This only works in 2008.
"""
frac_sec1 = 1.0; frac_sec2 = 0.0
while (frac_sec1 > frac_sec2) :
ftime1 = tm.time()
gmt = tm.gmtime()
ftime2 = tm.time()
frac_sec1 = ftime1 - int(ftime1)
frac_sec2 = ftime2 - int(ftime2)
yr = gmt[0]
mjd = 54465 + gmt[7] + (gmt[3] + (gmt[4] + gmt[5] / 60.0) / \
60.0) / 24.0 + (frac_sec1 + frac_sec2) / 172800.0
for i in range(2009, 2052) :
if i > yr :
break
if (i - 1) % 4 == 0 :
mjd += 366
else :
mjd += 365
return mjd

#print get_current_mjd()

def get_object_posn ( object ) :
if object == 'CasA' :
object = ''
elif object == 'CygA' :
object = ''
pos = st.split(object, ',')
if len(pos) != 2 :
print 'make_track(): Object position error', object
return False
return j2000_to_epoch(pos[0], pos[1], get_current_mjd())

#print get_object_posn('11:22:33,-5:33:77')

def map_azel_to_azel ( map_az_deg, map_el_deg, src_az_deg, src_el_deg ) :
"""
map_azel_to_azel ( map_az_deg, map_el_deg, src_az_deg, src_el_deg ) :
map_az, map_el -> Az, El
positive map_az drives the telescope to smaller azimuth
positive map_el drives the telescope to lower elevation
all arguments in degrees
"""
map_az_rad = np.pi * map_az_deg / 180.0
map_el_rad = np.pi * map_el_deg / 180.0
src_az_rad = np.pi * src_az_deg / 180.0
src_el_rad = np.pi * src_el_deg / 180.0
vec = sla_dcs2c(map_az_rad, map_el_rad)
#print vec
elg = src_el_rad - map_el_rad
azg = src_az_rad - map_az_rad / np.cos(elg)
err_tol = 2.0 * np.pi / (180.0 * 3600.0)
num_iter = 10
for i in range(num_iter) :
rmat = sla_deuler('yzx', elg, -azg, 0.0)
new_vec = sla_dmxv(rmat, vec)
#print new_vec
az, el = sla_dcc2s(new_vec)
#print i, '%10.5f %10.5f %10.5f %10.5f' % \
# (azg * 180 / np.pi, elg * 180 / np.pi,
# az * 180 / np.pi, el * 180 / np.pi)
az_diff = src_az_rad - az
if az_diff > np.pi :
az_diff -= 2.0 * np.pi
el_diff = src_el_rad - el
#print az_diff, el_diff
if abs(az_diff) < err_tol and abs(el_diff) < err_tol :
break
azg += az_diff
elg += el_diff
if i >= num_iter - 1 :
print 'map_azel_to_azel(%4.2f, %4.2f, %4.2f, %4.2f) did not converge' % (map_az_deg, map_el_deg, src_az_deg, src_el_deg)
if azg < 0.0 :
azg += 2.0 * np.pi
if azg > 2.0 * np.pi :
azg -= 2.0 * np.pi
return {'tel_az':180.0 * azg / np.pi, 'tel_el':180.0 * elg / np.pi}

#for az in np.arange(0.0, 361, 5.0) :
# for el in np.arange(0.0, 79.0, 5.0) :
# map_azel_to_azel(10.0, -10.0, az, el)

#print map_azel_to_azel(5.0, 5.0, 0.0, 5.0), '<- result'

DUT1 = -0.43 # UT1 - UTC in seconds UT1 = UTC + DUT1

def get_last_str ( mjd, utc_string ) :
utc_parsed = st.split(utc_string, ':')
utc_frac = str2frac(utc_string) / 24.0
mjd_frac = mjd + utc_frac
#desig, name, gbt_long, lat, elev = sla_obs(0, 'GBT')
desig, name, of_long, lat, elev = sla_obs(0, 'GBVA140')
lon_frac = of_long / (2.0 * np.pi)
gast_rad = sla_gmst(mjd_frac + DUT1 / 86400.0) + sla_eqeqx(mjd_frac)
return (gast_rad - of_long) / (2.0 * np.pi)

#print (get_last_str(54630.0, '0:0:0') - 0.50511288448786651) * 86400.0
#print (get_last_str(54630.0, '12:12:12') - 0.014977253263429011) * 86400.0

def get_last ( mjd_frac ) :
#desig, name, gbt_long, lat, elev = sla_obs(0, 'GBT')
desig, name, of_long, lat, elev = sla_obs(0, 'GBVA140')
lon_frac = of_long / (2.0 * np.pi)
gast_rad = sla_gmst(mjd_frac + DUT1 / 86400.0) + sla_eqeqx(mjd_frac)
return (gast_rad - of_long) / (2.0 * np.pi)

#print (get_last(54630.5084722222222) - 0.014977253263429011) * 86400.0


def make_track ( object, start_mjd, stop_mjd, cross_el_offset, el_offset ) :
"""
make_track ( object, start_mjd, stop_mjd, cross_el_offset, el_offset )
object = 'CasA', CygA', or 'hh:mm:ss.s,sdd:mm:ss' J2000 coordinates
start_mjd, stop_mjd include UTC, e.g., 54630.3456
positive cross_el_offset drives the telescope to smaller azimuth
positive el_offset drives the telescope to lower elevation
both arguments in degrees
"""
posn = get_object_posn(object)
if posn == False :
return
start_last = get_last(start_mjd)
start_ha = 24.0 * start_last - posn['ra_hrs']
if start_ha < -12.0 :
start_ha += 24.0
if start_ha > 12.0 :
start_ha -= 24.0
stop_last = get_last(stop_mjd)
stop_ha = 24.0 * stop_last - posn['ra_hrs']
if stop_ha < -12.0 :
stop_ha += 24.0
if stop_ha > 12.0 :
stop_ha -= 24.0
if stop_ha < start_ha :
stop_ha += 24.0

if stop_last < start_last :
stop_last += 1.0
num_iter = 3
max_num_pts = 2**num_iter
pt_mjd = np.zeros(max_num_pts + 1, dtype=np.float64)
pt_az = np.zeros(max_num_pts + 1, dtype=np.float64)
pt_az_vel = np.zeros(max_num_pts, dtype=np.float64)
pt_el = np.zeros(max_num_pts + 1, dtype=np.float64)
pt_el_vel = np.zeros(max_num_pts, dtype=np.float64)
posn_tol = 2.0 / 3600.0 # degrees
for iter in range(num_iter) :
num_pts = 2**(iter + 1)
ha_step = (stop_ha - start_ha) / num_pts
for pt in range(num_pts + 1) :
azel = hadec_to_azel(start_ha + ha_step * pt, posn['dec_deg'])
tel_azel = map_azel_to_azel(cross_el_offset, el_offset,
azel['az'], azel['el'])
pt_mjd[pt] = start_mjd + pt * ha_step * 0.997269 / 24.0
pt_az[pt] = tel_azel['tel_az']
pt_el[pt] = tel_azel['tel_el']
tol_flag = True
for pt in range(1, num_pts + 1, 2) :
az_err = pt_az[pt] - (pt_az[pt - 1] + pt_az[pt + 1]) / 2.0
az_err *= np.cos(pt_el[pt] * np.pi / 180.0)
el_err = pt_el[pt] - (pt_el[pt - 1] + pt_el[pt + 1]) / 2.0
err = np.sqrt(az_err**2 + el_err**2)
if err > posn_tol :
tol_flag = False
if tol_flag :
break
if not tol_flag :
print 'WARNING: make_track() did not converge!'
num_pts_out = num_pts / 2
step_sec = 2.0 * ha_step * 0.997269 * 3600.0 # seconds
for pt in range(num_pts_out) :
pt_mjd[pt] = pt_mjd[pt * 2]
pt_az[pt] = pt_az[pt * 2]
pt_el[pt] = pt_el[pt * 2]
pt_az_vel[pt] = (pt_az[pt * 2 + 2] - pt_az[pt * 2]) / step_sec
pt_el_vel[pt] = (pt_el[pt * 2 + 2] - pt_el[pt * 2]) / step_sec
return {'mjd':pt_mjd[:num_pts_out],
'az':pt_az[:num_pts_out], 'az_vel':pt_az_vel[:num_pts_out],
'el':pt_el[:num_pts_out], 'el_vel':pt_el_vel[:num_pts_out]}

print make_track('19:22:33,65:33:77', 54642.7803896, 54642.7803896 + 30.0/86400, 5.0, 5.0)