From db97908e4eecd708999e1e74e4424232bb6981ac Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Wed, 5 Oct 2016 17:10:48 +0200 Subject: [PATCH] Implementation of visibility computation using the Sun constraint (Moon constraint is still missing) --- cscripts/csviscube.par | 4 +++- cscripts/csviscube.py | 426 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------------------------------------------------------------------------------------------------- test/test_csviscube.py | 3 ++- 3 files changed, 330 insertions(+), 103 deletions(-) diff --git a/cscripts/csviscube.par b/cscripts/csviscube.par index ce09b9a..56a2116 100644 --- a/cscripts/csviscube.par +++ b/cscripts/csviscube.par @@ -28,7 +28,9 @@ outfile, f, a, viscube.fits,,, "Visibility cube output file" # # 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 diff --git a/cscripts/csviscube.py b/cscripts/csviscube.py index c7b3ad3..0b6a649 100755 --- a/cscripts/csviscube.py +++ b/cscripts/csviscube.py @@ -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()) diff --git a/test/test_csviscube.py b/test/test_csviscube.py index 6c6bc87..ffe7b61 100755 --- a/test/test_csviscube.py +++ b/test/test_csviscube.py @@ -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 -- libgit2 0.21.4