Документ взят из кэша поисковой машины. Адрес оригинального документа : http://www.astro.louisville.edu/software/sbig/archive/xmccd-4.1/xmccd-4.1e/scripts/dome_track.py Дата изменения: Thu Nov 29 07:36:11 2012 Дата индексирования: Thu Feb 27 22:05:48 2014 Кодировка:

#!/usr/bin/python

# * dome_track
#
# * Look up the right ascension and declination of the telescope field
# * Calculate the optimal azimuth of the dome shutter
# * Look up the current azimuth of the dome shutter
# * Drive the dome to the optimal azimuth

# Copyright (c) 2010-2012 John Kielkopf, Jeff Hay and Karen Collins
# kielkopf@louisville.edu
#
#
# Date: April 19, 2012
# Version: 1.1
#
# This file is part of PyDome
#
# PyDome is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# PyDome is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.



# Define values appropriate for this dome


# Switch settings for the DLI device

SWITCHIP="192.168.0.100"
USER=""
PASSWORD=""

# Rotation rate in degrees per second

dps = 2.5

# Azimuth tolerance for setting should be greater than tag increment

tolazimuth = 5.

# Lower limit of azimuth allowed

minazimuth = 0.

# Upper limit of azimuth allowed

maxazimuth = 270.


# Longitude is + west of prime meridian

# Longitude and latitude for Moore Observatory
# site_longitude = 85.5288888888
# site_latitude = 38.3334

# Longitude and latitude for Mt. Kent Observatory
# site_longitude = -151.855528
# site.latitude = -27.797778

site_longitude = 85.5288888888
site_latitude = 38.3334


# Radius of dome

# CDK20N and CDK20S
# dome_radius = 1.75

# RC24
# dome_radius = 2.83

dome_radius = 1.75


# Offset to dec axis from center of dome

# CDK20N
# mount_dec_offset = 0.24

# CDK20S
# mount_dec_offset = 0.24

# RC24
# mount_dec_offset =0.66

mount_dec_offset = 0.20


# Height of intersection of declination and polar axes above dome equator

# CDK20N
# mount_dec_height = 0.0

# CDK20S
# mount_dec_height = 0.0

# RC24
# mount_dec_height = -0.25

mount_dec_height = 0.0


# Length from polar to optical axes
# CDK20N
# mount_dec_length = 0.6

# CDK20S
# mount_dec_length = 0.6

# RC24
# mount_dec_length = 0.

mount_dec_length = 0.65


# Mount type altaz = 0, fork = 1, german equatorial = 2

telmount = 2


################################################################################
################################################################################

import sys
import argparse
import time
from datetime import datetime
from math import sin, cos, asin, acos, atan2, sqrt, pi

# The dli module needs to be in the python path.
# Add current directory to the path so we can find the dli module

sys.path.append('.')
import dli


parser= argparse.ArgumentParser(description='Track telescope with dome azimuth')

if len(sys.argv) > 1:
sys.exit("Usage: dome_track")
exit()


# Map a time in hours to the range 0 to 24

def map24(hour):

if (hour < 0.0):
n = int(hour / 24.0) - 1
hour24 = hour - n * 24.0
return (hour24)

elif (hour >= 24.0):
n = int(hour / 24.0)
hour24= hour - n * 24.0
return (hour24)

else:
hour24 = hour
return (hour24)


#Map an hourangle in hours to -12 <= ha < +12

def map12(hour):

hour12=map24(hour)

if (hour12 >= 12.0):
hour12 = hour12 - 24.
return (hour12)

else:
return (hour12)




# Map an angle in degrees to 0 <= angle < 360

def map360(angle):

if (angle < 0.0):
n = int(angle / 360.0) - 1
return (angle - float(n) * 360.0)

elif (angle >= 360.0):
n = int(angle / 360.0)
return (angle - float(n) * 360.0)

else:
return (angle)


# Map an angle in degrees to -180 <= angle < 180

def map180(angle):

angle360 = map360(angle)

if (angle360 >= 180.0):
return (angle360 - 360.0)

else:
return (angle360)


# Find the fractional part handling negative values properly

def frac(x):

newx = x - int(x)
if (newx < 0.):
newx = newx + 1.

return (newx)


# Calculate the current universal time in hours with millisecond resolution

def utnow():

dt = datetime.utcnow()

hours = dt.hour
minutes = dt.minute
seconds = dt.second
microseconds = dt.microsecond
ut = hours + minutes/60. + seconds/3600. + microseconds/3600000000.

return (ut)


# Calculate the Julian day number for a given year, month, day, and universal time

def calcjd(year,month,day,ut):

day = day + ut / 24.0


if int(month)== 1 or int(month)== 2:
year = year - 1.
month = month + 12.

