Implementation of visibility computation using the Sun constraint (Moon constraint is still missing)

Jürgen Knödlseder
1 parent 28cac9e0
Showing 3 changed files with 330 additions and 103 deletions
cscripts/csviscube.par
cscripts/csviscube.py
test/test_csviscube.py
@@ -28,7 +28,9 @@ outfile,    f, a, viscube.fits,,, &quot;Visibility cube output file&quot;
 #
 # Visibility cube parameters
 #===========================
-binsz, r, h, 1.0,,, "Image scale (in degrees/pixel)"
+binsz, r, h, 1.0,,, "Image scale (deg/pixel)"
+dz,    r, h, 1.0,,, "Zenith angle bin size (deg)"
+zmax,  r, h, 60.0,,, "Maximum zenith angle (deg)"
  
 #
 # Standard parameters
@@ -37,7 +37,7 @@ class csviscube(ctools.cscript):
     # Constructor
     def __init__(self, *argv):
         """
-        Constructor.
+        Constructor
  
         Parameters
         ----------
@@ -54,8 +54,10 @@ class csviscube(ctools.cscript):
         self._name    = 'csviscube'
         self._version = '1.0.0'
  
-        # Initialise application by calling the appropriate class
-        # constructor.
+        # Initialise members
+        self._cube = gammalib.GSkyMap()
+
+        # Initialise application by calling the appropriate class constructor
         self._init_cscript(argv)
  
         # Return
@@ -72,10 +74,7 @@ class csviscube(ctools.cscript):
         self['tmax'].real()
         self['geolon'].real()
         self['geolat'].real()
-        self['sunzenith'].real()
-        self['moonzenith'].real()
         self['outfile'].filename()
-        self['binsz'].real()
  
         #  Write input parameters into logger
         self._log_parameters(gammalib.TERSE)
@@ -83,38 +82,301 @@ class csviscube(ctools.cscript):
         # Return
         return
  
-    def _hour_angle_weight(self):
+    def _zenith_angle_map(self, nx, ny, dx, dy):
+        """
+        Compute zenith angle map
+
+        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
+        -------
+        map : `~gammalib.GSkyMap`
+            Allsky map comprising the zenith angle for an hour angle of 0.
+
+        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.
         """
-        Compute hour angle weight vector
+        # Initialise zenith angle map
+        map = 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
+                map[index] = zenith
+                index     += 1
+
+        # Log zenith angle map
+        self._log_header2(gammalib.EXPLICIT, 'Zenith angle map')
+        self._log_string(gammalib.EXPLICIT, str(map))
+
+        # Return zenith angle map
+        return map
+
+    def _sun_radec(self, time):
+        """
+        Compute Right Ascension and Declination of the Sun
+
+        Parameters
+        ----------
+        time : `~gammalib.GTime()`
+            Time for which declination is to be computed
+
+        Returns
+        -------
+        ra, dec : tuple of float
+            Right Ascension and Declination of the Sun in degrees
+
+        The formulae are inspired from
+        https://en.wikipedia.org/wiki/Position_of_the_Sun
+        """
+        # Compute number of days since Greenwich noon, Terrestrial Time, on
+        # 1 January 2000
+        n = time.jd() - 2451545.0
+
+        # Compute mean longitude of the Sun in degrees, corrected for the
+        # aberration of light
+        L = 280.460 + 0.9856474 * n
+        while L < 0.0:
+            L += 360.0
+        while L >= 360.0:
+            L -= 360.0
+
+        # Compute the mean anomaly of the Sun in radians
+        g = (357.528 + 0.9856003 * n) * gammalib.deg2rad
+
+        # Compute the ecliptic longitude of the Sun in degrees
+        lam = L + 1.915 * math.sin(g) + 0.020 * math.sin(2.0*g)
+
+        # Compute Right Ascension and Declination of the Sun in degrees
+        ra  = math.atan(math.cos(23.43711*gammalib.deg2rad) *
+                        math.tan(lam*gammalib.deg2rad))*gammalib.rad2deg
+        dec = math.asin(math.sin(23.43711*gammalib.deg2rad) *
+                        math.sin(lam*gammalib.deg2rad))*gammalib.rad2deg
+
+        # Put ra in the right quadrant (same quadrant as lam)
+        while lam < 0.0:
+            lam += 360.0
+        while lam >= 360.0:
+            lam -= 360.0
+        while ra < 0.0:
+            ra += 360.0
+        while ra >= 360.0:
+            ra -= 360.0
+        q_lam = int(lam / 90.0)
+        q_ra  = int(ra / 90.0)
+        ra   += (q_lam-q_ra)*90.0
+
+        # Return declination of the Sun
+        return ra, dec
+
+    def _sun_ra_exclusion(self, time):
+        """
+        Compute the Right Ascension difference in degrees before and after
+        noon that are excluded due to the Sun constraint
+
+        Parameters
+        ----------
+        time : `~gammalib.GTime()`
+            Time for which the Sun exclusion is to be computed
+
+        Returns
+        -------
+        dra : float
+            Positive Right Ascension difference in 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
+        sunra, sundec = self._sun_radec(time)
+        sundec *= gammalib.deg2rad
+
+        # 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 = math.acos((cos_sunzenith - sin_geolat * sin_sundec) /
+                        (cos_geolat * cos_sundec)) * gammalib.rad2deg
+
+        # Return
+        return dra
+
+    def _hour_angle_weight(self):
+        """
+        Compute hour angle weights
+
+        Computes an array specifying how many hours the array was observing
+        for a given hour angle during the time interval [tmin,tmax]. The hour
+        angle runs from 0 to 360 degrees.
+
         Compute during which time the array was observing due to hour angle
         constraints. We take here a very simple model where the night lasts
         between 6 and 10 hours (or 90-150 degrees), with a sinusoidal variation
         along the year. Every day the hour angle of the Sun set is advancing by
         4 minutes.
         """
