diff --git a/obstools.py b/obstools.py new file mode 100644 index 0000000..9e6153f --- /dev/null +++ b/obstools.py @@ -0,0 +1,2583 @@ +#Copyright 2008 Erik Tollerud +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + +""" + +=============================================== +obstools -- miscellaneous tools for observation +=============================================== + +The :mod:`obstools` module stores tools for oberving (pre- and post-) as well as +functioning as a module for various corrections and simple calculations that +don't have a better place to live. + +Most implementations are for optical astronomy, as that is what the primary +author does. + +Note that throughout this module it is assumed that UTC as used in +:mod:`datetime` is the same as UT1 (the UT used for calculations). This is only +an issue if leapseconds stop being updated or if something such as daylight +savings time calculations in :mod:`datetime` are important to be correct to less +than DUT1 - otherwise, the "UTC" in :mod:`datetime` can simply by replaced with +UT1 values when precision better than DUT1 is needed. + +Also note that some of these functions require the `dateutil +`_ package (it is installed by default with +matplotlib) + +.. todo:: examples/tutorials + +.. seealso:: + + :func:`astropysics.coords.earth_rotation_angle` + :func:`astropysics.coords.greenwich_sidereal_time` + + +Classes and Inheritance Structure +--------------------------------- + +.. inheritance-diagram:: astropysics.obstools + :parts: 1 + +Module API +---------- + +""" +#TODO: exposure time calculator (maybe in phot instead?) +#TODO: make Extinction classes spec.HasSpecUnits and follow models framework? +#TODO: return Pipeline support to Extinction +#TODO: Instruments and telescopes for sites +#TODO: get leap second from IERS web site + +from __future__ import division,with_statement +from .constants import pi +from .utils import DataObjectRegistry +from .pipeline import PipelineElement +import numpy as np + +from datetime import timedelta + +try: + from dateutil import tz as tzmod + tzoffset = tzmod.tzoffset +except ImportError: + import datetime + tzmod = None + class tzoffset(datetime.tzinfo): + """ + Backup class to do basic fixed-offset time zones if :mod:`dateutil.tz` + is missing. + """ + + def __init__(self,name,offset): + if not isinstance(name,basestring): + raise TypeError('name must be a string') + self._name = name + + self._hoffset = offset + + def dst(self): + return False + + def tzname(self): + return self._name + + def utcoffset(self): + from datetime import timedelta + return timedelta(hours=self._hoffset) + + +#<----------------------Time and calendar functions----------------------------> +jd2000 = 2451545.0 +mjdoffset = 2400000.5 +""" +Offset between Julian Date and Modified Julian Date - e.g. mjd = jd - mjdoffset +""" + +def jd_to_calendar(jd,rounding=1000000,output='datetime',gregorian=None,mjd=False): + """ + Converts a julian date to a calendar date and time. + + :param jd: + The Julian Date at which to compute the calendar date/time, a sequence + of JDs, or None for the current date/time at the moment the function is + called. + :type jd: scalar, array-like, or None + :param rounding: + If non-0, Performs a fix for floating-point errors. It specifies the + number of milliseconds by which to round the result to the nearest + second. If 1000000 (one second), no milliseconds are recorded. If + larger, a ValueError is raised. + :type rounding: scalar + :param output: + Determines the format of the returned object and can be: + + * 'datetime' + A list of :class:`datetime.datetime` objects in UTC will be + returned. If the input is a scalar, a single object will be + returned. + * 'array' + A Nx7 array will be returned of the form + [(year,month,day,hr,min,sec,msec),...] unless the input was a + scalar, in which case it will be a length-7 array. + * 'fracarray' + An Nx3 array (year,month,day) where day includes the decimal + portion. + + :param gregorian: + If True, the output will be in the Gregorian calendar. Otherwise, it + will be Julian. If None, it will be assumed to switch over on October + 4/15 1582. + :type gregorian: bool or None + :param bool mjd: + If True, the input is interpreted as a modified julian date instead of a + standard julian date. + + :returns: + The calendar date and time in a format determined by the `output` + parameter (see above). + + :except ValueError: + If `rounding` is larger than one second, or `output` is invalid. + + + **Examples** + + >>> jd_to_calendar(2451545) + datetime.datetime(2000, 1, 1, 12, 0, tzinfo=tzutc()) + >>> jd_to_calendar(2305812.5) + datetime.datetime(1600, 12, 31, 0, 0, tzinfo=tzutc()) + >>> jd_to_calendar([2415020.5,2305447.5],output='array') + array([[1900, 1, 1, 0, 0, 0, 0], + [1600, 1, 1, 0, 0, 0, 0]]) + >>> jd_to_calendar(0.0,output='fracarray') + array([[ -4.71200000e+03, 1.00000000e+00, 1.50000000e+00]]) + + """ + import datetime + from dateutil import tz + + if jd is None: + jd = calendar_to_jd(datetime.datetime.now(tz.tzlocal())) + + jd = np.array(jd,copy=True,dtype=float) + scalar = jd.shape == () + jd = jd.ravel() + + if mjd: + jd += mjdoffset + + if rounding > 1000000: + raise ValueError('rounding cannot exceed a second') + elif rounding <= 0: + jd += .5 + else: + rounding = int(rounding) + roundingfrac = rounding/86400000000 + jd += .5 + roundingfrac + + z = np.floor(jd).astype(int) + dec = jd - z #fractional piece + + #fix slight floating-point errors if they hapepn TOOD:check + dgtr1 = dec>=1.0 + dec[dgtr1] -= 1.0 + z[dgtr1] += 1 + + + if gregorian is None: + gregorian = 2299161 + + if gregorian is True: + alpha = ((z-1867216.25)/36524.25).astype(int) + z += 1 + alpha - alpha//4 + elif gregorian is False: + pass + else: + gmask = z >= gregorian + alpha = ((z[gmask]-1867216.25)/36524.25).astype(int) + z[gmask] += 1 + alpha - alpha//4 + + b = z + 1524 + c = ((b-122.1)/365.25).astype(int) + d = (365.25*c).astype(int) + e = ((b-d)/30.6001).astype(int) + + day = b - d - (30.6001*e).astype(int) + + mmask = e<14 + month = e + month[mmask] -= 1 + month[~mmask] -= 13 + year = c + year[month>2] -= 4716 + year[month<=2] -= 4715 + + if output == 'fracarray': + dec = dec-roundingfrac + dec[dec<0]=0 + return np.array((year,month,day+dec)).T + + if rounding == 1000000: + secdec = dec*86400 + sec = secdec.astype(int) + min = sec//60 + sec -= 60*min + hr = min//60 + min -= 60*hr + #sec[sec==secdec] -= 1 + msec = None + else: + msec = (dec*86400000000.).astype('int64') + if rounding > 0: + div = (msec//1000000)*1000000 + toround = (msec - div)<(2*rounding) + msec[toround] = div + rounding + msec -= rounding + + sec = msec//1000000 + msec -= 1000000*sec + + min = sec//60 + sec -= 60*min + hr = min//60 + min -= 60*hr + + if output == 'datetime': + tzi = tz.tzutc() + if msec is None: + ts = (year,month,day,hr%24,min%60,sec%60) + else: + ts = (year,month,day,hr%24,min%60,sec%60,msec%1000000) + res = [datetime.datetime(*t,**dict(tzinfo=tzi)) for t in zip(*ts)] + elif output == 'array': + msec = np.zeros_like(sec) if msec is None else msec + res = np.array([year,month,day,hr%24,min%60,sec%60,msec]).T + else: + raise ValueError('invlid output form '+str(output)) + if scalar: + return res[0] + else: + return res + + + + +def calendar_to_jd(caltime,tz=None,gregorian=True,mjd=False): + + """ + Convert a calendar date and time to julian date. + + :param caltime: + The date and time to compute the JD. Can be in one of these forms: + + * A sequence of floats in the order (yr,month,day,[hr,min,sec]). + * A sequence in the order (yr,month,day,[hr,min,sec]) where at least + one of the elements is a sequence (a sequence will be returned). + * A :class:`datetime.datetime` or :class:`datetime.date` object + * A sequence of :class:`datetime.datetime` or :class:`datetime.date` + objects (a sequence will be returned). + * None : returns the JD at the moment the function is called. + + If the time is unspecified, it is taken to be noon (i.e. Julian Date = + Julian Day Number) + + :param tz: + Sets the time zone to assume for the inputs for conversion to UTC. Can + be any of the following: + + * None + No time zone conversion will occur unless `caltime` is given as + :class:`datetime.datetime` or :class:`datetime.date` objects + with `tzinfo`, in which case they will be converted to UTC using + their own `tzinfo`. + * a string + Specifies a timezone name (resolved into a timezone using the + :func:`dateutil.tz.gettz` function). + * a scalar + The hour offset of the timezone. + * a :class:`datetime.tzinfo` object, + This object will be used for timezone information. + + :param gregorian: + If True, the input will be interpreted as in the Gregorian calendar. + Otherwise, it will be Julian. If None, it will be assumed to switch over + on October 4/15, 1582. + :type gregorian: bool or None + :param bool mjd: + If True, a modified julian date is returned instead of the standard + julian date. + + :returns: JD as a float, or a sequence of JDs if sequences were input. + + + **Examples** + + >>> import datetime,dateutil + >>> calendar_to_jd((2010,1,1)) + 2455198.0 + >>> calendar_to_jd(datetime.datetime(2000,12,21,3,0,0)) + 2451899.625 + >>> calendar_to_jd([2004,3,(5,6)]) + array([ 2453070., 2453071.]) + >>> dates = [datetime.datetime(2004,3,5),datetime.datetime(2004,3,9)] + >>> calendar_to_jd(dates) + array([ 2453069.5, 2453073.5]) + >>> tz = dateutil.tz.tzoffset('2',3*3600) + >>> calendar_to_jd((2010,1,1),tz) + 2455197.875 + + + """ + #Adapted from xidl jdcnv.pro + from datetime import datetime,date,tzinfo + + if caltime is None: + from dateutil.tz import tzlocal + datetimes = [datetime.now(tzlocal())] + scalarout = True + elif isinstance(caltime,datetime) or isinstance(caltime,date): + datetimes = [caltime] + scalarout = True + elif all([isinstance(ct,datetime) or isinstance(ct,date) for ct in caltime]): + datetimes = caltime + scalarout = False + else: + datetimes = None + caltime = list(caltime) + if not (3 <= len(caltime) < 8): + raise ValueError('caltime input sequence is invalid size') + while len(caltime) < 7: + if len(caltime) == 3: + #make hours 12 + caltime.append(12*np.ones_like(caltime[-1])) + else: + caltime.append(np.zeros_like(caltime[-1])) + yr,month,day,hr,min,sec,msec = caltime + scalarout = all([np.shape(v) is tuple() for v in caltime]) + + #if input objects are datetime objects, generate arrays + if datetimes is not None: + yr,month,day,hr,min,sec,msec = [],[],[],[],[],[],[] + for dt in datetimes: + if not hasattr(dt,'hour'): + dt = datetime(dt.year,dt.month,dt.day,12) + + if tz is None: + off = dt.utcoffset() + if off is not None: + dt = dt - off + + yr.append(dt.year) + month.append(dt.month) + day.append(dt.day) + hr.append(dt.hour) + min.append(dt.minute) + sec.append(dt.second) + msec.append(dt.microsecond) + + + + yr = np.array(yr,dtype='int64',copy=False).ravel() + month = np.array(month,dtype='int64',copy=False).ravel() + day = np.array(day,dtype='int64',copy=False).ravel() + hr = np.array(hr,dtype=float,copy=False).ravel() + min = np.array(min,dtype=float,copy=False).ravel() + sec = np.array(sec,dtype=float,copy=False).ravel() + msec = np.array(msec,dtype=float,copy=False).ravel() + + #do tz conversion if tz is provided + if isinstance(tz,basestring) or isinstance(tz,tzinfo): + if isinstance(tz,basestring): + from dateutil import tz + tzi = tz.gettz(tz) + else: + tzi = tz + + utcoffset = [] + for t in zip(yr,month,day,hr,min,sec,msec): + #microsecond from float component of seconds + + dt = datetime(*[int(ti) for ti in t],**dict(tzinfo=tzi)) + utcdt = dt.utcoffset() + if utcdt is None: + utcoffset.append(0) + else: + utcoffset.append(utcdt.days*24 + (utcdt.seconds + utcdt.microseconds*1e-6)/3600) + else: + utcoffset = tz + +# ly = ((month-14)/12).astype(int) #In leap years, -1 for Jan, Feb, else 0 +# jdn = day - 32075l + 1461l*(yr+4800l+ly)//4 + +# jdn += 367l*(month - 2-ly*12)//12 - 3*((yr+4900l+ly)//100)//4 + +# res = jdn + (hr/24.0) + min/1440.0 + sec/86400.0 - 0.5 + + #this algorithm from meeus 2ed + m3 = month < 3 + yr[m3] -= 1 + month[m3] += 12 + + cen = yr//100 + + if gregorian is None: + gregorian = (1582,10,4) + if gregorian is True: + gregoffset = 2 - cen + cen//4 + elif gregorian is False: + gregoffset = 0 + else: + gregoffset = 2 - cen + cen//4 + gmask = (yr>gregorian[0])&(month>gregorian[1])&(day>gregorian[2]) + gregoffset[~gmask] = 0 + + + jdn = (365.25*(yr+4716)).astype(int) + \ + (30.6001*(month + 1)).astype(int) + \ + day + gregoffset - 1524.5 + res = jdn + hr/24.0 + min/1440.0 + sec/86400.0 + + if mjd: + res -= mjdoffset + + if np.any(utcoffset): + res -= np.array(utcoffset)/24.0 + + if scalarout: + return res[0] + else: + return res + + +def jd_to_epoch(jd,julian=True,asstring=False,mjd=False): + """ + Converts a Julian Date to a Julian or Besselian Epoch expressed in decimal + years. + + :param jd: Julian Date for computing the epoch, or None for current epoch. + :type jd: scalar, array-like, or None + :param bool julian: + If True, a Julian Epoch will be used (the year is exactly 365.25 days + long). Otherwise, the epoch will be Besselian (assuming a tropical year + of 365.242198781 days). + :param bool asstring: + If True, a string of the form 'J2000.0' will be returned. If it is an + integer, the number sets the number of significant figures in the output + string Otherwise, scalars are returned (an int if a whole year, float if + not). + :param bool mjd: + If True, a modified julian date is returned instead of the standard + julian date. + + :returns: + The epoch as a string (or list of strings if `jd` was array-like) if + `asstring` is True. If not, an int if a whole year, or a float (or array + if `jd` was array-like). + + :Reference: http://www.iau-sofa.rl.ac.uk/2003_0429/sofa/epj.html + + """ + if jd is None: + jd = calendar_to_jd(None) + else: + jd = np.array(jd,copy=False) + + if mjd: + jd += mjdoffset + + if julian: + epoch = 2000.0 + (jd - 2451545.0)/365.25 + else: + epoch = 1900 + (jd - 2415020.31352)/365.242198781 + + + + if asstring: + if asstring is not True: + fmt = ('J' if julian else 'B')+'{0:.'+str(int(asstring))+'}' + else: + fmt = 'J{0}' if julian else 'B{0}' + + if len(epoch.shape) == 0: + return fmt.format(epoch) + else: + return [fmt.format(e) for e in epoch] + else: + if len(epoch.shape) == 0: + if round(epoch)==epoch: + return int(epoch) + else: + return float(epoch) + else: + return epoch + +def epoch_to_jd(epoch,julian=True,mjd=False): + """ + Converts a Julian or Besselian Epoch to a Julian Day. + + :param epoch: The epoch as a decimal year. + :type epoch: string, scalar, or array-like + :param bool julian: + If True, a Julian Epoch will be used (the year is exactly 365.25 days + long). Otherwise, the epoch will be Besselian (assuming a tropical year + of 365.242198781 days). If `epoch` is a string and starts with 'B' or + 'J', this parameter will be ignored, and the 'B' or 'J' specifies the + epoch type. + :param bool mjd: + If True, a modified julian date is returned instead of the standard + julian date. + + :returns: The Julian Day as a float or array (if `epoch` is array-like) + + :Reference: http://www.iau-sofa.rl.ac.uk/2003_0429/sofa/epj.html + """ + if isinstance(epoch,basestring): + if epoch[0]=='J': + julian = True + epoch = epoch[1:] + elif epoch[0]=='B': + julian = False + epoch = epoch[1:] + epoch = float(epoch) + else: + epoch = np.array(epoch,copy=False) + if epoch.dtype.kind == 'S': + return np.array([epoch_to_jd(e,julian) for e in epoch]) + + + if julian: + res = (epoch - 2000)*365.25 + 2451545.0 + else: + res = (epoch - 1900)*365.242198781 + 2415020.31352 + + if mjd: + return res - mjdoffset + else: + return res + + +def delta_AT(jdutc,usett=False): + """ + Computes the difference between International Atomic Time (TAI) and + UTC, known as delta(AT). + + Note that this is not valid before UTC (Jan 1,1960) beganand it is not + correct for future dates, as leap seconds are not predictable. Hence, + warnings are issued if before UTC or >5 years from the date of the + algorithm. + + Implementation adapted from the matching `SOFA `_ + algorithm (dat.c). + + :param utc: + UTC time as a Julian Date (use :func:`calendar_to_jd` for calendar form + inputs.) + :type utc: float + :param usett: + If True, the return value will be the difference between UTC and + Terrestrial Time (TT) instead of TAI (TT - TAI = 32.184 s). + :type usett: bool + + :returns: TAI - UTC in seconds as a float (or TT - UTC if `usett` is True) + + """ + from warnings import warn + + + dt = jd_to_calendar(jdutc) + + #get data from arrays defined below + drift = __dat_drift + cyear,cmonth,cdelat = __dat_changes + + if dt.year > __dat_valid_year+5: + warn('delta(AT) requested for more than 5 years after current leap seconds (%i)'%dt.year) + + m = 12*dt.year + dt.month + i = np.sum(__dat_m + + +class Site(object): + """ + This class represents a location on Earth from which skies are observable. + + lat/long are coords.AngularCoordinate objects, altitude in meters. + `latitude` here is the geographic/geodetic latitude. + """ +# tznames = {'EST':-5,'CST':-6,'MST':-7,'PST':-8, +# 'EDT':-4,'CDT':-5,'MDT':-6,'PDT':-7} + def __init__(self,lat,long,alt=0,tz=None,name=None): + """ + generate a site by specifying latitude and longitude (as + coords.AngularCoordinate objects or initializers for one), optionally + providing altitude (in meters), time zone (either as a timezone name + provided by the system or as an offset from UTC), and/or a site name. + """ + self.latitude = lat + self.latitude.range = (-90,90) + self.longitude = long + self.longitude.range = (-180,180) + self.altitude = alt + if tz is None: + self.tz = self._tzFromLong(self._long) + elif isinstance(tz,basestring) and tzmod is not None: + self.tz = tzmod.gettz(tz) + if self.tz is None: + raise ValueError('unrecognized time zone string '+tz) + else: + self.tz = tzoffset(str(tz),int(tz*60*60)) + if name is None: + name = 'Default Site' + self.name = name + self._currjd = None + + @staticmethod + def _tzFromLong(long): + offset = int(np.round(long.degrees/15)) + return tzoffset(str(offset),offset*60*60) + + def _getLatitude(self): + return self._lat + def _setLatitude(self,val): + from .coords import AngularCoordinate + from operator import isSequenceType + if isinstance(val,AngularCoordinate): + self._lat = val + elif isSequenceType(val) and not isinstance(val,basestring): + self._lat = AngularCoordinate(*val) + else: + self._lat = AngularCoordinate(val) + latitude = property(_getLatitude,_setLatitude,doc='Geographic/Geodetic Latitude of the site as an :class:`AngularCoordinate` object') + + def _getGeocentriclat(self): + from .coords import geographic_to_geocentric_latitude + return geographic_to_geocentric_latitude(self._lat) + def _setGeocentriclat(self,val): + from .coords import geocentric_to_geographic_latitude + self._lat = geocentric_to_geographic_latitude(val) + geocentriclat = property(_getGeocentriclat,_setGeocentriclat,doc=None) + + + def _getLongitude(self): + return self._long + def _setLongitude(self,val): + from .coords import AngularCoordinate + from operator import isSequenceType + if isinstance(val,AngularCoordinate): + self._long = val + elif isSequenceType(val) and not isinstance(val,basestring): + self._long = AngularCoordinate(*val) + else: + self._long = AngularCoordinate(val) + longitude = property(_getLongitude,_setLongitude,doc='Longitude of the site as an :class:`AngularCoordinate` object') + + def _getAltitude(self): + return self._alt + def _setAltitude(self,val): + if val is None: + self._alt = None + else: + self._alt = float(val) + altitude = property(_getAltitude,_setAltitude,doc='Altitude of the site in meters') + + def _getCurrentobsjd(self): + if self._currjd is None: + from datetime import datetime + return calendar_to_jd(datetime.utcnow(),tz=None) + else: + return self._currjd + def _setCurrentobsjd(self,val): + if val is None: + self._currjd = None + else: + if np.isscalar(val): + self._currjd = val + else: + self._currjd = calendar_to_jd(val) + currentobsjd = property(_getCurrentobsjd,_setCurrentobsjd,doc=""" + Date and time to use for computing time-dependent values. If set to None, + the jd at the instant of calling will be used. It can also be set as + datetime objects or (yr,mon,day,hr,min,sec) tuples. + """) + + def localSiderialTime(self,*args,**kwargs): + """ + Compute the local siderial time given an input civil time. The various + input forms are used to determine the interpretation of the civil time: + + * localSiderialTime() + current local siderial time for this Site or uses the value of the + :attr:`currentobsjd` property. + * localSiderialTime(JD) + input argument is julian date UT1 + * localSdierialTime(:class:`datetime.date`) + compute the local siderial time for midnight on the given date + * localSiderialTime(:class:`datetime.datetime`) + the datetime object specifies the local time. If it has tzinfo, the + object's time zone will be used, otherwise the :class:`Site's` + * localSiderialTime(time,year,month,day) + input arguments determine local time - time is in hours + * localSiderialTime(year,month,day,hr,min,sec) + local time - hours and minutes will be interpreted as integers + + *keywords* + + * `apparent` + if True (default) the returned time will be local apparent sidereal + time (i.e. nutation terms included), otherwise it will be local mean + sidereal time + * `returntype` + a string that determines the form of the returned LST as described + below + + returns the local siderial time in a format that depends on the + `returntype` keyword. It can be: + + * None/'hours' (default) + LST in decimal hours + * 'string' + LST as a hh:mm:ss.s + * 'datetime' + a :class:`datetime.time` object + + """ + #guts of calculation adapted from xidl + import datetime + from .coords import greenwich_sidereal_time + + rettype = kwargs.pop('returntype',None) + apparent = kwargs.pop('apparent',True) + if len(kwargs)>0: + raise TypeError('got unexpected argument '+kwargs.keys()[0]) + if len(args)==0: + jd = self.currentobsjd + elif len(args)==1: + if hasattr(args[0],'year'): + if hasattr(args[0],'hour'): + if args[0].tzinfo is None: + dtobj = datetime.datetime(args[0].year,args[0].month, + args[0].day,args[0].hour,args[0].minute, + args[0].second,args[0].microsecond,self.tz) + else: + dtobj = args[0] + else: #only date provided + dtobj = datetime.datetime(args[0].year,args[0].month, + args[0].day,tzinfo=self.tz) + jd = calendar_to_jd(dtobj,tz=None) + else: + jd = args[0] + elif len(args) == 4: + time,year,month,day = args + hr = int(np.floor(time)) + min = int(np.floor(60*(time - hr))) + sec = int(np.floor(60*(60*(time-hr) - min))) + msec = int(np.floor(1e6*(60*(60*(time-hr) - min) - sec))) + jd = calendar_to_jd(datetime.datetime(year,month,day,hr,min,sec,msec,self.tz),tz=None) + elif len(args) == 6: + year,month,day,hr,min,sec = args + msec = int(1e6*(sec - np.floor(sec))) + sec = int(np.floor(sec)) + jd = calendar_to_jd(datetime.datetime(year,month,day,hr,min,sec,msec,self.tz),tz=None) + else: + raise TypeError('invalid number of input arguments') + + lst = (greenwich_sidereal_time(jd,apparent) + self._long.d/15)%24.0 + +# #from idl astro ct2lst.pro +# jd2000 = 2451545.0 +# t0 = jd - jd2000 +# t = t0//36525 #TODO:check if true-div + +# #Compute GMST in seconds. constants from Meeus 1ed, pg 84 +# c1,c2,c3,c4 = 280.46061837,360.98564736629,0.000387933,38710000.0 +# theta = c1 + (c2 * t0) + t**2*(c3 - t/ c4 ) + +# #TODO: Add in mean->apparent corrections + +# #Compute LST in hours. +# lst = np.array((theta + self._long.d)/15.0) % 24.0 + + if lst.shape == tuple(): + lst = lst.ravel()[0] + + if rettype is None or rettype == 'hours': + return lst + elif rettype == 'string': + hr = int(lst) + min = 60*(lst - hr) + sec = 60*(min - int(min)) + min = int(min) + return '%02i:%02i:%f'%(hr,min,sec) + elif rettype == 'datetime': + hr = int(lst) + min = 60*(lst - hr) + sec = 60*(min - int(min)) + min = int(min) + msec = int(1e6*(sec-int(sec))) + sec = int(sec) + return datetime.time(hr,min,sec,msec) + else: + raise ValueError('invalid returntype argument') + + def localTime(self,lsts,date=None,apparent=True,returntype=None,utc=False): + """ + Computes the local civil time given a particular local siderial time. + + `lsts` are the input local times, and may be either a scalar or array, + and must be in decimal hours. Althernatively, it can be a + :class:`datetime.datetime` object, in which case the date will be + inferred from this object and the `date` argument will be ignored. + + `date` should be a :class:`datetime.date` object or a (year,month,day) + tuple. If None, the current date will be assumed as inferred from + :attr:`Site.currentobsjd` + + if `lsts` is a :class:`datetime.datetime` object, the object will be + interpreted as the local siderial time with the corresponding date (any + timezone info will be ignored), and the `date` argument will be ignored. + + if `apparent` is True, the inputs are assumed to be in local apparent + siderial time (i.e. nutation terms included), otherwise local mean + siderial time. + + + returns the local time in a format that depends on the `returntype` + keyword. It can be: + + * None/'hours' (default) + local time in decimal hours + * 'string' + local time as a hh:mm:ss.s + * 'datetime' + a :class:`datetime.time` object with the appropriate tzinfo + + If `utc` is True, the time will be converted to UTC before being + returned. + + """ + import datetime + from operator import isSequenceType + + if isinstance(lsts,datetime.datetime): + utcoffset = lsts.replace(tzinfo=self.tz).utcoffset() + date = lsts.date() + lsts = datetime.time() + lsts = lsts.hour+lsts.minute/60+(lsts.second+lsts.microsecond*1e6)/3600 + else: + jds,dt = self._processDate(date) + utcoffset = dt.replace(tzinfo=self.tz).utcoffset() + if isSequenceType(dt): + raise ValueError('must provide only one date for localTime') + date = dt.date() + + lsts = np.array(lsts,copy=False) + scalarout = False + if lsts.shape == tuple(): + scalarout = True + lsts = lsts.ravel() + + + lst0 = self.localSiderialTime(date) + lthrs = (lsts - lst0)%24 + #NB floor rounds towards -ve infinity. This is undesirable because: + #Sometimes current lst - lst0 < 0 as the LST has looped since daybreak. + #So floor will round a negative fraction to a whole day offset. + dayoffs = [ int((lst - lst0) / 24.0) for lst in lsts ] + + lthrs /= 1.0027378507871321 #correct for siderial day != civil day + + if utc: + lthrs = (lthrs - utcoffset.days*24 - utcoffset.seconds/3600)%24 + + if returntype is None or returntype == 'hours': + res = lthrs + elif returntype == 'string': + res = [] + for lthr in lthrs: + hr = int(lthr) + min = 60*(lthr - hr) + sec = 60*(min - int(min)) + min = int(min) + res.append('%02i:%02i:%f'%(hr,min,sec)) + elif returntype == 'datetime': + res = [] + for lthr,dayoff in zip(lthrs,dayoffs): + hr = int(lthr) + min = 60*(lthr - hr) + sec = 60*(min - int(min)) + min = int(min) + msec = int(1e6*(sec-int(sec))) + sec = int(sec) + + dtobj = datetime.datetime(date.year,date.month,date.day,hr,min,sec,msec,self.tz) + res.append(dtobj + datetime.timedelta(dayoff)) + else: + raise ValueError('invalid returntype argument') + + if scalarout: + return res[0] + else: + return res + + + def equatorialToHorizontal(self,eqpos,lsts,epoch=None): + """ + Generates a list of horizontal positions (or just one) from a provided + equatorial position and local siderial time (or sequence of LSTs) in + decimal hours and a latitude in degrees. If `epoch` is not None, it + will be used to set the epoch in the equatorial system. + + Note that this generally should *not* be used to compute the observed + horizontal coordinates directly - :meth:`apparentCoordinates` includes all + the corrections and formatting. This method is purely for the + coordinate conversion. + """ + from .coords import EquatorialCoordinatesEquinox,HorizontalCoordinates + + latitude = self.latitude.d + + if epoch is not None: + #make copy so as not to change the epoch of the input + eqpos = EquatorialCoordinatesEquinox(eqpos) + eqpos.epoch = epoch + + lsts = np.array(lsts,copy=False) + singleout = lsts.shape == tuple() + + HA = lsts.ravel() - eqpos.ra.hours + sHA = np.sin(pi*HA/12) + cHA = np.cos(pi*HA/12) + + sdec = np.sin(eqpos.dec.radians) + cdec = np.cos(eqpos.dec.radians) + slat = np.sin(np.radians(latitude)) + clat = np.cos(np.radians(latitude)) + + alts = np.arcsin(slat*sdec+clat*cdec*cHA) + calts = np.cos(alts) + azs = np.arctan2(-cdec*sHA,clat*sdec-slat*cdec*cHA)%(2*pi) + + if eqpos.decerr is not None or eqpos.raerr is not None: + decerr = eqpos.decerr.radians if eqpos.decerr is not None else 0 + raerr = eqpos.raerr.radians if eqpos.raerr is not None else 0 + + dcosalt = np.cos(alts) + daltdH = -clat*cdec*sHA/dcosalt + daltddec = (slat*cdec-clat*sdec*cHA)/dcosalt + + dalts = ((daltdH*raerr)**2 + (daltddec*decerr)**2)**0.5 + + #error propogation computed with sympy following standard rules + dtanaz = 1+np.tan(azs)**2 + dazdH = (cHA*cdec/(cHA*cdec*slat-clat*sdec) \ + + cdec*cdec*sHA*sHA*slat*(cHA*cdec*slat-clat*sdec)**-2) \ + /dtanaz + dazddec = ((sHA*sdec)/(clat*sdec - cHA*cdec*slat) \ + + ((cdec*clat+cHA*sdec*slat)*cdec*sHA)*(cHA*cdec*slat-clat*sdec)**-2) \ + /dtanaz + + dazs = ((dazdH*raerr)**2 + (dazddec*decerr)**2)**0.5 + else: + dazs = dalts = None + + if singleout: + if dazs is None: + return HorizontalCoordinates(*np.degrees((alts[0],azs[0]))) + else: + return HorizontalCoordinates(*np.degrees((alts[0],azs[0],dazs[0],dalts[0]))) + else: + if dazs is None: + return [HorizontalCoordinates(alt,az) for alt,az in np.degrees((alts,azs)).T] + else: + return [HorizontalCoordinates(alt,az,daz,dalt) for alt,az,daz,dalt in np.degrees((alts,azs,dazs,dalts)).T] + + def riseSetTransit(self,eqpos,date=None,alt=-.5667,timeobj=False,utc=False): + """ + Computes the rise, set, and transit times of a provided equatorial + position in local time. + + `alt` determines the altitude to be considered as risen or set in + degrees. Default is for approximate rise/set including refraction. + + `date` determines the date at which to do the computation. See + :func:`calendar_to_jd` for acceptable formats. Note that this date is + the date of the *transit*, so the rise and set times may be on the + previous/following day. + + *returns* + (rise,set,transit) as :class:datetime.datetime objects if `timeobj` + is True or if it is False, they are in decimal hours. If the object is + circumpolar, rise and set are both None. If it is never visible, + rise,set, and transit are all None + """ + import datetime + from dateutil import tz + from math import sin,cos,radians,acos + + alt = radians(alt) + + transit = self.localTime(eqpos.ra.hours,date,utc=utc) + + #TODO:iterative algorithm for rise/set + lat = self.latitude.radians + dec = eqpos.dec.radians + #local hour angle for alt + coslha = (sin(alt) - sin(lat)*sin(dec))/(cos(lat)*cos(dec)) + if coslha<-1: + #circumpolar + rise=set=None + elif coslha>1: + #never visible + rise=set=transit=None + else: + lha = acos(coslha)*12/pi + lha /= 1.0027378507871321 #correction for siderial day != calendar day + rise = (transit - lha)%24 + set = (transit + lha)%24 #TODO:iterative algorithm + + if timeobj: + transitdate = self.localTime(eqpos.ra.hours,date,utc=utc,returntype='datetime').date() + def hrtodatetime(t): + if t is None: + return None + else: + hr = np.floor(t) + t = (t-hr)*60 + min = np.floor(t) + t = (t-min)*60 + sec = np.floor(t) + t = (t-sec)*1e6 + msec = np.floor(t) + timeobj = datetime.time(int(hr),int(min),int(sec),int(msec), + tzinfo=tz.tzutc() if utc else tz.tzlocal()) + return datetime.datetime.combine(transitdate,timeobj) + + + if rise is not None: + #now do the conversion + if rise < transit: + rise = hrtodatetime(rise) + else: + rise = hrtodatetime(rise) + datetime.timedelta(-1) + if set > transit: + set = hrtodatetime(set) + else: + set = hrtodatetime(set) + datetime.timedelta(1) + transit = hrtodatetime(transit) + + return rise,set,transit + + def nextRiseSetTransit(self, eqpos, dtime=None, alt= -.5667): + """ + Get 'next' rise,set,transit for an equatorial position at a given + datetime (see details). + + If on-sky at the specified datetime, returns the rise / set times + bracketing current transit, + i.e. returns rise, set, transit such that rise < now < set; + so rise and possibly transit will be in the past at specified time. + + If the target is off-sky currently, then rise/transit/set will all be + in the future, so now < rise < transit < set. + + If circumpolar (always visible) rise/set are None, + and next transit is returned. + If never visible, rise / set / transit are all None. + + NB The only dates checked for visibility are 'today' and 'tomorrow' + (relative to given datetime). + If the target is an edge case that's below the horizon all day today, + but is visible in 6 months, you will still get a None result. + + See also: `riseSetTransit` + + *args* + alt` determines the altitude to be considered as risen or set in + degrees. Default is for approximate rise/set including refraction. + + `dtime` determines the datetime at which to do the computation. See + :func:`calendar_to_jd` for acceptable formats. If None, the current + datetime will be assumed as inferred from :attr:`Site.currentobsjd` + + *returns* + (rise,set,transit) as :class:datetime.datetime objects (UTC). + If the object is circumpolar, rise and set are both None. + If it is never visible, rise, set, and transit are all None. + + """ + if dtime is None: + jds, dtime = self._processDate(dtime) + + today = dtime.date() + tomorrow = today + timedelta(days=1) + + #Today's transit + rise1, set1, transit1 = self.riseSetTransit(eqpos, today, + alt, + timeobj=True, utc=True) + + #Tomorrow's transit (may rise before midnight today) + rise2, set2, transit2 = self.riseSetTransit(eqpos, tomorrow, + alt, + timeobj=True, utc=True) + + #Circumpolar case, have we passed today's transit? + if rise1 is None and transit1 is not None: + if dtime < transit1: + return (rise1, set1, transit1) #i.e. (None,None,t1) + else: + return (rise2, set2, transit2) #i.e. (None,None,t2) + elif transit1 is None: #Always below horizon + return (rise1, set1, transit1) #(None,None,None) + + + #Yesterday's transit (might have transited just before midnight + # and still be on-sky). + yesterday = today - timedelta(days=1) + rise0, set0, transit0 = self.riseSetTransit(eqpos, yesterday, + alt, + timeobj=True, utc=True) + + #Yesterday's transit has not set: + if dtime < set0: + return (rise0, set0, transit0) + #Today's transit has not yet risen or is still up: + elif dtime < set1: + return (rise1, set1, transit1) + else: + return (rise2, set2, transit2) #Tomorrow's transit. + + def onSky(self, eqpos, dtime=None, alt= -.5667): + """ + Checks if a target is on-sky at specified datetime. + + NB on-sky means above minimum elevation as specifed by 'alt' parameter. + + *args* + `eqpos` Ra,Dec in equatorial co-ordinates. + + `dtime` determines the datetime at which to do the computation. See + :func:`calendar_to_jd` for acceptable formats. If None, the current + datetime will be assumed as inferred from :attr:`Site.currentobsjd` + + `alt` determines the altitude to be considered as risen or set in + degrees. Default is for approximate rise/set including refraction. + + *returns* + Boolean. + """ + if dtime is None: + jds, dtime = self._processDate(dtime) + + r, s, t = self.nextRiseSetTransit(eqpos, dtime, alt) + + #Wrong hemisphere + if t is None: + return False + + #Circumpolar, correct hemisphere. + if r is None: + return True + + #Equatorial + if r <= dtime <= s: + return True + else: + return False + + def apparentCoordinates(self,coords,datetime=None,precess=True,refraction=True): + """ + computes the positions in horizontal coordinates of an object with the + provided fixed coordinates at the requested time(s). + + `datetime` can be a sequence of datetime objects or similar representations, + or it can be a sequnce of JD's. If None, it uses the current time or + the value of the :attr:`currentobsjd` property. + + If `precess` is True, the position will be precessed to the epoch of + the observation (almost always the right thing to do) + + If `refraction` is True, an added correction to the altitude due to + atmospheric refraction at STP (formula from Meeus ch 16) is included. If + `refraction` is a (non-0) float, it will be taken as the temperature at + which to perform the refraction calculation. If it evaluates to False, + no refraction correction is performed + + + *returns* + a sequence of :class:`astropysics.coords.HorizontalCoordinates` objects + """ + from .coords import HorizontalCoordinates + from operator import isSequenceType + + if datetime is None: + return self.equatorialToHorizontal(coords,self.localSiderialTime()) + elif hasattr(datetime,'year') or (isSequenceType(datetime) and hasattr(datetime[0],'year')): + jd = calendar_to_jd(datetime,self.tz).ravel() + else: + jd = np.array(datetime,copy=False,ndmin=1) + if len(jd.shape)>1: + jd = np.array([calendar_to_jd(v,self.tz) for v in jd]) + lsts = [self.localSiderialTime(d) for d in jd] + + if precess: + res = self.equatorialToHorizontal(coords,lsts,epoch=jd_to_epoch(jd[0])) + else: + res = self.equatorialToHorizontal(coords,lsts) + + if refraction: + from math import tan,radians + if not hasattr(self,'_cached_pressure') or self._cached_pressure is None: + from math import exp + from.constants import g0,Rb + t0 = 273 #K + M = 28.9644 #g/mol + self._cached_pressure = P = exp(g0*M*self.altitude/(Rb*t0)) + else: + P = self._cached_pressure + + if refraction is True: + T = 273 + else: + T = float(refraction) + if isinstance(res,list): + res_list = res + else: + res_list = [res] + for this_res in res_list: + h = this_res.alt.radians + R = 1.02/tan(h+(10.3/(h+5.11))) #additive correction in arcmin + #for inverse problem of apparent h->true/airless h, use: + #R = 1/tan(h0+(7.31/(h0+4.4))) + + this_res.alt._decval += radians(R/60.) + + return res + + def _processDate(self,date): + """ + utitily function to convert a date in a variety of formats to a + tuple jds,datetimes + """ + if date is None: + jd = self.currentobsjd + elif np.isscalar(date): + jd = date + else: + jd = calendar_to_jd(date) + + return jd,jd_to_calendar(jd) + + def nightTable(self,coord,date=None,strtablename=None,localtime=True, + hrrange=(18,6,25)): + """ + tabulates altitude, azimuth, and airmass values for the provided fixed + position on a particular date, specified either as a datetime.date + object, a (yr,mon,day) tuple, or a julain date (will be rounded to + nearest) + + if `date` is None, the current date is used + + If `localtime` is True, the hours output (and input) will be in local + time for this Site. Otherwise, it is UTC. + + `hrrange` determines the size of the table - it should be a 3-tuple + (starthr,endhr,n) where starthr is on the day specified in date and + endhr is date + 1 day + + if `strtablename` is True, a string is returned with a printable table + of observing data. If `strtable` is a string, it will be used as the + title for the table. Otherwise, a record array is returned with the + hour(UTC), alt, az, and airmass + """ + from .astropysics import coords + import datetime + + #for objects that can get a position with no argument + if hasattr(coord,'equatorialCoordinates'): + coord = coord.equatorialCoordinates() + + if date is None: + jd = self.currentobsjd + elif np.isscalar(date): + jd = date + else: + jd = calendar_to_jd(date) + + date = jd_to_calendar(jd).date() + jd0 = np.floor(jd) + + if localtime: + dt = datetime.datetime.combine(date,datetime.time(12,tzinfo=self.tz)) + offs = dt.utcoffset() + utcoffset = offs.days*24+offs.seconds/3600 + else: + utcoffset = 0 + + starthr,endhr,n = hrrange + starthr -= utcoffset #local->UTC + endhr -= utcoffset #local->UTC + + startjd = jd0 + (starthr - 12)/24 + endjd = jd0 + (endhr + 12)/24 + jds = np.linspace(startjd,endjd,n) + + timehr = (jds-np.round(np.mean(jds))+.5)*24+utcoffset #UTC hr + + + alt,az = coords.objects_to_coordinate_arrays(self.apparentCoordinates(coord,jds,precess=False)) + + z = 90 - alt + airmass = 1/np.cos(np.radians(z)) + + ra = np.rec.fromarrays((timehr,alt,az,airmass),names = 'hour,alt,az,airmass') + + + if strtablename is not None: + rise,set,transit = self.riseSetTransit(coord,date=date,timeobj=True) + + lines = [] + + if isinstance(strtablename,basestring): + lines.append('Object:'+str(strtablename)) + lines.append('{0} {1}'.format(coord.ra.getHmsStr(),coord.dec.getDmsStr())) + lines.append(str(date)) + if transit is None: + lines.append('Not visible from this site!') + elif rise is None: + lines.append('Circumpolar') + lines.append('Transit : {0}'.format(transit)) + else: + lines.append('Rise : {0}'.format(rise)) + lines.append('Set : {0}'.format(set)) + lines.append('Transit : {0}'.format(transit)) + lines.append('') + lines.append('time\talt\taz\tairmass') + lines.append('-'*(len(lines[-1])+lines[-1].count('\t')*4)) + + currdate = date+datetime.timedelta(int(np.floor(starthr/24))) + lines.append('\t'+str(currdate)) + + oldhr24 = int(np.floor(ra[0].hour)) + for r in ra: + hr = int(np.floor(r.hour)) + if hr%24 < oldhr24: + currdate = currdate+datetime.timedelta(1) + lines.append('\t'+str(currdate)) + min = int(np.round((r.hour-hr)*60)) + if min==60: + hr+=1 + min = 0 + alt = r.alt + az = r.az + am = r.airmass if r.airmass>0 else '...' + line = '{0:2}:{1:02}\t{2:.3}\t{3:.4}\t{4:.3}'.format(hr%24,min,alt,az,am) + lines.append(line) + oldhr24 = hr%24 + return '\n'.join(lines) + else: + rise,set,transit = self.riseSetTransit(coord,date=date,timeobj=False) + ra.sitedate = date + ra.utcoffset = utcoffset + ra.rise = rise + ra.set = set + ra.transit = transit + return ra + + def nightPlot(self,coords,date=None,plottype='altam',onlynight=False, + clf=True,utc=False,colors=None, + plotkwargs=None): + + """ + Generates plots of important observability quantities for the provided + coordinate objects for a single night. + + `coords` should be a :class:`astropysics.coords.LatLongCoordinates`, a + sequence of such objects, or a dictionary mapping the name of an object + to the object itself. These names will be used to label the plot. + + `plottype` can be one of: + + * 'altam' : a plot of time vs. altitude with a secondary axis for sec(z) + * 'am' : a plot of time vs. airmass/sec(z) + * 'altaz': a plot of azimuth vs. altitude + * 'sky': polar projection of the path on the sky + + If `onlynight` is True, the plot will be shrunk to only show the times + between sunset and sunrise. Otherwise, an entire day/night will be + plotted. + + If `clf` is True, the figure is cleared before the observing plot is + made. + + If `utc` is True, the times will be in UTC instead of local time. + + `colors` should be a sequence of :mod:`matplotlib` color specifiers, or + None to use the default color cycle. + + `plotkwargs` will be provided as a keyword dictionary to + :func:`matplotlib.pyplot.plot`, unless it is None + """ + from operator import isMappingType,isSequenceType + from .coords import LatLongCoordinates + from .plotting import add_mapped_axis,mpl_context + + sun = moon = False # Sun/Moon classes removed for now. May return later + + if isMappingType(coords): + names = coords.keys() + coords = coords.values() + nonames = False + elif not isSequenceType(coords): + coords = [coords] + names = ['' for c in coords] + nonames = True + else: + names = ['' for c in coords] + nonames = True + + if colors: + plt.gca().set_color_cycle(colors) + + if isMappingType(plotkwargs): + plotkwargs = [plotkwargs for c in coords] + elif plotkwargs is None: + plotkwargs = [None for c in coords] + + + with mpl_context(clf=clf) as plt: + oldright = None + if plottype == 'altam' or plottype == 'alt' or plottype == 'am': + if plottype == 'altam': + oldright = plt.gcf().subplotpars.right + plt.subplots_adjust(right=0.86) + + for n,c,kw in zip(names,coords,plotkwargs): + ra = self.nightTable(c,date,hrrange=(12,12,100),localtime=True) + if kw is None: + kw = {} + kw.setdefault('label',n) + kw.setdefault('zorder',3) + kw.setdefault('lw',2) + if plottype == 'am': + ammask = ra.airmass>0 + x = ra.hour[ammask] + y = ra.airmass[ammask] + else: + x = ra.hour + y = ra.alt + plt.plot(x,y,**kw) + + plt.title(str(ra.sitedate)) + + if 'alt' in plottype: + if plt.ylim()[0] < 0: + lowery = 0 + else: + lowery = plt.ylim()[0] + uppery = plt.ylim()[1] if plt.ylim()[1]<90 else 90 + plt.ylim(lowery,uppery) + plt.ylabel(r'${\rm altitude} [{\rm degrees}]$') + elif 'am' == plottype: + if plt.ylim()[0] < 1: + lowery = 1 + else: + lowery = plt.ylim()[0] + uppery = plt.ylim()[1] + plt.ylim(lowery,uppery) + if lowery==1: + yticks = list(plt.yticks()[0]) + yticks.insert(0,1) + plt.yticks(yticks) + plt.ylim(uppery,lowery) + plt.ylabel(r'${\rm airmass} / \sec(z)$') + + if not nonames: + plt.legend(loc=0) + + if sun: + seqp = Sun(ra.sitedate).equatorialCoordinates() + rise,set,t = self.riseSetTransit(seqp,date,0) + rise12,set12,t12 = self.riseSetTransit(seqp,date,-12) + rise18,set18,t18 = self.riseSetTransit(seqp,date,-18) + #need to correct them back to the previous day + set -= 24 + set12 -= 24 + set18 -= 24 + + xl,xu = plt.xlim() + yl,yu = plt.ylim() + + grey1 = (0.5,0.5,0.5) + grey2 = (0.65,0.65,0.65) + grey3 = (0.8,0.8,0.8) + plt.fill_between((set,set12),lowery,uppery,lw=0,color=grey3,zorder=1) + plt.fill_between((set12,set18),lowery,uppery,lw=0,color=grey2,zorder=1) + plt.fill_between((set18,rise18),lowery,uppery,lw=0,color=grey1,zorder=1) + plt.fill_between((rise18,rise12),lowery,uppery,lw=0,color=grey2,zorder=1) + plt.fill_between((rise12,rise),lowery,uppery,lw=0,color=grey3,zorder=1) + + + + plt.xlim(xl,xu) + plt.ylim(yl,yu) + + if moon: + xls = plt.xlim() + yls = plt.ylim() + + m = Moon(ra.sitedate) + ram = self.nightTable(m,date,hrrange=(12,12,100),localtime=True) + if not isMappingType(moon): + moon = {} + moon.setdefault('label','Moon') + moon.setdefault('zorder',2) + moon.setdefault('ls','--') + moon.setdefault('lw',1) + moon.setdefault('c','k') + if plottype == 'am': + ammask = ram.airmass>0 + plt.plot(ram.hour[ammask],ram.airmass[ammask],**moon) + else: + plt.plot(ram.hour,ram.alt,**moon) + + plt.xlim(*xls) + plt.ylim(*yls) + + if plottype == 'altam': + firstax = plt.gca() + yticks = plt.yticks()[0] + newticks = list(yticks[1:]) + newticks.insert(0,(yticks[0]+newticks[0])/2) + newticks.insert(0,(yticks[0]+newticks[0])/2) + plt.twinx() + plt.ylim(lowery,uppery) + plt.yticks(newticks,['%2.2f'%(1/np.cos(np.radians(90-yt))) for yt in newticks]) + plt.ylabel(r'$\sec(z)$') + plt.axes(firstax) + + xtcks = np.round(np.arange(13)*(x[-1]-x[0])/12+x[0]).astype(int) + if utc: + plt.xticks(xtcks,['%i'%np.round(xt%24) for xt in xtcks]) + plt.xlabel(r'$\rm UTC (hours)$') + else: + plt.xticks(xtcks,[xt%24 for xt in xtcks]) + plt.xlabel(r'$\rm Local Time (hours)$') + + if onlynight: + plt.xlim(xtcks[0]/2,xtcks[-1]/2) + else: + plt.xlim(xtcks[0],xtcks[-1]) + + + elif plottype == 'altaz': + for n,c,kw in zip(names,coords,plotkwargs): + ra = self.nightTable(c,date,hrrange=(0,0,100)) + + if kw is None: + plt.plot(ra.az,ra.alt,label=n) + else: + kw['label'] = n + plt.plot(ra.az,ra.alt,**kw) + + plt.xlim(0,360) + plt.xticks(np.arange(9)*360/8) + plt.ylim(0,90) + + plt.title(str(ra.sitedate)) + plt.xlabel(r'${\rm azimuth} [{\rm degrees}]$') + plt.ylabel(r'${\rm altitude} [{\rm degrees}]$') + + if not nonames: + plt.legend(loc=0) + + elif plottype == 'sky': + for n,c,kw in zip(names,coords,plotkwargs): + + ra = self.nightTable(c,date,hrrange=(0,0,100)) + if kw is None: + plt.polar(np.radians(ra.az),90-ra.alt,label=n) + else: + kw['label'] = n + plt.polar(np.radians(ra.az),90-ra.alt,**kw) + + xticks = [0,45,90,135,180,225,270,315] + xtlabs = ['N',r'$45^\circ$','E',r'$135^\circ$','S',r'$225^\circ$','W',r'$315^\circ$'] + plt.xticks(np.radians(xticks),xtlabs) + plt.ylim(0,90) + yticks = [15,30,45,60,75] + plt.yticks(yticks,[r'${0}^\circ$'.format(int(90-yt)) for yt in yticks]) + + plt.title(str(ra.sitedate)) + + if not nonames: + plt.legend(loc=0) + else: + raise ValueError('unrecognized plottype {0}'.format(plottype)) + if oldright is not None: + plt.gcf().subplotpars.right = oldright + + def yearPlot(self,coords,startdate=None,n=13,months=12, + utc=False,clf=True,colors=None): + """ + Plots the location transit and rise/set times of object(s) over ~year + timescales. + + `coords` should be a :class:`astropysics.coords.LatLongCoordinates`, a + sequence of such objects, or a dictionary mapping the name of an object + to the object itself. These names will be used to label the plot. + + `startdate` is a :class:`datetime.date` object or a date tuple + (year,month,day) that is used as the start of the plot. + + `n` specifies the number of points to include of the objects + + `months` is the number of months to draw the plot for. + + If `utc` is True, the times will be in UTC instead of local time. + + If `clf` is True, the figure is cleared before the observing plot is + made. + + `colors` should be a sequence of :mod:`matplotlib` color specifiers, or + None to use the default color cycle. + """ + import datetime + from matplotlib.dates import MonthLocator,WeekdayLocator,DateFormatter, \ + MONDAY,DayLocator,YearLocator + from operator import isMappingType,isSequenceType + from .plotting import mpl_context + + sun = moon = False # Sun/Moon classes removed for now. May return later + + if isMappingType(coords): + names = coords.keys() + coords = coords.values() + elif not isSequenceType(coords): + coords = [coords] + names = ['' for c in coords] + else: + names = ['' for c in coords] + + jdstart,startdt = self._processDate(startdate) + jds = jdstart+np.linspace(0,365.25*months/12,n) + jd1 = jds - calendar_to_jd((1,1,1)) + + if utc: + utco = startdt.replace(tzinfo=self.tz).utcoffset() + center = -utco.days*24 - utco.seconds/3600 + else: + center = 0 + cp12 = center + 12 + + with mpl_context(clf=clf) as plt: + if colors: + plt.gca().set_color_cycle(colors) + for c,nm in zip(coords,names): + rst = [self.riseSetTransit(c,jd_to_calendar(jd),0,utc=utc) for jd in jds] + rise,set,transit = np.array(rst).T + transit[transit>cp12] -= 24 #do this to get the plot to cross over night time + + c = plt.gca()._get_lines.color_cycle.next() + + tline = plt.plot_date(jd1,transit,label=nm,color=c)[0] + rise[rise>transit]-=24 + set[setcp12] -= 24 + + grey1 = (0.5,0.5,0.5) + grey2 = (0.65,0.65,0.65) + grey3 = (0.8,0.8,0.8) + plt.fill_between(jd1sun,rise18,set18,color=grey1) + plt.fill_between(jd1sun,rise18,rise12,color=grey2) + plt.fill_between(jd1sun,rise12,rise,color=grey3) + plt.fill_between(jd1sun,set12,set18,color=grey2) + plt.fill_between(jd1sun,set,set12,color=grey3) + + if moon: + jdmoon = jdstart+np.linspace(0,365.25*months/12,30*months) + jd1moon = jdmoon - calendar_to_jd((1,1,1)) + moon = Moon() + + rst = [] + for jd in jdmoon: + moon.jd = jd + rst.append(self.riseSetTransit(moon.equatorialCoordinates(),jd_to_calendar(jd),0,utc=utc)) + rise,set,t = np.array(rst).T + + t[t>cp12] -= 24 + transitionmask = np.roll(t,1)>t + transitionmask[0] = transitionmask[-1] = False + + lasti = 0 + for i in np.where(transitionmask)[0]: + lastplot = plt.plot(jd1moon[lasti:i],t[lasti:i],'--k')[0] + lasti = i + if lasti == len(t): + lastplot._label = 'Moon' + else: + plt.plot(jd1moon[lasti:],t[lasti:],'--k',label='Moon') + if utc: + plt.ylabel('UTC (hours)') + else: + plt.ylabel('Local Time (hours)') + plt.xlabel('Month (small ticks on Mondays)') + + + + if months <=3: + plt.gca().xaxis.set_major_locator(MonthLocator()) + plt.gca().xaxis.set_major_formatter(DateFormatter("%b")) + plt.gca().xaxis.set_minor_locator(DayLocator((6,10,15,20,25))) + plt.gca().xaxis.set_minor_formatter(DateFormatter('$%d$')) + elif 36 >= months > 12: + plt.gca().xaxis.set_major_locator(MonthLocator(interval=3)) + plt.gca().xaxis.set_major_formatter(DateFormatter("%b '%y")) + plt.gcf().autofmt_xdate() + elif months > 36: + plt.gca().xaxis.set_major_locator(YearLocator()) + plt.gca().xaxis.set_major_formatter(DateFormatter("%Y")) + plt.gcf().autofmt_xdate() + else: + plt.gca().xaxis.set_major_locator(MonthLocator()) + plt.gca().xaxis.set_major_formatter(DateFormatter("%b '%y")) + plt.gca().xaxis.set_minor_locator(WeekdayLocator(MONDAY)) + plt.gcf().autofmt_xdate() + + plt.xlim(jd1[0],jd1[-1]) + plt.ylim(center-12,center+12) + ytcks = np.round(np.arange(13)*2-12+center).astype(int) + plt.yticks(ytcks,ytcks%24) + + if any([nm!='' for nm in names]): + plt.legend(loc=0) + + +def __loadobsdb(sitereg): + from .utils.io import get_package_data + obsdb = get_package_data('obsdb.dat') + from .coords import AngularCoordinate + + obs = None + for l in obsdb.split('\n'): + ls = l.strip() + if len(ls)==0 or ls[0]=='#': + continue + k,v = [ss.strip() for ss in l.split('=')] + if k == 'observatory': + if obs is not None: + o = Site(lat,long,alt,tz,name) + sitereg[obs] = o + obs = v.replace('"','') + name = long = lat = alt = tz = None + elif k == 'name': + name = v.replace('"','') + elif k == 'longitude': + vs = v.split(':') + dec = float(vs[0]) + if len(vs)>1: + dec += float(vs[1])/60 + #longitudes here are deg W instead of (-180, 180) + dec*=-1 + if dec <=-180: + dec += 360 + long = AngularCoordinate(dec) + elif k == 'latitude': + vs = v.split(':') + dec = float(vs[0]) + if len(vs)>1: + dec += float(vs[1])/60 + lat = AngularCoordinate(dec) + elif k == 'altitude': + alt = float(v) + elif k == 'timezone': + #time zones are also flipped + tz = -1*float(v) + if obs is not None: + o = Site(lat,long,alt,tz,name) + sitereg[obs] = o + + +sites = DataObjectRegistry('sites',Site) +try: + sites['uciobs'] = Site(33.63614044191056,-117.83079922199249,80,'US/Pacific','UC Irvine Observatory') +except ValueError: #in case US/Pacific is not present for some reason + sites['uciobs'] = Site(33.63614044191056,-117.83079922199249,80,-8,'UC Irvine Observatory') +sites['greenwich'] = Site('51d28m38s',0,7,0,'Royal Observatory,Greenwich') +__loadobsdb(sites) + + +#<-----------------Attenuation/Reddening and dust-related----------------------> + +class Extinction(object): + """ + This is the base class for extinction-law objects. Extinction laws can be + passed in as a function to the initializer, or subclasses should override + the function f with the preferred law as f(self,lambda), and their + initializers should set the zero point + + Note that functions are interpreted as magnitude extinction laws ... if + optical depth is desired, f should return (-2.5/log(10))*tau(lambda) + + A0 is the normalization factor that gets multiplied into the reddening law. + """ + def __init__(self,f=None,A0=1): + if f is not None: + if not callable(f): + self.f = f + else: + raise ValueError('function must be a callable') + self.A0 = A0 + + self._plbuffer = None + + def f(self,lamb): + raise NotImplementedError('must specify an extinction function ') + + def __call__(self,*args,**kwargs): + return self.A0*self.f(*args,**kwargs) + + def correctPhotometry(self,mags,band): + """ + Uses the extinction law to correct a magnitude (or array of magnitudes) + + bands is either a string specifying the band, or a wavelength to use + """ + from .phot import bandwl + #TODO:better band support + if isinstance(band,basestring): + wl = bandwl[band] + else: + wl = band + + return mags-self(wl) + + def Alambda(self,band): + """ + determines the extinction for this extinction law in a given band or bands + + band can be a wavelength or a string specifying a band + """ + from operator import isSequenceType + if isSequenceType(band) and not isinstance(band,basestring): + return -1*np.array([self.correctPhotometry(0,b) for b in band]) + else: + return -1*self.correctPhotometry(0,band) + + def correctColor(self,colors,bands): + """ + Uses the supplied extinction law to correct a color (or array of colors) + where the color is in the specified bands + + bands is a length-2 sequence with either a band name or a wavelength for + the band, or of the form 'bandname1-bandname2' or 'E(band1-band2)' + """ + + if isinstance(bands,basestring): + b1,b2=bands.replace('E(','').replace(')','').split(',') + else: + b1,b2=bands + + from .phot import bandwl + #TODO:better band support + wl1 = bandwl[b1] if isinstance(b1,basestring) else b1 + wl2 = bandwl[b2] if isinstance(b1,basestring) else b2 + + return colors-self(wl1)+self(wl2) + + def correctSpectrum(self,spec,newspec=True): + """ + Uses the supplied extinction law to correct a spectrum for extinction. + + if newspec is True, a copy of the supplied spectrum will have the + extinction correction applied + + returns the corrected spectrum + """ + + if newspec: + spec = spec.copy() + + oldunit = spec.unit + spec.unit = 'wavelength-angstrom' + corr = 10**(self(spec.x)/2.5) + spec.flux *= corr + spec.err *= corr + + spec.unit = oldunit + return spec + + __builtinlines ={ + 'Ha':6562.82, + 'Hb':4861.33, + 'Hg':4340.46, + 'Hd':4101.74, + 'He':3970.07 + } + def correctFlux(self,flux,lamb): + if isinstance(lamb,basestring) and lamb in Extinction.__builtinlines: + lamb = Extinction.__builtinlines[lamb] + return flux*10**(self(lamb)/2.5) + + + __balmerratios={ + 'Hab':(2.86,6562.82,4861.33), + 'Hag':(6.16,6562.82,4340.46), + 'Had':(11.21,6562.82,4101.74), + 'Hae':(18.16,6562.82,3970.07), + 'Hbg':(2.15,4861.33,4340.46), + 'Hbd':(3.91,4861.33,4101.74), + 'Hbe':(6.33,4861.33,3970.07), + 'Hgd':(1.82,4340.46,4101.74), + 'Hge':(2.95,4340.46,3970.07), + 'Hde':(1.62,4101.74,3970.07) + } + def computeA0FromFluxRatio(self,measured,expected,lambda1=None,lambda2=None,filterfunc=None): + """ + This derives the normalization of the Extinction function from provided + ratios for theoretically expected fluxes. If multiple measurements are + provided, the mean will be used + + measured is the measured line ratio (possibly an array), while expected + is either the expected line ratios, or a string specifying the + appropriate balmer flux ratio as "Hab","Hde",etc. (for Halpha/Hbeta or + Hdelta/Hepsilon) to assume case B recombination fluxes (see e.g. + Osterbrock 2007). + + lambda1 and lambda2 are the wavelengths for the ratios F1/F2, or None if + a string is provided + + filterfunc is a function to be applied as the last step - it can + either be used to contract an array to (e.g. np.mean), or filter out + invalid values (e.g. lambda x:x[np.isfinite(x)]). Default + does nothing + + returns A0,standard deviation of measurements + """ + if isinstance(expected,basestring): + R = measured + balmertuple = Extinction.__balmerratios[expected] + R0=balmertuple[0] + lambda1,lambda2=balmertuple[1],balmertuple[2] + elif isinstance(expected[0],basestring): + if not np.all([isinstance(e,basestring) for e in expected]): + raise ValueError('expected must be only transitions or only numerical') + R0,lambda1,lambda2=[],[],[] + for e in expected: + balmertuple = Extinction.__balmerratios[e] + R0.append(balmertuple[0]) + lambda1.append(balmertuple[1]) + lambda2.append(balmertuple[2]) + R0,lambda1,lambda2=np.array(R0),np.array(lambda1),np.array(lambda2) + else: + if lambda1 and lambda2: + R,R0 = np.array(measured,copy=False),np.array(expected,copy=False) + lambda1,lambda2 = np.array(lambda1,copy=False),np.array(lambda2,copy=False) + else: + raise ValueError('need to provide wavelengths if not specifying transitions') + + A0 = -2.5*np.log10(R/R0)/(self.f(lambda1)-self.f(lambda2)) + + if filterfunc is None: + filterfunc = lambda x:x + self.A0 = A0 = filterfunc(A0) + return A0,np.std(A0) + + #PipelineElement methods +# def _plFeed(self,data,src): +# from .pipeline import PipelineError +# from .spec import Spectrum +# if self._plbuffer is None: +# self._plbuffer = {'in':[],'out':[]} + +# if isinstance(data,Spectrum): +# self._plbuffer['in'].append(('spec',data)) +# else: +# raise PipelineError('unrecognized Extinction correction input data') + +# def _plProcess(self): +# from .pipeline import PipelineError +# if self._plbuffer is None: +# self._plbuffer = {'in':[],'out':[]} + +# type,data = self._plbuffer['in'].pop(0) +# try: +# if type=='spec': +# newspec = self.correctSpectrum(data) +# self._plbuffer['out'].append(newspec) +# else: +# assert False,'Impossible point - code error in Extinction pipeline' +# except: +# #TODO:check this +# self.insert(0,spec(type,data)) + +# def _plExtract(self): +# if self._plbuffer is None: +# self._plbuffer = {'in':[],'out':[]} + +# if len(self._plbuffer['out']<1): +# return None +# else: +# return self._plbuffer['out'].pop(0) + +# def plClear(self): +# self._plbuffer = None + +class CalzettiExtinction(Extinction): + """ + This is the average extinction law derived in Calzetti et al. 1994: + x=1/lambda in mu^-1 + Q(x)=-2.156+1.509*x-0.198*x**2+0.011*x**3 + """ + _poly=np.poly1d((0.011,-0.198,1.509,-2.156)) + def __init__(self,A0=1): + super(CalzettiExtinction,self).__init__(A0=A0) + + def f(self,lamb): + return (-2.5/np.log(10))*self._poly(1e4/lamb) + +class _EBmVExtinction(Extinction): + """ + Base class for Extinction classes that get normalization from E(B-V) + """ + + def __init__(self,EBmV=1,Rv=3.1): + super(_EBmVExtinction,self).__init__(f=None,A0=1) + self.Rv=Rv + self.EBmV=EBmV + + def _getEBmV(self): + from .phot import bandwl + return self.A0*self.f(bandwl['V'])/self.Rv + def _setEBmV(self,val): + from .phot import bandwl + Av=self.Rv*val + self.A0=Av/self.f(bandwl['V']) + EBmV = property(_getEBmV,_setEBmV) + + def f(self,lamb): + raise NotImplementedError + +class FMExtinction(_EBmVExtinction): + """ + Base class for Extinction classes that use the form from + Fitzpatrick & Massa 90 + """ + + def __init__(self,C1,C2,C3,C4,x0,gamma,EBmV=1,Rv=3.1): + self.x0 = x0 + self.gamma = gamma + self.C1 = C1 + self.C2 = C2 + self.C3 = C3 + self.C4 = C4 + + super(FMExtinction,self).__init__(EBmV=EBmV,Rv=Rv) + + def f(self,lamb): + x=1e4/np.array(lamb,copy=False) + C1,C2,C3,C4 = self.C1,self.C2,self.C3,self.C4 + gamma,x0 = self.gamma,self.x0 + + xsq=x*x + D=xsq*((xsq-x0*x0)**2+xsq*gamma*gamma)**-1 + FMf=C1+C2*x+C3*D + + if np.isscalar(FMf): + if x>=5.9: + FMf+=C4*(0.5392*(x-5.9)**2+0.05644*(x-5.9)**3) + else: + C4m=x>=5.9 + FMf[C4m]+=C4*(0.5392*(x[C4m]-5.9)**2+0.05644*(x[C4m]-5.9)**3) + + return FMf+self.Rv #EBmV is the normalization and is multiplied in at the end + +class CardelliExtinction(_EBmVExtinction): + """ + Milky Way Extinction law from Cardelli et al. 1989 + """ + def f(self,lamb): + scalar=np.isscalar(lamb) + x=1e4/np.array(lamb,ndmin=1) #CCM x is 1/microns + a,b=np.ndarray(x.shape,x.dtype),np.ndarray(x.shape,x.dtype) + + if any((x<0.3)|(10 0) > 0: + print 'using ngp file when lat < 0 present... put %s wherever "ngp" or "sgp" should go in filename' + elif polename=='sgp': + n=[-1] + if sum(b < 0) > 0: + print 'using sgp file when lat > 0 present... put %s wherever "ngp" or "sgp" should go in filename' + else: + raise ValueError("couldn't determine South/North from filename - should have 'sgp' or 'ngp in it somewhere") + masks = [ones_like(b).astype(bool)] + else: #need to do things seperately for north and south files + nmask = b >= 0 + smask = ~nmask + + masks = [nmask,smask] + ns = [1,-1] + + mapds=[] + f=open(dustmapfn%'ngp') + try: + mapds.append(f[0].data) + finally: + f.close() + assert mapds[-1].shape[0] == mapds[-1].shape[1],'map dimensions not equal - incorrect map file?' + f=open(dustmapfn%'sgp') + try: + mapds.append(f[0].data) + finally: + f.close() + assert mapds[-1].shape[0] == mapds[-1].shape[1],'map dimensions not equal - incorrect map file?' + + retvals=[] + for n,mapd,m in zip(ns,mapds,masks): + #project from galactic longitude/latitude to lambert pixels (see SFD98) + npix=mapd.shape[0] + + x=npix/2*cos(l[m])*(1-n*sin(b[m]))**0.5+npix/2-0.5 + y=-npix/2*n*sin(l[m])*(1-n*sin(b[m]))**0.5+npix/2-0.5 + #now remap indecies - numpy arrays have y and x convention switched from SFD98 appendix + x,y=y,x + + if interpolate: + from scipy.ndimage import map_coordinates + if type(interpolate) is int: + retvals.append(map_coordinates(mapd,[x,y],order=interpolate)) + else: + retvals.append(map_coordinates(mapd,[x,y])) + else: + x=round(x).astype(int) + y=round(y).astype(int) + retvals.append(mapd[x,y]) + + + + + if isscalar(long) or isscalar(lat): + for r in retvals: + if len(r)>0: + return r[0] + assert False,'None of the return value arrays were populated - incorrect inputs?' + else: + #now recombine the possibly two arrays from above into one that looks like the original + retval=ndarray(l.shape) + for m,val in zip(masks,retvals): + retval[m] = val + return retval + + +def get_dust_radec(ra,dec,dustmap,interpolate=True): + from .coords import equatorial_to_galactic + l,b = equatorial_to_galactic(ra,dec) + return get_SFD_dust(l,b,dustmap,interpolate) + + + +#DEPRECATED! +def extinction_correction(lineflux,linewl,EBmV,Rv=3.1,exttype='MW'): + """ + + Extinction correct a la Cardelli et al 89 from the supplied line data and + a given E(B-V) along the sightline + + inputs may be numpy arrays + + linewl is in angstroms, lineflux in erg s^-1 cm^-2 + + if the input lineflux is None (or NaN but NOT False or 0) , Alambda is returned instead + """ + from numpy import array,logical_and,logical_or,polyval,log10,isscalar,where + from .phot import bandwl + + from warnings import warn + warn('extinction_correction function is deprecated - use Extinction class instead',DeprecationWarning) + + if exttype=='LMC': + eo = LMCExtinction(EBmV=EBmV,Rv=Rv) + return eo.correctFlux(lineflux,linewl) + elif exttype == 'SMC': + eo = SMCExtinction(EBmV=EBmV,Rv=Rv) + return eo.correctFlux(lineflux,linewl) + + if isinstance(linewl,basestring): + linewl=bandwl[linewl] + + if lineflux is None or lineflux is False: + lineflux=np.nan + + + lineflux=np.array(lineflux,ndmin=1,dtype=float) + linewl=np.array(linewl,ndmin=1) + EBmV=np.array(EBmV,ndmin=1) + n=np.max((lineflux.size,linewl.size,EBmV.size)) + + if n!=1: + if lineflux.size == 1: + lineflux = np.ones(n)*lineflux[0] + elif lineflux.size != n: + raise ValueError("lineflux is incorrect length") + if linewl.size == 1: + linewl = np.ones(n)*linewl[0] + elif linewl.size != n: + raise ValueError("linewl is incorrect length") + if EBmV.size == 1: + EBmV = np.ones(n)*EBmV[0] + elif EBmV.size != n: + raise ValueError("EBmV is incorrect length") + + x=1e4/linewl #CCM x is 1/microns + + a=array(x) + b=array(x) + + if any(logical_or(x<0.3,10