#! /usr/bin/env python # ========================================================================== # Compute a visibility cube # # Copyright (C) 2016-2019 Juergen Knoedlseder # # This program 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. # # This program 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. # # You should have received a copy of the GNU General Public License # along with this program. If not, see . # # ========================================================================== import sys import math import gammalib import ctools # =============== # # csviscube class # # =============== # class csviscube(ctools.cscript): """ Compute a visibility cube """ # Constructor def __init__(self, *argv): """ Constructor Parameters ---------- argv : list of str List of IRAF command line parameter strings of the form ``parameter=3``. Raises ------ TypeError An invalid number of command line arguments was provided. """ # Initialise application by calling the base class constructor self._init_cscript(self.__class__.__name__, ctools.__version__, argv) # Initialise members self._cube = gammalib.GSkyMap() self._results = [] # Return return # Private methods def _get_parameters(self): """ Get all parameters """ # Query parameters self['mjdref'].real() self['tmin'].time() self['tmax'].time() self['geolon'].real() self['geolat'].real() self['outfile'].filename() # Write input parameters into logger self._log_parameters(gammalib.TERSE) # Return return def _zenith_angle_map(self, nx, ny, dx, dy): """ Compute zenith angle map The zenith angle of a position (ra,dec) depends on the declination and the hour angle h and is given by zenith(h,dec) = arccos( sin(lat) * sin(dec) + cos(lat) * cos(dec) * cos(h) ) The hour angle h (or local hour angle, LHA) is defined as the difference between local siderial time (LST) and the Right Ascension h = LST - ra The map is computed for h=-ra which is equivalent to GST=lon (or LST=0). In other words, the map corresponds to the time when ra=0 passes through the local meridian. Parameters ---------- nx : int Number of Right Ascension pixels. ny : int Number of Declination pixels. dx : float Right Ascension pixel size in degrees. dy : float Declination pixel size in degrees. Returns ------- zmap : `~gammalib.GSkyMap` Allsky map comprising the zenith angle for an hour angle of 0. """ # Initialise zenith angle map zmap = gammalib.GSkyMap('CAR','CEL',0.0,0.0,-dx,dy,nx,ny) # Set hour angle and declination vectors hours = [(float(i)+0.5)*dx for i in range(nx)] decs = [(float(i)+0.5)*dy-90.0 for i in range(ny)] # Get array geographic latitude geolat = self['geolat'].real() # Precompute latitude terms cos_lat = math.cos(geolat*gammalib.deg2rad) sin_lat = math.sin(geolat*gammalib.deg2rad) # Loop over all declinations and hour angles and compute the zenith # angle. Store the zenith angle in the map index = 0 for dec in decs: cos_dec = math.cos(dec*gammalib.deg2rad) sin_dec = math.sin(dec*gammalib.deg2rad) for h in hours: cos_h = math.cos(h*gammalib.deg2rad) zenith = math.acos(sin_lat*sin_dec + cos_lat*cos_dec*cos_h)*gammalib.rad2deg zmap[index] = zenith index += 1 # Log zenith angle map self._log_header2(gammalib.EXPLICIT, 'Zenith angle map') self._log_string(gammalib.EXPLICIT, str(zmap)) # Return zenith angle map return zmap def _sun_ra_exclusion(self, time): """ Compute the half length of the Right Ascension interval to exclude due to the Sun constraint The Sun zenith angle constraint implies an interval in Right Ascension that is to be excluded (that's the interval during which it is day). This method computes half the length of this interval in degrees. Parameters ---------- time : `~gammalib.GTime()` Time for which the Sun exclusion is to be computed Returns ------- dra : float Right Ascension interval half length (degrees) """ # Get array geographic latitude and minimum Sun zenith angle in radians geolat = self['geolat'].real() * gammalib.deg2rad sunzenith = self['sunzenith'].real() * gammalib.deg2rad # Compute Right Ascension and Declination of Sun in degrees and # convert the Declination into radians sun = gammalib.GSkyDir() sun.sun(time) sundec = sun.dec() # Compute some sines and cosines cos_sunzenith = math.cos(sunzenith) sin_geolat = math.sin(geolat) cos_geolat = math.cos(geolat) sin_sundec = math.sin(sundec) cos_sundec = math.cos(sundec) # Compute Right Ascension difference when Sun is below the mimimum # zenith angle in degrees dra = 0.0 nom = cos_sunzenith - sin_geolat * sin_sundec denom = cos_geolat * cos_sundec if denom != 0.0: arg = nom/denom if arg >= -1.0 and arg <= 1.0: dra = math.acos(arg) * gammalib.rad2deg # Return return dra def _moon_ra_exclusion(self, time): """ Compute the half length of the Right Ascension interval to exclude due to the Moon constraint The Moon zenith angle constraint implies an interval in Right Ascension that is to be excluded. This method computes half the length of this interval. Parameters ---------- time : `~gammalib.GTime()` Time for which the Moon exclusion is to be computed Returns ------- dra : float Right Ascension interval half length (degrees) """ # Get array geographic latitude and minimum Moon zenith angle in radians geolat = self['geolat'].real() * gammalib.deg2rad moonzenith = self['moonzenith'].real() * gammalib.deg2rad # Compute Right Ascension and Declination of Moon in degrees and # convert the Declination into radians moon = gammalib.GSkyDir() moon.moon(time) moondec = moon.dec() # Compute some sines and cosines cos_moonzenith = math.cos(moonzenith) sin_geolat = math.sin(geolat) cos_geolat = math.cos(geolat) sin_moondec = math.sin(moondec) cos_moondec = math.cos(moondec) # Compute Right Ascension difference when Moon is below the mimimum # zenith angle in degrees dra = 0.0 nom = cos_moonzenith - sin_geolat * sin_moondec denom = cos_geolat * cos_moondec if denom != 0.0: arg = nom/denom if arg >= -1.0 and arg <= 1.0: dra = math.acos(arg) * gammalib.rad2deg # Return return dra def _adjust_ra_interval(self, ra_start, ra_stop): """ Adjust Right Ascension interval so that it overlaps with the [0,360[ interval Parameters ---------- ra_start : float Start of interval to exclude (degrees) ra_stop : float Stop of interval to exclude (degrees) Returns ------- ra_start : float Adjusted start of interval to exclude (degrees) ra_stop : float Adjusted stop of interval to exclude (degrees) """ # Adjust interval if both boundaries are negative while ra_start < 0.0 and ra_stop < 0.0: ra_start += 360.0 ra_stop += 360.0 # Adjust interval if both boundaries are equal or larger than 360 # degrees while ra_start >= 360.0 and ra_stop >= 360.0: ra_start -= 360.0 ra_stop -= 360.0 # Return interval return ra_start, ra_stop def _exclude_ra_interval(self, hours, ra_start, ra_stop): """ Exclude Right Ascension interval from array of hour angles Parameters ---------- hours : list of floats Array of hours ra_start : float Start of interval to exclude (degrees) ra_stop : float Stop of interval to exclude (degrees) Returns ------- hours : list of floats Array of hours with excluded Right Ascension interval """ # Conversion from RA (degrees) to index nhours = len(hours) ra2inx = float(nhours)/360.0 # Adjust interval ra_start, ra_stop = self._adjust_ra_interval(ra_start, ra_stop) # Case 1: The RA interval is fully comprised within the [0,360] # interval, hence we simply exclude this interval. if ra_start >= 0.0 and ra_stop <= 360.0: inx_start = int(ra_start * ra2inx + 0.5) inx_stop = int(ra_stop * ra2inx + 0.5) for i in range(inx_start,inx_stop): hours[i] = 0.0 # Case 2: The RA interval is starting at negative RA. In that case # there are two intervals to exclude: # [ra_start+360,360] and [0,ra_stop] elif ra_start < 0.0: inx_start = int((ra_start+360.0) * ra2inx + 0.5) inx_stop = int(ra_stop * ra2inx + 0.5) for i in range(inx_start,nhours): hours[i] = 0.0 for i in range(0,inx_stop): hours[i] = 0.0 # Case 3: The RA interval is stopping at RA>360. In that case there # are two intervals to exclude: # [0,ra_stop-360] and [0,ra_start] elif ra_stop > 360.0: inx_start = int(ra_start * ra2inx + 0.5) inx_stop = int((ra_stop-360.0) * ra2inx + 0.5) for i in range(0,inx_stop): hours[i] = 0.0 for i in range(inx_start,nhours): hours[i] = 0.0 # Return hours return hours def _hour_angle_weight(self): """ Compute hour angle weights Computes an array that specifies the time during which a given hour angle was observed during the observing time interval [tmin,tmax]. The hour angle runs from 0 to 360 degrees, array values are in units of hours. The method loops over the days in the time interval [tmin,tmax], with time = tmin + i*86400 seconds. At each time step the position in Right Ascension and Declination of the Sun is computed (the Sun's Declination is in fact not used here). An interval of [ra_sun-dra_sun, ra_sun+dra_sun] is assumed to be the day. The method _sun_ra_exclusion() is used to compute dra_sun, which is half of the length of the day in degrees (recall that 15 degrees is one hour). The length of the day is computed using the sunzenith constraint. """ # Write header self._log_header2(gammalib.NORMAL, 'Hour angle weights') # Get MET time reference tref = gammalib.GTimeReference(self['mjdref'].real(),'s','TT','LOCAL') # Get time interval tmin = self['tmin'].time(tref) tmax = self['tmax'].time(tref) # Initialise hour angle list and results hour_angles = [] self._results = [] # Set number of hour angle bins and compute conversion factor and weight nsteps = 1000 ra2inx = float(nsteps)/360.0 # Conversion from RA (deg) to index weight = 24.0/float(nsteps) # Weight per hour angle step # Initialise hour angle array and weight hours = [0.0 for i in range(nsteps)] # Initialise time and loop until the end time is reached time = tmin while time <= tmax: # Initialise hours array for this time step hours_time = [weight for i in range(nsteps)] # Compute half the length of the exclusion intervals in degrees sun_dra = self._sun_ra_exclusion(time) moon_dra = self._moon_ra_exclusion(time) # Compute Right Ascension and Declination of Sun in degrees and # derive Sun exclusion interval sun = gammalib.GSkyDir() sun.sun(time) sun_ra = sun.ra_deg() sun_dec = sun.dec_deg() sun_ra_start = sun_ra - sun_dra sun_ra_stop = sun_ra + sun_dra # Adjust interval sun_ra_start, sun_ra_stop = self._adjust_ra_interval(sun_ra_start, sun_ra_stop) # Compute Right Ascension and Declination of Moon in degrees and # derive Moon exclusion interval moon = gammalib.GSkyDir() moon.moon(time) moon_ra = moon.ra_deg() moon_dec = moon.dec_deg() moon_ra_start = moon_ra - moon_dra moon_ra_stop = moon_ra + moon_dra # Adjust interval moon_ra_start, moon_ra_stop = self._adjust_ra_interval(moon_ra_start, moon_ra_stop) # Compute Moon elongation and illumation fraction (Moon phase) elongation = sun.dist_deg(moon) fli = (1.0 - math.cos(elongation * gammalib.deg2rad))/2.0 # Exclude hours due to Sun constraint (this defines the night) hours_time = self._exclude_ra_interval(hours_time, sun_ra_start, sun_ra_stop) # Exclude hours due to Moon constraint if the illumination fraction # is equal to or above the maximim fraction of illumination if fli >= self['maxfli'].real(): hours_time = self._exclude_ra_interval(hours_time, moon_ra_start, moon_ra_stop) # Add hours dark_time = 0.0 for i in range(nsteps): hours[i] += hours_time[i] dark_time += hours_time[i] # Set result record result = {'time': time.copy(), 'sun_ra': sun_ra, 'sun_dec': sun_dec, 'moon_ra': moon_ra, 'moon_dec': moon_dec, 'elongation': elongation, 'fli': fli, 'sun_ra_start': sun_ra_start, 'sun_ra_stop': sun_ra_stop, 'moon_ra_start': moon_ra_start, 'moon_ra_stop': moon_ra_stop, 'dark_time' : dark_time} # Append result record to results self._results.append(result) # Log results logs = 'Sun=(%8.3f,%7.3f) Moon=(%8.3f,%7.3f) '\ 'Sun_RA_excl=[%8.3f,%8.3f] Moon_RA_excl=[%8.3f,%8.3f] '\ 'FLI=%4.2f Dark time=%5.2f h' % \ (sun_ra, sun_dec, moon_ra, moon_dec, sun_ra_start, sun_ra_stop, moon_ra_start, moon_ra_stop, fli, dark_time) self._log_value(gammalib.VERBOSE, time.utc(), logs) # Add seconds of one day and start with next day time += 86400.0 # Build hour angle list dh = 360.0/float(nsteps) for i in range(nsteps): hour_angle = {'angle': dh*float(i), 'hours': hours[i]} hour_angles.append(hour_angle) # Log hour angle weights total = 0.0 for hour_angle in hour_angles: total += hour_angle['hours'] self._log_value(gammalib.EXPLICIT, 'Observing time', str(total)+' h') self._log_value(gammalib.EXPLICIT, 'Number of hour angle bins', len(hour_angles)) # Return hour angle weights return hour_angles def _visibility_cube(self): """ Compute visibility cube Compute visibility cube by displacing the zenith angle map for all hour angles. The visibility cube contains the number of hours a given celestial position is visible under a given zenith angle interval. Summing over all zenith angle intervals specifies for how long a given celestial position will be visible. """ # Write header self._log_header1(gammalib.TERSE, 'Compute visibility cube') # Set visibility cube and zenith angle map dimensions and bin size binsz = self['binsz'].real() nx = int(360.0/binsz+0.5) ny = int(180.0/binsz+0.5) dx = 360.0/float(nx) dy = 180.0/float(ny) # Set zenith angle dimension and bin size dz = self['dz'].real() zmax = self['zmax'].real() nz = int(zmax/dz+0.5) dz = zmax/float(nz) # Log dimensions self._log_header2(gammalib.NORMAL, 'Visibility cube dimensions') self._log_value(gammalib.NORMAL, 'Number of RA bins', nx) self._log_value(gammalib.NORMAL, 'Number of Dec bins', ny) self._log_value(gammalib.NORMAL, 'Number of zenith angles', nz) self._log_value(gammalib.NORMAL, 'RA pixel size', str(dx)+' deg') self._log_value(gammalib.NORMAL, 'Dec pixel size', str(dy)+' deg') self._log_value(gammalib.NORMAL, 'Zenith angle bin size', str(dz)+' deg') self._log_value(gammalib.NORMAL, 'Maximum zenith angle', str(zmax)+' deg') # Compute zenith angle map zmap = self._zenith_angle_map(nx,ny,dx,dy) # Setup hour angle weights. They specify for how many hours a given # hour angle is observed during the covered time period. hour_angles = self._hour_angle_weight() # Initialise visibility cube self._cube = gammalib.GSkyMap('CAR','CEL',0.0,0.0,-dx,dy,nx,ny,nz) # Loop over all hour angle weights for hour_angle in hour_angles: # Compute temporal shift in degrees shift = hour_angle['angle'] # Initialise shifted zenith angle map zmap_shift = gammalib.GSkyMap('CAR','CEL',shift,0.0,-dx,dy,nx,ny) # Merge zenith angle map into shifted map zmap_shift += zmap # Loop over all pixels of the shifted map and add the hours during # which the shifted map occurs to the relevant zenith angle bin for i, zenith in enumerate(zmap_shift): iz = int(zenith/dz) if iz < nz: self._cube[i,iz] += hour_angle['hours'] # Return return def _save_results(self, outfile, clobber): """ Save results in VISIBILITY FITS table Parameters ---------- outfile : str Result FITS file name clobber : bool Overwrite existing file? """ # Create FITS table columns nrows = len(self._results) time = gammalib.GFitsTableStringCol('Time', nrows, 20) mjd = gammalib.GFitsTableDoubleCol('MJD', nrows) sun_ra = gammalib.GFitsTableDoubleCol('RA_sun', nrows) sun_dec = gammalib.GFitsTableDoubleCol('DEC_sun', nrows) moon_ra = gammalib.GFitsTableDoubleCol('RA_moon', nrows) moon_dec = gammalib.GFitsTableDoubleCol('DEC_moon', nrows) sun_ra_start = gammalib.GFitsTableDoubleCol('RA_sun_start', nrows) sun_ra_stop = gammalib.GFitsTableDoubleCol('RA_sun_stop', nrows) moon_ra_start = gammalib.GFitsTableDoubleCol('RA_moon_start', nrows) moon_ra_stop = gammalib.GFitsTableDoubleCol('RA_moon_stop', nrows) elongation = gammalib.GFitsTableDoubleCol('Elongation', nrows) fli = gammalib.GFitsTableDoubleCol('FLI', nrows) dark_time = gammalib.GFitsTableDoubleCol('Darktime', nrows) # Set units of table columns mjd.unit('days') sun_ra.unit('deg') sun_dec.unit('deg') moon_ra.unit('deg') moon_dec.unit('deg') sun_ra_start.unit('deg') sun_ra_stop.unit('deg') moon_ra_start.unit('deg') moon_ra_stop.unit('deg') elongation.unit('deg') dark_time.unit('hours') # File FITS table columns for i, result in enumerate(self._results): time[i] = result['time'].utc() mjd[i] = result['time'].mjd() sun_ra[i] = result['sun_ra'] sun_dec[i] = result['sun_dec'] moon_ra[i] = result['moon_ra'] moon_dec[i] = result['moon_dec'] sun_ra_start[i] = result['sun_ra_start'] sun_ra_stop[i] = result['sun_ra_stop'] moon_ra_start[i] = result['moon_ra_start'] moon_ra_stop[i] = result['moon_ra_stop'] elongation[i] = result['elongation'] fli[i] = result['fli'] dark_time[i] = result['dark_time'] # Initialise FITS Table with extension "VISIBILITY" table = gammalib.GFitsBinTable(nrows) table.extname('VISIBILITY') # Add Header for compatibility with gammalib.GMWLSpectrum table.card('INSTRUME', 'CTA', 'Name of Instrument') table.card('TELESCOP', 'CTA', 'Name of Telescope') # Append filled columns to fits table table.append(time) table.append(mjd) table.append(sun_ra) table.append(sun_dec) table.append(moon_ra) table.append(moon_dec) table.append(sun_ra_start) table.append(sun_ra_stop) table.append(moon_ra_start) table.append(moon_ra_stop) table.append(elongation) table.append(fli) table.append(dark_time) # Append table to result FITS file fits = gammalib.GFits(outfile) fits.append(table) fits.save(clobber) # Return return # Public methods def run(self): """ Run the script """ # Switch screen logging on in debug mode if self._logDebug(): self._log.cout(True) # Get parameters self._get_parameters() # Compute visibility cube self._visibility_cube() # Optionally publish map if self['publish'].boolean(): self.publish() # Return return def save(self): """ Save the visibility cube """ # Write header self._log_header1(gammalib.TERSE, 'Save visibility cube') # Get outfile parameter outfile = self['outfile'].filename() # Log file name self._log_value(gammalib.NORMAL, 'Visibility cube file', outfile.url()) # Save the visibility cube self._cube.save(outfile, self['clobber'].boolean()) # Save the results in VISIBILITY table extension self._save_results(outfile, self['clobber'].boolean()) # Return return def publish(self, name=''): """ Publish visibility cube Parameters ---------- name : str, optional Name of visibility cube """ # Write header self._log_header1(gammalib.TERSE, 'Publish visibility cube') # Set default name is user name is empty if not name: user_name = self._name() else: user_name = name # Log cube name self._log_value(gammalib.TERSE, 'Visibility cube name', user_name) # Publish cube self._cube.publish(user_name) # Return return # ======================== # # Main routine entry point # # ======================== # if __name__ == '__main__': # Create instance of application app = csviscube(sys.argv) # Execute application app.execute()