-        # Set number of hour angle bins
+        # Write header
+        self._log_header2(gammalib.NORMAL, 'Hour angle weights')
+
+        # Get array geographic longitude
+        geolon = self['geolon'].real()
+
+        # Get time interval
+        tmin = gammalib.GTime(self['tmin'].real(), 's')
+        tmax = gammalib.GTime(self['tmax'].real(), 's')
+
+        # Initialise hour angle list
+        hour_angles = []
+        
+        # Set number of hour angle bins and compute conversion factor and weight
         nsteps = 1000
-        d_min  =  6.0/24.0*nsteps
-        d_max  = 10.0/24.0*nsteps
+        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)]
-        weight = 24.0/float(nsteps)  # Weight per hour angle step
  
-        # Loop over all days in one year
-        for day in range(nsteps):
-            h_set  = day
-            h_rise = day + \
-                     int(d_min + (d_max-d_min)*math.sin(day/float(nsteps)*gammalib.pi))
-            for i in range(h_set,h_rise):
-                if i >= nsteps:
-                    hours[i-nsteps] += weight
-                else:
+        # Initialise time and loop until the end time is reached
+        time = tmin
+        while time <= tmax:
+
+            # Compute Right Ascension and Declination of Sun in degrees
+            sunra, sundec = self._sun_radec(time)
+
+            # Compute by how much the zenith angle map needs to be shifted
+            # to correspond to the actual time
+            #
+            # NOTE: The local apparent siderial time only is relevant if
+            # the time interval is shorter than a day. Only in that case
+            # we have to set the corresponding hour angles to zero. We need
+            # to think how to properly implement that.
+            #last = time.last(geolon) * 15.0
+            #print(last)
+
+            # Compute the time from now when the Sun will culminate. The time
+            # is here expressed in degrees
+            dra_sun = self._sun_ra_exclusion(time)
+            
+            # Set [0,ra_start] and [ra_stop,360.0]
+            ra_start = sunra - dra_sun
+            ra_stop  = sunra + dra_sun
+            
+            # Case 1: The RA interval during which the Sun is above the zenith
+            # angle constraint is fully comprised within the [0,360] interval.
+            # In that case the dark time is comprised of two intervals:
+            # [0,ra_start] and [ra_stop,360]
+            if ra_start >= 0.0 and ra_stop <= 360.0:
+                inx_stop  = int(ra_start * ra2inx + 0.5)
+                inx_start = int(ra_stop  * ra2inx + 0.5)
+                for i in range(0,inx_stop):
+                    hours[i] += weight
+                for i in range(inx_start,nsteps):
                     hours[i] += weight
  
