diff --git a/gwemopt/io/skymap.py b/gwemopt/io/skymap.py index df3e3aa..9bf59e9 100644 --- a/gwemopt/io/skymap.py +++ b/gwemopt/io/skymap.py @@ -267,7 +267,7 @@ def read_skymap(params, map_struct=None): if params["do_3d"]: try: healpix_data, header = hp.read_map( - filename, field=(0, 1, 2, 3), verbose=False, h=True + filename, field=(0, 1, 2, 3), h=True ) except: skymap = read_sky_map(filename, moc=True, distances=True) diff --git a/gwemopt/moc.py b/gwemopt/moc.py index 219257c..796435f 100644 --- a/gwemopt/moc.py +++ b/gwemopt/moc.py @@ -55,7 +55,7 @@ def create_moc(params, map_struct=None): delayed(Fov2Moc)( params, config_struct, telescope, tess[1], tess[2], nside ) - for tess in tesselation + for tess in tqdm(tesselation) ) for ii, tess in tqdm(enumerate(tesselation), total=len(tesselation)): index, ra, dec = tess[0], tess[1], tess[2] diff --git a/gwemopt/plotting/style.py b/gwemopt/plotting/style.py index 26eded0..3275c14 100644 --- a/gwemopt/plotting/style.py +++ b/gwemopt/plotting/style.py @@ -3,7 +3,7 @@ import matplotlib import matplotlib.pyplot as plt import numpy as np -from astropy.coordinates import get_moon, get_sun +from astropy.coordinates import get_body from astropy.time import Time cmap = "cylon" @@ -21,7 +21,7 @@ def add_edges(): """ Add longitude and latitude lines to a healpix map. """ - hp.graticule(verbose=False) + hp.graticule() plt.grid(True) lons = np.arange(-150.0, 180, 30.0) lats = np.zeros(lons.shape) @@ -41,8 +41,8 @@ def add_sun_moon(params): start_time = params["gpstime"] start_time = Time(start_time, format="gps", scale="utc") - sun_position = get_sun(start_time) - moon_position = get_moon(start_time) + sun_position = get_body("sun", start_time) + moon_position = get_body("moon", start_time) hp.visufunc.projscatter( sun_position.ra, sun_position.dec, lonlat=True, color="yellow" @@ -61,7 +61,7 @@ def add_sun_moon(params): new_time = Time(new_time, format="gps", scale="utc") - new_moon_position = get_moon(new_time) + new_moon_position = get_body("moon", new_time) hp.visufunc.projscatter( new_moon_position.ra, new_moon_position.dec, @@ -73,7 +73,7 @@ def add_sun_moon(params): if not i % 8: # only plot point for the sun every two days in order to avoid overlap - new_sun_position = get_sun(new_time) + new_sun_position = get_body("sun", new_time) hp.visufunc.projscatter( new_sun_position.ra, new_sun_position.dec, diff --git a/gwemopt/scheduler.py b/gwemopt/scheduler.py index 4f7e4d5..b974495 100644 --- a/gwemopt/scheduler.py +++ b/gwemopt/scheduler.py @@ -2,6 +2,7 @@ import astropy.coordinates import astropy.units as u +import ephem import ligo.segments as segments import numpy as np from astropy.time import Time @@ -11,19 +12,18 @@ from gwemopt.utils import angular_distance -def get_altaz_tiles(ras, decs, observatory, obstime): - # Convert to RA, Dec. - radecs = astropy.coordinates.SkyCoord( - ra=np.array(ras) * u.degree, dec=np.array(decs) * u.degree, frame="icrs" - ) +def get_altaz_tiles(ra, dec, observer, obstime): - # Alt/az reference frame at observatory, now - frame = astropy.coordinates.AltAz(obstime=obstime, location=observatory) + observer.date = ephem.Date(obstime.iso) - # Transform grid to alt/az coordinates at observatory, now - altaz = radecs.transform_to(frame) + fxdbdy = ephem.FixedBody() + fxdbdy._ra = ephem.degrees(str(ra)) + fxdbdy._dec = ephem.degrees(str(dec)) + fxdbdy.compute(observer) - return altaz + return float(repr(fxdbdy.alt)) * (360 / (2 * np.pi)), float(repr(fxdbdy.az)) * ( + 360 / (2 * np.pi) + ) def find_tile( @@ -89,7 +89,7 @@ def find_tile( def get_order( - params, tile_struct, tilesegmentlists, exposurelist, observatory, config_struct + params, tile_struct, tilesegmentlists, exposurelist, observer, config_struct ): """ tile_struct: dictionary. key -> struct info. @@ -219,8 +219,9 @@ def get_order( for ii in np.arange(len(exposurelist)): # first, create an array of airmass-weighted probabilities t = Time(exposurelist[ii][0], format="mjd") - altaz = get_altaz_tiles(ras, decs, observatory, t) - alts = altaz.alt.degree + alts, azs = [ + get_altaz_tile(ra, dec, observer, t) for ra, dec in zip(ras, decs) + ] horizon = config_struct["horizon"] horizon_mask = alts <= horizon airmass = 1 / np.cos((90.0 - alts) * np.pi / 180.0) @@ -442,10 +443,13 @@ def scheduler(params, config_struct, tile_struct): if params["tilesType"] == "galaxy": coverage_struct["galaxies"] = [] - observatory = astropy.coordinates.EarthLocation( - lat=config_struct["latitude"] * u.deg, - lon=config_struct["longitude"] * u.deg, - height=config_struct["elevation"] * u.m, + observer = ephem.Observer() + observer.lat = str(config_struct["latitude"]) + observer.lon = str(config_struct["longitude"]) + observer.horizon = str(config_struct["horizon"]) + observer.elevation = config_struct["elevation"] + observer.horizon = ephem.degrees( + str(90 - np.arccos(1 / params["airmass"]) * 180 / np.pi) ) exposurelist = config_struct["exposurelist"] @@ -455,9 +459,8 @@ def scheduler(params, config_struct, tile_struct): for key in keys: # segments.py: tile_struct[key]["segmentlist"] is a list of segments when the tile is available for observation tilesegmentlists.append(tile_struct[key]["segmentlist"]) - print("Generating schedule order...") keys, filts = get_order( - params, tile_struct, tilesegmentlists, exposurelist, observatory, config_struct + params, tile_struct, tilesegmentlists, exposurelist, observer, config_struct ) if params["doPlots"]: @@ -504,10 +507,9 @@ def scheduler(params, config_struct, tile_struct): # calculate airmass for each tile at the start of its exposure: t = Time(mjd_exposure_start, format="mjd") - altaz = get_altaz_tiles( - tile_struct_hold["ra"], tile_struct_hold["dec"], observatory, t + alt, az = get_altaz_tiles( + tile_struct_hold["ra"], tile_struct_hold["dec"], observer, t ) - alt = altaz.alt.degree airmass = 1 / np.cos((90.0 - alt) * np.pi / 180) # total duration of the observation (?) diff --git a/gwemopt/segments.py b/gwemopt/segments.py index 635e0b7..2237aff 100644 --- a/gwemopt/segments.py +++ b/gwemopt/segments.py @@ -7,9 +7,11 @@ import numpy as np from astropy.time import Time from joblib import Parallel, delayed +from tqdm import tqdm from gwemopt.utils.geometry import angular_distance from gwemopt.utils.misc import get_exposures +from gwemopt.utils.sidereal_time import greenwich_sidereal_time def get_telescope_segments(params): @@ -38,13 +40,7 @@ def get_telescope_segments(params): return params -def get_moon_segments(config_struct, segmentlist, observer, fxdbdy, radec): - if "moon_constraint" in config_struct: - moon_constraint = float(config_struct["moon_constraint"]) - else: - moon_constraint = 20.0 - - moonsegmentlist = segments.segmentlist() +def get_moon_radecs(segmentlist, observer): dt = 1.0 / 24.0 tt = np.arange(segmentlist[0][0], segmentlist[-1][1] + dt, dt) conv = ( @@ -52,30 +48,41 @@ def get_moon_segments(config_struct, segmentlist, observer, fxdbdy, radec): ) # conversion between MJD (tt) and DJD (what ephem uses) tt_DJD = tt - conv - ra2 = radec.ra.radian - d2 = radec.dec.radian - + moon_radecs = [] # Where is the moon? moon = ephem.Moon() for ii in range(len(tt) - 1): observer.date = ephem.Date(tt_DJD[ii]) moon.compute(observer) - fxdbdy.compute(observer) - alt_target = float(repr(fxdbdy.alt)) * (360 / (2 * np.pi)) - az_target = float(repr(fxdbdy.az)) * (360 / (2 * np.pi)) - # print("Altitude / Azimuth of target: %.5f / %.5f"%(alt_target,az_target)) + # Coverting both target and moon ra and dec to radians + ra1 = float(repr(moon.ra)) + d1 = float(repr(moon.dec)) + + moon_radecs.append([ra1, d1]) - alt_moon = float(repr(moon.alt)) * (360 / (2 * np.pi)) - az_moon = float(repr(moon.az)) * (360 / (2 * np.pi)) - # print("Altitude / Azimuth of moon: %.5f / %.5f"%(alt_moon,az_moon)) + return moon_radecs - ra_moon = (180 / np.pi) * float(repr(moon.ra)) - dec_moon = (180 / np.pi) * float(repr(moon.dec)) +def get_moon_segments(config_struct, segmentlist, moon_radecs, radec): + if "moon_constraint" in config_struct: + moon_constraint = float(config_struct["moon_constraint"]) + else: + moon_constraint = 20.0 + + moonsegments = [] + moonsegmentlist = segments.segmentlist() + dt = 1.0 / 24.0 + tt = np.arange(segmentlist[0][0], segmentlist[-1][1] + dt, dt) + + ra2 = radec.ra.radian + d2 = radec.dec.radian + + # Where is the moon? + for ii in range(len(tt) - 1): # Coverting both target and moon ra and dec to radians - ra1 = float(repr(moon.ra)) - d1 = float(repr(moon.dec)) + ra1 = moon_radecs[ii][0] + d1 = moon_radecs[ii][1] # Calculate angle between target and moon cosA = np.sin(d1) * np.sin(d2) + np.cos(d1) * np.cos(d2) * np.cos(ra1 - ra2) @@ -84,9 +91,12 @@ def get_moon_segments(config_struct, segmentlist, observer, fxdbdy, radec): # if angle >= 50.0*moon.moon_phase**2: if angle >= moon_constraint: - segment = segments.segment(tt[ii], tt[ii + 1]) - moonsegmentlist = moonsegmentlist + segments.segmentlist([segment]) - moonsegmentlist.coalesce() + moonsegments.append([tt[ii], tt[ii + 1]]) + + moonsegmentlist = segments.segmentlist() + for seg in moonsegments: + moonsegmentlist.append(segments.segment(seg)) + moonsegmentlist.coalesce() moonsegmentlistdic = segments.segmentlistdict() moonsegmentlistdic["observations"] = segmentlist @@ -97,7 +107,7 @@ def get_moon_segments(config_struct, segmentlist, observer, fxdbdy, radec): return moonsegmentlist -def get_ha_segments(config_struct, segmentlist, observer, fxdbdy, radec): +def get_ha_segments(config_struct, segmentlist, observer, radec): if "ha_constraint" in config_struct: ha_constraint = config_struct["ha_constraint"].split(",") ha_min = float(ha_constraint[0]) @@ -113,17 +123,14 @@ def get_ha_segments(config_struct, segmentlist, observer, fxdbdy, radec): 35.0 - radec.dec.deg ), 0.644981 * np.sqrt(35.0 - radec.dec.deg) - location = astropy.coordinates.EarthLocation( - config_struct["longitude"], - config_struct["latitude"], - config_struct["elevation"], - ) - halist = segments.segmentlist() for seg in segmentlist: mjds = np.linspace(seg[0], seg[1], 100) - tt = Time(mjds, format="mjd", scale="utc", location=location) - lst = tt.sidereal_time("mean") + tt = Time(mjds, format="mjd", scale="utc") + lst = astropy.coordinates.Longitude( + np.deg2rad(config_struct["longitude"]) + greenwich_sidereal_time(tt), + u.radian, + ) ha = (lst - radec.ra).hour idx = np.where((ha >= ha_min) & (ha <= ha_max))[0] if len(idx) >= 2: @@ -193,7 +200,7 @@ def get_segments(params, config_struct): return segmentlist -def get_segments_tile(config_struct, observatory, radec, segmentlist, airmass): +def get_segments_tile(config_struct, radec, segmentlist, moon_radecs, airmass): # check for empty segmentlist and immediately return if len(segmentlist) == 0: @@ -205,12 +212,11 @@ def get_segments_tile(config_struct, observatory, radec, segmentlist, airmass): observer.horizon = str(config_struct["horizon"]) observer.elevation = config_struct["elevation"] observer.horizon = ephem.degrees(str(90 - np.arccos(1 / airmass) * 180 / np.pi)) + observer.date = ephem.Date(Time(segmentlist[0][0], format="mjd", scale="utc").iso) fxdbdy = ephem.FixedBody() fxdbdy._ra = ephem.degrees(str(radec.ra.degree)) fxdbdy._dec = ephem.degrees(str(radec.dec.degree)) - - observer.date = ephem.Date(Time(segmentlist[0][0], format="mjd", scale="utc").iso) fxdbdy.compute(observer) date_start = ephem.Date(Time(segmentlist[0][0], format="mjd", scale="utc").iso) @@ -252,11 +258,9 @@ def get_segments_tile(config_struct, observatory, radec, segmentlist, airmass): # moonsegmentlist = get_skybrightness(\ # config_struct,segmentlist,observer,fxdbdy,radec) - halist = get_ha_segments(config_struct, segmentlist, observer, fxdbdy, radec) + halist = get_ha_segments(config_struct, segmentlist, observer, radec) - moonsegmentlist = get_moon_segments( - config_struct, segmentlist, observer, fxdbdy, radec - ) + moonsegmentlist = get_moon_segments(config_struct, segmentlist, moon_radecs, radec) tilesegmentlistdic = segments.segmentlistdict() tilesegmentlistdic["observations"] = segmentlist @@ -273,12 +277,6 @@ def get_segments_tile(config_struct, observatory, radec, segmentlist, airmass): def get_segments_tiles(params, config_struct, tile_struct): - observatory = astropy.coordinates.EarthLocation( - lat=config_struct["latitude"] * u.deg, - lon=config_struct["longitude"] * u.deg, - height=config_struct["elevation"] * u.m, - ) - segmentlist = config_struct["segmentlist"] print("Generating segments for tiles...") @@ -300,6 +298,17 @@ def get_segments_tiles(params, config_struct, tile_struct): ra=np.array(ras) * u.degree, dec=np.array(decs) * u.degree, frame="icrs" ) + observer = ephem.Observer() + observer.lat = str(config_struct["latitude"]) + observer.lon = str(config_struct["longitude"]) + observer.horizon = str(config_struct["horizon"]) + observer.elevation = config_struct["elevation"] + observer.horizon = ephem.degrees( + str(90 - np.arccos(1 / params["airmass"]) * 180 / np.pi) + ) + + moon_radecs = get_moon_radecs(segmentlist, observer) + if params["doParallel"]: tilesegmentlists = Parallel( n_jobs=params["Ncores"], @@ -307,16 +316,14 @@ def get_segments_tiles(params, config_struct, tile_struct): batch_size=int(len(radecs) / params["Ncores"]) + 1, )( delayed(get_segments_tile)( - config_struct, observatory, radec, segmentlist, params["airmass"] + config_struct, radec, segmentlist, moon_radecs, params["airmass"] ) - for radec in radecs + for radec in tqdm(radecs) ) for ii, key in enumerate(keys): tile_struct[key]["segmentlist"] = tilesegmentlists[ii] else: - for ii, key in enumerate(keys): - # if np.mod(ii,100) == 0: - # print("Generating segments for tile %d/%d"%(ii+1,len(radecs))) + for ii, key in tqdm(enumerate(keys), total=len(keys)): radec = radecs[ii] if params["doMinimalTiling"]: @@ -325,9 +332,9 @@ def get_segments_tiles(params, config_struct, tile_struct): radecs_computed = np.atleast_2d([radec.ra.value, radec.dec.value]) tilesegmentlist = get_segments_tile( config_struct, - observatory, radec, segmentlist, + moon_radecs, params["airmass"], ) tile_struct[key]["segmentlist"] = tilesegmentlist @@ -352,16 +359,16 @@ def get_segments_tiles(params, config_struct, tile_struct): ) tilesegmentlist = get_segments_tile( config_struct, - observatory, radec, segmentlist, + moon_radecs, params["airmass"], ) tile_struct[key]["segmentlist"] = tilesegmentlist else: tilesegmentlist = get_segments_tile( - config_struct, observatory, radec, segmentlist, params["airmass"] + config_struct, radec, segmentlist, moon_radecs, params["airmass"] ) tile_struct[key]["segmentlist"] = tilesegmentlist diff --git a/gwemopt/tiles.py b/gwemopt/tiles.py index f4280e6..a9aaa7c 100644 --- a/gwemopt/tiles.py +++ b/gwemopt/tiles.py @@ -9,6 +9,7 @@ from astropy.time import Time from scipy.stats import norm from shapely.geometry import MultiPoint +from tqdm import tqdm import gwemopt import gwemopt.moc @@ -589,7 +590,7 @@ def schedule_alternating( coverage_structs = [] maxidx = 0 - for i in range(len(exposuretimes)): + for i in tqdm(range(len(exposuretimes))): params["filters"] = [filters[i]] params["exposuretimes"] = [exposuretimes[i]] config_struct["exposurelist"] = segments.segmentlist( diff --git a/gwemopt/utils/__init__.py b/gwemopt/utils/__init__.py index 6417440..9fbf2df 100644 --- a/gwemopt/utils/__init__.py +++ b/gwemopt/utils/__init__.py @@ -8,4 +8,5 @@ getRectanglePixels, getSquarePixels, ) +from gwemopt.utils.sidereal_time import greenwich_sidereal_time from gwemopt.utils.treasuremap import get_treasuremap_pointings diff --git a/gwemopt/utils/sidereal_time.py b/gwemopt/utils/sidereal_time.py new file mode 100644 index 0000000..f1289c3 --- /dev/null +++ b/gwemopt/utils/sidereal_time.py @@ -0,0 +1,57 @@ +from datetime import datetime, timedelta +from math import pi + +import numpy as np + +GPS_EPOCH = datetime(1980, 1, 6, 0, 0, 0) +EPOCH_J2000_0_JD = 2451545.0 +_LEAP_SECONDS = np.asarray( + [ + 46828800, + 78364801, + 109900802, + 173059203, + 252028804, + 315187205, + 346723206, + 393984007, + 425520008, + 457056009, + 504489610, + 551750411, + 599184012, + 820108813, + 914803214, + 1025136015, + 1119744016, + 1167264017, + ] +) + + +def greenwich_sidereal_time(tt, equation_of_equinoxes=0): + """ + Compute the Greenwich mean sidereal time from the GPS time and equation of + equinoxes. + + Based on XLALGreenwichSiderealTime in lalsuite/lal/lib/XLALSiderealTime.c. + + Parameters + ---------- + tt: astropy.time.Time + The astropy time to convert + """ + t_hi = (tt.jd - EPOCH_J2000_0_JD) / 36525.0 + t_lo = (tt.gps % 1) / (36525.0 * 86400.0) + + t = t_hi + t_lo + + sidereal_time = ( + equation_of_equinoxes + (-6.2e-6 * t + 0.093104) * t**2 + 67310.54841 + ) + sidereal_time += 8640184.812866 * t_lo + sidereal_time += 3155760000.0 * t_lo + sidereal_time += 8640184.812866 * t_hi + sidereal_time += 3155760000.0 * t_hi + + return sidereal_time * pi / 43200.0