if (year + month/12.0 + day/365.25) >= (1582.0 + 10.0 / 12.0 + 15.0 / 365.25):
na = int(year/100.)
a = float(na)
b = 2.0 - a + int(a/4.0)
else:
b = 0.

if year < 0.:
nc = int(float(int(365.25 * year)) - 0.75)
c = float(nc)
else:
nc = int(365.25 * year)
c = float(nc)

nd = int(30.6001 * (month + 1.))
d = float(nd)
jd = b + c + d + day + 1720994.5

return (jd)


# Calculate the Julian day number at this moment

def jdnow():
dt = datetime.utcnow()

year = dt.year
month = dt.month
day = dt.day
hours = dt.hour
minutes = dt.minute
seconds = dt.second
microseconds = dt.microsecond
ut = hours + minutes/60. + seconds/3600. + microseconds/3600000000.
jd = calcjd(year, month, day, ut)

return (jd)

# Calculate the sidereal time for a year, month, day and universal time at the define site latitude

def calclst(year, month, day, ut, glong):
jdeod = calcjd(year, month, day, 0.0)
tu = (jdeod - 2451545.0) / 36525.0
t0 = (24110.54841 / 3600.0) + (8640184.812866 / 3600.0) * tu + (0.093104 / 3600.0) * tu * tu - (6.2e-6 / 3600.0) * tu * tu * tu
t0 = map24(t0)
gmst = map24(t0 + ut * 1.002737909)
lmst = 24.0 * frac( (gmst - (glong / 15.0) ) / 24.0)

return (lmst)


# Calculate the sidereal time at this moment for the defined site latitude

def lstnow():

dt = datetime.utcnow()

year = dt.year
month = dt.month
day = dt.day
hours = dt.hour
minutes = dt.minute
seconds = dt.second
microseconds = dt.microsecond
ut = hours + minutes/60. + seconds/3600. + microseconds/3600000000.
glong = site_longitude
lst = calclst(year,month,day,ut,glong)

return (lst)


# Convert local ha and dec to local az and alt
# Geographic azimuth convention is followed:
# Due north is zero and azimuth increases from north to east

def equatorial_to_horizontal(ha, dec):

ha = ha*pi/12.
phi = site_latitude*pi/180.
dec = dec*pi/180.
altitude = asin(sin(phi)*sin(dec)+cos(phi)*cos(dec)*cos(ha))
altitude = altitude*180.0/pi
azimuth = atan2(-cos(dec)*sin(ha),
sin(dec)*cos(phi)-sin(phi)*cos(dec)*cos(ha))
azimuth = azimuth*180.0/pi
azimuth = map360(azimuth)

return (azimuth, altitude)


# Iteratively solve for the azimuth of the dome given the telescope RA and Dec

def solve_dome_azimuth(ra, dec, nloops):

# Initialize the dome and mount parameters

R = dome_radius
H = mount_dec_height
L = mount_dec_length
P = mount_dec_offset

# Find the current hour angle for the telescope pointing

ha = lstnow()
ha = ha - ra
ha = map24(ha)
ha0 = map12(ha)

# Find the altitude and azimuth of the current pointing
# Function requires ha in 0 to 24 basis
# This should be valid in either hemisphere

telaz, telalt = equatorial_to_horizontal(ha, dec)

# Find the reference point on the optical axis in dome coordinates
# z: vertical
# y: north - south
# x: east - west
# theta: altitude of the polar axis from the horizontal plane
# phi: rotation of dec axis about polar axis

# Meaning of x and y compared to sky depends on the hemisphere!
# z: always + up
# y: always + toward N for +lat or S for -lat
# x: maintains right-handed coordinate system with z and y
# thus x is + to E for +lat and to W for -lat
# theta: always positive from the horizontal plane to the pole
# phi: rotation about polar axis
# and is +90 for a horizontal axis with OTA toward +x
# and is 0 for OTA over mount with dec axis counterweight down
# and is -90 for a horizontal axis with OTA toward -x


# Set default parameters

if (telmount != 2) :

# We have an altaz or fork mount and the origin is simple

x0 = 0.
y0 = P
z0 = H

else:

# We have a German equatorial and the origin changes with ha

# Assign phi based on an assumption about the basis derived from the ha

if (site_latitude >= 0.):

if (ha0 > 0.):

# Looking west with OTA east of pier
# HA is between 0 and +12 hours

phi = (6. - ha0)*pi/12.


else:

# Looking east with OTA west of pier
# HA is between -12 and 0 hours

phi = -(6. + ha0)*pi/12.



else:

if (ha0 > 0.):

# Looking west with OTA east of pier
# HA is between 0 and +12 hours

phi = -(6. - ha0)*pi/12.


else:

# Looking east with OTA west of pier
# HA is between -12 and 0 hours

phi = (6. + ha0)*pi/12.


# Assign theta