-        # Return
-        return hours
+                # Log setting
+                logs = 'Sun=(%.3f,%.3f) dRA=%.3f RA_excl=[%.3f,%.3f] '\
+                       'Indices=[0-%d] & [%d-%d] Dark time=%.2f h %s' % \
+                       (sunra, sundec, dra_sun, ra_start, ra_stop, inx_stop, inx_start,
+                        nsteps, float(inx_stop+(nsteps-inx_start))*weight,
+                        '(Case 1)')
+                self._log_value(gammalib.VERBOSE, time.utc(), logs)
+
+            # Case 2: The RA interval during which the Sun is above the zenith
+            # angle constraint is starting at negative RA. In that case the
+            # dark time is comprised of a single interval
+            # [ra_stop,ra_start+360]
+            elif ra_start < 0.0:
+                inx_start = int(ra_stop * ra2inx + 0.5)
+                inx_stop  = int((ra_start+360.0) * ra2inx + 0.5)
+                for i in range(inx_start,inx_stop):
+                    hours[i] += weight
+
+                # Log setting
+                logs = 'Sun=(%.3f,%.3f) dRA=%.3f RA_excl=[%.3f,%.3f] '\
+                       'Indices=[%d-%d] Dark time=%.2f h %s' % \
+                       (sunra, sundec, dra_sun, ra_start, ra_stop, inx_start, inx_stop,
+                        float(inx_stop-inx_start)*weight, '(Case 2)')
+                self._log_value(gammalib.VERBOSE, time.utc(), logs)
+
+            # Case 3: The RA interval during which the Sun is above the zenith
+            # angle constraint is stopping at RA>360. In that case the dark
+            # time is comprised of a single interval [ra_stop-360,ra_start]
+            elif ra_stop > 360.0:
+                inx_start = int((ra_stop-360.0) * ra2inx + 0.5)
+                inx_stop  = int(ra_start * ra2inx + 0.5)
+                for i in range(inx_start,inx_stop):
+                    hours[i] += weight
+
+                # Log setting
+                logs = 'Sun=(%.3f,%.3f) dRA=%.3f RA_excl=[%.3f,%.3f] '\
+                       'Indices=[%d-%d] Dark time=%.2f h %s' % \
+                       (sunra, sundec, dra_sun, ra_start, ra_stop, inx_start, inx_stop,
+                        float(inx_stop-inx_start)*weight, '(Case 3)')
+                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):
         """
@@ -126,102 +388,64 @@ class csviscube(ctools.cscript):
         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
-        map = self._zenith_angle_map()
-        
-        # Set zenith angle range and bin size
-        dz = 1.0
-        nz = int(60.0/dz)
-    
-        # Initialise visibility cube
-        self._cube = gammalib.GSkyMap('CAR','CEL',0.0,0.0,-1.0,1.0,360,180,nz)
+        map = self._zenith_angle_map(nx,ny,dx,dy)
  
         # Setup hour angle weights. They specify for how many hours a given
-        # hour angle will be observed during one year.
-        hours = self._hour_angle_weight()
+        # hour angle is observed during the covered time period.
+        hour_angles = self._hour_angle_weight()
  
-        # Compute normalisation factor
-        norm = 365.0/float(len(hours))
+        # 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 h, weight in enumerate(hours):
+        for hour_angle in hour_angles:
  
-            # Compute temporal shift
-            shift = 360.0/float(len(hours)) * float(h)
+            # Compute temporal shift in degrees
+            shift = hour_angle['angle']
  
             # Initialise shifted zenith angle map
-            map_shift = gammalib.GSkyMap('CAR','CEL',shift,0.0,-1.0,1.0,360,180)
+            map_shift = gammalib.GSkyMap('CAR','CEL',shift,0.0,-dx,dy,nx,ny)
  
             # Merge zenith angle map into shifted map
             map_shift += map
  
-            # Loop over all pixels of shifted map
+            # 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(map_shift):
-                iz     = int(zenith/dz)
+                iz = int(zenith/dz)
                 if iz < nz:
-                    self._cube[i,iz] += weight * norm
+                    self._cube[i,iz] += hour_angle['hours']
  
         # Return
         return
  
-    def _zenith_angle_map(self):
-        """
-        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
-
-        LST is the difference between the Greenwich siderial time (GST) and the
-        geographic longitude of the observer
-    
-        LST = GST - lon
-    
-        where lon is counted positive west from the meridian. Hence
-    
-        h = GST - lon - 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 the
-        local meridian.
-        """
-        # Initialise zenith angle map
-        map = gammalib.GSkyMap('CAR','CEL',0.0,0.0,-1.0,1.0,360,180)
-
-        # Set hour angle and declination vectors
-        hours = [float(i) for i in range(360)]
-        decs  = [float(i)-89.5 for i in range(180)]
-
-        # Get array geographic longitude and latitude
-        geolon = self['geolon'].real()
-        geolat = self['geolat'].real()
-
-        # Precompute latitude terms
-        cos_lat = math.cos(geolat*gammalib.deg2rad)
-        sin_lat = math.sin(geolat*gammalib.deg2rad)
-
-        # Loop over all declination 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
-                map[index] = zenith
-                index     += 1
-
-        # Return zenith angle map
-        return map
-
  
     # Public methods
     def run(self):
@@ -252,7 +476,7 @@ class csviscube(ctools.cscript):
         outfile = self['outfile'].filename()
  
         # Log file name
-        self._log_value(gammalib.TERSE, 'Visibility cube file', outfile.url())
+        self._log_value(gammalib.NORMAL, 'Visibility cube file', outfile.url())
  
         # Save the visibility cube
         self._cube.save(outfile, self['clobber'].boolean())
@@ -67,7 +67,7 @@ class Test(test):
  
         # Setup csviscube command
         cmd = csviscube+' tmin=0.0 tmax=31557600.0 geolon=79.4041'+\
-                        ' geolat=-24.6272 outfile=csviscube_cmd1.fits'+\
+                        ' geolat=-24.6272 binsz=3.0 outfile=csviscube_cmd1.fits'+\
                         ' logfile="csviscube_cmd1.log" chatter=1'
  
         # Check if execution of wrong command fails
@@ -92,6 +92,7 @@ class Test(test):
         viscube['tmax']    =  31557600.0
         viscube['geolon']  =  79.4041
         viscube['geolat']  = -24.6272
+        viscube['binsz']   = 3.0
         viscube['outfile'] = 'csviscube_py1.fits'
         viscube['logfile'] = 'csviscube_py1.log'
         viscube['chatter'] = 2