'''A library of general astronomical functions, many of which interface with PyNovas.''' from math import * import novas import sys, os # get OBSPREPDIR environment variable obsprepDir=os.getenv('OBSPREPDIR') if obsprepDir == None: raise EnvironmentError('*** environment variable OBSPREPDIR not defined!') DEBUG_ASTRO = False # True -> print out extra debug messages # False -> do not print out extra debug messages # Initialize ephemeris file # use next line if not using the --prefix option for the install # novas.EphemerisSet('%s/local/novas/JPLEPH' % sys.exec_prefix) # novas.EphemerisSet('/home/scr/bzavala/software/ObsPrep/noiPlan/pnovas-2.0.1b/JPLEPH') novas.EphemerisSet(obsprepDir+'/pnovas-2.0.1b/JPLEPH') def npoi(): '''Return a NOVAS Location structure for the NPOI site.''' site = novas.Location() site.lat = 35.09666 site.lon = -111.535 site.h = 2200.66 return site def nofs(): '''Return a NOVAS Location structure for the NOFS site.''' site = novas.Location() site.lat = 35.18333 site.lon = -111.7405278 site.h = 2200.66 return site def body(bodyName): '''Return a NOVAS Body structure for the specified solar system body. ARGS: bodyName (str): Sun, Moon, or one of the planets ''' bodies = {'Mercury':1, 'Venus':2, 'Earth':3, 'Mars':4, 'Jupiter':5, 'Saturn':6, 'Uranus':7, 'Neptune':8, 'Pluto':9, 'Sun':10, 'Moon':11} bodyName = bodyName.capitalize() if bodyName not in bodies: raise ValueError, 'Unrecognized body "%s"' % bodyName b = novas.Body() b.type = 0 b.number = bodies[bodyName] b.name = bodyName return b def inRange(x,min,max): '''Return the input number with the specified range. For example, inRange(26.2,0,24) = 2.2 ARGS: x (double): Input number min (double): minimum of the range max (double): maximum of the range ''' assert max > min dx = max - min y = x while y < min: y += dx while y > max: y -= dx return y def jd2mjd(jd): '''Return the Modified Julian Date(MJD) for the specified Julian Date(JD). ''' return jd - 2400000.5 def mjd2jd(mjd): '''Return the Julian Date(JD) for the specified Modified Julian Date(MJD). ''' return mjd + 2400000.5 def tai2tt(): '''Return the increment to be applied to International Atomic Time (TAI) to give Terrestial Time (TT), in seconds. TT = TAI + 32.184 sec, always. ''' return 32.184 def utc2tai(mjd): '''Return the increment to be applied to Coordinated Universal Time (UTC) to give International Atomic Time (TAI), in seconds. NOTE: This routine needs to be modified everytime a leap second is added. Code was taken from the slalib function slaDat(). ARGS: mjd (double): Modified Julian Date. ''' # 1999 January 1 if mjd >= 51179.0: return 32.0 # 1997 July 1 if mjd >= 50630.0: return 31.0 # 1996 January 1 if mjd >= 50083.0: return 30.0 # 1994 July 1 if mjd >= 49534.0: return 29.0 # 1993 July 1 if mjd >= 49169.0: return 28.0 # 1992 July 1 if mjd >= 48804.0: return 27.0 # 1991 January 1 if mjd >= 48257.0: return 26.0 # 1990 January 1 if mjd >= 47892.0: return 25.0 # 1988 January 1 if mjd >= 47161.0: return 24.0 # 1985 July 1 if mjd >= 46247.0: return 23.0 # 1983 July 1 if mjd >= 45516.0: return 22.0 # 1982 July 1 if mjd >= 45151.0: return 21.0 # 1981 July 1 if mjd >= 44786.0: return 20.0 # 1980 January 1 if mjd >= 44239.0: return 19.0 # 1979 January 1 if mjd >= 43874.0: return 18.0 # 1978 January 1 if mjd >= 43509.0: return 17.0 # 1977 January 1 if mjd >= 43144.0: return 16.0 # 1976 January 1 if mjd >= 42778.0: return 15.0 # 1975 January 1 if mjd >= 42413.0: return 14.0 # 1974 January 1 if mjd >= 42048.0: return 13.0 # 1973 January 1 if mjd >= 41683.0: return 12.0 # 1972 July 1 if mjd >= 41499.0: return 11.0 # 1972 January 1 if mjd >= 41371.0: return 10.0 # 1968 February 1 if mjd >= 39887.0: return 4.2131700 + (mjd - 39126.0) * 0.002592 # 1966 January 1 if mjd >= 39126.0: return 4.3131700 + (mjd - 39126.0) * 0.002592 # 1965 September 1 if mjd >= 39004.0: return 3.8401300 + (mjd - 38761.0) * 0.001296 # 1965 July 1 if mjd >= 38942.0: return 3.7401300 + (mjd - 38761.0) * 0.001296 # 1965 March 1 if mjd >= 38820.0: return 3.6401300 + (mjd - 38761.0) * 0.001296 # 1965 January 1 if mjd >= 38761.0: return 3.5401300 + (mjd - 38761.0) * 0.001296 # 1964 September 1 if mjd >= 38639.0: return 3.4401300 + (mjd - 38761.0) * 0.001296 # 1964 April 1 if mjd >= 38486.0: return 3.3401300 + (mjd - 38761.0) * 0.001296 # 1964 January 1 if mjd >= 38395.0: return 3.2401300 + (mjd - 38761.0) * 0.001296 # 1963 November 1 if mjd >= 38334.0: return 1.9458580 + (mjd - 37665.0) * 0.0011232 # 1962 January 1 if mjd >= 37665.0: return 1.8458580 + (mjd - 37665.0) * 0.0011232 # 1961 August 1 if mjd >= 37512.0: return 1.3728180 + (mjd - 37300.0) * 0.001296 # 1961 January 1 if mjd >= 37300.0: return 1.4228180 + (mjd - 37300.0) * 0.001296 # Before that. return 1.4178180 + (mjd - 37300.0) * 0.001296 def eqOfEquinox(mjd): '''Return the equation of the equinox (in seconds) for the specified Modified Julian Date. ARGS: mjd (double): TDB Modified Julian Date. ''' out = novas.EarthTilt(mjd2jd(mjd)) return out[2] def gast(mjd): '''Return the Greenwich apparent sidereal time, in hours, for the specified Modified Julian Date. ARGS: mjd (double): UT1 Modified Julian Date (same as UTC within 1 sec) ''' # We're passed in MJD(UT1). We also need MJD(TDB). # TDB and TT are within a mas of each other, so we'll use TT. # UTC and UT1 are within one second of each other; that difference # won't matter when calculating the equation of the equinoxes. mjdTDB = mjd + (utc2tai(mjd) + tai2tt()) / 86400. ee = eqOfEquinox(mjdTDB) jd = mjd2jd(mjd) return novas.SiderealTime(floor(jd), jd - floor(jd), ee) def gmst(mjd): '''Return the Greenwich mean sidereal time, in hours, for the specified Modified Julian Date. ARGS: mjd (double): UT1 Modified Julian Date (same as UTC within 1 sec) ''' jd = mjd2jd(mjd) return novas.SiderealTime(floor(jd), jd - floor(jd), 0.0) def last(mjd, site): '''Return the local apparent sidereal time, in hours, for the specified Modified Julian Date at the specified site. ARGS: mjd (double): UT1 Modified Julian Date (same as UTC within 1 sec) site (Location): NOVAS Location structure specifying the site ''' last = gast(mjd) + site.lon / 15. last = inRange(last,0,24) return last def lmst(mjd, site): '''Return the local mean sidereal time, in hours, for the specified Modified Julian Date at the specified site. ARGS: mjd (double): UT1 Modified Julian Date (same as UTC within 1 sec) site (Location): NOVAS Location structure specifying the site ''' lmst = gmst(mjd) + site.lon / 15. lmst = inRange(lmst,0,24) return lmst def observedPlanet(bodyName, mjd, site, norefract=False): '''Return the observed place (OP) of a solar system body at the specified Modified Julian Date at the specified site. If "norefract" is set to True, then suppress refraction, returning the topocentric place. ARGS: bodyName (string): Name of the solar system body (Sun, Moon, or one of the planets) mjd (double): UT1 Modified Julian Date (within 1 sec of UTC) site (Location): NOVAS Location structure specifying the site norefract (bool): If True, dont apply refraction (return topocentric plate rather than observed place). RETURNS: (OP zenith distance (deg), OP azimuth (deg), OP right ascension (hour), OP declination (deg)) ''' import sys # Validate inputs assert norefract == True or norefract == False assert isinstance(site,novas.Location) # Fetch the object and earth earth = body('Earth') object = body(bodyName) # We need TT - UT1. We'll settle for TT - UTC, which is within 1 sec # of correct. deltat = utc2tai(mjd) + tai2tt() # We need JD(TT) jdTT = mjd2jd(mjd) + deltat / 86400. # Fetch the topocentric place (everything but refraction). # Must change to the novas installation directory to read the ephemeris. saveDir = os.getcwd() #@@@ #@@@os.chdir('%s/local/novas' % sys.exec_prefix) #@@@os.chdir('/home/scr/bzavala/software/ObsPrep/noiPlan/pnovas-2.0.1b/local/novas') o = novas.TopoPlanet(jdTT, object, earth, deltat, site) os.chdir(saveDir) if o[0] > 0: raise ValueError, 'TopoPlanet() returned an error' # Convert to azimuth and elevation if norefract: refract = 0 else: refract = 1 a = novas.Equ2Hor(jdTT, deltat, 0., 0., site, o[1], o[2], refract) return a def observedStar(star, mjd, site, norefract=False): '''Return the observed place (OP) of a star at the specified Modified Julian Date at the specified site. If "norefract" is set to True, then suppress refraction, returning the topocentric place. ARGS: star (Cat): NOVAS "Cat" structure for the star. mjd (double): UT1 Modified Julian Date (within 1 sec of UTC) site (Location): NOVAS Location structure specifying the site norefract (bool): If True, dont apply refraction (return topocentric plate rather than observed place). RETURNS: (OP zenith distance (deg), OP azimuth (deg), OP right ascension (hour), OP declination (deg)) ''' import sys # Validate inputs assert norefract == True or norefract == False assert isinstance(site,novas.Location) assert isinstance(star,novas.Cat) # Fetch the earth earth = body('Earth') # We need TT - UT1. We'll settle for TT - UTC, which is within 1 sec # of correct. deltat = utc2tai(mjd) + tai2tt() # We need JD(TT) jdTT = mjd2jd(mjd) + deltat / 86400. # Fetch the topocentric place (everything but refraction). # Must change to the novas installation directory to read the ephemeris. saveDir = os.getcwd() #@@@os.chdir('%s/local/novas' % sys.exec_prefix) os.chdir('/home/scr/bzavala/software/ObsPrep/obsprep.1.5.2/pnovas-2.0.1b/local/novas') o = novas.TopoStar(jdTT, earth, deltat, star, site) os.chdir(saveDir) if o[0] > 0: raise ValueError, 'TopoStar() returned an error' # Convert to azimuth and elevation if norefract: refract = 0 else: refract = 1 a = novas.Equ2Hor(jdTT, deltat, 0., 0., site, o[1], o[2], refract) return a def apparentStar(star, mjd): '''Return the apparent place (AP) of a star at the specified Modified Julian Date at the specified site. ARGS: star (Cat): NOVAS "Cat" structure for the star. mjd (double): UT1 Modified Julian Date (within 1 sec of UTC) RETURNS: (status (int, 0=ok, 1=error), AP right ascension (hour), AP declination (deg)) ''' import sys # Validate inputs assert isinstance(star,novas.Cat) # Fetch the earth earth = body('Earth') # We need TT - UT1. We'll settle for TT - UTC, which is within 1 sec # of correct. deltat = utc2tai(mjd) + tai2tt() # We need JD(TT) jdTT = mjd2jd(mjd) + deltat / 86400. # Fetch the apparent place. # Must change to the novas installation directory to read the ephemeris. saveDir = os.getcwd() #@@@ if DEBUG_ASTRO: #@@@ print 'novas exec is: %s/local/novas' % sys.exec_prefix if DEBUG_ASTRO: print 'Ready to call novas.AppStar' o = novas.AppStar(jdTT, earth, star) os.chdir(saveDir) if o[0] > 0: raise ValueError, 'AppStar() returned an error' return o def airmass2zd(airmass): '''Return the zenith distance, in degrees, for the specified airmass. ARGS: airmass (double): airmass for which we want the zenith distance ''' scale = 750. # Approximate atmospheric scale height coszd = (2 * scale + 1 - airmass ** 2) / (2 * scale * airmass) return novas.RADTODEG * acos(coszd) def hourAngle(zd, dec, site): '''Returns the hour angle (in hours) for the specified declination, earth latitude, and zenith distance. Returns None if no solution (object never reaches that zenith distance). NOTE: Formula from Astronomical Formulae for Calculators, 3rd edition, Jean Meeus, p. 47. ARGS: zd (double): Zenith distance (degrees) dec (double): Declination (degrees) site (Location): NOVAS Location structure for the observing site ''' # Compute the hour angle cosha = (sin(novas.DEG2RAD * (90. - zd)) - \ sin(novas.DEG2RAD * site.lat) * sin(novas.DEG2RAD * dec)) / \ (cos(novas.DEG2RAD * site.lat) * cos(novas.DEG2RAD * dec)) if cosha >= 1. or cosha <= -1.: return None ha = novas.RAD2DEG * acos(cosha) / 15. return fabs(ha) def zenithDistance(ha, dec, site): '''Returns the zenith distance (in degrees) for the specified hour angle, declination, and earth latitude. Returns None if no solution (object never reaches that hour angle). NOTE: Formula from Astronomical Formulae for Calculators, 3rd edition, Jean Meeus, p. 47. ARGS: ha (hours): Hour angle (hours) dec (double): Declination (degrees) site (Location): NOVAS Location structure for the observing site ''' # Compute the hour angle sinElev = cos(ha * 15. * novas.DEG2RAD) * cos(novas.DEG2RAD * site.lat) * \ cos(novas.DEG2RAD * dec) + \ sin(novas.DEG2RAD * site.lat) * sin(novas.DEG2RAD * dec) if sinElev >= 1. or sinElev <= -1.: return None zd = 90. - asin(sinElev) / novas.DEG2RAD return zd def riseSet(ra,dec,site,timeZone,mjd,zd): '''Returns the local time of rise and set for an object with the specified ra and dec, at the specified site, for the given Modified Julian Date. Rising and setting are defined with respect to the specified zenith distance. Returns None if object doesnt rise or set. It attempts to keep either the rise or set time within -8 to 8 hours local time. NOTE: Formula from Astronomical Formulae for Calculators, 3rd edition, Jean Meeus, p. 47 and Astronomical Almanac 1989, p. A12). ARGS: ra (double): Right ascension (hours) dec (double): Declination (degrees) site (Location): Observing site timeZone (int) Time zone (hr) mjd (int): Modified Julian Date zd (double) Zenith distance defining rising and setting (deg) RETURNS: (rise (hr), set (hr)) ''' assert isinstance(mjd,int) # Compute the hour angle of rise and set ha = hourAngle(zd, dec, site) if ha == None: return None # Return the appropriate times ut2sidereal = 1.002737811906 # Convert UT interval to a sidereal interval rise = ra - ha - last(mjd+timeZone/24., site) rise = inRange(rise,-12,12) / ut2sidereal set = rise + 2 * ha / ut2sidereal return (rise,set) def twilight(site,timeZone,mjd,zd=96.): '''Returns dawn and dusk for the specified site for the specified Modified Julian Date. Dawn and dusk are defined with respect to the specified zenith distance. ARGS: site (Location): Observing site timeZone (int) Time zone (hr) mjd (int): Modified Julian Date zd (double) Zenith distance defining rising and setting (deg) defaults to civil twilight RETURNS: (dusk (hr, zone time; will be negative), dawn(hr, zone time)) ''' assert isinstance(mjd,int) # Iterate 4 times, so we can base dawn and dusk on the Sun's actual # location at those times. dusk = 0 dawn = 0 nIters = 4 for i in xrange(0,nIters): mjd2 = mjd + (timeZone + dawn) / 24. sun = observedPlanet('Sun', mjd2, site) dawn = riseSet(sun[2], sun[3], site, timeZone, mjd, zd)[0] mjd2 = mjd - 1 + (timeZone + dusk) / 24. sun = observedPlanet('Sun', mjd2, site) dusk = riseSet(sun[2], sun[3], site, timeZone, mjd, zd)[1] return (dusk,dawn) def astroTest(): jd=novas.Cal2JD(2004,5,28,0) mjd = jd2mjd(jd) print mjd site=npoi() a=observedPlanet('sun',mjd,site) print a zd = 100. print hourAngle(zd,a[3],site) b=riseSet(a[2],a[3],site,7,mjd,zd) print b b=twilight(site,7,mjd,zd) duskHr = floor(b[0]) duskMin =floor(60.*(b[0]-duskHr)) dawnHr = floor(b[1]) dawnMin =floor(60.*(b[1]-dawnHr)) print '%d:%02d %d:%02d' % (duskHr, duskMin, dawnHr, dawnMin)