#-----------------------------------------------------
# gblst() returns the Green Bank LST in hours.
# this returns the mean sidereal time.
#-----------------------------------------------------
# get the time stuff and slalib
from datetime import datetime
import slalib
import math
def gblst() :
gbtlong = 79.839840 # GBT west longitude in degrees.
now = datetime.utcnow() # get now date and time
# get MJD from slalib routine
mjd = slalib.sla_cldj(now.year, now.month, now.day)
# get time as fraction of day
tfrac = ( now.hour + ((now.minute + (now.second/60.0))/60.0)) /24.0
# get the GMST in degrees
gmst = (180.0/math.pi)*slalib.sla_gmsta(mjd, tfrac)
# rotate to Green Bank and convert to hours
gblst = (gmst-gbtlong)/15.0
# avoid slightly negative problems
if -0.0001 <= gblst < 0.0001 : gblst = 0 # within 0.5 second means zero
# normalize to be within 0 to 24 hours
if gblst >= 24 : gblst = gblst-24
if gblst < -0.0001 : gblst = gblst+24
return gblst
#----------------------------------------------------------
# lst2ut( hr, min, sec, year, month, day)
# (hr,min,sec) is the input LST time.
# (year, month, day) is the input UT date.
# Returns the UT (in hours) for the given date and LST.
# if year=0, use the current UT date.
#----------------------------------------------------------
from datetime import datetime
import slalib
import math
lst2ut( hr, min, sec, year, month, day) :
gbtlong = 79.839840 # GBT west longitude in degrees.
now = datetime.utcnow() # get now date and time
# get MJD from slalib routine
mjdnow = slalib.sla_cldj(now.year, now.month, now.day)
# if year < 10, use present date
if year < 10 : mjd = mjdnow
else : mjd = slalib.sla_cldj(year, month, day)
# get time as fraction of day
tfracn = ( now.hour + ((now.minute + (now.second/60.0))/60.0)) /24.0
ddnow = mjdnow[0] + tfracn # current date as MJD+Fraction
# convert requested LST to degrees
lstin = 15.0*( hr + (min + (sec/60.0))/60.0)
# get the local LMST for 0h UT in degrees
lst0 = (180.0/math.pi)*slalib.sla_gmst(mjd) - gbtlong
# normalize lst0 to be in (0,360) :
if -0.004 <= lst0 <= 0.004 : lst0 = 0 # less than one second means zero
if lst0 >= 360 : lst0 = lst0 - 360
if lst0 < -0.004 : lst0 = lst0 + 360
# get LST difference between 0h UT and requested LST
tdiff = lstin - lst0
if -0.004 <= tdiff <= 0.004 : tdiff = 0 # within one second means zero
if tdiff < -0.004 : tdiff = tdiff + 360 # make it be between 0 and 360
if tdiff >= 360 : tdiff = tdiff - 360
rr = (365.25/366.25) # ratio of solar to sidereal day
ddif = 24.0*(1.0 - rr) # difference in length of sid day and sol day
# convert to solar time in hours
tsol = (tdiff/15.0)*rr
# check for ambiguous times
if tsol < ddif : # if true there is an ambiguity in the UT time
tsolx = ((tdiff+360)/15.0)*rr # alternate UT time
d1 = mjd[0] + tsol/24.0
d2 = mjd[0] + tsolx/24.0
# if present time is earlier than d1, use d1
if ddnow <= d1 : return tsol
# if present time is more than 3 hours past d1, use d2
if ddnow > d1+0.125 : return tsolx
return tsol
#------------------------------------------------------ # rcorr(elev) returns refraction correction in degrees # input is elevation in degrees # Using Ron's handy approximate formula #------------------------------------------------------ import math dtr = math.pi/180.0 def rcorr( elev) : aa = elev + (4.56934/(2.19001 + elev)) del_el = (62.2/3600.0)/( math.tan( aa*dtr)) return del_el
#----------------------------------------------------------
# hahor ( elev, dec ) returns hour angle at which a source
# at declination dec intersects the elevation = elev
# Inputs elev and dec are in degrees.
# the input elevation is Eenc (encoder elevation)
# The elevation is corrected for refraction.
# returned value is in hours.
#----------------------------------------------------------
import math
dtr = pi/180.0 # convert degrees to radians
gbtlat = 38.4331294 # GBT latitude in degrees
slat = math.sin(gbtlat*dtr)
clat = math.cos(gbtlat*dtr)
hahor( elev, dec) :
delev = rcorr(elev) # get refraction correction
e2 = elev - delev
sinelv = math.sin(e2 * dtr)
sdec = math.sin(dec * dtr)
cdec = math.cos(dec * dtr)
cosH = (sinelv - sdec*slat)/(cdec*clat)
if cosH > 1.0 : # source never rises
return 0.0
if cosH < -1.0 : # source is circumpolar
return 12.0
dh = math.acos( cosH)/dtr # result in degrees
DH = math.abs(dh)/15 # result in hours
return DH
#----------------------------------------------------------
# riset( ra, dec, elev) get next rise and set times in LST
# for a source at apparent position (ra,dec),
# and for a horizon at elevation=elev.
# Input angles : ra in hours, dec in degrees.
# Output times in hours of LST
# Return value is ( LSTrise, LSTset )
#----------------------------------------------------------
riset( ra, dec, elev) :
sha = gblst() - ra # current HA of source
if sha < -12 : sha = sha + 24.0
if sha > 12 : sha = sha - 24.0 # HA must be between -12 and +12 hours.
dha = hahor( elev, dec) # get HA at horizon
if dha >= 12.0 : # its circumpolar
return (gblst(), None)
if dha <= 0 : # never rises
return (None, gblst())
# rise LST is when the HA = -dha
# set time is when HA = +dha
lstrise = ra -dha
lstset = ra + dha
return (lstrise, lstset)
#----------------------------------------------------------
# cconvert.py -- convert coordinates to "GAPPT",
# i.e., Apparent RA/DEC.
#--------------------------------------------------
#
# The above routines, "hahor" and "riset" require apparent
# coordinates.
# This routine uses SLALIB to convert from input coordinates
# to Apparent.
#
# Input: lon and lat are longitude-like and latitude-like
# coordinates, i.e., (ra,dec), or (glon, glat).
# input units are degrees.
# ctype is one of "J2000", "B1950", "JMEAN", "GAPPT",
# or "GALACTIC"
# if ctype="JMEAN", then equinox must be specified.
#
# Output: apparent ra and dec in radians for the present time.
#--------------------------------------------------
#
# (if not in astrid, source /home/gbt7/sparrow/release/sparrow.bash )
#
import math
import slalib
from datetime import datetime
def cconvert( lon, lat, ctype, equinox) :
dtr = math.pi/180.0
rr = lon*dtr
dd = lat*dtr
if ctype == "J2000" :
eq = 2000.0
elif ctype == "B1950" :
eq = 1950.0
elif ctype == "JMEAN" :
eq = equinox
elif ctype == "GALACTIC" :
# convert galactic to J2000
glc = slalib.sla_galeq( lon*dtr, lat*dtr)
rr = glc[0]
dd = glc[1]
eq = 2000.0
elif ctype == "GAPPT" :
return(lon,lat) # no conversion necessary!
else :
print "Unknown ctype"
return (0, 0)
# get todays date in MJD
now = datetime.utcnow()
mjd1 = slalib.sla_cldj(now.year, now.month, now.day)
tfrac = (now.hour + ((now.minute + (now.second/60.0))/60.0))/24.0
mjd = mjd1[0] + tfrac
# print "eq,mjd=", eq, mjd, "rr,dd=", rr,dd
# get conversion parameters
amprms = slalib.sla_mappa( eq, mjd)
# convert to apparent place
rada = slalib.sla_mapqkz( rr, dd, amprms) # result in radians
# print "cconvert:", rada
rr = slalib.sla_cr2tf( 2, rrdd[0]) # convert to hh:mm:ss
dd = slalib.sla_cr2af( 2, rrdd[1]) # convert to dd:mm:ss
print rr, dd # print results as hh mm ss .ss dd mm ss .ss
return rada
#----------------------------------------------------------
-- FrankGhigo - 4 May 2006
| Topic HorizonCodePrototypes . { Edit | Attach | Ref-By | Printable | Diffs | r1.7 | > | r1.6 | > | r1.5 | More } |
|
Revision r1.7 - 05 May 2006 - 13:56 GMT - FrankGhigo Parents: PlanOfRecordC32006 > ModificationRequest14C306 |
Content copyright © 1999-2007 by the contributing authors. All material on this collaboration platform is the property of the contributing authors. |