if (site_latitude >= 0.):

theta = site_latitude*pi/180.

else:

theta = -1.*site_latitude*pi/180.


# Find the dome coordinates of the OTA reference point for a German equatorial
# This works in either hemisphere

x0 = L*sin(phi)
y0 = -L*cos(phi)*sin(theta) + P
z0 = L*cos(phi)*cos(theta) + H


# (x,y,z) is on the optical axis
# Iterate to make this point also lie on the dome surface
# Begin iteration assuming the zero point is at the center of the dome
# Telescope azimuth is measured from the direction to the pole

d = 0.
r = R

if (site_latitude >= 0.):

telaz2 = telaz * pi / 180.
telalt2 = telalt * pi / 180.

else:

telaz2 = telaz - 180.
telaz2 = map360(telaz2)
telaz2 = telaz2 * pi / 180.
telalt2 = telalt * pi / 180.


# Iterate for convergence

n = 0
while (n < nloops):
d = d - (r - R)
rp = R + d
x = x0 + rp * cos(telalt2) * sin(telaz2)
y = y0 + rp * cos(telalt2) * cos(telaz2)
z = z0 + rp * sin(telalt2)
r = sqrt(x*x + y*y + z*z)

# Repeat the calculation correcting the assumed OTA path length each time
# This converges quickly so just do a few loops and assume it is ok

# print "n, r, rp: ", n, r, rp

n = n + 1

# Use (x,y,0) from the interation to find the azimuth of the dome
# Azimuth is N (0), E (90), S (180), W (270) in both hemispheres
# However x and y are different in the hemispheres so we fix that here

zeta = atan2(x,y)
if ( (zeta > - 2.*pi) and (zeta < 2.*pi)):

if ( site_latitude > 0. ):

zeta = (180./pi)*zeta
zeta = map360(zeta)

else :

zeta = (180./pi)*zeta
zeta = zeta + 180.
zeta = map360(zeta)


else:

zeta = telaz

return (zeta)

# Look up the current telescope coordinates

telcoordsstatus = open('/usr/local/observatory/status/telcoords','r')
telcoordsstr = telcoordsstatus.read(40)
telcoordsstatus.close()

# print ("Coordinates: %s " % telcoordsstr)

rastr, decstr = telcoordsstr.split()
ra = float(rastr)
dec = float(decstr)

# Iteratively solve for the optimal dome azimuth

nloops = 5
newazimuth = solve_dome_azimuth(ra, dec, nloops)

print ("New azimuth: " '{0:.2f}'.format(newazimuth))

# Look up the current dome azimuth

nowazimuthstatus = open('/usr/local/observatory/status/domeazimuth','r')
nowazimuthstr = nowazimuthstatus.read(40)
nowazimuthstatus.close()
nowazimuth = float(nowazimuthstr)

print ("Current azimuth: " '{0:.2f}'.format(nowazimuth))

# Limit motion if needed

#if minazimuth > newazimuth:
# newazimuth = minazimuth
# print ("Limiting dome track to minimum azimuth ...")

#if maxazimuth < newazimuth:
# newazimuth = maxazimuth
# print (" Limiting dome track to maximum azimuth ...")


# Calculate the required change in azimuth

# Use this if crossing the zero point is allowed

#delazimuth = map180(newazimuth - nowazimuth)

# Use this if zero crossing is forbidden

delazimuth = newazimuth - nowazimuth

if abs(delazimuth) > tolazimuth:

# Use the DLI switch to rotate the dome

# Open (connect) to the power switch, if the
# parameters are not passed it defaults to the
# IP address, user and password the switch is
# set with from the factory. Only changed items
# actually need to be passed.
switch=dli.powerswitch(hostname=SWITCHIP,userid=USER,password=PASSWORD)

# Verify the switch
if not switch.verify():
print "Cannot talk to the switch"
sys.exit(1)

# Print the current state of all the outlets
# switch.printstatus()

if (delazimuth > 0.):

# Power increase azimuth on
switch.on(2)

# print 'Increase dome azimuth switch is ',switch.status(2)

# Wait the dome to rotate the requested angle
dps = 2.5
delay = delazimuth/dps
time.sleep(delay)

# Power increase dome azimuth off
switch.off(2)

# print 'Increase dome azimuth switch is ',switch.status(2)

elif (delazimuth < 0.):

# Power decrease azimuth on
switch.on(1)

# print 'Decrease dome azimuth switch is ',switch.status(1)

# Wait the dome to rotate the requested angle
delay = -1.*delazimuth/dps
time.sleep(delay)

# Power decrease dome azimuth off
switch.off(1)

# print 'Decrease dome azimuth switch is ',switch.status(1)

else:
pass
# print 'No rotation requested.'

else:
pass
# print "Dome azimuth within tolerance"

exit()