From 46c5fea95265b194dec8997e8a51c5cd0c668eab Mon Sep 17 00:00:00 2001 From: Michael Coughlin Date: Tue, 28 May 2024 18:38:16 +0200 Subject: [PATCH] MOC calculation (#154) This PR uses MOCs. --- gwemopt/args.py | 12 +- gwemopt/catalogs/get.py | 45 +- gwemopt/chipgaps/__init__.py | 4 +- gwemopt/chipgaps/decam.py | 17 +- gwemopt/chipgaps/ztf.py | 18 +- gwemopt/coverage.py | 430 +----------- gwemopt/efficiency.py | 72 +- gwemopt/io/__init__.py | 5 +- gwemopt/io/export_tiles.py | 10 - gwemopt/io/schedule.py | 57 +- gwemopt/io/skymap.py | 225 +++--- gwemopt/moc.py | 82 +-- gwemopt/params.py | 11 +- gwemopt/plotting/__init__.py | 21 - gwemopt/plotting/observability.py | 70 -- gwemopt/plotting/plot_coverage.py | 309 +++----- gwemopt/plotting/plot_efficiency.py | 12 - gwemopt/plotting/plot_skymap.py | 75 +- gwemopt/plotting/plot_tiles.py | 82 +-- gwemopt/plotting/style.py | 47 +- gwemopt/run.py | 40 -- gwemopt/scheduler.py | 184 +---- gwemopt/skyportal.py | 69 +- .../data/expected_results/coverage_DECam.dat | 10 +- .../data/expected_results/coverage_KPED.dat | 332 ++++----- .../data/expected_results/coverage_ZTF.dat | 4 +- gwemopt/tests/data/expected_results/map.dat | 1 - .../data/expected_results/schedule_DECam.dat | 14 +- .../data/expected_results/schedule_KPED.dat | 662 ++++++++---------- .../data/expected_results/schedule_ZTF.dat | 8 +- .../tests/data/expected_results/summary.dat | 6 +- .../expected_results/summary_coverage.dat | 6 +- .../tiles_coverage_int_DECam.txt | 14 +- .../tiles_coverage_int_KPED.txt | 662 ++++++++---------- .../tiles_coverage_int_ZTF.txt | 8 +- gwemopt/tests/test_schedule.py | 4 +- gwemopt/tiles.py | 454 ++---------- gwemopt/utils/pixels.py | 134 +--- gwemopt/utils/treasuremap.py | 16 +- pyproject.toml | 2 +- 40 files changed, 1337 insertions(+), 2897 deletions(-) delete mode 100644 gwemopt/plotting/observability.py delete mode 100644 gwemopt/tests/data/expected_results/map.dat diff --git a/gwemopt/args.py b/gwemopt/args.py index 96d531fd..8f547dab 100644 --- a/gwemopt/args.py +++ b/gwemopt/args.py @@ -37,7 +37,6 @@ def parse_args(args): parser.add_argument("--doPerturbativeTiling", action="store_true", default=False) parser.add_argument("--doOrderByObservability", action="store_true", default=False) - parser.add_argument("--doUseCatalog", action="store_true", default=False) parser.add_argument("--catalogDir", help="catalog directory", default=CATALOG_DIR) parser.add_argument( "--catalog", help="Galaxy catalog name (e.g GLADE)", default=None @@ -54,7 +53,6 @@ def parse_args(args): type=int, default=2000, ) - parser.add_argument("--doObservability", action="store_true", default=False) parser.add_argument("--doEfficiency", action="store_true", default=False) parser.add_argument( @@ -83,7 +81,6 @@ def parse_args(args): parser.add_argument("--Ndet", default=1, type=int) parser.add_argument("--nside", default=256, type=int) - parser.add_argument("--DScale", default=1.0, type=float) parser.add_argument("--Tobs", default="0.0,1.0") parser.add_argument("--mindiff", default=0.0, type=float) @@ -111,8 +108,6 @@ def parse_args(args): parser.add_argument("--doSingleExposure", action="store_true", default=False) parser.add_argument("--filters", default="r,g,r") parser.add_argument("--doAlternatingFilters", action="store_true", default=False) - parser.add_argument("--doRASlices", action="store_true", default=False) - parser.add_argument("--nside_down", default=2, type=int) parser.add_argument("--max_filter_sets", default=4, type=int) parser.add_argument("--iterativeOverlap", default=0.0, type=float) parser.add_argument("--maximumOverlap", default=1.0, type=float) @@ -137,10 +132,8 @@ def parse_args(args): parser.add_argument("--doBlocks", action="store_true", default=False) parser.add_argument("--Nblocks", default=4, type=int) - parser.add_argument("--doRASlice", action="store_true", default=False) - parser.add_argument("--raslice", default="0.0,24.0") - parser.add_argument("--absmag", default=-15.0, type=float) + parser.add_argument("--raslice", default=None) parser.add_argument( "--galactic_limit", help="Galactic limit.", default=0.0, type=float @@ -155,10 +148,9 @@ def parse_args(args): parser.add_argument("--true_inclination", default=0.0, type=float) parser.add_argument("--doTrueLocation", action="store_true", default=False) - parser.add_argument("--observability_thresh", default=0.05, type=float) - parser.add_argument("--doObservabilityExit", action="store_true", default=False) parser.add_argument("--inclination", action="store_true", default=False) + parser.add_argument("--projection", default="astro mollweide") parser.add_argument("--solverType", default="heuristic") parser.add_argument("--milpSolver", default="PULP_CBC_CMD") diff --git a/gwemopt/catalogs/get.py b/gwemopt/catalogs/get.py index 8f5cc145..6cd39f9c 100644 --- a/gwemopt/catalogs/get.py +++ b/gwemopt/catalogs/get.py @@ -2,8 +2,10 @@ import healpy as hp import numpy as np +from astropy.table import Table from astroquery.vizier import Vizier from ligo.skymap import distance +from ligo.skymap.bayestar import derasterize from scipy.stats import norm from gwemopt.catalogs.clu import CluCatalog @@ -63,7 +65,7 @@ def get_catalog(params, map_struct, export_catalog: bool = True): n, cl = params["powerlaw_n"], params["powerlaw_cl"] - prob_scaled = copy.deepcopy(map_struct["prob"]) + prob_scaled = copy.deepcopy(map_struct["skymap_raster"]["PROB"]) prob_sorted = np.sort(prob_scaled)[::-1] prob_indexes = np.argsort(prob_scaled)[::-1] prob_cumsum = np.cumsum(prob_sorted) @@ -72,34 +74,42 @@ def get_catalog(params, map_struct, export_catalog: bool = True): prob_scaled = prob_scaled**n ipix = hp.ang2pix( - map_struct["nside"], + params["nside"], np.array(cat_df["ra"]), np.array(cat_df["dec"]), lonlat=True, ) - if "distnorm" in map_struct: - if map_struct["distnorm"] is not None: + if "DISTNORM" in map_struct: + if map_struct["skymap_raster"]["DISTNORM"] is not None: # creat an mask to cut at 3 sigma in distance mask = np.zeros(len(cat_df["distmpc"])) - # calculate the moments from distmu, distsigma and distnorm - mom_mean, mom_std, mom_norm = distance.parameters_to_moments( - map_struct["distmu"], map_struct["distsigma"] - ) - condition_indexer = np.where( - (cat_df["distmpc"] < (mom_mean[ipix] + (3 * mom_std[ipix]))) - & (cat_df["distmpc"] > (mom_mean[ipix] - (3 * mom_std[ipix]))) + ( + cat_df["distmpc"] + < ( + map_struct["skymap_raster"]["DISTMEAN"][ipix] + + (3 * map_struct["skymap_raster"]["DISTSTD"][ipix]) + ) + ) + & ( + cat_df["distmpc"] + > ( + map_struct["skymap_raster"]["DISTMEAN"][ipix] + - (3 * map_struct["skymap_raster"]["DISTSTD"][ipix]) + ) + ) ) mask[condition_indexer] = 1 s_loc = ( prob_scaled[ipix] * ( - map_struct["distnorm"][ipix] + map_struct["skymap_raster"]["DISTNORM"][ipix] * norm( - map_struct["distmu"][ipix], map_struct["distsigma"][ipix] + map_struct["skymap_raster"]["DISTMU"][ipix], + map_struct["skymap_raster"]["DISTSIGMA"][ipix], ).pdf(cat_df["distmpc"]) ) ** params["powerlaw_dist_exp"] @@ -202,7 +212,7 @@ def get_catalog(params, map_struct, export_catalog: bool = True): s_mass = s_loc * (1 + (alpha_mass * beta_mass * s_mass)) s = np.array(s_loc * Slum * sdet) - prob = np.zeros(map_struct["prob"].shape) + prob = np.zeros(map_struct["skymap_raster"]["PROB"].shape) if params["galaxy_grade"] == "Sloc": for j in range(len(ipix)): prob[ipix[j]] += s_loc[j] @@ -220,11 +230,12 @@ def get_catalog(params, map_struct, export_catalog: bool = True): "You are trying to use a galaxy grade that is not implemented yet." ) + prob[np.isnan(prob)] = 0.0 prob = prob / np.sum(prob) - map_struct["prob_catalog"] = prob - if params["doUseCatalog"]: - map_struct["prob"] = prob + map_struct["skymap_raster"]["PROB"] = prob + skymap = map_struct["skymap_raster"].copy() + map_struct["skymap"] = derasterize(skymap) cat_df["grade"] = grade cat_df["S"] = s diff --git a/gwemopt/chipgaps/__init__.py b/gwemopt/chipgaps/__init__.py index 29fe3d2f..c5fac626 100644 --- a/gwemopt/chipgaps/__init__.py +++ b/gwemopt/chipgaps/__init__.py @@ -2,5 +2,5 @@ Module for chip gaps. """ -from gwemopt.chipgaps.decam import get_decam_quadrant_ipix -from gwemopt.chipgaps.ztf import get_ztf_quadrant_ipix +from gwemopt.chipgaps.decam import get_decam_quadrant_moc +from gwemopt.chipgaps.ztf import get_ztf_quadrant_moc diff --git a/gwemopt/chipgaps/decam.py b/gwemopt/chipgaps/decam.py index 2ad5a2a4..62d16a1c 100644 --- a/gwemopt/chipgaps/decam.py +++ b/gwemopt/chipgaps/decam.py @@ -24,6 +24,7 @@ from astropy import units as u from astropy.coordinates import ICRS, SkyCoord from astropy.wcs import WCS +from mocpy import MOC class DECamtile: @@ -210,7 +211,7 @@ def get_decam_ccds(ra, dec, save_footprint=False): return np.transpose(offsets, (2, 0, 1)) -def get_decam_quadrant_ipix(nside, ra, dec): +def get_decam_quadrant_moc(ra, dec): ccd_coords = get_decam_ccds(0, 0) skyoffset_frames = SkyCoord(ra, dec, unit=u.deg).skyoffset_frame() @@ -219,13 +220,15 @@ def get_decam_quadrant_ipix(nside, ra, dec): unit=u.deg, frame=skyoffset_frames[:, np.newaxis, np.newaxis], ).transform_to(ICRS) - ccd_xyz = np.moveaxis(ccd_coords_icrs.cartesian.xyz.value, 0, -1)[0] - ipixs = [] - for xyz in ccd_xyz: - ipix = hp.query_polygon(nside, xyz) - ipixs.append(ipix.tolist()) - return ipixs + moc = None + for subfield_id, coords in enumerate(ccd_coords_icrs[0]): + if moc is None: + moc = MOC.from_polygon_skycoord(coords) + else: + moc = moc + MOC.from_polygon_skycoord(coords) + + return moc def ccd_xy_to_radec(alpha_p, delta_p, ccd_centers_xy): diff --git a/gwemopt/chipgaps/ztf.py b/gwemopt/chipgaps/ztf.py index 047413e8..8f3bd6e9 100644 --- a/gwemopt/chipgaps/ztf.py +++ b/gwemopt/chipgaps/ztf.py @@ -25,6 +25,7 @@ from astropy import units as u from astropy.coordinates import ICRS, SkyCoord from astropy.wcs import WCS +from mocpy import MOC class ZTFtile: @@ -273,7 +274,7 @@ def get_ztf_quadrants(): return np.transpose(offsets, (2, 0, 1)) -def get_ztf_quadrant_ipix(nside, ra, dec, subfield_ids=None): +def get_ztf_quadrant_moc(ra, dec, subfield_ids=None): quadrant_coords = get_ztf_quadrants() skyoffset_frames = SkyCoord(ra, dec, unit=u.deg).skyoffset_frame() @@ -282,13 +283,16 @@ def get_ztf_quadrant_ipix(nside, ra, dec, subfield_ids=None): unit=u.deg, frame=skyoffset_frames[:, np.newaxis, np.newaxis], ).transform_to(ICRS) - quadrant_xyz = np.moveaxis(quadrant_coords_icrs.cartesian.xyz.value, 0, -1)[0] - ipixs = [] - for subfield_id, xyz in enumerate(quadrant_xyz): + moc = None + for subfield_id, coords in enumerate(quadrant_coords_icrs[0]): if subfield_ids is not None: if subfield_id not in subfield_ids: continue - ipix = hp.query_polygon(nside, xyz) - ipixs.append(ipix.tolist()) - return ipixs + + if moc is None: + moc = MOC.from_polygon_skycoord(coords) + else: + moc = moc + MOC.from_polygon_skycoord(coords) + + return moc diff --git a/gwemopt/coverage.py b/gwemopt/coverage.py index d3fa0b59..b533fa7f 100644 --- a/gwemopt/coverage.py +++ b/gwemopt/coverage.py @@ -30,10 +30,7 @@ def combine_coverage_structs(coverage_structs): coverage_struct_combined = {} coverage_struct_combined["data"] = np.empty((0, 8)) coverage_struct_combined["filters"] = np.empty((0, 1)) - coverage_struct_combined["ipix"] = [] - coverage_struct_combined["patch"] = [] - coverage_struct_combined["FOV"] = np.empty((0, 1)) - coverage_struct_combined["area"] = np.empty((0, 1)) + coverage_struct_combined["moc"] = [] coverage_struct_combined["telescope"] = np.empty((0, 1)) coverage_struct_combined["galaxies"] = [] coverage_struct_combined["exposureused"] = [] @@ -45,17 +42,8 @@ def combine_coverage_structs(coverage_structs): coverage_struct_combined["filters"] = np.append( coverage_struct_combined["filters"], coverage_struct["filters"] ) - coverage_struct_combined["ipix"] = ( - coverage_struct_combined["ipix"] + coverage_struct["ipix"] - ) - coverage_struct_combined["patch"] = ( - coverage_struct_combined["patch"] + coverage_struct["patch"] - ) - coverage_struct_combined["FOV"] = np.append( - coverage_struct_combined["FOV"], coverage_struct["FOV"] - ) - coverage_struct_combined["area"] = np.append( - coverage_struct_combined["area"], coverage_struct["area"] + coverage_struct_combined["moc"] = ( + coverage_struct_combined["moc"] + coverage_struct["moc"] ) coverage_struct_combined["telescope"] = np.append( coverage_struct_combined["telescope"], coverage_struct["telescope"] @@ -79,9 +67,7 @@ def read_coverage(params, telescope, filename, moc_struct=None): coverage_struct = {} coverage_struct["data"] = np.empty((0, 8)) coverage_struct["filters"] = [] - coverage_struct["ipix"] = [] - coverage_struct["patch"] = [] - coverage_struct["area"] = [] + coverage_struct["moc"] = [] for ii, row in schedule_table.iterrows(): ra, dec = row["ra"], row["dec"] @@ -112,34 +98,27 @@ def read_coverage(params, telescope, filename, moc_struct=None): color = "#6c71c4" if config_struct["FOV_coverage_type"] == "square": - ipix, radecs, patch, area = getSquarePixels( + moc = getSquarePixels( ra, dec, config_struct["FOV_coverage"], - nside, alpha=alpha, color=color, ) elif config_struct["FOV_coverage_type"] == "circle": - ipix, radecs, patch, area = getCirclePixels( + moc = getCirclePixels( ra, dec, config_struct["FOV_coverage"], - nside, alpha=alpha, color=color, ) else: - ipix = moc_struct[field]["ipix"] - patch = moc_struct[field]["patch"] - area = moc_struct[field]["area"] + moc = moc_struct[field]["moc"] - coverage_struct["patch"].append(patch) - coverage_struct["ipix"].append(ipix) - coverage_struct["area"].append(area) + coverage_struct["moc"].append(moc) coverage_struct["filters"] = np.array(coverage_struct["filters"]) - coverage_struct["area"] = np.array(coverage_struct["area"]) coverage_struct["FOV"] = config_struct["FOV_coverage"] * np.ones( (len(coverage_struct["filters"]),) ) @@ -173,7 +152,7 @@ def powerlaw(params, map_struct, tile_structs, previous_coverage_struct=None): coverage_structs = [] n_scope = 0 - full_prob_map = map_struct["prob"] + full_prob_map = map_struct["skymap"] for jj, telescope in enumerate(params["telescopes"]): if params["splitType"] is not None: @@ -233,16 +212,7 @@ def powerlaw(params, map_struct, tile_structs, previous_coverage_struct=None): previous_coverage_struct, ) - if params["doRASlices"]: - coverage_struct = gwemopt.scheduler.schedule_ra_splits( - params, - config_struct, - map_struct_hold, - tile_struct, - telescope, - previous_coverage_struct, - ) - elif params["doBalanceExposure"]: + if params["doBalanceExposure"]: if ( params["max_nb_tiles"] is None ): # optimize max tiles (iff max tiles not already specified) @@ -393,11 +363,6 @@ def powerlaw(params, map_struct, tile_structs, previous_coverage_struct=None): "need to specify tiles that have been observed using --observedTiles" ) - if params["treasuremap_token"] is not None and previous_coverage_struct: - tile_struct = update_observed_tiles( - params, tile_struct, previous_coverage_struct - ) # coverage_struct of the previous round - coverage_struct = gwemopt.scheduler.scheduler( params, config_struct, tile_struct ) @@ -435,327 +400,16 @@ def powerlaw(params, map_struct, tile_structs, previous_coverage_struct=None): if params["doIterativeTiling"]: map_struct_hold = slice_map_tiles(params, map_struct_hold, coverage_struct) - map_struct["prob"] = full_prob_map - - return tile_structs, combine_coverage_structs(coverage_structs) - - -def absmag(params, map_struct, tile_structs, previous_coverage_struct=None): - if "distmu" not in map_struct: - print("timeallocationType absmag requires 3d geometry") - exit(0) - - map_struct_hold = copy.deepcopy(map_struct) - - coverage_structs = [] - n_scope = 0 - full_prob_map = map_struct["prob"] - - for jj, telescope in enumerate(params["telescopes"]): - if params["splitType"] is not None: - if "observability" in map_struct: - map_struct["observability"][telescope]["prob"] = map_struct["groups"][ - n_scope - ] - else: - map_struct["prob"] = map_struct["groups"][n_scope] - if n_scope < len(map_struct["groups"]) - 1: - n_scope += 1 - else: - n_scope = 0 - - config_struct = params["config"][telescope] - tile_struct = tile_structs[telescope] - - # Try to load the minimum duration of time from telescope config file - # Otherwise set it to zero - try: - min_obs_duration = config_struct["min_observability_duration"] / 24 - except: - min_obs_duration = 0.0 - - if params["doIterativeTiling"] and (params["tilesType"] == "galaxy"): - tile_struct = slice_galaxy_tiles( - params, tile_struct, combine_coverage_structs(coverage_structs) - ) - - if ( - params["doPerturbativeTiling"] - and (jj > 0) - and (not params["tilesType"] == "galaxy") - ): - tile_struct = perturb_tiles( - params, config_struct, telescope, map_struct_hold, tile_struct - ) - - if params["doOverlappingScheduling"]: - tile_struct = check_overlapping_tiles( - params, tile_struct, combine_coverage_structs(coverage_structs) - ) - - if params["doAlternatingFilters"]: - if params["doBlocks"]: - tile_struct = eject_tiles(params, telescope, tile_struct) - - params_hold = copy.copy(params) - config_struct_hold = copy.copy(config_struct) - - coverage_struct, tile_struct = gwemopt.scheduler.schedule_alternating( - params_hold, - config_struct_hold, - telescope, - map_struct_hold, - tile_struct, - previous_coverage_struct, - ) - if params["doRASlices"]: - coverage_struct = gwemopt.scheduler.schedule_ra_splits( - params, - config_struct, - map_struct_hold, - tile_struct, - telescope, - previous_coverage_struct, - ) - elif params["doBalanceExposure"]: - if params["max_nb_tiles"] is None: - # optimize max tiles (iff max tiles not already specified) - optimized_max, coverage_struct, tile_struct = optimize_max_tiles( - params, - tile_struct, - coverage_struct, - config_struct, - telescope, - map_struct_hold, - ) - params["max_nb_tiles"] = np.array([optimized_max], dtype=float) - - else: - params_hold = copy.copy(params) - config_struct_hold = copy.copy(config_struct) - - ( - coverage_struct, - tile_struct, - ) = gwemopt.scheduler.schedule_alternating( - params_hold, - config_struct_hold, - telescope, - map_struct_hold, - tile_struct, - previous_coverage_struct, - ) - do_reschedule, balanced_fields = balance_tiles( - params_hold, tile_struct, coverage_struct - ) - - if do_reschedule: - config_struct_hold = copy.copy(config_struct) - ( - coverage_struct, - tile_struct, - ) = gwemopt.scheduler.schedule_alternating( - params_hold, - config_struct_hold, - telescope, - map_struct_hold, - tile_struct, - previous_coverage_struct, - ) - - # coverage_struct = gwemopt.utils.erase_unbalanced_tiles(params_hold,coverage_struct) - else: - # load the sun retriction for a satelite - try: - sat_sun_restriction = config_struct["sat_sun_restriction"] - except: - sat_sun_restriction = 0.0 - - if not params["tilesType"] == "galaxy": - tile_struct = gwemopt.tiles.absmag_tiles_struct( - params, config_struct, telescope, map_struct_hold, tile_struct - ) - - elif sat_sun_restriction == 0.0: - for key in tile_struct.keys(): - # Check that a given tile is observable a minimum amount of time - # If not set the proba associated to the tile to zero - if ( - "segmentlist" - and "prob" in tile_struct[key] - and tile_struct[key]["segmentlist"] - and min_obs_duration > 0.0 - ): - observability_duration = 0.0 - for counter in range(len(tile_struct[key]["segmentlist"])): - observability_duration += ( - tile_struct[key]["segmentlist"][counter][1] - - tile_struct[key]["segmentlist"][counter][0] - ) - if ( - tile_struct[key]["prob"] > 0.0 - and observability_duration < min_obs_duration - ): - tile_struct[key]["prob"] = 0.0 - - else: - # Check that a given tile is not to close to the sun for the satelite - # If it's to close set the proba associated to the tile to zero - - time = map_struct["trigtime"] - time = Time(time, format="isot", scale="utc") - sun_position = get_sun(time) - - # astropy don't like the output of get sun in the following separator function, need here to redefine the skycoord - sun_position = SkyCoord(sun_position.ra, sun_position.dec, frame="gcrs") - - for key in tile_struct.keys(): - if ( - "segmentlist" - and "prob" in tile_struct[key] - and tile_struct[key]["segmentlist"] - ): - for counter in range(len(tile_struct[key]["segmentlist"])): - tile_position = SkyCoord( - ra=tile_struct[key]["ra"] * u.degree, - dec=tile_struct[key]["dec"] * u.degree, - frame="icrs", - ) - ang_dist = sun_position.separation(tile_position).deg - - if ang_dist < sat_sun_restriction: - tile_struct[key]["prob"] = 0.0 - - if ( - params["timeallocationType"] == "manual" - ): # only works if using same telescope - try: - for field_id in params["observedTiles"]: - field_id = int(field_id) - if field_id in tile_struct: - tile_struct[field_id]["prob"] = 0.0 - - except: - raise ValueError( - "need to specify tiles that have been observed using --observedTiles" - ) - - if params["treasuremap_token"] is not None and previous_coverage_struct: - tile_struct = update_observed_tiles( - params, tile_struct, previous_coverage_struct - ) # coverage_struct of the previous round - - coverage_struct = gwemopt.scheduler.scheduler( - params, config_struct, tile_struct - ) - - if params["doBalanceExposure"]: - ntrials = 10 - for _ in tqdm(range(ntrials)): - do_reschedule, balanced_fields = balance_tiles( - params, tile_struct, coverage_struct - ) - if do_reschedule: - for key in params["unbalanced_tiles"]: - tile_struct[key]["prob"] = 0.0 - coverage_struct = gwemopt.scheduler.scheduler( - params, config_struct, tile_struct - ) - - if params["max_nb_tiles"] is not None: - tile_struct, do_reschedule = slice_number_tiles( - params, telescope, tile_struct, coverage_struct - ) - if do_reschedule: - coverage_struct = gwemopt.scheduler.scheduler( - params, config_struct, tile_struct - ) - - tile_structs[telescope] = tile_struct - coverage_structs.append(coverage_struct) - - if params["doIterativeTiling"]: - map_struct_hold = slice_map_tiles(params, map_struct_hold, coverage_struct) - - map_struct["prob"] = full_prob_map + map_struct["skymap"] = full_prob_map return tile_structs, combine_coverage_structs(coverage_structs) -def pem(params, map_struct, tile_structs): - map_struct_hold = copy.deepcopy(map_struct) - - coverage_structs = [] - for telescope in params["telescopes"]: - config_struct = params["config"][telescope] - tile_struct = tile_structs[telescope] - - if params["doAlternatingFilters"]: - filters, exposuretimes = params["filters"], params["exposuretimes"] - tile_struct_hold = copy.copy(tile_struct) - coverage_structs_hold = [] - for filt, exposuretime in zip(filters, exposuretimes): - params["filters"] = [filt] - params["exposuretimes"] = [exposuretime] - tile_struct_hold = gwemopt.tiles.pem_tiles_struct( - params, config_struct, telescope, map_struct_hold, tile_struct_hold - ) - coverage_struct_hold = gwemopt.scheduler.scheduler( - params, config_struct, tile_struct_hold - ) - coverage_structs_hold.append(coverage_struct_hold) - coverage_struct = combine_coverage_structs(coverage_structs_hold) - else: - tile_struct = gwemopt.tiles.pem_tiles_struct( - params, config_struct, telescope, map_struct_hold, tile_struct - ) - coverage_struct = gwemopt.scheduler.scheduler( - params, config_struct, tile_struct - ) - coverage_structs.append(coverage_struct) - - if params["doIterativeTiling"]: - map_struct_hold = slice_map_tiles(map_struct_hold, coverage_struct) - - return combine_coverage_structs(coverage_structs) - - -def update_observed_tiles(params, tile_struct, previous_coverage_struct): - if not params["doAlternatingFilters"]: - tile_struct = check_overlapping_tiles( - params, tile_struct, previous_coverage_struct - ) # maps field ids to tile_struct - - for key in tile_struct.keys(): # sets tile to 0 if previously observed - if "epochs" not in tile_struct[key]: - continue - ipix = tile_struct[key]["ipix"] - - tot_overlap = sum( - tile_struct[key]["epochs_overlap"] - ) # sums over list of overlapping ipix lengths - - if params["doAlternatingFilters"]: - # only takes into account fields with same filters for total overlap - for ii, filt in enumerate(tile_struct[key]["epochs_filters"]): - if filt != params["filters"][0]: - tot_overlap -= tile_struct[key]["epochs_overlap"][ii] - - rat = tot_overlap / len(ipix) - - if rat > 0.3: - tile_struct[key]["prob"] = 0.0 - - return tile_struct - - def timeallocation(params, map_struct, tile_structs, previous_coverage_struct=None): if len(params["telescopes"]) > 1 and params["doOrderByObservability"]: order_by_observability(params, tile_structs) - if (params["timeallocationType"] == "powerlaw") or ( - params["timeallocationType"] == "absmag" - ): + if params["timeallocationType"] == "powerlaw": print("Generating powerlaw schedule...") if params["treasuremap_token"] is not None: @@ -802,44 +456,23 @@ def timeallocation(params, map_struct, tile_structs, previous_coverage_struct=No exposurelist.append(segments.segment(seg[0], seg[1])) params_hold["config"][telescope]["exposurelist"] = exposurelist - if params["timeallocationType"] == "absmag": - tile_structs_hold[telescope] = ( - gwemopt.tiles.absmag_tiles_struct( - params_hold, - config_struct, - telescope, - map_struct, - tile_structs_hold[telescope], - ) - ) - elif params["timeallocationType"] == "powerlaw": - tile_structs_hold[telescope] = ( - gwemopt.tiles.powerlaw_tiles_struct( - params_hold, - config_struct, - telescope, - map_struct, - tile_structs_hold[telescope], - ) - ) - - if params["timeallocationType"] == "absmag": - tile_structs_hold, coverage_struct = gwemopt.coverage.absmag( - params_hold, - map_struct, - tile_structs_hold, - previous_coverage_struct, - ) - elif params["timeallocationType"] == "powerlaw": - tile_structs_hold, coverage_struct = gwemopt.coverage.powerlaw( + tile_structs_hold[telescope] = gwemopt.tiles.powerlaw_tiles_struct( params_hold, + config_struct, + telescope, map_struct, - tile_structs_hold, - previous_coverage_struct, + tile_structs_hold[telescope], ) + tile_structs_hold, coverage_struct = gwemopt.coverage.powerlaw( + params_hold, + map_struct, + tile_structs_hold, + previous_coverage_struct, + ) + coverage_structs.append(coverage_struct) - for ii in range(len(coverage_struct["ipix"])): + for ii in range(len(coverage_struct["moc"])): telescope = coverage_struct["telescope"][ii] scheduled_fields[telescope].append( coverage_struct["data"][ii, 5] @@ -847,14 +480,9 @@ def timeallocation(params, map_struct, tile_structs, previous_coverage_struct=No coverage_struct = combine_coverage_structs(coverage_structs) else: - if params["timeallocationType"] == "absmag": - tile_structs, coverage_struct = gwemopt.coverage.absmag( - params, map_struct, tile_structs, previous_coverage_struct - ) - elif params["timeallocationType"] == "powerlaw": - tile_structs, coverage_struct = gwemopt.coverage.powerlaw( - params, map_struct, tile_structs, previous_coverage_struct - ) + tile_structs, coverage_struct = gwemopt.coverage.powerlaw( + params, map_struct, tile_structs, previous_coverage_struct + ) elif params["timeallocationType"] == "manual": print("Generating manual schedule...") @@ -862,8 +490,4 @@ def timeallocation(params, map_struct, tile_structs, previous_coverage_struct=No params, map_struct, tile_structs ) - elif params["timeallocationType"] == "pem": - print("Generating PEM schedule...") - coverage_struct = gwemopt.coverage.pem(params, map_struct, tile_structs) - return tile_structs, coverage_struct diff --git a/gwemopt/efficiency.py b/gwemopt/efficiency.py index 1ead4010..5ab70403 100644 --- a/gwemopt/efficiency.py +++ b/gwemopt/efficiency.py @@ -1,11 +1,13 @@ import os +import astropy.units as u import healpy as hp import numpy as np import scipy.stats +from astropy.coordinates import SkyCoord from astropy.time import Time from ligo.skymap import distance -from scipy.interpolate import interpolate as interp +from scipy.interpolate import interp1d from gwemopt.io.export_efficiency import ( export_efficiency_data, @@ -23,24 +25,22 @@ def compute_efficiency(params, map_struct, lightcurve_struct, coverage_struct): gpstime = params["gpstime"] mjd_inj = Time(gpstime, format="gps", scale="utc").mjd - if params["catalog"] is not None: - distn = scipy.stats.rv_discrete( - values=(np.arange(npix), map_struct["prob_catalog"]) - ) - else: - distn = scipy.stats.rv_discrete(values=(np.arange(npix), map_struct["prob"])) + distn = scipy.stats.rv_discrete( + values=(np.arange(npix), map_struct["skymap_raster"]["PROB"]) + ) ipix = distn.rvs(size=Ninj) ras, decs = hp.pix2ang(nside, ipix, lonlat=True) dists = np.logspace(-1, 3, 1000) ndetections = np.zeros((len(dists),)) - for pinpoint in ipix: + for ra, dec in zip(ras, decs): + coords = SkyCoord(ra * u.deg, dec * u.deg) idxs = [] - for jj in range(len(coverage_struct["ipix"])): - expPixels = coverage_struct["ipix"][jj] + for jj in range(len(coverage_struct["moc"])): + moc = coverage_struct["moc"][jj] - if pinpoint in expPixels: + if moc.contains_skycoords(coords): idxs.append(jj) if len(idxs) == 0: continue @@ -56,7 +56,7 @@ def compute_efficiency(params, map_struct, lightcurve_struct, coverage_struct): lightcurve_mag = lightcurve_struct[filt] idx = np.where(np.isfinite(lightcurve_mag))[0] - f = interp.interp1d( + f = interp1d( lightcurve_t[idx], lightcurve_mag[idx], fill_value="extrapolate" ) lightcurve_mag_interp = f(mjd) @@ -113,16 +113,13 @@ def compute_true_efficiency(params, map_struct, lightcurve_struct, coverage_stru gpstime = params["gpstime"] mjd_inj = Time(gpstime, format="gps", scale="utc").mjd - ipix = hp.ang2pix(nside, params["true_ra"], params["true_dec"], lonlat=True) - # ipix = np.argmax(map_struct["prob"]) - # lat, lon = hp.pix2ang(nside, ipix, lonlat=True) dist = params["true_distance"] - + coords = SkyCoord(params["true_ra"] * u.deg, params["true_dec"] * u.deg) idxs = [] - for jj in range(len(coverage_struct["ipix"])): - expPixels = coverage_struct["ipix"][jj] + for jj in range(len(coverage_struct["moc"])): + moc = coverage_struct["moc"][jj] - if ipix in expPixels: + if moc.contains_skycoords(coords): idxs.append(jj) mjds = coverage_struct["data"][idxs, 2] @@ -154,38 +151,39 @@ def compute_3d_efficiency(params, map_struct, lightcurve_struct, coverage_struct gpstime = params["gpstime"] mjd_inj = Time(gpstime, format="gps", scale="utc").mjd - if params["catalog"] is not None: - distn = scipy.stats.rv_discrete( - values=(np.arange(npix), map_struct["prob_catalog"]) - ) - else: - distn = scipy.stats.rv_discrete(values=(np.arange(npix), map_struct["prob"])) - ipix = distn.rvs(size=Ninj) - - # calculate the moments from distmu, distsigma and distnorm - mom_mean, mom_std, mom_norm = distance.parameters_to_moments( - map_struct["distmu"], map_struct["distsigma"] + distn = scipy.stats.rv_discrete( + values=(np.arange(npix), map_struct["skymap_raster"]["PROB"]) ) + ipix = distn.rvs(size=Ninj) + ras, decs = hp.pix2ang(nside, ipix, lonlat=True) detections = 0 dists_inj = {} dists_inj["recovered"], dists_inj["tot"] = [], [] - for pinpoint in ipix: + for pinpoint, ra, dec in zip(ipix, ras, decs): + coords = SkyCoord(ra * u.deg, dec * u.deg) + dist = -1 while dist < 0: - if np.isinf(mom_mean[pinpoint]) or np.isinf(mom_std[pinpoint]): + if np.isinf(map_struct["skymap_raster"]["DISTMEAN"][pinpoint]) or np.isinf( + map_struct["skymap_raster"]["DISTSTD"][pinpoint] + ): dist = np.inf else: - dist = mom_mean[pinpoint] + mom_std[pinpoint] * np.random.normal() + dist = ( + map_struct["skymap_raster"]["DISTMEAN"][pinpoint] + + map_struct["skymap_raster"]["DISTSTD"][pinpoint] + * np.random.normal() + ) if dist != np.inf: dists_inj["tot"].append(dist) idxs = [] - for jj in range(len(coverage_struct["ipix"])): - expPixels = coverage_struct["ipix"][jj] + for jj in range(len(coverage_struct["moc"])): + moc = coverage_struct["moc"][jj] - if pinpoint in expPixels: + if moc.contains_skycoords(coords): idxs.append(jj) if len(idxs) == 0: @@ -200,7 +198,7 @@ def compute_3d_efficiency(params, map_struct, lightcurve_struct, coverage_struct lightcurve_mag = lightcurve_struct[filt] idx = np.where(np.isfinite(lightcurve_mag))[0] - f = interp.interp1d( + f = interp1d( lightcurve_t[idx], lightcurve_mag[idx], fill_value="extrapolate" ) lightcurve_mag_interp = f(mjd) diff --git a/gwemopt/io/__init__.py b/gwemopt/io/__init__.py index 6a59ac8e..758e542d 100644 --- a/gwemopt/io/__init__.py +++ b/gwemopt/io/__init__.py @@ -1,6 +1,3 @@ -from gwemopt.io.export_tiles import ( - export_tiles_coverage_hist, - export_tiles_coverage_int, -) +from gwemopt.io.export_tiles import export_tiles_coverage_int from gwemopt.io.schedule import export_schedule_xml, read_schedule, summary from gwemopt.io.skymap import get_skymap, read_skymap diff --git a/gwemopt/io/export_tiles.py b/gwemopt/io/export_tiles.py index e6926c03..50112696 100644 --- a/gwemopt/io/export_tiles.py +++ b/gwemopt/io/export_tiles.py @@ -1,13 +1,3 @@ -def export_tiles_coverage_hist(params, diffs): - """ - Export the histogram of the tiles coverage. - """ - filename = params["outputDir"].joinpath("tiles_coverage_hist.dat") - with open(filename, "w") as fid: - for ii in range(len(diffs)): - fid.write("%.10f\n" % diffs[ii]) - - def export_tiles_coverage_int(filename, tts, cum_probs, cum_areas): with open(filename, "w") as fid: for ii in range(len(tts)): diff --git a/gwemopt/io/schedule.py b/gwemopt/io/schedule.py index a6745d75..3f5365ea 100644 --- a/gwemopt/io/schedule.py +++ b/gwemopt/io/schedule.py @@ -151,13 +151,10 @@ def export_schedule_xml(xmlfile, map_struct, coverage_struct, config_struct): Field(name="priority", ucd="", unit="", dataType="int", Description=[""]) ) table_field = utilityTable(table) - table_field.blankTable(len(coverage_struct["ipix"])) + table_field.blankTable(len(coverage_struct["moc"])) - for ii in range(len(coverage_struct["ipix"])): + for ii in range(len(coverage_struct["moc"])): data = coverage_struct["data"][ii, :] - ipix = coverage_struct["ipix"][ii] - - prob = np.sum(map_struct["prob"][ipix]) ra, dec = data[0], data[1] observ_time, exposure_time, field_id, prob, airmass = ( @@ -198,15 +195,6 @@ def export_schedule_xml(xmlfile, map_struct, coverage_struct, config_struct): def summary(params, map_struct, coverage_struct, catalog_struct=None): - idx50 = len(map_struct["cumprob"]) - np.argmin(np.abs(map_struct["cumprob"] - 0.50)) - idx90 = len(map_struct["cumprob"]) - np.argmin(np.abs(map_struct["cumprob"] - 0.90)) - - mapfile = params["outputDir"].joinpath("map.dat") - with open(mapfile, "w") as fid: - fid.write( - "%.5f %.5f\n" - % (map_struct["pixarea_deg2"] * idx50, map_struct["pixarea_deg2"] * idx90) - ) filts = list(set(coverage_struct["filters"])) for jj, telescope in enumerate(params["telescopes"]): @@ -225,15 +213,12 @@ def summary(params, map_struct, coverage_struct, catalog_struct=None): totexp = 0 with open(schedulefile, "w") as fid: - for ii in range(len(coverage_struct["ipix"])): + for ii in range(len(coverage_struct["filters"])): if not telescope == coverage_struct["telescope"][ii]: continue data = coverage_struct["data"][ii, :] filt = coverage_struct["filters"][ii] - ipix = coverage_struct["ipix"][ii] - - prob = np.sum(map_struct["prob"][ipix]) ra, dec = data[0], data[1] observ_time, mag, exposure_time, field_id, prob, airmass = ( @@ -285,15 +270,7 @@ def summary(params, map_struct, coverage_struct, catalog_struct=None): fields_sum = np.sum(fields[:, 2:], axis=1) idx = np.where(fields_sum >= 2)[0] - print( - "%d/%d fields were observed at least twice\n" % (len(idx), len(fields_sum)) - ) - - print( - "Integrated probability, All: %.5f, 2+: %.5f" - % (np.sum(fields[:, 1]), np.sum(fields[idx, 1])) - ) - + print("%d/%d fields were observed at least twice" % (len(idx), len(fields_sum))) print(f"Expected time spent on exposures: {totexp / 3600:.1f} hr.") slew_readout_time = computeSlewReadoutTime(config_struct, coverage_struct) print(f"Expected time spent on slewing and readout: {slew_readout_time:.0f} s.") @@ -307,6 +284,8 @@ def summary(params, map_struct, coverage_struct, catalog_struct=None): fid.write("\n") summaryfile = params["outputDir"].joinpath("summary.dat") + cummoc = None + with open(summaryfile, "w") as fid: gpstime = params["gpstime"] event_mjd = Time(gpstime, format="gps", scale="utc").mjd @@ -322,26 +301,36 @@ def summary(params, map_struct, coverage_struct, catalog_struct=None): if params["tilesType"] == "galaxy": galaxies = np.empty((0, 2)) - for ii in range(len(coverage_struct["ipix"])): + for ii in range(len(coverage_struct["filters"])): data = coverage_struct["data"][ii, :] filt = coverage_struct["filters"][ii] - ipix = coverage_struct["ipix"][ii] + moc = coverage_struct["moc"][ii] - prob = np.sum(map_struct["prob"][ipix]) + ra, dec = data[0], data[1] + observ_time, mag, exposure_time, field_id, prob, airmass = ( + data[2], + data[3], + data[4], + data[5], + data[6], + data[7], + ) if data[2] > event_mjd + tt: continue - ipixs = np.append(ipixs, ipix) - ipixs = np.unique(ipixs).astype(int) - cum_prob = np.sum(map_struct["prob"][ipixs]) + if cummoc is None: + cummoc = moc + else: + cummoc = cummoc + moc + cum_prob = cummoc.probability_in_multiordermap(map_struct["skymap"]) if params["tilesType"] == "galaxy": galaxies = np.append(galaxies, coverage_struct["galaxies"][ii]) galaxies = np.unique(galaxies).astype(int) cum_prob = np.sum(catalog_struct[params["galaxy_grade"]][galaxies]) - cum_area = len(ipixs) * map_struct["pixarea_deg2"] + cum_area = cummoc.sky_fraction * 360**2 / np.pi mjds.append(data[2]) mjds_floor.append(int(np.floor(data[2]))) diff --git a/gwemopt/io/skymap.py b/gwemopt/io/skymap.py index 9bf59e98..3058a583 100644 --- a/gwemopt/io/skymap.py +++ b/gwemopt/io/skymap.py @@ -5,8 +5,10 @@ import os from pathlib import Path +import astropy_healpix as ah import healpy as hp import ligo.skymap.distance as ligodist +import ligo.skymap.plot import lxml.etree import numpy as np import requests @@ -15,9 +17,11 @@ from astropy.io import fits from astropy.time import Time from ligo.gracedb.rest import GraceDb +from ligo.skymap import moc from ligo.skymap.bayestar import rasterize from ligo.skymap.io import read_sky_map -from ligo.skymap.moc import uniq2nest +from matplotlib import pyplot as plt +from mocpy import MOC from scipy.interpolate import PchipInterpolator from scipy.stats import norm @@ -137,14 +141,14 @@ def read_inclination(skymap, params, map_struct): theta = 0.5 * np.pi - np.deg2rad(dec) # make use of the maP nside maP_idx = np.argmax(skymap["PROBDENSITY"]) - order, _ = uniq2nest(skymap[maP_idx]["UNIQ"]) + order, _ = moc.uniq2nest(skymap[maP_idx]["UNIQ"]) nside = hp.order2nside(order) # get the nested idx for the given sky location nest_idx = hp.ang2pix(nside, theta, phi, nest=True) # find the row with the closest nested index nest_idxs = [] for row in skymap: - order_per_row, nest_idx_per_row = uniq2nest(row["UNIQ"]) + order_per_row, nest_idx_per_row = moc.uniq2nest(row["UNIQ"]) if order_per_row == order: nest_idxs.append(nest_idx_per_row) else: @@ -156,7 +160,7 @@ def read_inclination(skymap, params, map_struct): maP_idx = np.argmax(skymap["PROBDENSITY"]) uniq_idx = skymap[maP_idx]["UNIQ"] # convert to nested indexing and find the location of that index - order, nest_idx = uniq2nest(uniq_idx) + order, nest_idx = moc.uniq2nest(uniq_idx) nside = hp.order2nside(order) theta, phi = hp.pix2ang(nside, int(nest_idx), nest=True) # convert theta and phi to ra and dec @@ -258,151 +262,108 @@ def read_skymap(params, map_struct=None): if "do_3d" not in params: params["do_3d"] = is_3d - header = [] - if map_struct is None: - map_struct = {} - - filename = params["skymap"] - - if params["do_3d"]: - try: - healpix_data, header = hp.read_map( - filename, field=(0, 1, 2, 3), h=True - ) - except: - skymap = read_sky_map(filename, moc=True, distances=True) - order = hp.nside2order(params["nside"]) - - # for colname in skymap.colnames: - # if colname.startswith('PROB'): - # newname = colname.replace('PROB', 'PROBDENSITY') - # skymap.rename_column(colname, newname) - # skymap[newname] *= len(skymap) / (4 * np.pi) - # skymap[newname].unit = u.steradian ** -1 - - if "PROBDENSITY_SAMPLES" in skymap.columns: - if params["inclination"]: - map_struct = read_inclination(skymap, params, map_struct) - - skymap.remove_columns( - [ - f"{name}_SAMPLES" - for name in [ - "PROBDENSITY", - "DISTMU", - "DISTSIGMA", - "DISTNORM", - ] - ] - ) - - t = rasterize(skymap) - result = t["PROB"], t["DISTMU"], t["DISTSIGMA"], t["DISTNORM"] - healpix_data = hp.reorder(result, "NESTED", "RING") - - distmu_data = healpix_data[1] - distsigma_data = healpix_data[2] - prob_data = healpix_data[0] - norm_data = healpix_data[3] - - map_struct["distmu"] = distmu_data / params["DScale"] - map_struct["distsigma"] = distsigma_data / params["DScale"] - map_struct["prob"] = prob_data - map_struct["distnorm"] = norm_data + map_struct = {} - else: - prob_data, header = hp.read_map( - filename, field=0, verbose=False, h=True - ) - prob_data = prob_data / np.sum(prob_data) - - map_struct["prob"] = prob_data - - for j in range(len(header)): - if header[j][0] == "DATE": - map_struct["trigtime"] = header[j][1] - - natural_nside = hp.pixelfunc.get_nside(map_struct["prob"]) - nside = params["nside"] - - print("natural_nside =", natural_nside) - print("nside =", nside) - - if not params["do_3d"]: - map_struct["prob"] = hp.ud_grade(map_struct["prob"], nside, power=-2) - - if params["do_3d"]: - if natural_nside != nside: - map_struct["prob"] = hp.pixelfunc.ud_grade( - map_struct["prob"], nside, power=-2 - ) - map_struct["distmu"] = hp.pixelfunc.ud_grade(map_struct["distmu"], nside) - map_struct["distsigma"] = hp.pixelfunc.ud_grade( - map_struct["distsigma"], nside - ) - map_struct["distnorm"] = hp.pixelfunc.ud_grade( - map_struct["distnorm"], nside - ) - - map_struct["distmu"][map_struct["distmu"] < -1e30] = np.inf - - nside_down = 32 - - distmu_down = hp.pixelfunc.ud_grade(map_struct["distmu"], nside_down) - - ( - map_struct["distmed"], - map_struct["diststd"], - mom_norm, - ) = ligodist.parameters_to_moments( - map_struct["distmu"], map_struct["distsigma"] - ) - - distmu_down[distmu_down < -1e30] = np.inf + filename = params["skymap"] - map_struct["distmed"] = hp.ud_grade(map_struct["distmed"], nside, power=-2) - map_struct["diststd"] = hp.ud_grade(map_struct["diststd"], nside, power=-2) + if params["do_3d"]: + skymap = read_sky_map(filename, moc=True, distances=True) - npix = hp.nside2npix(nside) - theta, phi = hp.pix2ang(nside, np.arange(npix)) - ra = np.rad2deg(phi) - dec = np.rad2deg(0.5 * np.pi - theta) + if "PROBDENSITY_SAMPLES" in skymap.columns: + if params["inclination"]: + map_struct = read_inclination(skymap, params, map_struct) - map_struct["ra"] = ra - map_struct["dec"] = dec + skymap.remove_columns( + [ + f"{name}_SAMPLES" + for name in [ + "PROBDENSITY", + "DISTMU", + "DISTSIGMA", + "DISTNORM", + ] + ] + ) - if params["doRASlice"]: + map_struct["skymap"] = skymap + else: + skymap = read_sky_map(filename, moc=True, distances=False) + map_struct["skymap"] = skymap + + level, ipix = ah.uniq_to_level_ipix(map_struct["skymap"]["UNIQ"]) + LEVEL = MOC.MAX_ORDER + shift = 2 * (LEVEL - level) + hpx = np.array(np.vstack([ipix << shift, (ipix + 1) << shift]), dtype=np.uint64).T + nside = ah.level_to_nside(level) + pixel_area = ah.nside_to_pixel_area(ah.level_to_nside(level)) + ra, dec = ah.healpix_to_lonlat(ipix, nside, order="nested") + map_struct["skymap"]["ra"] = ra.deg + map_struct["skymap"]["dec"] = dec.deg + + if params["raslice"] is not None: ra_low, ra_high = params["raslice"][0], params["raslice"][1] if ra_low <= ra_high: ipix = np.where( - (ra_high * 360.0 / 24.0 < ra) | (ra_low * 360.0 / 24.0 > ra) + (ra_high * 360.0 / 24.0 < ra.deg) | (ra_low * 360.0 / 24.0 > ra.deg) )[0] else: ipix = np.where( - (ra_high * 360.0 / 24.0 < ra) & (ra_low * 360.0 / 24.0 > ra) + (ra_high * 360.0 / 24.0 < ra.deg) & (ra_low * 360.0 / 24.0 > ra.deg) )[0] - map_struct["prob"][ipix] = 0.0 - map_struct["prob"] = map_struct["prob"] / np.sum(map_struct["prob"]) + map_struct["skymap"]["PROBDENSITY"][ipix] = 0.0 if params["galactic_limit"] > 0.0: - coords = SkyCoord(ra=ra * u.deg, dec=dec * u.deg) + coords = SkyCoord(ra=ra, dec=dec) ipix = np.where(np.abs(coords.galactic.b.deg) <= params["galactic_limit"])[0] - map_struct["prob"][ipix] = 0.0 - map_struct["prob"] = map_struct["prob"] / np.sum(map_struct["prob"]) - - sort_idx = np.argsort(map_struct["prob"])[::-1] - csm = np.empty(len(map_struct["prob"])) - csm[sort_idx] = np.cumsum(map_struct["prob"][sort_idx]) + map_struct["skymap"]["PROBDENSITY"][ipix] = 0.0 - map_struct["cumprob"] = csm - map_struct["ipix_keep"] = np.where(csm <= params["iterativeOverlap"])[0] + map_struct["skymap_raster"] = rasterize( + map_struct["skymap"], order=hp.nside2order(int(params["nside"])) + ) + peak = map_struct["skymap_raster"][ + map_struct["skymap_raster"]["PROB"] + == np.max(map_struct["skymap_raster"]["PROB"]) + ] + map_struct["center"] = SkyCoord(peak["ra"][0] * u.deg, peak["dec"][0] * u.deg) - pixarea = hp.nside2pixarea(nside) - pixarea_deg2 = hp.nside2pixarea(nside, degrees=True) + if "DISTMU" in map_struct["skymap_raster"].columns: + ( + map_struct["skymap_raster"]["DISTMEAN"], + map_struct["skymap_raster"]["DISTSTD"], + mom_norm, + ) = ligodist.parameters_to_moments( + map_struct["skymap_raster"]["DISTMU"], + map_struct["skymap_raster"]["DISTSIGMA"], + ) - map_struct["nside"] = nside - map_struct["npix"] = npix - map_struct["pixarea"] = pixarea - map_struct["pixarea_deg2"] = pixarea_deg2 + extra_header = [ + ("PIXTYPE", "HEALPIX", "HEALPIX pixelisation"), + ("ORDERING", "NESTED", "Pixel ordering scheme: RING, NESTED, or NUNIQ"), + ("COORDSYS", "C", "Ecliptic, Galactic or Celestial (equatorial)"), + ( + "MOCORDER", + moc.uniq2order(map_struct["skymap"]["UNIQ"].max()), + "MOC resolution (best order)", + ), + ("INDXSCHM", "EXPLICIT", "Indexing: IMPLICIT or EXPLICIT"), + ] + + hdu = fits.table_to_hdu(map_struct["skymap_raster"]) + hdu.header.extend(extra_header) + map_struct["hdu"] = hdu + + skymap_sorted = map_struct["skymap"].copy() + skymap_sorted.sort("PROBDENSITY") + skymap_sorted.reverse() + + level, ipix = ah.uniq_to_level_ipix(skymap_sorted["UNIQ"]) + pixel_area = ah.nside_to_pixel_area(ah.level_to_nside(level)) + + prob = pixel_area * skymap_sorted["PROBDENSITY"] + cumprob = np.cumsum(prob) + + i = cumprob.searchsorted(params["iterativeOverlap"]) + skymap_keep = skymap_sorted[:i] + map_struct["moc_keep"] = skymap_keep return params, map_struct diff --git a/gwemopt/moc.py b/gwemopt/moc.py index 796435f1..b75b03dc 100644 --- a/gwemopt/moc.py +++ b/gwemopt/moc.py @@ -8,7 +8,7 @@ from tqdm import tqdm import gwemopt.tiles -from gwemopt.chipgaps import get_decam_quadrant_ipix, get_ztf_quadrant_ipix +from gwemopt.chipgaps import get_decam_quadrant_moc, get_ztf_quadrant_moc from gwemopt.paths import CONFIG_DIR from gwemopt.utils.pixels import getCirclePixels, getRegionPixels, getSquarePixels @@ -16,36 +16,12 @@ def create_moc(params, map_struct=None): nside = params["nside"] - if params["doMinimalTiling"]: - prob = map_struct["prob"] - cl = params["powerlaw_cl"] - prob_scaled = copy.deepcopy(prob) - prob_sorted = np.sort(prob_scaled)[::-1] - prob_indexes = np.argsort(prob_scaled)[::-1] - prob_cumsum = np.cumsum(prob_sorted) - index = np.argmin(np.abs(prob_cumsum - cl)) + 1 - prob_indexes = prob_indexes[: index + 1] - moc_structs = {} for telescope in params["telescopes"]: config_struct = params["config"][telescope] tesselation = config_struct["tesselation"] moc_struct = {} - if params["doMinimalTiling"]: - if config_struct["FOV_type"] == "region" or config_struct["FOV"] < 1.0: - idxs = hp.pixelfunc.ang2pix( - map_struct["nside"], - tesselation[:, 1], - tesselation[:, 2], - lonlat=True, - ) - isin = np.isin(idxs, prob_indexes) - - idxs = [i for i, x in enumerate(isin) if x] - print("Keeping %d/%d tiles" % (len(idxs), len(tesselation))) - tesselation = tesselation[idxs, :] - if params["doParallel"]: moclists = Parallel( n_jobs=params["Ncores"], @@ -77,18 +53,27 @@ def create_moc(params, map_struct=None): ) if map_struct is not None: - ipix_keep = map_struct["ipix_keep"] + moc_keep = map_struct["moc_keep"] else: - ipix_keep = [] + moc_keep = None + if params["doMinimalTiling"]: moc_struct_new = copy.copy(moc_struct) if params["tilesType"] == "galaxy": tile_probs = gwemopt.tiles.compute_tiles_map( - params, moc_struct_new, prob, func="center", ipix_keep=ipix_keep + params, + moc_struct_new, + map_struct["skymap"], + func="center", + moc_keep=moc_keep, ) else: tile_probs = gwemopt.tiles.compute_tiles_map( - params, moc_struct_new, prob, func="np.sum(x)", ipix_keep=ipix_keep + params, + moc_struct_new, + map_struct["skymap"], + func="np.sum(x)", + moc_keep=moc_keep, ) keys = moc_struct.keys() @@ -96,7 +81,7 @@ def create_moc(params, map_struct=None): sort_idx = np.argsort(tile_probs)[::-1] csm = np.empty(len(tile_probs)) csm[sort_idx] = np.cumsum(tile_probs[sort_idx]) - ipix_keep = np.where(csm <= cl)[0] + ipix_keep = np.where(csm <= params["powerlaw_cl"])[0] moc_struct = {} cnt = 0 @@ -124,45 +109,40 @@ def Fov2Moc(params, config_struct, telescope, ra_pointing, dec_pointing, nside): moc_struct = {} - if "rotation" in params: - rotation = params["rotation"] - else: - rotation = None - if config_struct["FOV_type"] == "square": - ipix, radecs, patch, area = getSquarePixels( - ra_pointing, dec_pointing, config_struct["FOV"], nside, rotation=rotation + moc = getSquarePixels( + ra_pointing, + dec_pointing, + config_struct["FOV"], ) elif config_struct["FOV_type"] == "circle": - ipix, radecs, patch, area = getCirclePixels( - ra_pointing, dec_pointing, config_struct["FOV"], nside, rotation=rotation + moc = getCirclePixels( + ra_pointing, + dec_pointing, + config_struct["FOV"], ) elif config_struct["FOV_type"] == "region": region_file = os.path.join(CONFIG_DIR, config_struct["FOV"]) regions = Regions.read(region_file, format="ds9") - ipix, radecs, patch, area = getRegionPixels( - ra_pointing, dec_pointing, regions, nside, rotation=rotation + moc = getRegionPixels( + ra_pointing, + dec_pointing, + regions, + nside, ) else: raise ValueError("FOV_type must be square, circle or region") if params["doChipGaps"]: if telescope == "ZTF": - ipixs = get_ztf_quadrant_ipix(nside, ra_pointing, dec_pointing) - ipix = list({y for x in ipixs for y in x}) + moc = get_ztf_quadrant_moc(ra_pointing, dec_pointing) elif telescope == "DECam": - ipixs = get_decam_quadrant_ipix(nside, ra_pointing, dec_pointing) - ipix = list({y for x in ipixs for y in x}) + moc = get_decam_quadrant_moc(ra_pointing, dec_pointing) else: raise ValueError("Chip gaps only available for DECam and ZTF") moc_struct["ra"] = ra_pointing moc_struct["dec"] = dec_pointing - moc_struct["ipix"] = ipix - moc_struct["corners"] = radecs - moc_struct["patch"] = patch - moc_struct["area"] = area - - moc_struct["moc"] = [] + moc_struct["moc"] = moc return moc_struct diff --git a/gwemopt/params.py b/gwemopt/params.py index b0c9f9b1..07f04be4 100644 --- a/gwemopt/params.py +++ b/gwemopt/params.py @@ -99,7 +99,10 @@ def params_struct(opts): params["treasuremap_status"] = opts.treasuremap_status.split(",") - params["raslice"] = np.array(opts.raslice.split(","), dtype=float) + if opts.raslice is not None: + params["raslice"] = np.array(opts.raslice.split(","), dtype=float) + else: + params["raslice"] = None params["unbalanced_tiles"] = None params["filters"] = opts.filters.split(",") @@ -143,8 +146,10 @@ def params_struct(opts): else: params["end_time"] = time.Time(opts.end_time, format="isot", scale="utc") - params["inclination"] = ( - opts.inclination if hasattr(opts, "true_location") else False + params["inclination"] = opts.inclination if hasattr(opts, "inclination") else False + + params["projection"] = ( + opts.projection if hasattr(opts, "projection") else "astro mollweide" ) params["solverType"] = ( diff --git a/gwemopt/plotting/__init__.py b/gwemopt/plotting/__init__.py index 1a87c30f..e5f9b92a 100644 --- a/gwemopt/plotting/__init__.py +++ b/gwemopt/plotting/__init__.py @@ -1,28 +1,7 @@ -import os - -import matplotlib.pyplot as plt - from gwemopt.plotting.movie import make_movie -from gwemopt.plotting.observability import plot_observability from gwemopt.plotting.plot_coverage import make_coverage_plots from gwemopt.plotting.plot_efficiency import make_efficiency_plots from gwemopt.plotting.plot_inclination import plot_inclination from gwemopt.plotting.plot_schedule import make_schedule_plots from gwemopt.plotting.plot_skymap import plot_skymap from gwemopt.plotting.plot_tiles import make_tile_plots -from gwemopt.plotting.style import add_edges, cmap - - -def tauprob(params, tau, prob): - plotName = os.path.join(params["outputDir"], "tau_prob.pdf") - plt.figure() - plt.plot(tau, prob) - plt.grid() - plt.xscale("log") - plt.yscale("log") - plt.xlabel("Log of observing time $\tau$", fontsize=20) - plt.ylabel( - "Log of detection prob. given the target is at the observing field", fontsize=20 - ) - plt.savefig(plotName, dpi=200) - plt.close() diff --git a/gwemopt/plotting/observability.py b/gwemopt/plotting/observability.py deleted file mode 100644 index 70b7bf3a..00000000 --- a/gwemopt/plotting/observability.py +++ /dev/null @@ -1,70 +0,0 @@ -import healpy as hp -import matplotlib.pyplot as plt -import numpy as np -from tqdm import tqdm - -from gwemopt.plotting.movie import make_movie -from gwemopt.plotting.style import CBAR_BOOL, UNIT, add_edges, cmap - - -def plot_observability(params, map_struct): - """ - Function to plot the observability - """ - observability_struct = map_struct["observability"] - - for telescope in observability_struct.keys(): - plot_name = params["outputDir"].joinpath(f"observable_area_{telescope}.pdf") - - vals = map_struct["prob"] * observability_struct[telescope]["observability"] - vals[~observability_struct[telescope]["observability"].astype(bool)] = np.nan - - hp.mollview( - vals, - title=f"Observable Area - {telescope} (integrated)", - unit=UNIT, - cbar=CBAR_BOOL, - min=np.min(map_struct["prob"]), - max=np.max(map_struct["prob"]), - cmap=cmap, - ) - add_edges() - print(f"Saving to {plot_name}") - plt.savefig(plot_name, dpi=200) - plt.close() - - if params["doMovie"]: - moviedir = params["outputDir"].joinpath("movie") - moviedir.mkdir(parents=True, exist_ok=True) - - for telescope in observability_struct.keys(): - dts = list(observability_struct[telescope]["dts"].keys()) - dts = np.sort(dts) - - for ii, dt in tqdm(enumerate(dts), total=len(dts)): - plot_name = moviedir.joinpath(f"observability-{ii:04d}.png") - title = f"Observability Map: {dt:.2f} Days" - - vals = map_struct["prob"] * observability_struct[telescope]["dts"][dt] - vals[~observability_struct[telescope]["dts"][dt].astype(bool)] = np.nan - - hp.mollview( - vals, - title=title, - cbar=CBAR_BOOL, - unit=UNIT, - min=np.min(map_struct["prob"]), - max=np.max(map_struct["prob"]), - cmap=cmap, - ) - - add_edges() - plt.savefig(plot_name, dpi=200) - plt.close() - - moviefiles = moviedir.joinpath("observability-%04d.png") - filename = params["outputDir"].joinpath( - f"observability_timelapse_{telescope}.mpg" - ) - - make_movie(moviefiles, filename) diff --git a/gwemopt/plotting/plot_coverage.py b/gwemopt/plotting/plot_coverage.py index b888a1de..11cc495a 100644 --- a/gwemopt/plotting/plot_coverage.py +++ b/gwemopt/plotting/plot_coverage.py @@ -5,103 +5,48 @@ import numpy as np from astropy.time import Time from matplotlib.pyplot import cm +from tqdm import tqdm -from gwemopt.io import export_tiles_coverage_hist, export_tiles_coverage_int +from gwemopt.io import export_tiles_coverage_int from gwemopt.plotting.movie import make_movie -from gwemopt.plotting.style import CBAR_BOOL, UNIT, add_edges, add_sun_moon, cmap +from gwemopt.plotting.style import add_sun_moon from gwemopt.utils.geometry import angular_distance -def plot_mollweide_coverage(params, map_struct, coverage_struct): - """ - Plot the coverage in Mollweide projection. - """ - plot_name = params["outputDir"].joinpath("mollview_coverage.pdf") - plt.figure() - hp.mollview(map_struct["prob"], title="", unit=UNIT, cbar=CBAR_BOOL, cmap=cmap) - hp.projplot( - coverage_struct["data"][:, 0], - coverage_struct["data"][:, 1], - "wx", - lonlat=True, - coord="G", - ) - add_edges() - plt.savefig(plot_name, dpi=200) - plt.close() - - -def plot_coverage(params, coverage_struct): - """ - Plot the coverage - """ - plot_name = params["outputDir"].joinpath("coverage.pdf") - plt.figure(figsize=(10, 8)) - for ii in range(len(coverage_struct["ipix"])): - data = coverage_struct["data"][ii, :] - filt = coverage_struct["filters"][ii] - - if filt == "g": - color = "g" - elif filt == "r": - color = "r" - else: - color = "k" - - plt.scatter(data[2], data[5], s=20, color=color) - - plt.xlabel("Time [MJD]") - plt.ylabel("Tile Number") - plt.savefig(plot_name, dpi=200) - plt.close() - - def plot_tiles_coverage(params, map_struct, coverage_struct, plot_sun_moon=False): """ - Plot the tiles coverage in Mollweide projection. + Plot the tiles coverage """ plot_name = params["outputDir"].joinpath("tiles_coverage.pdf") - plt.figure() - hp.mollview(map_struct["prob"], title="", unit=UNIT, cbar=CBAR_BOOL, cmap=cmap) - add_edges() - ax = plt.gca() - for ii in range(len(coverage_struct["ipix"])): - patch = coverage_struct["patch"][ii] - if patch == []: - continue + hdu = map_struct["hdu"] + columns = [col.name for col in hdu.columns] + + fig = plt.figure(figsize=(8, 6), dpi=100) + ax = plt.axes( + [0.05, 0.05, 0.9, 0.9], + center=map_struct["center"], + projection=params["projection"], + ) + ax.imshow_hpx(hdu, field=columns.index("PROB"), cmap="cylon") + ax.grid() - if not type(patch) == list: - patch = [patch] + for ii in range(len(coverage_struct["moc"])): + moc = coverage_struct["moc"][ii] + data = coverage_struct["data"][ii] - for p in patch: - patch_cpy = copy.copy(p) - patch_cpy.axes = None - patch_cpy.figure = None - patch_cpy.set_transform(ax.transData) - hp.projaxes.HpxMollweideAxes.add_patch(ax, patch_cpy) + moc.fill( + ax=ax, wcs=ax.wcs, alpha=data[6], fill=True, color="black", linewidth=1 + ) + moc.border(ax=ax, wcs=ax.wcs, alpha=1, color="black") if plot_sun_moon: - add_sun_moon(params) + add_sun_moon(params, ax) plt.savefig(plot_name, dpi=200) plt.close() -def plot_tiles_coverage_hist(params, diffs): - """ - Plot the histogram of the tiles coverage. - """ - plot_name = params["outputDir"].joinpath("tiles_coverage_hist.pdf") - plt.figure(figsize=(12, 8)) - bins = np.linspace(0.0, 24.0, 25) - plt.hist(24.0 * np.array(diffs), bins=bins) - plt.xlabel("Difference Between Observations [hours]") - plt.ylabel("Number of Observations") - plt.savefig(plot_name, dpi=200) - plt.close("all") - - def plot_tiles_coverage_int( params, map_struct, catalog_struct, coverage_struct, plot_sun_moon=False ): @@ -110,19 +55,22 @@ def plot_tiles_coverage_int( colors = cm.rainbow(np.linspace(0, 1, len(params["telescopes"]))) + hdu = map_struct["hdu"] + columns = [col.name for col in hdu.columns] + plot_name = params["outputDir"].joinpath("tiles_coverage_int.pdf") fig = plt.figure(figsize=(12, 8)) gs = fig.add_gridspec(4, 1) - ax1 = fig.add_subplot(gs[0:3, 0], projection="astro hours mollweide") + ax1 = fig.add_subplot( + gs[0:3, 0], center=map_struct["center"], projection=params["projection"] + ) ax2 = fig.add_subplot(gs[3, 0]) ax3 = ax2.twinx() # mirror them plt.axes(ax1) - hp.mollview( - map_struct["prob"], title="", unit=UNIT, cbar=CBAR_BOOL, cmap=cmap, hold=True - ) - add_edges() + ax1.imshow_hpx(hdu, field=columns.index("PROB"), cmap="cylon") + ax1.grid() ax = plt.gca() if params["tilesType"] == "galaxy": @@ -136,30 +84,18 @@ def plot_tiles_coverage_int( color=color, ) else: - for ii in range(len(coverage_struct["ipix"])): - patch = coverage_struct["patch"][ii] - - idx = params["telescopes"].index(coverage_struct["telescope"][ii]) - - if patch == []: - continue - - if not type(patch) == list: - patch = [patch] - - for p in patch: - patch_cpy = copy.copy(p) - patch_cpy.axes = None - patch_cpy.figure = None - patch_cpy.set_transform(ax.transData) - patch_cpy.set_facecolor(colors[idx]) - - hp.projaxes.HpxMollweideAxes.add_patch(ax, patch_cpy) + for ii in range(len(coverage_struct["moc"])): + moc = coverage_struct["moc"][ii] + data = coverage_struct["data"][ii] + moc.fill( + ax=ax, wcs=ax.wcs, alpha=data[6], fill=True, color="black", linewidth=1 + ) + moc.border(ax=ax, wcs=ax.wcs, alpha=1, color="black") idxs = np.argsort(coverage_struct["data"][:, 2]) plt.axes(ax2) for telescope, color in zip(params["telescopes"], colors): - ipixs = np.empty((0, 2)) + moc = None cum_prob = 0.0 tts, cum_probs, cum_areas = [], [], [] @@ -171,7 +107,7 @@ def plot_tiles_coverage_int( print("%s: %d/%d" % (telescope, jj, len(idxs))) data = coverage_struct["data"][ii, :] - ipix = coverage_struct["ipix"][ii] + m = coverage_struct["moc"][ii] if params["tilesType"] == "galaxy": galaxies = coverage_struct["galaxies"][ii] @@ -196,11 +132,13 @@ def plot_tiles_coverage_int( cum_galaxies = np.unique(cum_galaxies).astype(int) cum_area = len(cum_galaxies) else: - ipixs = np.append(ipixs, ipix) - ipixs = np.unique(ipixs).astype(int) + if moc is None: + moc = m + else: + moc = moc.union(m) - cum_prob = np.sum(map_struct["prob"][ipixs]) - cum_area = len(ipixs) * map_struct["pixarea_deg2"] + cum_prob = moc.probability_in_multiordermap(map_struct["skymap"]) + cum_area = moc.sky_fraction * 360**2 / np.pi cum_probs.append(cum_prob) cum_areas.append(cum_area) @@ -223,7 +161,7 @@ def plot_tiles_coverage_int( else: ax3.set_ylabel("Sky area [sq. deg.]") - ipixs = np.empty((0, 2)) + moc = None cum_prob = 0.0 tts, cum_probs, cum_areas = [], [], [] @@ -232,7 +170,7 @@ def plot_tiles_coverage_int( for jj, ii in enumerate(idxs): data = coverage_struct["data"][ii, :] - ipix = coverage_struct["ipix"][ii] + m = coverage_struct["moc"][ii] if params["tilesType"] == "galaxy": galaxies = coverage_struct["galaxies"][ii] @@ -254,11 +192,13 @@ def plot_tiles_coverage_int( cum_galaxies = np.unique(cum_galaxies).astype(int) cum_area = len(cum_galaxies) else: - ipixs = np.append(ipixs, ipix) - ipixs = np.unique(ipixs).astype(int) + if moc is None: + moc = m + else: + moc = moc.union(m) - cum_prob = np.sum(map_struct["prob"][ipixs]) - cum_area = len(ipixs) * map_struct["pixarea_deg2"] + cum_prob = moc.probability_in_multiordermap(map_struct["skymap"]) + cum_area = moc.sky_fraction * 360**2 / np.pi tts.append(data[2] - event_mjd) cum_probs.append(cum_prob) @@ -287,7 +227,7 @@ def plot_tiles_coverage_int( ax2.legend(loc=1) if plot_sun_moon: - add_sun_moon(params) + add_sun_moon(params, ax) plt.savefig(plot_name, dpi=200) plt.close() @@ -295,41 +235,39 @@ def plot_tiles_coverage_int( def plot_coverage_scaled(params, map_struct, coverage_struct, plot_sun_moon, max_time): plot_name = params["outputDir"].joinpath("tiles_coverage_scaled.pdf") - plt.figure() - hp.mollview(map_struct["prob"], title="", unit=UNIT, cbar=CBAR_BOOL, cmap=cmap) - add_edges() - ax = plt.gca() - for ii in range(len(coverage_struct["ipix"])): - data = coverage_struct["data"][ii, :] - patch = coverage_struct["patch"][ii] - if patch == []: - continue + hdu = map_struct["hdu"] + columns = [col.name for col in hdu.columns] + + fig = plt.figure(figsize=(8, 6), dpi=100) + ax = plt.axes( + [0.05, 0.05, 0.9, 0.9], + center=map_struct["center"], + projection=params["projection"], + ) + ax.imshow_hpx(hdu, field=columns.index("PROB"), cmap="cylon") + ax.grid() - if not type(patch) == list: - patch = [patch] + for ii in range(len(coverage_struct["moc"])): + moc = coverage_struct["moc"][ii] + data = coverage_struct["data"][ii] - for p in patch: - patch_cpy = copy.copy(p) - patch_cpy.axes = None - patch_cpy.figure = None - patch_cpy.set_transform(ax.transData) - current_alpha = patch_cpy.get_alpha() + alpha = data[4] / max_time + if alpha > 1: + alpha = 1.0 - if current_alpha > 0.0: - alpha = data[4] / max_time - if alpha > 1: - alpha = 1.0 - patch_cpy.set_alpha(alpha) - hp.projaxes.HpxMollweideAxes.add_patch(ax, patch_cpy) + moc.fill(ax=ax, wcs=ax.wcs, alpha=alpha, fill=True, color="black", linewidth=1) + moc.border(ax=ax, wcs=ax.wcs, alpha=1, color="black") if plot_sun_moon: - add_sun_moon(params) + add_sun_moon(params, ax) plt.savefig(plot_name, dpi=200) plt.close() if params["doMovie"]: + print("Creating movie from schedule...") + idx = np.isfinite(coverage_struct["data"][:, 2]) mjd_min = np.min(coverage_struct["data"][idx, 2]) mjd_max = np.max(coverage_struct["data"][idx, 2]) @@ -339,34 +277,37 @@ def plot_coverage_scaled(params, map_struct, coverage_struct, plot_sun_moon, max moviedir = params["outputDir"].joinpath("movie") moviedir.mkdir(exist_ok=True, parents=True) - for jj in range(len(mjds)): + for jj in tqdm(range(len(mjds))): mjd = mjds[jj] plot_name = moviedir.joinpath(f"coverage-{jj:04d}.png") title = f"Coverage Map: {mjd:.2f}" - plt.figure() - hp.mollview( - map_struct["prob"], title=title, unit=UNIT, cbar=CBAR_BOOL, cmap=cmap + fig = plt.figure(figsize=(8, 6), dpi=100) + ax = plt.axes( + [0.05, 0.05, 0.9, 0.9], + center=map_struct["center"], + projection=params["projection"], ) - add_edges() - ax = plt.gca() + ax.imshow_hpx(hdu, field=columns.index("PROB"), cmap="cylon") + ax.grid() idx = np.where(coverage_struct["data"][:, 2] <= mjd)[0] for ii in idx: - patch = coverage_struct["patch"][ii] - - if patch == []: - continue + moc = coverage_struct["moc"][ii] + data = coverage_struct["data"][ii] - if not type(patch) == list: - patch = [patch] + alpha = data[4] / max_time + if alpha > 1: + alpha = 1.0 - for p in patch: - patch_cpy = copy.copy(p) - patch_cpy.axes = None - patch_cpy.figure = None - patch_cpy.set_transform(ax.transData) - hp.projaxes.HpxMollweideAxes.add_patch(ax, patch_cpy) + moc.fill( + ax=ax, + wcs=ax.wcs, + alpha=alpha, + fill=True, + color="black", + linewidth=1, + ) plt.savefig(plot_name, dpi=200) plt.close() @@ -380,7 +321,6 @@ def plot_coverage_scaled(params, map_struct, coverage_struct, plot_sun_moon, max def make_coverage_plots( params, map_struct, coverage_struct, catalog_struct=None, plot_sun_moon: bool = True ): - plot_mollweide_coverage(params, map_struct, coverage_struct) idx = np.isfinite(coverage_struct["data"][:, 4]) if not idx.size: @@ -388,59 +328,6 @@ def make_coverage_plots( max_time = np.max(coverage_struct["data"][idx, 4]) - plot_coverage(params, coverage_struct) - - plot_tiles_coverage( - params, map_struct, coverage_struct, plot_sun_moon=plot_sun_moon - ) - - diffs = [] - if params["tilesType"] == "galaxy": - coverage_ras = coverage_struct["data"][:, 0] - coverage_decs = coverage_struct["data"][:, 1] - - for ii in range(len(coverage_ras) - 1): - current_ra, current_dec = coverage_ras[ii], coverage_decs[ii] - - dist = angular_distance( - current_ra, current_dec, coverage_ras[ii + 1 :], coverage_decs[ii + 1 :] - ) - idx = np.where(dist <= 1 / 3600.0)[0] - if len(idx) > 0: - jj = idx[0] - diffs.append( - np.abs( - coverage_struct["data"][ii, 2] - coverage_struct["data"][jj, 2] - ) - ) - else: - for ii in range(len(coverage_struct["ipix"])): - ipix = coverage_struct["ipix"][ii] - for jj in range(len(coverage_struct["ipix"])): - if ii >= jj: - continue - if coverage_struct["telescope"][ii] == coverage_struct["telescope"][jj]: - continue - ipix2 = coverage_struct["ipix"][jj] - overlap = np.intersect1d(ipix, ipix2) - rat = np.array( - [ - float(len(overlap)) / float(len(ipix)), - float(len(overlap)) / float(len(ipix2)), - ] - ) - if np.any(rat > 0.5): - diffs.append( - np.abs( - coverage_struct["data"][ii, 2] - - coverage_struct["data"][jj, 2] - ) - ) - - export_tiles_coverage_hist(params, diffs) - - plot_tiles_coverage_hist(params, diffs) - plot_tiles_coverage_int( params, map_struct, catalog_struct, coverage_struct, plot_sun_moon ) diff --git a/gwemopt/plotting/plot_efficiency.py b/gwemopt/plotting/plot_efficiency.py index efc87c0f..614051ab 100644 --- a/gwemopt/plotting/plot_efficiency.py +++ b/gwemopt/plotting/plot_efficiency.py @@ -2,8 +2,6 @@ import matplotlib.pyplot as plt import numpy as np -from gwemopt.plotting.style import CBAR_BOOL, UNIT, add_edges - def make_efficiency_plots(params, map_struct, efficiency_structs): plot_name = params["outputDir"].joinpath("efficiency.pdf") @@ -30,16 +28,6 @@ def make_efficiency_plots(params, map_struct, efficiency_structs): plt.savefig(plot_name, dpi=200) plt.close() - plot_name = params["outputDir"].joinpath("mollview_injs.pdf") - plt.figure() - hp.mollview(map_struct["prob"], unit=UNIT, cbar=CBAR_BOOL) - hp.projplot( - efficiency_struct["ra"], efficiency_struct["dec"], "wx", lonlat=True, coord="G" - ) - add_edges() - plt.savefig(plot_name, dpi=200) - plt.close() - if params["do_3d"]: for key in efficiency_structs: efficiency_struct = efficiency_structs[key] diff --git a/gwemopt/plotting/plot_skymap.py b/gwemopt/plotting/plot_skymap.py index 14b56308..e33c11ec 100644 --- a/gwemopt/plotting/plot_skymap.py +++ b/gwemopt/plotting/plot_skymap.py @@ -1,64 +1,33 @@ import healpy as hp +import ligo.skymap.plot import matplotlib.pyplot as plt import numpy as np +from astropy.io import fits +from ligo.skymap import moc +from matplotlib import pyplot as plt -from gwemopt.plotting.style import CBAR_BOOL, UNIT, add_edges, cmap - -def plot_skymap(params, map_struct): +def plot_skymap(params, map_struct, colnames=["PROB", "DISTMEAN", "DISTSTD"]): """ Function to plot the skymap """ plot_name = params["outputDir"].joinpath("prob.pdf") - if np.percentile(map_struct["prob"], 99) > 0: - hp.mollview( - map_struct["prob"], - title="", - unit=UNIT, - cbar=CBAR_BOOL, - min=np.percentile(map_struct["prob"], 1), - max=np.percentile(map_struct["prob"], 99), - cmap=cmap, - ) - else: - hp.mollview( - map_struct["prob"], - title="", - unit=UNIT, - cbar=CBAR_BOOL, - min=np.percentile(map_struct["prob"], 1), - cmap=cmap, - ) - - add_edges() - plt.savefig(plot_name, dpi=200) - plt.close() - - if "distmu" in map_struct: - fin = np.copy(map_struct["distmu"]) - fin[~np.isfinite(fin)] = np.nan - plot_name = params["outputDir"].joinpath("dist.pdf") - hp.mollview( - map_struct["distmu"], - unit="Distance [Mpc]", - min=np.nanpercentile(fin, 10), - max=np.nanpercentile(fin, 90), - ) - add_edges() - plt.savefig(plot_name, dpi=200) - plt.close() - - fin = np.copy(map_struct["distmed"]) - fin[~np.isfinite(fin)] = np.nan - plot_name = params["outputDir"].joinpath("dist_median.pdf") - hp.mollview( - map_struct["distmed"], - unit="Distance [Mpc]", - min=np.nanpercentile(fin, 10), - max=np.nanpercentile(fin, 90), - ) - add_edges() - plt.savefig(plot_name, dpi=200) - plt.close() + hdu = map_struct["hdu"] + columns = [col.name for col in hdu.columns] + + for col in colnames: + if col in columns: + + fig = plt.figure(figsize=(8, 6), dpi=100) + ax = plt.axes( + [0.05, 0.05, 0.9, 0.9], + center=map_struct["center"], + projection=params["projection"], + ) + ax.imshow_hpx(hdu, field=columns.index(col), cmap="cylon") + ax.grid() + plot_name = params["outputDir"].joinpath(f"{col.lower()}.pdf") + plt.savefig(plot_name, dpi=200) + plt.close() diff --git a/gwemopt/plotting/plot_tiles.py b/gwemopt/plotting/plot_tiles.py index 7cc3736e..7fe6568c 100644 --- a/gwemopt/plotting/plot_tiles.py +++ b/gwemopt/plotting/plot_tiles.py @@ -3,8 +3,9 @@ import healpy as hp import matplotlib.pyplot as plt import numpy as np +from tqdm import tqdm -from gwemopt.plotting.style import CBAR_BOOL, UNIT, add_edges, add_sun_moon, cmap +from gwemopt.plotting.style import add_sun_moon def make_tile_plots(params, map_struct, tiles_structs, plot_sun_moon=True): @@ -13,62 +14,37 @@ def make_tile_plots(params, map_struct, tiles_structs, plot_sun_moon=True): """ plot_name = params["outputDir"].joinpath("tiles.pdf") - plt.figure() - hp.mollview(map_struct["prob"], title="", unit=UNIT, cbar=CBAR_BOOL, cmap=cmap) - ax = plt.gca() + + hdu = map_struct["hdu"] + columns = [col.name for col in hdu.columns] + + fig = plt.figure(figsize=(8, 6), dpi=100) + ax = plt.axes( + [0.05, 0.05, 0.9, 0.9], + center=map_struct["center"], + projection=params["projection"], + ) + ax.imshow_hpx(hdu, field=columns.index("PROB"), cmap="cylon") for telescope in tiles_structs: tiles_struct = tiles_structs[telescope] - for index in tiles_struct.keys(): - patch = tiles_struct[index]["patch"] - if not patch: - continue - if type(patch) == list: - for p in patch: - hp.projaxes.HpxMollweideAxes.add_patch(ax, p) - else: - hp.projaxes.HpxMollweideAxes.add_patch(ax, patch) + keys = list(tiles_struct.keys()) + probs = np.array([tiles_struct[key]["prob"] for key in keys]) + alphas = probs / np.max(probs) + + for ii, index in tqdm(enumerate(keys), total=len(keys)): + moc = tiles_struct[index]["moc"] + moc.fill( + ax=ax, + wcs=ax.wcs, + alpha=alphas[ii], + fill=True, + color="black", + linewidth=1, + ) + moc.border(ax=ax, wcs=ax.wcs, alpha=1, color="black") if plot_sun_moon: - add_sun_moon(params) + add_sun_moon(params, ax) - add_edges() plt.savefig(plot_name, dpi=200) plt.close() - - if params["doReferences"]: - rot = (90, 0, 0) - plot_name = params["outputDir"].joinpath("tiles_ref.pdf") - plt.figure() - hp.mollview( - np.zeros(map_struct["prob"].shape), - title="", - unit=UNIT, - cbar=CBAR_BOOL, - cmap=cmap, - ) - ax = plt.gca() - for telescope in tiles_structs: - config_struct = params["config"][telescope] - tiles_struct = tiles_structs[telescope] - if "reference_images" not in config_struct: - continue - for index in tiles_struct.keys(): - if index not in config_struct["reference_images"]: - continue - if len(params["filters"]) == 1: - if ( - not params["filters"][0] - in config_struct["reference_images"][index] - ): - continue - patch = tiles_struct[index]["patch"] - if not patch: - continue - patch_cpy = copy.copy(patch) - patch_cpy.axes = None - patch_cpy.figure = None - patch_cpy.set_transform(ax.transData) - hp.projaxes.HpxMollweideAxes.add_patch(ax, patch_cpy) - add_edges() - plt.savefig(plot_name, dpi=200) - plt.close() diff --git a/gwemopt/plotting/style.py b/gwemopt/plotting/style.py index 3275c146..743c5ae6 100644 --- a/gwemopt/plotting/style.py +++ b/gwemopt/plotting/style.py @@ -13,27 +13,8 @@ matplotlib.rcParams.update({"font.size": 16}) matplotlib.rcParams["contour.negative_linestyle"] = "solid" -UNIT = "Gravitational-wave probability" -CBAR_BOOL = False - -def add_edges(): - """ - Add longitude and latitude lines to a healpix map. - """ - hp.graticule() - plt.grid(True) - lons = np.arange(-150.0, 180, 30.0) - lats = np.zeros(lons.shape) - for lon, lat in zip(lons, lats): - hp.projtext(lon, lat, "%.0f" % lon, lonlat=True) - lats = np.arange(-60.0, 90, 30.0) - lons = np.zeros(lons.shape) - for lon, lat in zip(lons, lats): - hp.projtext(lon, lat, "%.0f" % lat, lonlat=True) - - -def add_sun_moon(params): +def add_sun_moon(params, ax): """ Add sun and moon position to the skymap """ @@ -44,11 +25,17 @@ def add_sun_moon(params): 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" + ax.plot( + sun_position.ra, + sun_position.dec, + transform=ax.get_transform("world"), + color="yellow", ) - hp.visufunc.projscatter( - moon_position.ra, moon_position.dec, lonlat=True, color="grey" + ax.plot( + moon_position.ra, + moon_position.dec, + transform=ax.get_transform("world"), + color="grey", ) # also plot (in smaller scale) sun and moon position for the 7 following days @@ -62,23 +49,23 @@ def add_sun_moon(params): new_time = Time(new_time, format="gps", scale="utc") new_moon_position = get_body("moon", new_time) - hp.visufunc.projscatter( + ax.plot( new_moon_position.ra, new_moon_position.dec, - lonlat=True, + transform=ax.get_transform("world"), color="black", marker=".", - s=1, + markersize=1, ) if not i % 8: # only plot point for the sun every two days in order to avoid overlap new_sun_position = get_body("sun", new_time) - hp.visufunc.projscatter( + ax.plot( new_sun_position.ra, new_sun_position.dec, - lonlat=True, + transform=ax.get_transform("world"), color="black", marker=".", - s=1, + markersize=1, ) diff --git a/gwemopt/run.py b/gwemopt/run.py index a389be7f..d0bc03cc 100644 --- a/gwemopt/run.py +++ b/gwemopt/run.py @@ -19,7 +19,6 @@ make_efficiency_plots, make_tile_plots, plot_inclination, - plot_observability, plot_skymap, ) from gwemopt.utils import calculate_observability @@ -71,45 +70,6 @@ def run(args=None): if args.inclination: plot_inclination(params, map_struct) - if args.doObservability: - print("Calculating observability") - observability_struct = calculate_observability(params, map_struct) - map_struct["observability"] = observability_struct - if args.doPlots: - print("Plotting observability...") - plot_observability(params, map_struct) - - if args.doObservabilityExit: - for telescope in params["telescopes"]: - if ( - np.sum(observability_struct[telescope]["prob"]) - < args.observability_thresh - ): - print( - "Observability for %s: %.5f < %.5f... exiting." - % ( - telescope, - np.sum(observability_struct[telescope]["prob"]), - args.observability_thresh, - ) - ) - - if params["doTrueLocation"]: - lightcurve_structs = gwemopt.lightcurve.read_files( - params, params["lightcurveFiles"] - ) - for key in lightcurve_structs.keys(): - filename = os.path.join( - params["outputDir"], - "efficiency_true_" - + lightcurve_structs[key]["name"] - + ".txt", - ) - fid = open(filename, "w") - fid.write("0") - fid.close() - exit(0) - if params["splitType"] is not None: print("Splitting skymap...") map_struct["groups"] = gwemopt.mapsplit.similar_range(params, map_struct) diff --git a/gwemopt/scheduler.py b/gwemopt/scheduler.py index 3ec9efa6..b60ee2e5 100644 --- a/gwemopt/scheduler.py +++ b/gwemopt/scheduler.py @@ -118,13 +118,13 @@ def get_order_heuristic( for jj, key in enumerate(keys): tileprobs[jj] = tile_struct[key]["prob"] tilenexps[jj] = tile_struct[key]["nexposures"] - try: + + if type(tile_struct[key]["exposureTime"]) in [float, np.float64]: tileexpdur[jj] = tile_struct[key]["exposureTime"] - except: - try: - tileexpdur[jj] = tile_struct[key]["exposureTime"][0] - except: - tileexpdur[jj] = 0.0 + elif type(tile_struct[key]["exposureTime"]) in [list, np.ndarray]: + tileexpdur[jj] = tile_struct[key]["exposureTime"][0] + else: + tileexpdur[jj] = 0.0 tilefilts[key] = copy.deepcopy(tile_struct[key]["filt"]) tileavailable_tiles[jj] = [] @@ -603,9 +603,7 @@ def scheduler(params, config_struct, tile_struct): coverage_struct = {} coverage_struct["data"] = np.empty((0, 8)) coverage_struct["filters"] = [] - coverage_struct["ipix"] = [] - coverage_struct["patch"] = [] - coverage_struct["area"] = [] + coverage_struct["moc"] = [] if params["tilesType"] == "galaxy": coverage_struct["galaxies"] = [] @@ -714,13 +712,10 @@ def scheduler(params, config_struct, tile_struct): ) coverage_struct["filters"].append(filt) - coverage_struct["patch"].append(tile_struct_hold["patch"]) - coverage_struct["ipix"].append(tile_struct_hold["ipix"]) - coverage_struct["area"].append(tile_struct_hold["area"]) + coverage_struct["moc"].append(tile_struct_hold["moc"]) if params["tilesType"] == "galaxy": coverage_struct["galaxies"].append(tile_struct_hold["galaxies"]) - coverage_struct["area"] = np.array(coverage_struct["area"]) coverage_struct["filters"] = np.array(coverage_struct["filters"]) coverage_struct["FOV"] = [config_struct["FOV"]] * len(coverage_struct["filters"]) coverage_struct["telescope"] = [config_struct["telescope"]] * len( @@ -743,166 +738,3 @@ def computeSlewReadoutTime(config_struct, coverage_struct): prev_dec = dat[0] prev_ra = dat[1] return acc_time - - -def schedule_ra_splits( - params, - config_struct, - map_struct_hold, - tile_struct, - telescope, - previous_coverage_struct, -): - location = astropy.coordinates.EarthLocation( - config_struct["longitude"], - config_struct["latitude"], - config_struct["elevation"], - ) - - raslices = gwemopt.utils.utils.auto_rasplit( - params, map_struct_hold, params["nside_down"] - ) - - maxidx = 0 - coverage_structs = [] - skip = False - while len(raslices) != 0: - params["unbalanced_tiles"] = [] - config_struct["exposurelist"] = segments.segmentlist( - config_struct["exposurelist"][maxidx:] - ) - if len(config_struct["exposurelist"]) < 2: - break - - map_struct_slice = copy.deepcopy(map_struct_hold) - - exposurelist = np.array_split(config_struct["exposurelist"], len(raslices))[0] - minhas = [] - minhas_late = [] - try_end = False - if len(raslices) == 1: - raslice = raslices[0] - del raslices[0] - else: - for raslice in raslices: - has = [] - has_late = [] - for seg in exposurelist: - mjds = np.linspace(seg[0], seg[1], 100) - tt = Time(mjds, format="mjd", scale="utc", location=location) - lst = tt.sidereal_time("mean") / u.hourangle - ha = np.abs(lst - raslice[0]) - ha_late = np.abs(lst - raslice[1]) - - idx = np.where(ha > 12.0)[0] - ha[idx] = 24.0 - ha[idx] - idx_late = np.where(ha_late > 12.0)[0] - ha_late[idx_late] = 24.0 - ha_late[idx_late] - has += list(ha) - has_late += list(ha_late) - if len(has) > 0: - minhas.append(np.min(has)) - if len(has_late) > 0: - minhas_late.append(np.min(has_late)) - - if (len(minhas_late) > 0) and (len(has_late) > 0): - # conditions for trying to schedule end of slice - if np.min(minhas_late) <= 5.0 and np.min(has) > 4.0 and not skip: - try_end = True - min = np.argmin(minhas_late) - raslice = raslices[min] - else: - min = np.argmin(minhas) - raslice = raslices[min] - del raslices[min] - else: - min = np.argmin(minhas) - raslice = raslices[min] - del raslices[min] - - # do RA slicing - ra_low, ra_high = raslice[0], raslice[1] - ra = map_struct_slice["ra"] - if ra_low <= ra_high: - ipix = np.where( - (ra_high * 360.0 / 24.0 < ra) | (ra_low * 360.0 / 24.0 > ra) - )[0] - else: - ipix = np.where( - (ra_high * 360.0 / 24.0 < ra) & (ra_low * 360.0 / 24.0 > ra) - )[0] - - map_struct_slice["prob"][ipix] = 0.0 - - if params["timeallocationType"] == "absmag": - tile_struct = gwemopt.tiles.absmag_tiles_struct( - params, config_struct, telescope, map_struct_slice, tile_struct - ) - else: - tile_struct = gwemopt.tiles.powerlaw_tiles_struct( - params, config_struct, telescope, map_struct_slice, tile_struct - ) - - config_struct_hold = copy.copy(config_struct) - coverage_struct, tile_struct = schedule_alternating( - params, - config_struct_hold, - telescope, - map_struct_slice, - tile_struct, - previous_coverage_struct, - ) - if len(coverage_struct["ipix"]) == 0: - continue - optimized_max, coverage_struct, tile_struct = optimize_max_tiles( - params, - tile_struct, - coverage_struct, - config_struct, - telescope, - map_struct_slice, - ) - params["max_nb_tiles"] = np.array([optimized_max], dtype=float) - balanced_fields = 0 - coverage_struct, tile_struct = schedule_alternating( - params, - config_struct, - telescope, - map_struct_slice, - tile_struct, - previous_coverage_struct, - ) - - doReschedule, balanced_fields = balance_tiles( - params, tile_struct, coverage_struct - ) - config_struct_hold = copy.copy(config_struct) - - if balanced_fields == 0: - if try_end: - skip = True - continue - elif try_end: - del raslices[min] - skip = False - - if len(coverage_struct["exposureused"]) > 0: - maxidx = int(coverage_struct["exposureused"][-1]) - - coverage_struct = gwemopt.utils.utils.erase_unbalanced_tiles( - params, coverage_struct - ) - - # limit to max number of filter sets - if len(coverage_structs) < params["max_filter_sets"]: - coverage_structs.append(coverage_struct) - else: - prob_structs = [ - np.sum(prev_struct["data"][:, 6]) for prev_struct in coverage_structs - ] - if np.any(np.array(prob_structs) < np.sum(coverage_struct["data"][:, 6])): - argmin = np.argmin(prob_structs) - del coverage_structs[argmin] - coverage_structs.append(coverage_struct) - - return gwemopt.coverage.combine_coverage_structs(coverage_structs) diff --git a/gwemopt/skyportal.py b/gwemopt/skyportal.py index 372ba1df..b4ccca77 100644 --- a/gwemopt/skyportal.py +++ b/gwemopt/skyportal.py @@ -7,12 +7,7 @@ from regions import CircleSkyRegion, PolygonSkyRegion, RectangleSkyRegion import gwemopt.tiles -from gwemopt.tiles import ( - absmag_tiles_struct, - angular_distance, - get_rectangle, - powerlaw_tiles_struct, -) +from gwemopt.tiles import angular_distance, get_rectangle, powerlaw_tiles_struct def create_galaxy_from_skyportal(params, map_struct, catalog_struct, regions=None): @@ -194,19 +189,14 @@ def create_galaxy_from_skyportal(params, map_struct, catalog_struct, regions=Non moc_struct[int(cnt)]["galaxies"] = galaxies cnt = cnt + 1 - if params["timeallocationType"] == "absmag": - tile_struct = absmag_tiles_struct( - params, config_struct, telescope, map_struct, moc_struct - ) - elif params["timeallocationType"] == "powerlaw": - tile_struct = powerlaw_tiles_struct( - params, - config_struct, - telescope, - map_struct, - moc_struct, - catalog_struct=catalog_struct, - ) + tile_struct = powerlaw_tiles_struct( + params, + config_struct, + telescope, + map_struct, + moc_struct, + catalog_struct=catalog_struct, + ) tile_struct = gwemopt.segments.get_segments_tiles( params, config_struct, tile_struct @@ -283,13 +273,13 @@ def create_moc_from_skyportal(params, map_struct=None, field_ids=None): config_struct = params["config"][telescope] tesselation = config_struct["tesselation"] moc_struct = {} - ipixs = [] + mocs = [] for ii, tess in enumerate(tesselation): if field_ids is not None: if tess.field_id not in field_ids[telescope]: - ipixs.append([]) + mocs.append(MOC.new_empty(29)) continue - ipixs.append(skyportal2FOV(tess, nside)) + mocs.append(skyportal2FOV(tess, nside)) for ii, tess in enumerate(tesselation): index = tess.field_id @@ -298,38 +288,24 @@ def create_moc_from_skyportal(params, map_struct=None, field_ids=None): if (telescope == "ZTF") and doUseSecondary and (index < 1000): continue - ipix = ipixs[ii] - if len(ipix) == 0: + moc = mocs[ii] + if moc.empty(): continue moc_struct[index] = {} moc_struct[index]["ra"] = tess.ra moc_struct[index]["dec"] = tess.dec - - moc_struct[index]["ipix"] = ipix - moc_struct[index]["corners"] = [ - [np.min(ra[ipix]), np.min(dec[ipix])], - [np.min(ra[ipix]), np.max(dec[ipix])], - [np.max(ra[ipix]), np.max(dec[ipix])], - [np.max(ra[ipix]), np.min(dec[ipix])], - ] - moc_struct[index]["patch"] = [] - moc_struct[index]["area"] = len(ipix) * pixarea - - if map_struct is not None: - ipix_keep = map_struct["ipix_keep"] - else: - ipix_keep = [] + moc_struct[index]["moc"] = moc if params["doMinimalTiling"]: moc_struct_new = copy.copy(moc_struct) if params["tilesType"] == "galaxy": tile_probs = gwemopt.tiles.compute_tiles_map( - params, moc_struct_new, prob, func="center", ipix_keep=ipix_keep + params, moc_struct_new, prob, func="center" ) else: tile_probs = gwemopt.tiles.compute_tiles_map( - params, moc_struct_new, prob, func="np.sum(x)", ipix_keep=ipix_keep + params, moc_struct_new, prob, func="np.sum(x)" ) keys = moc_struct.keys() @@ -337,13 +313,13 @@ def create_moc_from_skyportal(params, map_struct=None, field_ids=None): sort_idx = np.argsort(tile_probs)[::-1] csm = np.empty(len(tile_probs)) csm[sort_idx] = np.cumsum(tile_probs[sort_idx]) - ipix_keep = np.where(csm <= cl)[0] + moc_keep = np.where(csm <= cl)[0] probs = [] moc_struct = {} cnt = 0 for ii, key in enumerate(keys): - if ii in ipix_keep: + if ii in moc_keep: moc_struct[key] = moc_struct_new[key] cnt = cnt + 1 @@ -354,13 +330,8 @@ def create_moc_from_skyportal(params, map_struct=None, field_ids=None): def skyportal2FOV(tess, nside): moc = moc_from_tiles([tile.healpix for tile in tess.tiles], 2**29) - pix_id = moc.degrade_to_order(int(np.log2(nside))).flatten() - if len(pix_id) > 0: - ipix = hp.nest2ring(int(nside), pix_id.tolist()) - else: - ipix = [] - return ipix + return moc def moc_from_tiles(rangeSet, nside): diff --git a/gwemopt/tests/data/expected_results/coverage_DECam.dat b/gwemopt/tests/data/expected_results/coverage_DECam.dat index e81308fa..f5baab5c 100644 --- a/gwemopt/tests/data/expected_results/coverage_DECam.dat +++ b/gwemopt/tests/data/expected_results/coverage_DECam.dat @@ -1,5 +1,5 @@ -1000246 0.2923655835 1 1 -1000823 0.1599655995 1 1 -1000381 0.1358876298 1 1 -1002413 0.0474212709 1 1 -1000165 0.0428240474 1 1 +1000246 0.2923700000 1 1 +1000823 0.1599700000 1 1 +1000381 0.1358900000 1 1 +1002413 0.0474200000 1 1 +1000165 0.0428200000 1 1 diff --git a/gwemopt/tests/data/expected_results/coverage_KPED.dat b/gwemopt/tests/data/expected_results/coverage_KPED.dat index 12c933f4..2b45d9e5 100644 --- a/gwemopt/tests/data/expected_results/coverage_KPED.dat +++ b/gwemopt/tests/data/expected_results/coverage_KPED.dat @@ -1,178 +1,154 @@ -0 0.0465378779 1 1 -1 0.0402247012 1 1 -5 0.0356842995 1 1 -2 0.0325496266 1 1 -3 0.0312539172 1 1 -8 0.0258489281 1 1 -4 0.0242680407 1 1 -6 0.0224916673 1 1 -7 0.0222180370 1 1 -12 0.0207573647 1 1 -24 0.0181768065 1 1 -9 0.0165753933 1 1 -36 0.0152967863 1 1 -23 0.0148210156 1 1 -10 0.0147876985 1 1 -11 0.0147330794 1 1 -13 0.0141169770 1 1 -14 0.0139952222 1 1 -15 0.0129470393 1 1 -16 0.0127856048 1 1 -17 0.0123939346 1 1 -41 0.0119388785 1 1 -18 0.0119253888 1 1 -19 0.0114672691 1 1 -43 0.0113705133 1 1 -20 0.0106183452 1 1 -21 0.0105591114 1 1 -32 0.0103824061 1 1 -22 0.0102762503 1 1 -25 0.0096948466 1 1 -40 0.0094337458 1 1 -37 0.0091500256 1 1 -26 0.0090819922 1 1 -27 0.0090310087 1 1 -28 0.0089964061 1 1 -29 0.0088269543 1 1 -30 0.0087371235 1 1 -33 0.0086458275 1 1 -45 0.0084772760 1 1 -48 0.0083434382 1 1 -34 0.0083355762 1 1 -35 0.0082674117 1 1 -46 0.0077663218 1 1 -51 0.0075851678 1 1 -38 0.0075122846 1 1 -39 0.0071812794 1 1 -56 0.0069248293 1 1 -42 0.0064726443 1 1 -72 0.0059788282 1 1 -44 0.0059045362 1 1 -63 0.0055916638 1 1 -47 0.0053229565 1 1 -54 0.0049203479 1 1 -49 0.0047591335 1 1 -70 0.0047122402 1 1 -50 0.0045860802 1 1 -52 0.0043533822 1 1 -53 0.0043188309 1 1 -55 0.0040491821 1 1 -61 0.0037364099 1 1 -117 0.0037145856 1 1 -64 0.0035426903 1 1 -57 0.0033823788 1 1 -59 0.0033717646 1 1 -58 0.0033684299 1 1 -73 0.0033623236 1 1 -81 0.0033384889 1 1 -84 0.0033015882 1 1 -60 0.0032635323 1 1 -93 0.0032330228 1 1 -62 0.0031477612 1 1 -75 0.0031280182 1 1 -82 0.0031159666 1 1 -65 0.0030052426 1 1 -107 0.0030019255 1 1 -113 0.0029395392 1 1 -66 0.0029058988 1 1 -67 0.0029047043 1 1 -68 0.0028976559 1 1 -69 0.0028366368 1 1 -76 0.0027227700 1 1 -31 0.0026626768 1 1 -88 0.0026566029 1 1 -112 0.0026344918 1 1 -71 0.0026091897 1 1 -74 0.0025718578 1 1 -118 0.0025250640 1 1 -78 0.0024887992 1 1 -97 0.0024401767 1 1 -77 0.0024394263 1 1 -79 0.0024118421 1 1 -86 0.0023867562 1 1 -80 0.0023800642 1 1 -115 0.0023007850 1 1 -111 0.0023000581 1 1 -148 0.0022800913 1 1 -83 0.0022166857 1 1 -85 0.0021326208 1 1 -116 0.0021077219 1 1 -87 0.0020811817 1 1 -125 0.0020811038 1 1 -96 0.0020126155 1 1 -89 0.0019433047 1 1 -90 0.0019383143 1 1 -91 0.0018575720 1 1 -92 0.0018411230 1 1 -99 0.0018106583 1 1 -94 0.0018083837 1 1 -95 0.0017901713 1 1 -108 0.0017392912 1 1 -98 0.0017373564 1 1 -100 0.0017140495 1 1 -101 0.0017137062 1 1 -102 0.0016839795 1 1 -103 0.0016392850 1 1 -142 0.0016386415 1 1 -104 0.0016321287 1 1 -105 0.0016316701 1 1 -106 0.0016237516 1 1 -127 0.0015730792 1 1 -114 0.0014979805 1 1 -137 0.0014667016 1 1 -109 0.0014605493 1 1 -110 0.0014507755 1 1 -151 0.0014349776 1 1 -134 0.0014052586 1 1 -164 0.0013861600 1 1 -157 0.0013640005 1 1 -120 0.0013403700 1 1 -143 0.0013114654 1 1 -156 0.0012757036 1 1 -119 0.0012651160 1 1 -149 0.0012594342 1 1 -121 0.0012282309 1 1 -155 0.0012212640 1 1 -122 0.0011819022 1 1 -123 0.0011596560 1 1 -173 0.0011288507 1 1 -124 0.0011196148 1 1 -126 0.0010999163 1 1 -128 0.0010911400 1 1 -129 0.0010890244 1 1 -146 0.0010783008 1 1 -130 0.0010751829 1 1 -167 0.0010585201 1 1 -131 0.0010475244 1 1 -132 0.0010471228 1 1 -133 0.0010363011 1 1 -135 0.0010183305 1 1 -136 0.0010152518 1 1 -138 0.0009979338 1 1 -139 0.0009854470 1 1 -140 0.0009449427 1 1 -141 0.0009439848 1 1 -192 0.0009300975 1 1 -144 0.0009043588 1 1 -145 0.0008923953 1 1 -147 0.0008798907 1 1 -150 0.0008646305 1 1 -198 0.0008553574 1 1 -201 0.0008520322 1 1 -152 0.0008277560 1 1 -153 0.0008241527 1 1 -154 0.0008015946 1 1 -170 0.0007958760 1 1 -158 0.0007756951 1 1 -159 0.0007724625 1 1 -160 0.0007708099 1 1 -161 0.0007540851 1 1 -162 0.0007476359 1 1 -180 0.0007444053 1 1 -163 0.0007274333 1 1 -165 0.0007011670 1 1 -166 0.0006897825 1 1 -197 0.0006890637 1 1 -168 0.0006812669 1 1 -176 0.0006333556 0 1 -195 0.0004927698 1 1 +0 0.0463680878 1 1 +1 0.0391439191 1 1 +2 0.0335266518 1 1 +3 0.0318553683 1 1 +4 0.0216635411 1 1 +18 0.0210519227 1 1 +14 0.0206739791 1 1 +5 0.0182095969 1 1 +23 0.0181958498 1 1 +8 0.0169505164 1 1 +24 0.0167605942 1 1 +17 0.0166000739 1 1 +6 0.0164201046 1 1 +7 0.0160856937 1 1 +9 0.0153050047 1 1 +10 0.0139981005 1 1 +13 0.0135998191 1 1 +11 0.0132702620 1 1 +22 0.0132211539 1 1 +12 0.0127169078 1 1 +16 0.0126497515 1 1 +15 0.0122651616 1 1 +30 0.0121237589 1 1 +29 0.0118554136 1 1 +19 0.0113732211 1 1 +39 0.0110561434 1 1 +20 0.0109698887 1 1 +21 0.0107861451 1 1 +31 0.0103469565 1 1 +25 0.0092939917 1 1 +28 0.0091815715 1 1 +34 0.0087700311 1 1 +26 0.0087170796 1 1 +38 0.0085575584 1 1 +33 0.0079226560 1 1 +44 0.0078138061 1 1 +32 0.0075021864 1 1 +40 0.0072273928 1 1 +53 0.0069471751 1 1 +64 0.0065839211 1 1 +35 0.0065742588 1 1 +47 0.0062338367 1 1 +79 0.0061106185 1 1 +37 0.0060770949 1 1 +41 0.0055869137 1 1 +42 0.0054746478 1 1 +55 0.0052835698 1 1 +43 0.0051617313 1 1 +45 0.0050315377 1 1 +46 0.0049791763 1 1 +48 0.0049188216 1 1 +49 0.0049163691 1 1 +51 0.0047761559 1 1 +50 0.0046946655 1 1 +52 0.0043790409 1 1 +57 0.0041968743 1 1 +90 0.0040080302 1 1 +56 0.0038637309 1 1 +58 0.0036963255 1 1 +70 0.0036947486 1 1 +59 0.0036805667 1 1 +60 0.0035963351 1 1 +61 0.0034618254 1 1 +63 0.0034193124 1 1 +104 0.0033819741 1 1 +62 0.0033703372 1 1 +65 0.0032443067 1 1 +66 0.0031583964 1 1 +67 0.0030551941 1 1 +68 0.0030530173 1 1 +69 0.0029211105 1 1 +72 0.0028488136 1 1 +118 0.0028303058 1 1 +74 0.0027857715 1 1 +109 0.0027766105 1 1 +77 0.0026888147 1 1 +120 0.0026794094 1 1 +78 0.0026748129 1 1 +143 0.0026473925 1 1 +80 0.0025811781 1 1 +81 0.0025711725 1 1 +82 0.0025561671 1 1 +83 0.0025013064 1 1 +85 0.0024600725 1 1 +88 0.0024255436 1 1 +91 0.0022399445 1 1 +92 0.0022165505 1 1 +101 0.0021823963 1 1 +94 0.0021619139 1 1 +96 0.0021268245 1 1 +99 0.0020665744 1 1 +100 0.0019938661 1 1 +102 0.0019413134 1 1 +122 0.0019331894 1 1 +113 0.0018610836 1 1 +105 0.0018272640 1 1 +107 0.0017973888 1 1 +139 0.0017909011 1 1 +112 0.0016738720 1 1 +133 0.0016411174 1 1 +115 0.0016307625 1 1 +116 0.0016190690 1 1 +124 0.0016138914 1 1 +123 0.0014776788 1 1 +159 0.0014367844 1 1 +127 0.0014120731 1 1 +128 0.0014103732 1 1 +136 0.0013848384 1 1 +129 0.0013644736 1 1 +131 0.0013572181 1 1 +132 0.0013446395 1 1 +134 0.0012919181 1 1 +135 0.0012907003 1 1 +189 0.0012504976 1 1 +137 0.0012476768 1 1 +138 0.0012466397 1 1 +140 0.0012370013 1 1 +142 0.0012265916 1 1 +145 0.0012087486 1 1 +181 0.0012080467 1 1 +146 0.0011971229 1 1 +148 0.0011388644 1 1 +151 0.0010684559 1 1 +157 0.0009707258 1 1 +158 0.0009689865 1 1 +164 0.0009486279 1 1 +161 0.0009376553 1 1 +241 0.0009323571 0 1 +163 0.0008969517 1 1 +223 0.0008860523 1 1 +214 0.0008767456 1 1 +168 0.0008382003 1 1 +173 0.0007782050 1 1 +207 0.0007671271 1 1 +220 0.0007380965 1 1 +182 0.0006875911 1 1 +27 0.0006750096 1 1 +184 0.0006696548 1 1 +188 0.0006454724 1 1 +226 0.0006382474 1 1 +191 0.0006267047 1 1 +193 0.0006100917 1 1 +195 0.0006038949 1 1 +196 0.0006019849 1 1 +197 0.0006008408 1 1 +199 0.0005746615 1 1 +202 0.0005713762 1 1 +203 0.0005619451 1 1 +205 0.0005589062 1 1 +206 0.0005556910 1 1 +208 0.0005388654 1 1 +209 0.0005321244 1 1 +215 0.0005075217 1 1 +217 0.0004825178 1 1 diff --git a/gwemopt/tests/data/expected_results/coverage_ZTF.dat b/gwemopt/tests/data/expected_results/coverage_ZTF.dat index 4eb72838..dd367a0b 100644 --- a/gwemopt/tests/data/expected_results/coverage_ZTF.dat +++ b/gwemopt/tests/data/expected_results/coverage_ZTF.dat @@ -1,2 +1,2 @@ -247 0.3402372379 1 1 -246 0.3099787805 1 1 +247 0.3402400000 1 1 +246 0.3099800000 1 1 diff --git a/gwemopt/tests/data/expected_results/map.dat b/gwemopt/tests/data/expected_results/map.dat deleted file mode 100644 index 1e1fd3de..00000000 --- a/gwemopt/tests/data/expected_results/map.dat +++ /dev/null @@ -1 +0,0 @@ -11734.95129 12003.05316 diff --git a/gwemopt/tests/data/expected_results/schedule_DECam.dat b/gwemopt/tests/data/expected_results/schedule_DECam.dat index 842da47b..62b0e11a 100644 --- a/gwemopt/tests/data/expected_results/schedule_DECam.dat +++ b/gwemopt/tests/data/expected_results/schedule_DECam.dat @@ -1,10 +1,4 @@ -1000246 13.23930 -24.74750 58710.30114 23.40360 30 0.29237 1.02371 g -1000823 12.39167 -26.30308 58710.30181 23.40360 30 0.15997 1.01795 g -1000381 12.39826 -24.27408 58710.30248 23.40360 30 0.13589 1.02076 g -1002413 11.49121 -23.73125 58710.30315 23.40360 30 0.04742 1.01887 g -1000165 13.23260 -22.71860 58710.30383 23.40360 30 0.04282 1.02529 g -1000246 13.23930 -24.74750 58710.32163 23.40360 30 0.29237 1.00778 r -1000823 12.39167 -26.30308 58710.32230 23.40360 30 0.15997 1.00429 r -1000381 12.39826 -24.27408 58710.32297 23.40360 30 0.13589 1.00722 r -1002413 11.49121 -23.73125 58710.32364 23.40360 30 0.04742 1.00738 r -1000165 13.23260 -22.71860 58710.32431 23.40360 30 0.04282 1.01084 r +1006132 12.63724 -25.57599 58710.30114 23.40360 30 0.39058 1.02043 g +1000381 12.39826 -24.27408 58710.30181 23.40360 30 0.36680 1.02144 g +1006132 12.63724 -25.57599 58710.32163 23.40360 30 0.39058 1.00577 r +1000381 12.39826 -24.27408 58710.32230 23.40360 30 0.36680 1.00744 r diff --git a/gwemopt/tests/data/expected_results/schedule_KPED.dat b/gwemopt/tests/data/expected_results/schedule_KPED.dat index 2f4ed454..ee780197 100644 --- a/gwemopt/tests/data/expected_results/schedule_KPED.dat +++ b/gwemopt/tests/data/expected_results/schedule_KPED.dat @@ -1,355 +1,307 @@ -96 11.01104 -22.36956 58710.33095 22.05000 30 0.00201 2.50529 g -96 11.01104 -22.36956 58710.33142 22.05000 30 0.00201 2.49607 r -195 10.64846 -22.68770 58710.33188 22.05000 30 0.00049 2.49266 g -195 10.64846 -22.68770 58710.33234 22.05000 30 0.00049 2.48368 r -176 10.28554 -23.29063 58710.33281 22.05000 30 0.00063 2.50390 g -146 10.31580 -23.40952 58710.33327 22.05000 30 0.00108 2.50630 g -146 10.31580 -23.40952 58710.33373 22.05000 30 0.00108 2.49743 r -101 11.50037 -22.81176 58710.33419 22.05000 30 0.00171 2.50351 g -101 11.50037 -22.81176 58710.33466 22.05000 30 0.00171 2.49445 r -68 11.21060 -23.15903 58710.33512 22.05000 30 0.00290 2.49797 g -68 11.21060 -23.15903 58710.33558 22.05000 30 0.00290 2.48910 r -49 11.86747 -23.02301 58710.33605 22.05000 30 0.00476 2.50426 g -49 11.86747 -23.02301 58710.33651 22.05000 30 0.00476 2.49527 r -73 10.78818 -23.87461 58710.33697 22.05000 30 0.00336 2.49891 g -73 10.78818 -23.87461 58710.33744 22.05000 30 0.00336 2.49029 r -47 10.83845 -24.08253 58710.33790 22.05000 30 0.00532 2.50133 g -48 11.47175 -23.78052 58710.33836 22.05000 30 0.00834 2.50066 g -48 11.47175 -23.78052 58710.33882 22.05000 30 0.00834 2.49198 r -1 12.18036 -23.56169 58710.33929 22.05000 30 0.04022 2.50237 g -1 12.18036 -23.56169 58710.33975 22.05000 30 0.04022 2.49360 r -36 11.05496 -24.32757 58710.34021 22.05000 30 0.01530 2.48998 g -5 12.25617 -23.81132 58710.34068 22.05000 30 0.03568 2.50050 g -5 12.25617 -23.81132 58710.34114 22.05000 30 0.03568 2.49184 r -0 11.78136 -24.37065 58710.34160 22.05000 30 0.04654 2.50493 g -0 11.78136 -24.37065 58710.34206 22.05000 30 0.04654 2.49644 r -36 11.05496 -24.32757 58710.34253 22.05000 30 0.01530 2.44908 r -23 12.72686 -23.63189 58710.34299 22.05000 30 0.01482 2.46728 g -23 12.72686 -23.63189 58710.34345 22.05000 30 0.01482 2.45890 r -16 12.54044 -23.28005 58710.34392 22.05000 30 0.01279 2.41425 g -16 12.54044 -23.28005 58710.34438 22.05000 30 0.01279 2.40630 r -43 11.27932 -25.02976 58710.34484 22.05000 30 0.01137 2.47682 g -43 11.27932 -25.02976 58710.34531 22.05000 30 0.01137 2.46888 r -10 13.80442 -24.04403 58710.34577 22.05000 30 0.01479 2.50456 g -3 12.09108 -25.12681 58710.34623 22.05000 30 0.03125 2.49989 g -3 12.09108 -25.12681 58710.34669 22.05000 30 0.03125 2.49175 r -9 11.87062 -25.44065 58710.34716 22.05000 30 0.01658 2.49932 g -9 11.87062 -25.44065 58710.34762 22.05000 30 0.01658 2.49131 r -10 13.80442 -24.04403 58710.34808 22.05000 30 0.01479 2.46239 r -11 12.10315 -25.59571 58710.34855 22.05000 30 0.01473 2.49950 g -11 12.10315 -25.59571 58710.34901 22.05000 30 0.01473 2.49155 r -13 12.42386 -25.05081 58710.34947 22.05000 30 0.01412 2.45376 g -13 12.42386 -25.05081 58710.34994 22.05000 30 0.01412 2.44607 r -15 12.87443 -24.64249 58710.35040 22.05000 30 0.01295 2.42644 g -15 12.87443 -24.64249 58710.35086 22.05000 30 0.01295 2.41887 r -24 12.55164 -25.97516 58710.35132 22.05000 30 0.01818 2.50540 g -24 12.55164 -25.97516 58710.35179 22.05000 30 0.01818 2.49755 r -2 13.17326 -25.73385 58710.35225 22.05000 30 0.03255 2.49847 g -2 13.17326 -25.73385 58710.35271 22.05000 30 0.03255 2.49059 r -4 12.79385 -25.95417 58710.35318 22.05000 30 0.02427 2.48377 g -4 12.79385 -25.95417 58710.35364 22.05000 30 0.02427 2.47613 r -8 12.82757 -26.34528 58710.35410 22.05000 30 0.02585 2.50317 g -8 12.82757 -26.34528 58710.35456 22.05000 30 0.02585 2.49550 r -6 12.32009 -26.21918 58710.35503 22.05000 30 0.02249 2.45487 g -6 12.32009 -26.21918 58710.35549 22.05000 30 0.02249 2.44765 r -12 12.33203 -26.47640 58710.35595 22.05000 30 0.02076 2.46231 g -12 12.33203 -26.47640 58710.35642 22.05000 30 0.02076 2.45512 r -7 13.80361 -26.33097 58710.35688 22.05000 30 0.02222 2.50070 g -7 13.80361 -26.33097 58710.35734 22.05000 30 0.02222 2.49305 r -14 13.60144 -25.46405 58710.35781 22.05000 30 0.01400 2.40578 g -14 13.60144 -25.46405 58710.35827 22.05000 30 0.01400 2.39875 r -17 13.49057 -24.54314 58710.35873 22.05000 30 0.01239 2.31794 g -17 13.49057 -24.54314 58710.35919 22.05000 30 0.01239 2.31144 r -41 12.38561 -26.53859 58710.35966 22.05000 30 0.01194 2.41402 g -41 12.38561 -26.53859 58710.36012 22.05000 30 0.01194 2.40736 r -18 12.82817 -26.16806 58710.36058 22.05000 30 0.01193 2.38873 g -18 12.82817 -26.16806 58710.36105 22.05000 30 0.01193 2.38217 r -19 12.44180 -26.44304 58710.36151 22.05000 30 0.01147 2.38239 g -19 12.44180 -26.44304 58710.36197 22.05000 30 0.01147 2.37602 r -20 12.22897 -25.06947 58710.36244 22.05000 30 0.01062 2.26026 g -20 12.22897 -25.06947 58710.36290 22.05000 30 0.01062 2.25455 r -21 13.78813 -25.45570 58710.36336 22.05000 30 0.01056 2.33239 g -21 13.78813 -25.45570 58710.36382 22.05000 30 0.01056 2.32611 r -22 13.06404 -24.69873 58710.36429 22.05000 30 0.01028 2.24019 g -22 13.06404 -24.69873 58710.36475 22.05000 30 0.01028 2.23453 r -25 11.54340 -24.65019 58710.36521 22.05000 30 0.00969 2.17730 g -25 11.54340 -24.65019 58710.36568 22.05000 30 0.00969 2.17225 r -40 13.19239 -22.97502 58710.36614 22.05000 30 0.00943 2.11205 g -40 13.19239 -22.97502 58710.36660 22.05000 30 0.00943 2.10703 r -37 13.26783 -26.17079 58710.36706 22.05000 30 0.00915 2.31817 g -37 13.26783 -26.17079 58710.36753 22.05000 30 0.00915 2.31233 r -26 12.83825 -26.98945 58710.36799 22.05000 30 0.00908 2.35392 g -26 12.83825 -26.98945 58710.36845 22.05000 30 0.00908 2.34808 r -27 12.37098 -25.99136 58710.36892 22.05000 30 0.00903 2.25290 g -27 12.37098 -25.99136 58710.36938 22.05000 30 0.00903 2.24765 r -28 12.64379 -23.61855 58710.36984 22.05000 30 0.00900 2.09606 g -28 12.64379 -23.61855 58710.37031 22.05000 30 0.00900 2.09143 r -29 11.77189 -24.23870 58710.37077 22.05000 30 0.00883 2.10069 g -29 11.77189 -24.23870 58710.37123 22.05000 30 0.00883 2.09625 r -30 12.47310 -24.70720 58710.37169 22.05000 30 0.00874 2.13987 g -30 12.47310 -24.70720 58710.37216 22.05000 30 0.00874 2.13522 r -33 12.66614 -26.81339 58710.37262 22.05000 30 0.00865 2.27970 g -33 12.66614 -26.81339 58710.37308 22.05000 30 0.00865 2.27454 r -45 13.57097 -23.55266 58710.37355 22.05000 30 0.00848 2.08088 g -45 13.57097 -23.55266 58710.37401 22.05000 30 0.00848 2.07638 r -34 13.35898 -26.59982 58710.37447 22.05000 30 0.00834 2.26486 g -34 13.35898 -26.59982 58710.37494 22.05000 30 0.00834 2.25975 r -35 13.65696 -25.06710 58710.37540 22.05000 30 0.00827 2.15867 g -35 13.65696 -25.06710 58710.37586 22.05000 30 0.00827 2.15397 r -46 11.45542 -25.53103 58710.37632 22.05000 30 0.00777 2.12132 g -46 11.45542 -25.53103 58710.37679 22.05000 30 0.00777 2.11720 r -51 12.81521 -23.29441 58710.37725 22.05000 30 0.00759 2.01311 g -51 12.81521 -23.29441 58710.37771 22.05000 30 0.00759 2.00919 r -38 13.70472 -26.37126 58710.37818 22.05000 30 0.00751 2.21926 g -38 13.70472 -26.37126 58710.37864 22.05000 30 0.00751 2.21452 r -39 12.17492 -23.36863 58710.37910 22.05000 30 0.00718 1.98746 g -39 12.17492 -23.36863 58710.37956 22.05000 30 0.00718 1.98382 r -56 11.65284 -26.16652 58710.38003 22.05000 30 0.00692 2.13443 g -56 11.65284 -26.16652 58710.38049 22.05000 30 0.00692 2.13047 r -42 13.26927 -24.70440 58710.38095 22.05000 30 0.00647 2.07371 g -42 13.26927 -24.70440 58710.38142 22.05000 30 0.00647 2.06973 r -72 10.75701 -24.65020 58710.38188 22.05000 30 0.00598 2.00812 g -72 10.75701 -24.65020 58710.38234 22.05000 30 0.00598 2.00479 r -44 12.88104 -26.07721 58710.38281 22.05000 30 0.00590 2.13404 g -44 12.88104 -26.07721 58710.38327 22.05000 30 0.00590 2.13004 r -63 12.68625 -23.33883 58710.38373 22.05000 30 0.00559 1.96138 g -63 12.68625 -23.33883 58710.38419 22.05000 30 0.00559 1.95800 r -47 10.83845 -24.08253 58710.38466 22.05000 30 0.00532 1.95904 r -54 12.72130 -23.08411 58710.38512 22.05000 30 0.00492 1.93869 g -54 12.72130 -23.08411 58710.38558 22.05000 30 0.00492 1.93544 r -50 12.66216 -26.40645 58710.38605 22.05000 30 0.00459 2.12309 g -50 12.66216 -26.40645 58710.38651 22.05000 30 0.00459 2.11936 r -52 13.52282 -23.19464 58710.38697 22.05000 30 0.00435 1.94708 g -52 13.52282 -23.19464 58710.38744 22.05000 30 0.00435 1.94379 r -53 11.73157 -24.81370 58710.38790 22.05000 30 0.00432 1.99423 g -53 11.73157 -24.81370 58710.38836 22.05000 30 0.00432 1.99113 r -55 12.95494 -27.07109 58710.38882 22.05000 30 0.00405 2.15097 g -55 12.95494 -27.07109 58710.38929 22.05000 30 0.00405 2.14729 r -117 13.96266 -23.29280 58710.38975 22.05000 30 0.00371 1.94125 g -117 13.96266 -23.29280 58710.39021 22.05000 30 0.00371 1.93806 r -64 13.90494 -23.82655 58710.39068 22.05000 30 0.00354 1.96204 g -64 13.90494 -23.82655 58710.39114 22.05000 30 0.00354 1.95885 r -57 11.55507 -25.95019 58710.39160 22.05000 30 0.00338 2.03240 g -57 11.55507 -25.95019 58710.39206 22.05000 30 0.00338 2.02944 r -58 13.32160 -25.32121 58710.39253 22.05000 30 0.00337 2.02178 g -58 13.32160 -25.32121 58710.39299 22.05000 30 0.00337 2.01862 r -81 14.47873 -24.43375 58710.39345 22.05000 30 0.00334 1.98702 g -81 14.47873 -24.43375 58710.39392 22.05000 30 0.00334 1.98382 r -60 12.31113 -23.85855 58710.39438 22.05000 30 0.00326 1.91247 g -60 12.31113 -23.85855 58710.39484 22.05000 30 0.00326 1.90983 r -93 11.34162 -25.53855 58710.39531 22.05000 30 0.00323 1.98259 g -93 11.34162 -25.53855 58710.39577 22.05000 30 0.00323 1.97999 r -62 11.27290 -25.76671 58710.39623 22.05000 30 0.00315 1.98930 g -62 11.27290 -25.76671 58710.39669 22.05000 30 0.00315 1.98675 r -75 13.10716 -22.73325 58710.39716 22.05000 30 0.00313 1.85346 g -75 13.10716 -22.73325 58710.39762 22.05000 30 0.00313 1.85099 r -82 13.43679 -27.60508 58710.39808 22.05000 30 0.00312 2.12800 g -82 13.43679 -27.60508 58710.39855 22.05000 30 0.00312 2.12488 r -65 12.80834 -26.46112 58710.39901 22.05000 30 0.00301 2.03860 g -65 12.80834 -26.46112 58710.39947 22.05000 30 0.00301 2.03585 r -66 14.01097 -27.50790 58710.39994 22.05000 30 0.00291 2.11982 g -66 14.01097 -27.50790 58710.40040 22.05000 30 0.00291 2.11674 r -67 10.88859 -23.95291 58710.40086 22.05000 30 0.00290 1.86477 g -67 10.88859 -23.95291 58710.40132 22.05000 30 0.00290 1.86276 r -69 14.08512 -24.34858 58710.40179 22.05000 30 0.00284 1.92398 g -69 14.08512 -24.34858 58710.40225 22.05000 30 0.00284 1.92145 r -76 13.62431 -23.59099 58710.40271 22.05000 30 0.00272 1.87377 g -76 13.62431 -23.59099 58710.40318 22.05000 30 0.00272 1.87146 r -88 13.50573 -24.13203 58710.40364 22.05000 30 0.00266 1.89467 g -88 13.50573 -24.13203 58710.40410 22.05000 30 0.00266 1.89238 r -112 11.18174 -23.85487 58710.40456 22.05000 30 0.00263 1.84788 g -112 11.18174 -23.85487 58710.40503 22.05000 30 0.00263 1.84604 r -71 11.18312 -24.82782 58710.40549 22.05000 30 0.00261 1.89285 g -71 11.18312 -24.82782 58710.40595 22.05000 30 0.00261 1.89098 r -74 13.14967 -26.75093 58710.40642 22.05000 30 0.00257 2.02026 g -74 13.14967 -26.75093 58710.40688 22.05000 30 0.00257 2.01792 r -78 14.24495 -23.86980 58710.40734 22.05000 30 0.00249 1.87346 g -78 14.24495 -23.86980 58710.40781 22.05000 30 0.00249 1.87129 r -77 13.76804 -26.77015 58710.40827 22.05000 30 0.00244 2.02071 g -77 13.76804 -26.77015 58710.40873 22.05000 30 0.00244 2.01838 r -79 12.24555 -24.51614 58710.40919 22.05000 30 0.00241 1.87394 g -79 12.24555 -24.51614 58710.40966 22.05000 30 0.00241 1.87214 r -86 11.49764 -26.50492 58710.41012 22.05000 30 0.00239 1.96872 g -86 11.49764 -26.50492 58710.41058 22.05000 30 0.00239 1.96693 r -80 14.44171 -26.34829 58710.41105 22.05000 30 0.00238 1.99138 g -80 14.44171 -26.34829 58710.41151 22.05000 30 0.00238 1.98915 r -115 13.06562 -27.34956 58710.41197 22.05000 30 0.00230 2.02884 g -115 13.06562 -27.34956 58710.41244 22.05000 30 0.00230 2.02681 r -85 10.85066 -25.29308 58710.41290 22.05000 30 0.00213 1.88687 g -85 10.85066 -25.29308 58710.41336 22.05000 30 0.00213 1.88543 r -116 14.32338 -24.00069 58710.41382 22.05000 30 0.00211 1.85280 g -116 14.32338 -24.00069 58710.41429 22.05000 30 0.00211 1.85098 r -87 10.54359 -23.59152 58710.41475 22.05000 30 0.00208 1.79604 g -87 10.54359 -23.59152 58710.41521 22.05000 30 0.00208 1.79483 r -89 13.37070 -25.03676 58710.41568 22.05000 30 0.00194 1.88784 g -89 13.37070 -25.03676 58710.41614 22.05000 30 0.00194 1.88618 r -90 12.26253 -26.61864 58710.41660 22.05000 30 0.00194 1.95931 g -90 12.26253 -26.61864 58710.41706 22.05000 30 0.00194 1.95776 r -91 14.80112 -26.47570 58710.41753 22.05000 30 0.00186 1.97411 g -91 14.80112 -26.47570 58710.41799 22.05000 30 0.00186 1.97221 r -92 11.84997 -26.81732 58710.41845 22.05000 30 0.00184 1.96098 g -92 11.84997 -26.81732 58710.41892 22.05000 30 0.00184 1.95959 r -94 10.07006 -24.45453 58710.41938 22.05000 30 0.00181 1.82261 g -94 10.07006 -24.45453 58710.41984 22.05000 30 0.00181 1.82165 r -95 11.72066 -23.03338 58710.42031 22.05000 30 0.00179 1.76519 g -95 11.72066 -23.03338 58710.42077 22.05000 30 0.00179 1.76413 r -108 13.83209 -23.53129 58710.42123 22.05000 30 0.00174 1.80054 g -108 13.83209 -23.53129 58710.42169 22.05000 30 0.00174 1.79921 r -84 22.14535 -32.29647 58710.42216 22.05000 30 0.00330 2.50255 g -84 22.14535 -32.29647 58710.42262 22.05000 30 0.00330 2.49843 r -98 13.68088 -23.87741 58710.42308 22.05000 30 0.00174 1.81031 g -98 13.68088 -23.87741 58710.42355 22.05000 30 0.00174 1.80907 r -100 13.72701 -23.18742 58710.42401 22.05000 30 0.00171 1.77657 g -100 13.72701 -23.18742 58710.42447 22.05000 30 0.00171 1.77541 r -103 10.58759 -23.68735 58710.42494 22.05000 30 0.00164 1.77923 g -103 10.58759 -23.68735 58710.42540 22.05000 30 0.00164 1.77850 r -83 22.93853 -32.41562 58710.42586 22.05000 30 0.00222 2.50070 g -83 22.93853 -32.41562 58710.42632 22.05000 30 0.00222 2.49671 r -102 22.75958 -32.51090 58710.42679 22.05000 30 0.00168 2.49767 g -102 22.75958 -32.51090 58710.42725 22.05000 30 0.00168 2.49381 r -104 21.91262 -32.11166 58710.42771 22.05000 30 0.00163 2.43443 g -104 21.91262 -32.11166 58710.42818 22.05000 30 0.00163 2.43098 r -107 22.23838 -32.82368 58710.42864 22.05000 30 0.00300 2.50093 g -107 22.23838 -32.82368 58710.42910 22.05000 30 0.00300 2.49733 r -99 23.38010 -32.61306 58710.42956 22.05000 30 0.00181 2.49868 g -99 23.38010 -32.61306 58710.43003 22.05000 30 0.00181 2.49490 r -70 22.53305 -32.91646 58710.43049 22.05000 30 0.00471 2.50188 g -70 22.53305 -32.91646 58710.43095 22.05000 30 0.00471 2.49836 r -61 23.14750 -32.75926 58710.43142 22.05000 30 0.00374 2.49276 g -61 23.14750 -32.75926 58710.43188 22.05000 30 0.00374 2.48919 r -31 23.51791 -32.80877 58710.43234 22.05000 30 0.00266 2.49831 g -32 23.61699 -32.84254 58710.43281 22.05000 30 0.01038 2.50010 g -32 23.61699 -32.84254 58710.43327 22.05000 30 0.01038 2.49652 r -31 23.51791 -32.80877 58710.43373 22.05000 30 0.00266 2.48769 r -97 23.64540 -32.78864 58710.43419 22.05000 30 0.00244 2.48495 g -97 23.64540 -32.78864 58710.43466 22.05000 30 0.00244 2.48152 r -59 23.87743 -32.97061 58710.43512 22.05000 30 0.00337 2.50032 g -59 23.87743 -32.97061 58710.43558 22.05000 30 0.00337 2.49687 r -113 23.76271 -33.00979 58710.43605 22.05000 30 0.00294 2.49493 g -113 23.76271 -33.00979 58710.43651 22.05000 30 0.00294 2.49160 r -111 23.92573 -32.86950 58710.43697 22.05000 30 0.00230 2.47811 g -118 23.92815 -33.13617 58710.43744 22.05000 30 0.00253 2.50048 g -118 23.92815 -33.13617 58710.43790 22.05000 30 0.00253 2.49721 r -111 23.92573 -32.86950 58710.43836 22.05000 30 0.00230 2.46847 r -148 24.74489 -32.96067 58710.43882 22.05000 30 0.00228 2.48984 g -148 24.74489 -32.96067 58710.43929 22.05000 30 0.00228 2.48652 r -125 23.32769 -32.89307 58710.43975 22.05000 30 0.00208 2.45088 g -125 23.32769 -32.89307 58710.44021 22.05000 30 0.00208 2.44804 r -142 24.10953 -33.19223 58710.44068 22.05000 30 0.00164 2.48718 g -142 24.10953 -33.19223 58710.44114 22.05000 30 0.00164 2.48417 r -105 11.25806 -24.68171 58710.44160 22.05000 30 0.00163 1.81236 g -105 11.25806 -24.68171 58710.44206 22.05000 30 0.00163 1.81229 r -106 13.28581 -22.96882 58710.44253 22.05000 30 0.00162 1.73658 g -106 13.28581 -22.96882 58710.44299 22.05000 30 0.00162 1.73632 r -127 11.11274 -25.38647 58710.44345 22.05000 30 0.00157 1.84668 g -127 11.11274 -25.38647 58710.44392 22.05000 30 0.00157 1.84673 r -114 11.56655 -27.10378 58710.44438 22.05000 30 0.00150 1.93760 g -114 11.56655 -27.10378 58710.44484 22.05000 30 0.00150 1.93763 r -137 22.34289 -32.84935 58710.44531 22.05000 30 0.00147 2.40210 g -137 22.34289 -32.84935 58710.44577 22.05000 30 0.00147 2.39997 r -109 10.47218 -24.08844 58710.44623 22.05000 30 0.00146 1.78575 g -109 10.47218 -24.08844 58710.44669 22.05000 30 0.00146 1.78599 r -110 22.29279 -31.57823 58710.44716 22.05000 30 0.00145 2.28509 g -110 22.29279 -31.57823 58710.44762 22.05000 30 0.00145 2.28328 r -151 23.34156 -32.79363 58710.44808 22.05000 30 0.00143 2.39708 g -151 23.34156 -32.79363 58710.44855 22.05000 30 0.00143 2.39497 r -134 24.33415 -33.33085 58710.44901 22.05000 30 0.00141 2.45544 g -134 24.33415 -33.33085 58710.44947 22.05000 30 0.00141 2.45309 r -164 13.66758 -22.90186 58710.44994 22.05000 30 0.00139 1.73188 g -164 13.66758 -22.90186 58710.45040 22.05000 30 0.00139 1.73190 r -157 22.44899 -33.02676 58710.45086 22.05000 30 0.00136 2.39602 g -157 22.44899 -33.02676 58710.45132 22.05000 30 0.00136 2.39429 r -120 23.42497 -32.22377 58710.45179 22.05000 30 0.00134 2.33271 g -120 23.42497 -32.22377 58710.45225 22.05000 30 0.00134 2.33094 r -143 13.22368 -22.43724 58710.45271 22.05000 30 0.00131 1.71307 g -143 13.22368 -22.43724 58710.45318 22.05000 30 0.00131 1.71326 r -156 23.34954 -32.95530 58710.45364 22.05000 30 0.00128 2.38860 g -156 23.34954 -32.95530 58710.45410 22.05000 30 0.00128 2.38690 r -119 11.78522 -26.42648 58710.45456 22.05000 30 0.00127 1.90585 g -119 11.78522 -26.42648 58710.45503 22.05000 30 0.00127 1.90636 r -149 24.19522 -33.39419 58710.45549 22.05000 30 0.00126 2.43051 g -149 24.19522 -33.39419 58710.45595 22.05000 30 0.00126 2.42873 r -121 11.44294 -25.92019 58710.45642 22.05000 30 0.00123 1.88245 g -121 11.44294 -25.92019 58710.45688 22.05000 30 0.00123 1.88309 r -155 23.08323 -31.79593 58710.45734 22.05000 30 0.00122 2.27592 g -155 23.08323 -31.79593 58710.45781 22.05000 30 0.00122 2.27467 r -122 14.60634 -25.19848 58710.45827 22.05000 30 0.00118 1.83913 g -122 14.60634 -25.19848 58710.45873 22.05000 30 0.00118 1.83942 r -123 22.15717 -32.34258 58710.45919 22.05000 30 0.00116 2.30991 g -123 22.15717 -32.34258 58710.45966 22.05000 30 0.00116 2.30894 r -173 24.39115 -32.22932 58710.46012 22.05000 30 0.00113 2.31397 g -173 24.39115 -32.22932 58710.46058 22.05000 30 0.00113 2.31262 r -124 13.35392 -25.82649 58710.46105 22.05000 30 0.00112 1.87674 g -124 13.35392 -25.82649 58710.46151 22.05000 30 0.00112 1.87735 r -126 11.56508 -25.36933 58710.46197 22.05000 30 0.00110 1.86261 g -126 11.56508 -25.36933 58710.46244 22.05000 30 0.00110 1.86349 r -128 22.24245 -32.06395 58710.46290 22.05000 30 0.00109 2.28051 g -128 22.24245 -32.06395 58710.46336 22.05000 30 0.00109 2.27979 r -129 12.76531 -27.68500 58710.46382 22.05000 30 0.00109 1.98667 g -129 12.76531 -27.68500 58710.46429 22.05000 30 0.00109 1.98758 r -130 11.04536 -23.09401 58710.46475 22.05000 30 0.00108 1.76335 g -130 11.04536 -23.09401 58710.46521 22.05000 30 0.00108 1.76435 r -167 23.63711 -31.78614 58710.46568 22.05000 30 0.00106 2.26003 g -167 23.63711 -31.78614 58710.46614 22.05000 30 0.00106 2.25925 r -131 10.47123 -23.38096 58710.46660 22.05000 30 0.00105 1.78415 g -131 10.47123 -23.38096 58710.46706 22.05000 30 0.00105 1.78533 r -132 10.93788 -25.49946 58710.46753 22.05000 30 0.00105 1.88602 g -132 10.93788 -25.49946 58710.46799 22.05000 30 0.00105 1.88729 r -133 10.16891 -24.26920 58710.46845 22.05000 30 0.00104 1.83298 g -133 10.16891 -24.26920 58710.46892 22.05000 30 0.00104 1.83433 r -135 10.44668 -24.38789 58710.46938 22.05000 30 0.00102 1.83920 g -135 10.44668 -24.38789 58710.46984 22.05000 30 0.00102 1.84057 r -136 12.67995 -22.88828 58710.47031 22.05000 30 0.00102 1.75657 g -136 12.67995 -22.88828 58710.47077 22.05000 30 0.00102 1.75761 r -138 13.26590 -22.74250 58710.47123 22.05000 30 0.00100 1.74873 g -138 13.26590 -22.74250 58710.47169 22.05000 30 0.00100 1.74973 r -139 11.09615 -25.79841 58710.47216 22.05000 30 0.00099 1.91428 g -139 11.09615 -25.79841 58710.47262 22.05000 30 0.00099 1.91579 r -140 11.47179 -27.00770 58710.47308 22.05000 30 0.00094 1.98138 g -140 11.47179 -27.00770 58710.47355 22.05000 30 0.00094 1.98298 r -141 23.57875 -33.59670 58710.47401 22.05000 30 0.00094 2.40488 g -141 23.57875 -33.59670 58710.47447 22.05000 30 0.00094 2.40465 r -192 23.54662 -33.25808 58710.47494 22.05000 30 0.00093 2.37384 g -192 23.54662 -33.25808 58710.47540 22.05000 30 0.00093 2.37369 r -144 11.70062 -25.59481 58710.47586 22.05000 30 0.00090 1.91024 g -144 11.70062 -25.59481 58710.47632 22.05000 30 0.00090 1.91185 r -145 13.54312 -22.52016 58710.47679 22.05000 30 0.00089 1.75041 g -145 13.54312 -22.52016 58710.47725 22.05000 30 0.00089 1.75162 r -147 11.79820 -23.75940 58710.47771 22.05000 30 0.00088 1.82348 g -147 11.79820 -23.75940 58710.47818 22.05000 30 0.00088 1.82505 r -150 12.86661 -27.64000 58710.47864 22.05000 30 0.00086 2.02517 g -150 12.86661 -27.64000 58710.47910 22.05000 30 0.00086 2.02693 r -198 23.35247 -31.86168 58710.47956 22.05000 30 0.00086 2.25656 g -198 23.35247 -31.86168 58710.48003 22.05000 30 0.00086 2.25677 r -201 11.33394 -26.95744 58710.48049 22.05000 30 0.00085 2.00907 g -201 11.33394 -26.95744 58710.48095 22.05000 30 0.00085 2.01115 r -152 11.33884 -25.61994 58710.48142 22.05000 30 0.00083 1.93695 g -152 11.33884 -25.61994 58710.48188 22.05000 30 0.00083 1.93895 r -153 23.53386 -32.60557 58710.48234 22.05000 30 0.00082 2.31894 g -153 23.53386 -32.60557 58710.48281 22.05000 30 0.00082 2.31932 r -154 11.10014 -25.47480 58710.48327 22.05000 30 0.00080 1.94019 g -154 11.10014 -25.47480 58710.48373 22.05000 30 0.00080 1.94232 r -170 23.18752 -33.12454 58710.48419 22.05000 30 0.00080 2.36664 g -170 23.18752 -33.12454 58710.48466 22.05000 30 0.00080 2.36725 r -158 23.18877 -32.61148 58710.48512 22.05000 30 0.00078 2.32344 g -158 23.18877 -32.61148 58710.48558 22.05000 30 0.00078 2.32409 r -159 14.16072 -26.64081 58710.48605 22.05000 30 0.00077 1.98093 g -159 14.16072 -26.64081 58710.48651 22.05000 30 0.00077 1.98283 r -160 14.33254 -26.82652 58710.48697 22.05000 30 0.00077 1.99358 g -160 14.33254 -26.82652 58710.48744 22.05000 30 0.00077 1.99553 r -161 13.00702 -22.53056 58710.48790 22.05000 30 0.00075 1.79177 g -161 13.00702 -22.53056 58710.48836 22.05000 30 0.00075 1.79360 r -162 22.05036 -32.16487 58710.48882 22.05000 30 0.00075 2.29885 g -162 22.05036 -32.16487 58710.48929 22.05000 30 0.00075 2.29996 r -180 23.75118 -33.34174 58710.48975 22.05000 30 0.00074 2.39246 g -180 23.75118 -33.34174 58710.49021 22.05000 30 0.00074 2.39337 r -163 23.51158 -32.70571 58710.49068 22.05000 30 0.00073 2.33959 g -163 23.51158 -32.70571 58710.49114 22.05000 30 0.00073 2.34058 r -165 11.12295 -23.29709 58710.49160 22.05000 30 0.00070 1.86819 g -165 11.12295 -23.29709 58710.49206 22.05000 30 0.00070 1.87060 r -166 22.38228 -32.95269 58710.49253 22.05000 30 0.00069 2.37361 g -166 22.38228 -32.95269 58710.49299 22.05000 30 0.00069 2.37499 r -197 23.50874 -32.29989 58710.49345 22.05000 30 0.00069 2.31159 g -197 23.50874 -32.29989 58710.49392 22.05000 30 0.00069 2.31276 r -168 13.10652 -27.07928 58710.49438 22.05000 30 0.00068 2.06241 g -168 13.10652 -27.07928 58710.49484 22.05000 30 0.00068 2.06511 r +0 319.01898 18.86289 58710.13142 22.05000 30 0.04637 1.75491 g +0 319.01898 18.86289 58710.13188 22.05000 30 0.04637 1.74869 r +1 318.72003 18.95082 58710.13234 22.05000 30 0.03914 1.72995 g +1 318.72003 18.95082 58710.13281 22.05000 30 0.03914 1.72394 r +2 329.89838 18.77292 58710.13327 22.05000 30 0.03353 2.27524 g +2 329.89838 18.77292 58710.13373 22.05000 30 0.03353 2.26389 r +3 331.01535 18.88890 58710.13419 22.05000 30 0.03186 2.32557 g +3 331.01535 18.88890 58710.13466 22.05000 30 0.03186 2.31367 r +18 224.18008 18.80293 58710.13512 22.05000 30 0.02105 1.21798 g +18 224.18008 18.80293 58710.13558 22.05000 30 0.02105 1.22002 r +14 224.96100 18.86283 58710.13605 22.05000 30 0.02067 1.21215 g +14 224.96100 18.86283 58710.13651 22.05000 30 0.02067 1.21416 r +5 282.39609 18.81629 58710.13697 22.05000 30 0.01821 1.09017 g +5 282.39609 18.81629 58710.13744 22.05000 30 0.01821 1.08914 r +23 212.78592 18.89713 58710.13790 22.05000 30 0.01820 1.41779 g +23 212.78592 18.89713 58710.13836 22.05000 30 0.01820 1.42132 r +8 222.75145 18.82670 58710.13882 22.05000 30 0.01695 1.25364 g +8 222.75145 18.82670 58710.13929 22.05000 30 0.01695 1.25594 r +4 335.19968 18.75072 58710.13975 22.05000 30 0.02166 2.50028 g +4 335.19968 18.75072 58710.14021 22.05000 30 0.02166 2.48635 r +24 221.25914 18.80288 58710.14068 22.05000 30 0.01676 1.28501 g +24 221.25914 18.80288 58710.14114 22.05000 30 0.01676 1.28755 r +17 323.91785 18.76300 58710.14160 22.05000 30 0.01660 1.80424 g +17 323.91785 18.76300 58710.14206 22.05000 30 0.01660 1.79758 r +7 336.23032 18.74458 58710.14253 22.05000 30 0.01609 2.50291 g +7 336.23032 18.74458 58710.14299 22.05000 30 0.01609 2.48896 r +9 319.72244 18.75153 58710.14345 22.05000 30 0.01531 1.63246 g +9 319.72244 18.75153 58710.14392 22.05000 30 0.01531 1.62728 r +6 336.95325 18.80861 58710.14438 22.05000 30 0.01642 2.50454 g +6 336.95325 18.80861 58710.14484 22.05000 30 0.01642 2.49057 r +10 323.87970 18.85852 58710.14531 22.05000 30 0.01400 1.74928 g +10 323.87970 18.85852 58710.14577 22.05000 30 0.01400 1.74310 r +13 255.07411 18.82923 58710.14623 22.05000 30 0.01360 1.03749 g +13 255.07411 18.82923 58710.14669 22.05000 30 0.01360 1.03789 r +11 255.45082 18.75540 58710.14716 22.05000 30 0.01327 1.03771 g +11 255.45082 18.75540 58710.14762 22.05000 30 0.01327 1.03810 r +22 320.04260 18.88901 58710.14808 22.05000 30 0.01322 1.58963 g +22 320.04260 18.88901 58710.14855 22.05000 30 0.01322 1.58479 r +12 320.75732 18.90317 58710.14901 22.05000 30 0.01272 1.60057 g +12 320.75732 18.90317 58710.14947 22.05000 30 0.01272 1.59565 r +16 223.54378 18.83519 58710.14994 22.05000 30 0.01265 1.30109 g +16 223.54378 18.83519 58710.15040 22.05000 30 0.01265 1.30375 r +15 321.89807 18.85492 58710.15086 22.05000 30 0.01227 1.61547 g +15 321.89807 18.85492 58710.15132 22.05000 30 0.01227 1.61042 r +30 327.17728 18.74647 58710.15179 22.05000 30 0.01212 1.78806 g +30 327.17728 18.74647 58710.15225 22.05000 30 0.01212 1.78154 r +29 268.83725 18.92122 58710.15271 22.05000 30 0.01186 1.02719 g +29 268.83725 18.92122 58710.15318 22.05000 30 0.01186 1.02709 r +39 336.96121 18.97250 58710.15364 22.05000 30 0.01106 2.24846 g +39 336.96121 18.97250 58710.15410 22.05000 30 0.01106 2.23741 r +20 326.63534 18.81716 58710.15456 22.05000 30 0.01097 1.72860 g +20 326.63534 18.81716 58710.15503 22.05000 30 0.01097 1.72261 r +21 334.27087 18.79971 58710.15549 22.05000 30 0.01079 2.05234 g +21 334.27087 18.79971 58710.15595 22.05000 30 0.01079 2.04333 r +25 221.93694 18.89068 58710.15642 22.05000 30 0.00929 1.36877 g +25 221.93694 18.89068 58710.15688 22.05000 30 0.00929 1.37192 r +38 218.49300 18.81443 58710.15734 22.05000 30 0.00856 1.44715 g +38 218.49300 18.81443 58710.15781 22.05000 30 0.00856 1.45090 r +33 269.78931 18.82039 58710.15827 22.05000 30 0.00792 1.02709 g +33 269.78931 18.82039 58710.15873 22.05000 30 0.00792 1.02704 r +32 266.93213 18.86146 58710.15919 22.05000 30 0.00750 1.02714 g +32 266.93213 18.86146 58710.15966 22.05000 30 0.00750 1.02722 r +64 223.42493 18.89317 58710.16012 22.05000 30 0.00658 1.36590 g +64 223.42493 18.89317 58710.16058 22.05000 30 0.00658 1.36904 r +35 285.14725 18.83654 58710.16105 22.05000 30 0.00657 1.05929 g +35 285.14725 18.83654 58710.16151 22.05000 30 0.00657 1.05858 r +47 213.31319 18.87345 58710.16197 22.05000 30 0.00623 1.62673 g +47 213.31319 18.87345 58710.16244 22.05000 30 0.00623 1.63192 r +79 334.13101 18.96335 58710.16290 22.05000 30 0.00611 1.90789 g +79 334.13101 18.96335 58710.16336 22.05000 30 0.00611 1.90028 r +37 254.31154 18.88439 58710.16382 22.05000 30 0.00608 1.06110 g +37 254.31154 18.88439 58710.16429 22.05000 30 0.00608 1.06184 r +41 324.31491 18.96594 58710.16475 22.05000 30 0.00559 1.53993 g +41 324.31491 18.96594 58710.16521 22.05000 30 0.00559 1.53549 r +42 338.35648 18.71071 58710.16568 22.05000 30 0.00547 2.07746 g +42 338.35648 18.71071 58710.16614 22.05000 30 0.00547 2.06820 r +43 324.91791 18.83660 58710.16660 22.05000 30 0.00516 1.53999 g +43 324.91791 18.83660 58710.16706 22.05000 30 0.00516 1.53555 r +45 333.38287 18.83717 58710.16753 22.05000 30 0.00503 1.80696 g +45 333.38287 18.83717 58710.16799 22.05000 30 0.00503 1.80027 r +46 319.12234 18.87075 58710.16845 22.05000 30 0.00498 1.39069 g +46 319.12234 18.87075 58710.16892 22.05000 30 0.00498 1.38739 r +48 329.22885 18.71823 58710.16938 22.05000 30 0.00492 1.63762 g +48 329.22885 18.71823 58710.16984 22.05000 30 0.00492 1.63239 r +49 315.51770 18.93398 58710.17031 22.05000 30 0.00492 1.31305 g +49 315.51770 18.93398 58710.17077 22.05000 30 0.00492 1.31033 r +51 212.41476 18.83974 58710.17123 22.05000 30 0.00478 1.77351 g +51 212.41476 18.83974 58710.17169 22.05000 30 0.00478 1.77996 r +52 269.15030 18.71571 58710.17216 22.05000 30 0.00438 1.02968 g +52 269.15030 18.71571 58710.17262 22.05000 30 0.00438 1.02987 r +57 211.11278 18.93243 58710.17308 22.05000 30 0.00420 1.85128 g +57 211.11278 18.93243 58710.17355 22.05000 30 0.00420 1.85844 r +90 342.35361 18.88591 58710.17401 22.05000 30 0.00401 2.12878 g +90 342.35361 18.88591 58710.17447 22.05000 30 0.00401 2.11900 r +56 327.47778 18.97655 58710.17494 22.05000 30 0.00386 1.52625 g +56 327.47778 18.97655 58710.17540 22.05000 30 0.00386 1.52191 r +58 329.97675 18.87427 58710.17586 22.05000 30 0.00370 1.58715 g +58 329.97675 18.87427 58710.17632 22.05000 30 0.00370 1.58233 r +70 219.85963 18.82935 58710.17679 22.05000 30 0.00369 1.59143 g +70 219.85963 18.82935 58710.17725 22.05000 30 0.00369 1.59632 r +19 348.91486 18.75476 58710.17771 22.05000 30 0.01137 2.50134 g +19 348.91486 18.75476 58710.17818 22.05000 30 0.01137 2.48740 r +59 333.13058 18.97282 58710.17864 22.05000 30 0.00368 1.65117 g +59 333.13058 18.97282 58710.17910 22.05000 30 0.00368 1.64583 r +60 326.20605 18.81969 58710.17956 22.05000 30 0.00360 1.45656 g +60 326.20605 18.81969 58710.18003 22.05000 30 0.00360 1.45276 r +61 224.74982 18.78513 58710.18049 22.05000 30 0.00346 1.49674 g +61 224.74982 18.78513 58710.18095 22.05000 30 0.00346 1.50087 r +63 339.67850 18.87129 58710.18142 22.05000 30 0.00342 1.85960 g +63 339.67850 18.87129 58710.18188 22.05000 30 0.00342 1.85243 r +104 220.64230 18.79817 58710.18234 22.05000 30 0.00338 1.62866 g +104 220.64230 18.79817 58710.18281 22.05000 30 0.00338 1.63386 r +65 315.28528 18.86793 58710.18327 22.05000 30 0.00324 1.24118 g +65 315.28528 18.86793 58710.18373 22.05000 30 0.00324 1.23898 r +66 340.19644 18.90464 58710.18419 22.05000 30 0.00316 1.83827 g +34 351.38257 18.78011 58710.18466 22.05000 30 0.00877 2.49699 g +34 351.38257 18.78011 58710.18512 22.05000 30 0.00877 2.48310 r +66 340.19644 18.90464 58710.18558 22.05000 30 0.00316 1.81755 r +67 324.42368 18.82575 58710.18605 22.05000 30 0.00306 1.37084 g +67 324.42368 18.82575 58710.18651 22.05000 30 0.00306 1.36769 r +68 284.51547 18.86341 58710.18697 22.05000 30 0.00305 1.03080 g +68 284.51547 18.86341 58710.18744 22.05000 30 0.00305 1.03055 r +62 352.78506 18.94286 58710.18790 22.05000 30 0.00337 2.50907 g +62 352.78506 18.94286 58710.18836 22.05000 30 0.00337 2.49505 r +69 258.12866 18.77626 58710.18882 22.05000 30 0.00292 1.08864 g +69 258.12866 18.77626 58710.18929 22.05000 30 0.00292 1.08966 r +72 266.79120 18.81468 58710.18975 22.05000 30 0.00285 1.04908 g +31 353.51883 18.83667 58710.19021 22.05000 30 0.01035 2.50543 g +31 353.51883 18.83667 58710.19068 22.05000 30 0.01035 2.49145 r +72 266.79120 18.81468 58710.19114 22.05000 30 0.00285 1.05085 r +118 223.06654 18.97986 58710.19160 22.05000 30 0.00283 1.65468 g +118 223.06654 18.97986 58710.19206 22.05000 30 0.00283 1.66009 r +74 342.46384 18.84121 58710.19253 22.05000 30 0.00279 1.80922 g +74 342.46384 18.84121 58710.19299 22.05000 30 0.00279 1.80251 r +109 220.25854 18.90368 58710.19345 22.05000 30 0.00278 1.77904 g +109 220.25854 18.90368 58710.19392 22.05000 30 0.00278 1.78554 r +77 278.71726 18.85578 58710.19438 22.05000 30 0.00269 1.02761 g +77 278.71726 18.85578 58710.19484 22.05000 30 0.00269 1.02773 r +120 223.69275 18.77448 58710.19531 22.05000 30 0.00268 1.68162 g +120 223.69275 18.77448 58710.19577 22.05000 30 0.00268 1.68727 r +78 224.49289 18.71036 58710.19623 22.05000 30 0.00267 1.66723 g +78 224.49289 18.71036 58710.19669 22.05000 30 0.00267 1.67275 r +143 349.32709 18.73942 58710.19716 22.05000 30 0.00265 2.05514 g +143 349.32709 18.73942 58710.19762 22.05000 30 0.00265 2.04610 r +80 221.62842 18.91681 58710.19808 22.05000 30 0.00258 1.79051 g +80 221.62842 18.91681 58710.19855 22.05000 30 0.00258 1.79711 r +82 318.34644 18.72255 58710.19901 22.05000 30 0.00256 1.20945 g +82 318.34644 18.72255 58710.19947 22.05000 30 0.00256 1.20748 r +83 274.65552 18.79406 58710.19994 22.05000 30 0.00250 1.03703 g +83 274.65552 18.79406 58710.20040 22.05000 30 0.00250 1.03742 r +85 342.21787 18.78295 58710.20086 22.05000 30 0.00246 1.68973 g +85 342.21787 18.78295 58710.20132 22.05000 30 0.00246 1.68407 r +88 329.55807 18.73863 58710.20179 22.05000 30 0.00243 1.36147 g +88 329.55807 18.73863 58710.20225 22.05000 30 0.00243 1.35839 r +92 337.25452 18.86579 58710.20271 22.05000 30 0.00222 1.52128 g +92 337.25452 18.86579 58710.20318 22.05000 30 0.00222 1.51699 r +101 349.93204 18.88200 58710.20364 22.05000 30 0.00218 1.96186 g +101 349.93204 18.88200 58710.20410 22.05000 30 0.00218 1.95374 r +94 278.58041 18.77123 58710.20456 22.05000 30 0.00216 1.03257 g +94 278.58041 18.77123 58710.20503 22.05000 30 0.00216 1.03286 r +96 320.25500 18.84582 58710.20549 22.05000 30 0.00213 1.20353 g +96 320.25500 18.84582 58710.20595 22.05000 30 0.00213 1.20161 r +99 257.27289 18.87008 58710.20642 22.05000 30 0.00207 1.14097 g +99 257.27289 18.87008 58710.20688 22.05000 30 0.00207 1.14243 r +102 317.09354 18.84816 58710.20734 22.05000 30 0.00194 1.16281 g +102 317.09354 18.84816 58710.20781 22.05000 30 0.00194 1.16119 r +122 322.55020 18.92177 58710.20827 22.05000 30 0.00193 1.21837 g +122 322.55020 18.92177 58710.20873 22.05000 30 0.00193 1.21633 r +113 274.03949 18.79164 58710.20919 22.05000 30 0.00186 1.04833 g +113 274.03949 18.79164 58710.20966 22.05000 30 0.00186 1.04891 r +105 218.50989 18.97303 58710.21012 22.05000 30 0.00183 2.15353 g +105 218.50989 18.97303 58710.21058 22.05000 30 0.00183 2.16368 r +107 353.28900 18.82143 58710.21105 22.05000 30 0.00180 1.99760 g +107 353.28900 18.82143 58710.21151 22.05000 30 0.00180 1.98913 r +139 349.64108 18.82275 58710.21197 22.05000 30 0.00179 1.81612 g +139 349.64108 18.82275 58710.21244 22.05000 30 0.00179 1.80935 r +133 283.18082 18.96308 58710.21290 22.05000 30 0.00164 1.02934 g +133 283.18082 18.96308 58710.21336 22.05000 30 0.00164 1.02956 r +115 215.80472 18.81650 58710.21382 22.05000 30 0.00163 2.43977 g +115 215.80472 18.81650 58710.21429 22.05000 30 0.00163 2.45314 r +116 282.18713 18.91864 58710.21475 22.05000 30 0.00162 1.03203 g +116 282.18713 18.91864 58710.21521 22.05000 30 0.00162 1.03233 r +123 336.26910 18.95957 58710.21568 22.05000 30 0.00148 1.39175 g +123 336.26910 18.95957 58710.21614 22.05000 30 0.00148 1.38844 r +127 342.57483 18.93278 58710.21660 22.05000 30 0.00141 1.52837 g +127 342.57483 18.93278 58710.21706 22.05000 30 0.00141 1.52402 r +128 281.43896 18.92883 58710.21753 22.05000 30 0.00141 1.03542 g +128 281.43896 18.92883 58710.21799 22.05000 30 0.00141 1.03579 r +136 351.01355 18.80262 58710.21845 22.05000 30 0.00138 1.77820 g +136 351.01355 18.80262 58710.21892 22.05000 30 0.00138 1.77177 r +129 353.41766 18.87845 58710.21938 22.05000 30 0.00136 1.86117 g +129 353.41766 18.87845 58710.21984 22.05000 30 0.00136 1.85399 r +131 353.80219 18.75913 58710.22031 22.05000 30 0.00136 1.86607 g +131 353.80219 18.75913 58710.22077 22.05000 30 0.00136 1.85885 r +132 315.02631 18.84650 58710.22123 22.05000 30 0.00134 1.10425 g +132 315.02631 18.84650 58710.22169 22.05000 30 0.00134 1.10310 r +134 342.63327 18.74779 58710.22216 22.05000 30 0.00129 1.48208 g +134 342.63327 18.74779 58710.22262 22.05000 30 0.00129 1.47809 r +189 283.79215 18.84751 58710.22308 22.05000 30 0.00125 1.03500 g +189 283.79215 18.84751 58710.22355 22.05000 30 0.00125 1.03536 r +137 284.31290 18.78095 58710.22401 22.05000 30 0.00125 1.03490 g +137 284.31290 18.78095 58710.22447 22.05000 30 0.00125 1.03524 r +140 281.05728 18.79243 58710.22494 22.05000 30 0.00124 1.04397 g +140 281.05728 18.79243 58710.22540 22.05000 30 0.00124 1.04448 r +142 340.79855 18.89365 58710.22586 22.05000 30 0.00123 1.40975 g +142 340.79855 18.89365 58710.22632 22.05000 30 0.00123 1.40631 r +145 254.82585 18.80448 58710.22679 22.05000 30 0.00121 1.24912 g +145 254.82585 18.80448 58710.22725 22.05000 30 0.00121 1.25140 r +146 351.68777 18.79003 58710.22771 22.05000 30 0.00120 1.68217 g +146 351.68777 18.79003 58710.22818 22.05000 30 0.00120 1.67657 r +151 265.63812 18.96231 58710.22864 22.05000 30 0.00107 1.13738 g +151 265.63812 18.96231 58710.22910 22.05000 30 0.00107 1.13881 r +157 325.80075 18.89421 58710.22956 22.05000 30 0.00097 1.16928 g +157 325.80075 18.89421 58710.23003 22.05000 30 0.00097 1.16760 r +158 340.54907 18.74728 58710.23049 22.05000 30 0.00097 1.37322 g +158 340.54907 18.74728 58710.23095 22.05000 30 0.00097 1.37005 r +164 222.11273 18.90058 58710.23142 22.05000 30 0.00095 2.43938 g +164 222.11273 18.90058 58710.23188 22.05000 30 0.00095 2.45274 r +241 221.80238 18.85918 58710.23234 22.05000 30 0.00093 2.49370 g +163 252.00386 18.85865 58710.23281 22.05000 30 0.00090 1.32436 g +163 252.00386 18.85865 58710.23327 22.05000 30 0.00090 1.32719 r +214 223.79947 18.89983 58710.23373 22.05000 30 0.00088 2.37369 g +214 223.79947 18.89983 58710.23419 22.05000 30 0.00088 2.38628 r +168 316.54068 18.75212 58710.23466 22.05000 30 0.00084 1.08363 g +168 316.54068 18.75212 58710.23512 22.05000 30 0.00084 1.08267 r +173 353.90216 18.94598 58710.23558 22.05000 30 0.00078 1.65880 g +173 353.90216 18.94598 58710.23605 22.05000 30 0.00078 1.65339 r +207 332.71875 18.81938 58710.23651 22.05000 30 0.00077 1.21879 g +207 332.71875 18.81938 58710.23697 22.05000 30 0.00077 1.21675 r +220 223.87431 18.83954 58710.23744 22.05000 30 0.00074 2.47522 g +220 223.87431 18.83954 58710.23790 22.05000 30 0.00074 2.48901 r +182 304.99896 18.75510 58710.23836 22.05000 30 0.00069 1.03449 g +182 304.99896 18.75510 58710.23882 22.05000 30 0.00069 1.03416 r +184 282.46033 18.96808 58710.23929 22.05000 30 0.00067 1.05662 g +184 282.46033 18.96808 58710.23975 22.05000 30 0.00067 1.05731 r +191 256.87500 18.77634 58710.24021 22.05000 30 0.00063 1.28968 g +191 256.87500 18.77634 58710.24068 22.05000 30 0.00063 1.29226 r +195 317.57877 18.78514 58710.24114 22.05000 30 0.00060 1.07623 g +195 317.57877 18.78514 58710.24160 22.05000 30 0.00060 1.07534 r +196 271.97540 18.76611 58710.24206 22.05000 30 0.00060 1.12611 g +196 271.97540 18.76611 58710.24253 22.05000 30 0.00060 1.12745 r +199 354.19092 18.73840 58710.24299 22.05000 30 0.00057 1.58887 g +199 354.19092 18.73840 58710.24345 22.05000 30 0.00057 1.58404 r +202 338.78482 18.73611 58710.24392 22.05000 30 0.00057 1.26424 g +202 338.78482 18.73611 58710.24438 22.05000 30 0.00057 1.26188 r +205 317.95248 18.87463 58710.24484 22.05000 30 0.00056 1.07079 g +205 317.95248 18.87463 58710.24531 22.05000 30 0.00056 1.06996 r +206 266.51761 18.82009 58710.24577 22.05000 30 0.00056 1.18977 g +206 266.51761 18.82009 58710.24623 22.05000 30 0.00056 1.19160 r +208 284.60147 18.95310 58710.24669 22.05000 30 0.00054 1.05890 g +208 284.60147 18.95310 58710.24716 22.05000 30 0.00054 1.05961 r +209 316.98190 18.74158 58710.24762 22.05000 30 0.00053 1.06208 g +209 316.98190 18.74158 58710.24808 22.05000 30 0.00053 1.06135 r +215 339.40707 18.76307 58710.24855 22.05000 30 0.00051 1.24945 g +215 339.40707 18.76307 58710.24901 22.05000 30 0.00051 1.24719 r +217 282.12311 18.81700 58710.24947 22.05000 30 0.00048 1.07625 g +217 282.12311 18.81700 58710.24994 22.05000 30 0.00048 1.07715 r +81 28.49257 18.55032 58710.28790 22.05000 30 0.00257 2.49657 g +81 28.49257 18.55032 58710.28836 22.05000 30 0.00257 2.48268 r +40 29.76173 18.66024 58710.29114 22.05000 30 0.00723 2.49994 g +40 29.76173 18.66024 58710.29160 22.05000 30 0.00723 2.48601 r +50 29.81568 18.67970 58710.29206 22.05000 30 0.00469 2.47581 g +53 30.23251 18.67836 58710.29253 22.05000 30 0.00695 2.49658 g +53 30.23251 18.67836 58710.29299 22.05000 30 0.00695 2.48270 r +55 30.51788 18.59102 58710.29345 22.05000 30 0.00528 2.49653 g +55 30.51788 18.59102 58710.29392 22.05000 30 0.00528 2.48264 r +50 29.81568 18.67970 58710.29438 22.05000 30 0.00469 2.40921 r +181 32.30503 18.70350 58710.29808 22.05000 30 0.00121 2.50122 g +181 32.30503 18.70350 58710.29855 22.05000 30 0.00121 2.48728 r +44 32.61850 18.75967 58710.29901 22.05000 30 0.00781 2.49692 g +44 32.61850 18.75967 58710.29947 22.05000 30 0.00781 2.48304 r +148 32.66304 18.59161 58710.29994 22.05000 30 0.00114 2.48053 g +148 32.66304 18.59161 58710.30040 22.05000 30 0.00114 2.46683 r +124 33.38802 18.77011 58710.30086 22.05000 30 0.00161 2.50498 g +124 33.38802 18.77011 58710.30132 22.05000 30 0.00161 2.49100 r +159 33.38687 18.82233 58710.30179 22.05000 30 0.00144 2.47475 g +159 33.38687 18.82233 58710.30225 22.05000 30 0.00144 2.46113 r +223 32.40858 18.75960 58710.30271 22.05000 30 0.00089 2.37446 g +223 32.40858 18.75960 58710.30318 22.05000 30 0.00089 2.36201 r +138 40.94326 18.68245 58710.32216 22.05000 30 0.00125 2.49856 g +138 40.94326 18.68245 58710.32262 22.05000 30 0.00125 2.48465 r +203 94.87891 18.58862 58710.47169 22.05000 30 0.00056 2.50386 g +203 94.87891 18.58862 58710.47216 22.05000 30 0.00056 2.48988 r +135 97.36169 18.76374 58710.47818 22.05000 30 0.00129 2.50801 g +135 97.36169 18.76374 58710.47864 22.05000 30 0.00129 2.49399 r +226 98.13447 18.62946 58710.48049 22.05000 30 0.00064 2.50898 g +161 98.24920 18.68051 58710.48095 22.05000 30 0.00094 2.50221 g +26 98.48325 18.65961 58710.48142 22.05000 30 0.00872 2.50883 g +26 98.48325 18.65961 58710.48188 22.05000 30 0.00872 2.49480 r +28 98.77230 18.75695 58710.48234 22.05000 30 0.00918 2.50056 g +28 98.77230 18.75695 58710.48281 22.05000 30 0.00918 2.48662 r +100 99.17927 18.76347 58710.48327 22.05000 30 0.00199 2.50640 g +100 99.17927 18.76347 58710.48373 22.05000 30 0.00199 2.49240 r +112 99.44604 18.66575 58710.48419 22.05000 30 0.00167 2.50524 g +112 99.44604 18.66575 58710.48466 22.05000 30 0.00167 2.49125 r +161 98.24920 18.68051 58710.48512 22.05000 30 0.00094 2.38238 r +27 98.77992 18.71250 58710.48558 22.05000 30 0.00068 2.40888 g +27 98.77992 18.71250 58710.48605 22.05000 30 0.00068 2.39603 r +226 98.13447 18.62946 58710.48651 22.05000 30 0.00064 2.33881 r +197 101.08033 18.68633 58710.48882 22.05000 30 0.00060 2.50130 g +197 101.08033 18.68633 58710.48929 22.05000 30 0.00060 2.48736 r +188 102.28793 18.80531 58710.49206 22.05000 30 0.00065 2.49908 g +188 102.28793 18.80531 58710.49253 22.05000 30 0.00065 2.48517 r +193 102.72141 18.59939 58710.49345 22.05000 30 0.00061 2.50289 g +193 102.72141 18.59939 58710.49392 22.05000 30 0.00061 2.48892 r +91 103.00915 18.56834 58710.49438 22.05000 30 0.00224 2.50043 g +91 103.00915 18.56834 58710.49484 22.05000 30 0.00224 2.48649 r diff --git a/gwemopt/tests/data/expected_results/schedule_ZTF.dat b/gwemopt/tests/data/expected_results/schedule_ZTF.dat index e063376d..d831a174 100644 --- a/gwemopt/tests/data/expected_results/schedule_ZTF.dat +++ b/gwemopt/tests/data/expected_results/schedule_ZTF.dat @@ -1,4 +1,4 @@ -246 8.55346 -24.25000 58710.35336 20.40000 30 0.30998 2.49913 g -247 15.94654 -24.25000 58710.37373 20.40000 30 0.34024 2.50095 g -247 15.94654 -24.25000 58710.37442 20.40000 30 0.34024 2.48909 r -246 8.55346 -24.25000 58710.37488 20.40000 30 0.30998 2.20549 r +246 8.55346 -24.25000 58710.35336 20.40000 30 0.39960 2.49913 g +247 15.94654 -24.25000 58710.37373 20.40000 30 0.65706 2.50095 g +247 15.94654 -24.25000 58710.37442 20.40000 30 0.65706 2.48909 r +246 8.55346 -24.25000 58710.37488 20.40000 30 0.39960 2.20549 r diff --git a/gwemopt/tests/data/expected_results/summary.dat b/gwemopt/tests/data/expected_results/summary.dat index d4abce8c..caf8025a 100644 --- a/gwemopt/tests/data/expected_results/summary.dat +++ b/gwemopt/tests/data/expected_results/summary.dat @@ -1,3 +1,3 @@ -1.0,10.04986,0.68596,12.64186,58710 -7.0,10.04986,0.68596,12.64186,58710 -60.0,10.04986,0.68596,12.64186,58710 +1.0,10.04986,0.60706,6.82582,58710 +7.0,10.04986,0.60706,6.82582,58710 +60.0,10.04986,0.60706,6.82582,58710 diff --git a/gwemopt/tests/data/expected_results/summary_coverage.dat b/gwemopt/tests/data/expected_results/summary_coverage.dat index 3bf93832..ad2b59d1 100644 --- a/gwemopt/tests/data/expected_results/summary_coverage.dat +++ b/gwemopt/tests/data/expected_results/summary_coverage.dat @@ -1,3 +1,3 @@ -1.0,10.04986,0.88843,93.58124,58710 -7.0,10.04986,0.88843,93.58124,58710 -60.0,10.04986,0.88843,93.58124,58710 +1.0,10.04986,0.88970,95.15492,58710 +7.0,10.04986,0.88970,95.15492,58710 +60.0,10.04986,0.88970,95.15492,58710 diff --git a/gwemopt/tests/data/expected_results/tiles_coverage_int_DECam.txt b/gwemopt/tests/data/expected_results/tiles_coverage_int_DECam.txt index 99d0e25e..685eccf1 100644 --- a/gwemopt/tests/data/expected_results/tiles_coverage_int_DECam.txt +++ b/gwemopt/tests/data/expected_results/tiles_coverage_int_DECam.txt @@ -1,10 +1,4 @@ -0.4187440827 2.9236558354e-01 3.1473511695 -0.4194153790 4.5801493250e-01 6.1897906334 -0.4200866753 6.0490832882e-01 8.0782013352 -0.4207579716 6.5886604363e-01 9.9666120369 -0.4214292679 6.8595736466e-01 12.6418605310 -0.4392301938 6.8595736466e-01 12.6418605310 -0.4399014901 6.8595736466e-01 12.6418605310 -0.4405727864 6.8595736466e-01 12.6418605310 -0.4412440827 6.8595736466e-01 12.6418605310 -0.4419153790 6.8595736466e-01 12.6418605310 +0.4187440827 3.9058285832e-01 4.0161512320 +0.4194153790 6.0705716806e-01 6.8258178489 +0.4392301938 6.0705716806e-01 6.8258178489 +0.4399014901 6.0705716806e-01 6.8258178489 diff --git a/gwemopt/tests/data/expected_results/tiles_coverage_int_KPED.txt b/gwemopt/tests/data/expected_results/tiles_coverage_int_KPED.txt index 7ffb96c4..3df7fc70 100644 --- a/gwemopt/tests/data/expected_results/tiles_coverage_int_KPED.txt +++ b/gwemopt/tests/data/expected_results/tiles_coverage_int_KPED.txt @@ -1,355 +1,307 @@ -0.4485578166 2.0126154540e-03 2.0000000000 -0.4490207795 2.0126154540e-03 2.0000000000 -0.4494837425 2.5053852190e-03 3.0000000000 -0.4499467054 2.5053852190e-03 3.0000000000 -0.4504096684 3.1387408580e-03 4.0000000000 -0.4508726314 4.2170416855e-03 6.0000000000 -0.4513355943 4.2170416855e-03 6.0000000000 -0.4517985573 5.9307479186e-03 7.0000000000 -0.4522615203 5.9307479186e-03 7.0000000000 -0.4527244832 8.8284037949e-03 8.0000000000 -0.4531874462 8.8284037949e-03 8.0000000000 -0.4536504091 1.3587537335e-02 9.0000000000 -0.4541133721 1.3587537335e-02 9.0000000000 -0.4545763351 1.6949860967e-02 11.0000000000 -0.4550392980 1.6949860967e-02 11.0000000000 -0.4555022610 2.2272817513e-02 12.0000000000 -0.4559652240 3.0616255671e-02 14.0000000000 -0.4564281869 3.0616255671e-02 14.0000000000 -0.4568911499 7.0840956838e-02 15.0000000000 -0.4573541128 7.0840956838e-02 15.0000000000 -0.4578170758 8.6137743180e-02 18.0000000000 -0.4582800388 1.2182204271e-01 20.0000000000 -0.4587430017 1.2182204271e-01 20.0000000000 -0.4592059647 1.6835992064e-01 21.0000000000 -0.4596689277 1.6835992064e-01 21.0000000000 -0.4601318906 1.6835992064e-01 21.0000000000 -0.4605948536 1.8318093626e-01 23.0000000000 -0.4610578165 1.8318093626e-01 23.0000000000 -0.4615207795 1.9596654110e-01 24.0000000000 -0.4619837425 1.9596654110e-01 24.0000000000 -0.4624467054 2.0733705438e-01 28.0000000000 -0.4629096684 2.0733705438e-01 28.0000000000 -0.4633726313 2.2212475292e-01 29.0000000000 -0.4638355943 2.5337867009e-01 30.0000000000 -0.4642985573 2.5337867009e-01 30.0000000000 -0.4647615202 2.6995406338e-01 31.0000000000 -0.4652244832 2.6995406338e-01 31.0000000000 -0.4656874462 2.6995406338e-01 31.0000000000 -0.4661504091 2.8468714278e-01 32.0000000000 -0.4666133721 2.8468714278e-01 32.0000000000 -0.4670763350 2.9880411981e-01 33.0000000000 -0.4675392980 2.9880411981e-01 33.0000000000 -0.4680022610 3.1175115909e-01 34.0000000000 -0.4684652239 3.1175115909e-01 34.0000000000 -0.4689281869 3.2992796562e-01 36.0000000000 -0.4693911499 3.2992796562e-01 36.0000000000 -0.4698541128 3.6247759222e-01 38.0000000000 -0.4703170758 3.6247759222e-01 38.0000000000 -0.4707800387 3.8674563293e-01 39.0000000000 -0.4712430017 3.8674563293e-01 39.0000000000 -0.4717059647 4.1259456103e-01 41.0000000000 -0.4721689276 4.1259456103e-01 41.0000000000 -0.4726318906 4.3508622829e-01 42.0000000000 -0.4730948536 4.3508622829e-01 42.0000000000 -0.4735578165 4.5584359294e-01 45.0000000000 -0.4740207795 4.5584359294e-01 45.0000000000 -0.4744837424 4.7806162994e-01 47.0000000000 -0.4749467054 4.7806162994e-01 47.0000000000 -0.4754096684 4.9205685215e-01 48.0000000000 -0.4758726313 4.9205685215e-01 48.0000000000 -0.4763355943 5.0445078676e-01 49.0000000000 -0.4767985572 5.0445078676e-01 49.0000000000 -0.4772615202 5.1638966530e-01 51.0000000000 -0.4777244832 5.1638966530e-01 51.0000000000 -0.4781874461 5.2831505414e-01 52.0000000000 -0.4786504091 5.2831505414e-01 52.0000000000 -0.4791133721 5.3978232324e-01 53.0000000000 -0.4795763350 5.3978232324e-01 53.0000000000 -0.4800392980 5.5040066841e-01 54.0000000000 -0.4805022609 5.5040066841e-01 54.0000000000 -0.4809652239 5.6095977985e-01 55.0000000000 -0.4814281869 5.6095977985e-01 55.0000000000 -0.4818911498 5.7123603012e-01 56.0000000000 -0.4823541128 5.7123603012e-01 56.0000000000 -0.4828170758 5.8093087673e-01 57.0000000000 -0.4832800387 5.8093087673e-01 57.0000000000 -0.4837430017 5.9036462250e-01 60.0000000000 -0.4842059646 5.9036462250e-01 60.0000000000 -0.4846689276 5.9951464811e-01 62.0000000000 -0.4851318906 5.9951464811e-01 62.0000000000 -0.4855948535 6.0859664028e-01 63.0000000000 -0.4860578165 6.0859664028e-01 63.0000000000 -0.4865207795 6.1762764899e-01 64.0000000000 -0.4869837424 6.1762764899e-01 64.0000000000 -0.4874467054 6.2662405508e-01 65.0000000000 -0.4879096683 6.2662405508e-01 65.0000000000 -0.4883726313 6.3545100942e-01 66.0000000000 -0.4888355943 6.3545100942e-01 66.0000000000 -0.4892985572 6.4418813288e-01 67.0000000000 -0.4897615202 6.4418813288e-01 67.0000000000 -0.4902244831 6.5283396042e-01 68.0000000000 -0.4906874461 6.5283396042e-01 68.0000000000 -0.4911504091 6.6131123640e-01 72.0000000000 -0.4916133720 6.6131123640e-01 72.0000000000 -0.4920763350 6.6964681260e-01 73.0000000000 -0.4925392980 6.6964681260e-01 73.0000000000 -0.4930022609 6.7791422432e-01 74.0000000000 -0.4934652239 6.7791422432e-01 74.0000000000 -0.4939281868 6.8568054607e-01 77.0000000000 -0.4943911498 6.8568054607e-01 77.0000000000 -0.4948541128 6.9326571385e-01 80.0000000000 -0.4953170757 6.9326571385e-01 80.0000000000 -0.4957800387 7.0077799846e-01 81.0000000000 -0.4962430017 7.0077799846e-01 81.0000000000 -0.4967059646 7.0795927788e-01 82.0000000000 -0.4971689276 7.0795927788e-01 82.0000000000 -0.4976318905 7.1488410713e-01 85.0000000000 -0.4980948535 7.1488410713e-01 85.0000000000 -0.4985578165 7.2135675148e-01 86.0000000000 -0.4990207794 7.2135675148e-01 86.0000000000 -0.4994837424 7.2733557965e-01 89.0000000000 -0.4999467053 7.2733557965e-01 89.0000000000 -0.5004096683 7.3324011584e-01 90.0000000000 -0.5008726313 7.3324011584e-01 90.0000000000 -0.5013355942 7.3883177968e-01 92.0000000000 -0.5017985572 7.3883177968e-01 92.0000000000 -0.5022615202 7.3883177968e-01 92.0000000000 -0.5027244831 7.4375212756e-01 94.0000000000 -0.5031874461 7.4375212756e-01 94.0000000000 -0.5036504090 7.4833820774e-01 95.0000000000 -0.5041133720 7.4833820774e-01 95.0000000000 -0.5045763350 7.5269158993e-01 96.0000000000 -0.5050392979 7.5269158993e-01 96.0000000000 -0.5055022609 7.5701042082e-01 97.0000000000 -0.5059652239 7.5701042082e-01 97.0000000000 -0.5064281868 7.6105960295e-01 98.0000000000 -0.5068911498 7.6105960295e-01 98.0000000000 -0.5073541127 7.6477418858e-01 101.0000000000 -0.5078170757 7.6477418858e-01 101.0000000000 -0.5082800387 7.6831687884e-01 103.0000000000 -0.5087430016 7.6831687884e-01 103.0000000000 -0.5092059646 7.7169925763e-01 104.0000000000 -0.5096689276 7.7169925763e-01 104.0000000000 -0.5101318905 7.7506768753e-01 105.0000000000 -0.5105948535 7.7506768753e-01 105.0000000000 -0.5110578164 7.7840617641e-01 107.0000000000 -0.5115207794 7.7840617641e-01 107.0000000000 -0.5119837424 7.8166970868e-01 108.0000000000 -0.5124467053 7.8166970868e-01 108.0000000000 -0.5129096683 7.8490273146e-01 110.0000000000 -0.5133726312 7.8490273146e-01 110.0000000000 -0.5138355942 7.8805049262e-01 111.0000000000 -0.5142985572 7.8805049262e-01 111.0000000000 -0.5147615201 7.9117851086e-01 113.0000000000 -0.5152244831 7.9117851086e-01 113.0000000000 -0.5156874461 7.9429447744e-01 115.0000000000 -0.5161504090 7.9429447744e-01 115.0000000000 -0.5166133720 7.9729972004e-01 116.0000000000 -0.5170763349 7.9729972004e-01 116.0000000000 -0.5175392979 8.0020561884e-01 117.0000000000 -0.5180022609 8.0020561884e-01 117.0000000000 -0.5184652238 8.0311032310e-01 118.0000000000 -0.5189281868 8.0311032310e-01 118.0000000000 -0.5193911498 8.0594695994e-01 119.0000000000 -0.5198541127 8.0594695994e-01 119.0000000000 -0.5203170757 8.0866972997e-01 121.0000000000 -0.5207800386 8.0866972997e-01 121.0000000000 -0.5212430016 8.1132633290e-01 123.0000000000 -0.5217059646 8.1132633290e-01 123.0000000000 -0.5221689275 8.1396082475e-01 125.0000000000 -0.5226318905 8.1396082475e-01 125.0000000000 -0.5230948535 8.1657001449e-01 126.0000000000 -0.5235578164 8.1657001449e-01 126.0000000000 -0.5240207794 8.1914187226e-01 127.0000000000 -0.5244837423 8.1914187226e-01 127.0000000000 -0.5249467053 8.2163067146e-01 129.0000000000 -0.5254096683 8.2163067146e-01 129.0000000000 -0.5258726312 8.2407009780e-01 130.0000000000 -0.5263355942 8.2407009780e-01 130.0000000000 -0.5267985571 8.2648193986e-01 131.0000000000 -0.5272615201 8.2648193986e-01 131.0000000000 -0.5277244831 8.2886869609e-01 133.0000000000 -0.5281874460 8.2886869609e-01 133.0000000000 -0.5286504090 8.3124876027e-01 134.0000000000 -0.5291133720 8.3124876027e-01 134.0000000000 -0.5295763349 8.3354954528e-01 137.0000000000 -0.5300392979 8.3354954528e-01 137.0000000000 -0.5305022608 8.3568216607e-01 138.0000000000 -0.5309652238 8.3568216607e-01 138.0000000000 -0.5314281868 8.3778988802e-01 140.0000000000 -0.5318911497 8.3778988802e-01 140.0000000000 -0.5323541127 8.3987106975e-01 141.0000000000 -0.5328170757 8.3987106975e-01 141.0000000000 -0.5332800386 8.4181437446e-01 142.0000000000 -0.5337430016 8.4181437446e-01 142.0000000000 -0.5342059645 8.4375268875e-01 143.0000000000 -0.5346689275 8.4375268875e-01 143.0000000000 -0.5351318905 8.4561026072e-01 144.0000000000 -0.5355948534 8.4561026072e-01 144.0000000000 -0.5360578164 8.4745138367e-01 145.0000000000 -0.5365207794 8.4745138367e-01 145.0000000000 -0.5369837423 8.4925976741e-01 146.0000000000 -0.5374467053 8.4925976741e-01 146.0000000000 -0.5379096682 8.5104993870e-01 147.0000000000 -0.5383726312 8.5104993870e-01 147.0000000000 -0.5388355942 8.5278922988e-01 150.0000000000 -0.5392985571 8.5278922988e-01 150.0000000000 -0.5397615201 8.5609081808e-01 153.0000000000 -0.5402244830 8.5609081808e-01 153.0000000000 -0.5406874460 8.5782817445e-01 154.0000000000 -0.5411504090 8.5782817445e-01 154.0000000000 -0.5416133719 8.5954222391e-01 155.0000000000 -0.5420763349 8.5954222391e-01 155.0000000000 -0.5425392979 8.6118150889e-01 156.0000000000 -0.5430022608 8.6118150889e-01 156.0000000000 -0.5434652238 8.6339819459e-01 157.0000000000 -0.5439281867 8.6339819459e-01 157.0000000000 -0.5443911497 8.6508217411e-01 158.0000000000 -0.5448541127 8.6508217411e-01 158.0000000000 -0.5453170756 8.6671430282e-01 159.0000000000 -0.5457800386 8.6671430282e-01 159.0000000000 -0.5462430016 8.6971622833e-01 161.0000000000 -0.5467059645 8.6971622833e-01 161.0000000000 -0.5471689275 8.7152688661e-01 163.0000000000 -0.5476318904 8.7152688661e-01 163.0000000000 -0.5480948534 8.7623912679e-01 167.0000000000 -0.5485578164 8.7623912679e-01 167.0000000000 -0.5490207793 8.7997553666e-01 169.0000000000 -0.5494837423 8.7997553666e-01 169.0000000000 -0.5499467053 8.8263821346e-01 172.0000000000 -0.5504096682 8.9302061952e-01 175.0000000000 -0.5508726312 8.9302061952e-01 175.0000000000 -0.5513355941 8.9302061952e-01 175.0000000000 -0.5517985571 8.9546079620e-01 179.0000000000 -0.5522615201 8.9546079620e-01 179.0000000000 -0.5527244830 8.9883256077e-01 181.0000000000 -0.5531874460 8.9883256077e-01 181.0000000000 -0.5536504089 9.0177209994e-01 184.0000000000 -0.5541133719 9.0177209994e-01 184.0000000000 -0.5545763349 9.0407215801e-01 186.0000000000 -0.5550392978 9.0659722198e-01 189.0000000000 -0.5555022608 9.0659722198e-01 189.0000000000 -0.5559652238 9.0659722198e-01 189.0000000000 -0.5564281867 9.0887731326e-01 193.0000000000 -0.5568911497 9.0887731326e-01 193.0000000000 -0.5573541126 9.1095841704e-01 195.0000000000 -0.5578170756 9.1095841704e-01 195.0000000000 -0.5582800386 9.1259705856e-01 198.0000000000 -0.5587430015 9.1259705856e-01 198.0000000000 -0.5592059645 9.1422872861e-01 199.0000000000 -0.5596689275 9.1422872861e-01 199.0000000000 -0.5601318904 9.1585248020e-01 200.0000000000 -0.5605948534 9.1585248020e-01 200.0000000000 -0.5610578163 9.1742555941e-01 202.0000000000 -0.5615207793 9.1742555941e-01 202.0000000000 -0.5619837423 9.1892353989e-01 205.0000000000 -0.5624467052 9.1892353989e-01 205.0000000000 -0.5629096682 9.2039024148e-01 208.0000000000 -0.5633726311 9.2039024148e-01 208.0000000000 -0.5638355941 9.2185079075e-01 209.0000000000 -0.5642985571 9.2185079075e-01 209.0000000000 -0.5647615200 9.2330156629e-01 210.0000000000 -0.5652244830 9.2330156629e-01 210.0000000000 -0.5656874460 9.2473654386e-01 212.0000000000 -0.5661504089 9.2473654386e-01 212.0000000000 -0.5666133719 9.2614180249e-01 214.0000000000 -0.5670763348 9.2614180249e-01 214.0000000000 -0.5675392978 9.2752796249e-01 216.0000000000 -0.5680022608 9.2752796249e-01 216.0000000000 -0.5684652237 9.2889196295e-01 218.0000000000 -0.5689281867 9.2889196295e-01 218.0000000000 -0.5693911497 9.3023233295e-01 220.0000000000 -0.5698541126 9.3023233295e-01 220.0000000000 -0.5703170756 9.3154379833e-01 223.0000000000 -0.5707800385 9.3154379833e-01 223.0000000000 -0.5712430015 9.3281950194e-01 225.0000000000 -0.5717059645 9.3281950194e-01 225.0000000000 -0.5721689274 9.3408461791e-01 226.0000000000 -0.5726318904 9.3408461791e-01 226.0000000000 -0.5730948534 9.3534405209e-01 228.0000000000 -0.5735578163 9.3534405209e-01 228.0000000000 -0.5740207793 9.3657228304e-01 229.0000000000 -0.5744837422 9.3657228304e-01 229.0000000000 -0.5749467052 9.3779354702e-01 231.0000000000 -0.5754096682 9.3779354702e-01 231.0000000000 -0.5758726311 9.3897544923e-01 232.0000000000 -0.5763355941 9.3897544923e-01 232.0000000000 -0.5767985570 9.4013510520e-01 233.0000000000 -0.5772615200 9.4013510520e-01 233.0000000000 -0.5777244830 9.4126395587e-01 235.0000000000 -0.5781874459 9.4126395587e-01 235.0000000000 -0.5786504089 9.4238357071e-01 236.0000000000 -0.5791133719 9.4238357071e-01 236.0000000000 -0.5795763348 9.4348348697e-01 237.0000000000 -0.5800392978 9.4348348697e-01 237.0000000000 -0.5805022607 9.4457462700e-01 238.0000000000 -0.5809652237 9.4457462700e-01 238.0000000000 -0.5814281867 9.4566365140e-01 239.0000000000 -0.5818911496 9.4566365140e-01 239.0000000000 -0.5823541126 9.4673883430e-01 240.0000000000 -0.5828170756 9.4673883430e-01 240.0000000000 -0.5832800385 9.4779735441e-01 243.0000000000 -0.5837430015 9.4779735441e-01 243.0000000000 -0.5842059644 9.4884487879e-01 244.0000000000 -0.5846689274 9.4884487879e-01 244.0000000000 -0.5851318904 9.4989200162e-01 245.0000000000 -0.5855948533 9.4989200162e-01 245.0000000000 -0.5860578163 9.5092830274e-01 246.0000000000 -0.5865207793 9.5092830274e-01 246.0000000000 -0.5869837422 9.5194663329e-01 247.0000000000 -0.5874467052 9.5194663329e-01 247.0000000000 -0.5879096681 9.5296188507e-01 248.0000000000 -0.5883726311 9.5296188507e-01 248.0000000000 -0.5888355941 9.5395981887e-01 249.0000000000 -0.5892985570 9.5395981887e-01 249.0000000000 -0.5897615200 9.5494526582e-01 250.0000000000 -0.5902244829 9.5494526582e-01 250.0000000000 -0.5906874459 9.5589020853e-01 251.0000000000 -0.5911504089 9.5589020853e-01 251.0000000000 -0.5916133718 9.5683419334e-01 252.0000000000 -0.5920763348 9.5683419334e-01 252.0000000000 -0.5925392978 9.5776429081e-01 254.0000000000 -0.5930022607 9.5776429081e-01 254.0000000000 -0.5934652237 9.5866864964e-01 255.0000000000 -0.5939281866 9.5866864964e-01 255.0000000000 -0.5943911496 9.5956104490e-01 256.0000000000 -0.5948541126 9.5956104490e-01 256.0000000000 -0.5953170755 9.6044093557e-01 257.0000000000 -0.5957800385 9.6044093557e-01 257.0000000000 -0.5962430015 9.6130556605e-01 258.0000000000 -0.5967059644 9.6130556605e-01 258.0000000000 -0.5971689274 9.6216092350e-01 261.0000000000 -0.5976318903 9.6216092350e-01 261.0000000000 -0.5980948533 9.6301295567e-01 263.0000000000 -0.5985578163 9.6301295567e-01 263.0000000000 -0.5990207792 9.6384071164e-01 264.0000000000 -0.5994837422 9.6384071164e-01 264.0000000000 -0.5999467052 9.6466486430e-01 265.0000000000 -0.6004096681 9.6466486430e-01 265.0000000000 -0.6008726311 9.6546645887e-01 266.0000000000 -0.6013355940 9.6546645887e-01 266.0000000000 -0.6017985570 9.6626233485e-01 268.0000000000 -0.6022615200 9.6626233485e-01 268.0000000000 -0.6027244829 9.6703802999e-01 269.0000000000 -0.6031874459 9.6703802999e-01 269.0000000000 -0.6036504088 9.6781049253e-01 270.0000000000 -0.6041133718 9.6781049253e-01 270.0000000000 -0.6045763348 9.6858130246e-01 271.0000000000 -0.6050392977 9.6858130246e-01 271.0000000000 -0.6055022607 9.6933538759e-01 272.0000000000 -0.6059652237 9.6933538759e-01 272.0000000000 -0.6064281866 9.7008302344e-01 273.0000000000 -0.6068911496 9.7008302344e-01 273.0000000000 -0.6073541125 9.7082742871e-01 275.0000000000 -0.6078170755 9.7082742871e-01 275.0000000000 -0.6082800385 9.7155486202e-01 276.0000000000 -0.6087430014 9.7155486202e-01 276.0000000000 -0.6092059644 9.7225602903e-01 277.0000000000 -0.6096689274 9.7225602903e-01 277.0000000000 -0.6101318903 9.7294581155e-01 278.0000000000 -0.6105948533 9.7294581155e-01 278.0000000000 -0.6110578162 9.7363487524e-01 280.0000000000 -0.6115207792 9.7363487524e-01 280.0000000000 -0.6119837422 9.7431614213e-01 281.0000000000 -0.6124467051 9.7431614213e-01 281.0000000000 +0.2490207799 4.6368087825e-02 2.0000000000 +0.2494837429 4.6368087825e-02 2.0000000000 +0.2499467058 8.5512006930e-02 3.0000000000 +0.2504096688 8.5512006930e-02 3.0000000000 +0.2508726318 1.1903865872e-01 4.0000000000 +0.2513355947 1.1903865872e-01 4.0000000000 +0.2517985577 1.5089402706e-01 5.0000000000 +0.2522615207 1.5089402706e-01 5.0000000000 +0.2527244836 1.7194594973e-01 7.0000000000 +0.2531874466 1.7194594973e-01 7.0000000000 +0.2536504095 1.9261992884e-01 9.0000000000 +0.2541133725 1.9261992884e-01 9.0000000000 +0.2545763355 2.1082952578e-01 11.0000000000 +0.2550392984 2.1082952578e-01 11.0000000000 +0.2555022614 2.2902537559e-01 13.0000000000 +0.2559652244 2.2902537559e-01 13.0000000000 +0.2564281873 2.4597589200e-01 17.0000000000 +0.2568911503 2.4597589200e-01 17.0000000000 +0.2573541132 2.6763943305e-01 18.0000000000 +0.2578170762 2.6763943305e-01 18.0000000000 +0.2582800392 2.8440002725e-01 20.0000000000 +0.2587430021 2.8440002725e-01 20.0000000000 +0.2592059651 3.0100010111e-01 22.0000000000 +0.2596689280 3.0100010111e-01 22.0000000000 +0.2601318910 3.1708579481e-01 23.0000000000 +0.2605948540 3.1708579481e-01 23.0000000000 +0.2610578169 3.3239079953e-01 24.0000000000 +0.2615207799 3.3239079953e-01 24.0000000000 +0.2619837429 3.4881090408e-01 25.0000000000 +0.2624467058 3.4881090408e-01 25.0000000000 +0.2629096688 3.6280900457e-01 26.0000000000 +0.2633726317 3.6280900457e-01 26.0000000000 +0.2638355947 3.7640882363e-01 29.0000000000 +0.2642985577 3.7640882363e-01 29.0000000000 +0.2647615206 3.8967908565e-01 30.0000000000 +0.2652244836 3.8967908565e-01 30.0000000000 +0.2656874466 4.0290023959e-01 32.0000000000 +0.2661504095 4.0290023959e-01 32.0000000000 +0.2666133725 4.1561714740e-01 33.0000000000 +0.2670763354 4.1561714740e-01 33.0000000000 +0.2675392984 4.2826689892e-01 35.0000000000 +0.2680022614 4.2826689892e-01 35.0000000000 +0.2684652243 4.4053206055e-01 36.0000000000 +0.2689281873 4.4053206055e-01 36.0000000000 +0.2693911503 4.5265581950e-01 38.0000000000 +0.2698541132 4.5265581950e-01 38.0000000000 +0.2703170762 4.6451123306e-01 40.0000000000 +0.2707800391 4.6451123306e-01 40.0000000000 +0.2712430021 4.7556737649e-01 42.0000000000 +0.2717059651 4.7556737649e-01 42.0000000000 +0.2721689280 4.8653726515e-01 43.0000000000 +0.2726318910 4.8653726515e-01 43.0000000000 +0.2730948539 4.9732341029e-01 44.0000000000 +0.2735578169 4.9732341029e-01 44.0000000000 +0.2740207799 5.0661740204e-01 45.0000000000 +0.2744837428 5.0661740204e-01 45.0000000000 +0.2749467058 5.1517496041e-01 53.0000000000 +0.2754096688 5.1517496041e-01 53.0000000000 +0.2758726317 5.2309761640e-01 55.0000000000 +0.2763355947 5.2309761640e-01 55.0000000000 +0.2767985576 5.3059980283e-01 56.0000000000 +0.2772615206 5.3059980283e-01 56.0000000000 +0.2777244836 5.3718372391e-01 59.0000000000 +0.2781874465 5.3718372391e-01 59.0000000000 +0.2786504095 5.4375798269e-01 60.0000000000 +0.2791133725 5.4375798269e-01 60.0000000000 +0.2795763354 5.4999181937e-01 63.0000000000 +0.2800392984 5.4999181937e-01 63.0000000000 +0.2805022613 5.5610243787e-01 66.0000000000 +0.2809652243 5.5610243787e-01 66.0000000000 +0.2814281873 5.6217953277e-01 67.0000000000 +0.2818911502 5.6217953277e-01 67.0000000000 +0.2823541132 5.6776644650e-01 68.0000000000 +0.2828170761 5.6776644650e-01 68.0000000000 +0.2832800391 5.7324109432e-01 69.0000000000 +0.2837430021 5.7324109432e-01 69.0000000000 +0.2842059650 5.7840282562e-01 70.0000000000 +0.2846689280 5.7840282562e-01 70.0000000000 +0.2851318910 5.8343436332e-01 71.0000000000 +0.2855948539 5.8343436332e-01 71.0000000000 +0.2860578169 5.8841353961e-01 72.0000000000 +0.2865207798 5.8841353961e-01 72.0000000000 +0.2869837428 5.9333236123e-01 73.0000000000 +0.2874467058 5.9333236123e-01 73.0000000000 +0.2879096687 5.9824873032e-01 74.0000000000 +0.2883726317 5.9824873032e-01 74.0000000000 +0.2888355947 6.0302488617e-01 76.0000000000 +0.2892985576 6.0302488617e-01 76.0000000000 +0.2897615206 6.0740392712e-01 77.0000000000 +0.2902244835 6.0740392712e-01 77.0000000000 +0.2906874465 6.1160080141e-01 79.0000000000 +0.2911504095 6.1160080141e-01 79.0000000000 +0.2916133724 6.1560883158e-01 81.0000000000 +0.2920763354 6.1560883158e-01 81.0000000000 +0.2925392984 6.1947256253e-01 82.0000000000 +0.2930022613 6.1947256253e-01 82.0000000000 +0.2934652243 6.2316888804e-01 83.0000000000 +0.2939281872 6.2316888804e-01 83.0000000000 +0.2943911502 6.2686363669e-01 85.0000000000 +0.2948541132 6.2686363669e-01 85.0000000000 +0.2953170761 6.3823685779e-01 87.0000000000 +0.2957800391 6.3823685779e-01 87.0000000000 +0.2962430020 6.4191742451e-01 88.0000000000 +0.2967059650 6.4191742451e-01 88.0000000000 +0.2971689280 6.4551375960e-01 89.0000000000 +0.2976318909 6.4551375960e-01 89.0000000000 +0.2980948539 6.4897558496e-01 90.0000000000 +0.2985578169 6.4897558496e-01 90.0000000000 +0.2990207798 6.5239489739e-01 92.0000000000 +0.2994837428 6.5239489739e-01 92.0000000000 +0.2999467057 6.5577687146e-01 95.0000000000 +0.3004096687 6.5577687146e-01 95.0000000000 +0.3008726317 6.5902117817e-01 96.0000000000 +0.3013355946 6.5902117817e-01 96.0000000000 +0.3017985576 6.6217957455e-01 97.0000000000 +0.3022615206 6.7094960561e-01 99.0000000000 +0.3027244835 6.7094960561e-01 99.0000000000 +0.3031874465 6.7094960561e-01 99.0000000000 +0.3036504094 6.7400479976e-01 100.0000000000 +0.3041133724 6.7400479976e-01 100.0000000000 +0.3045763354 6.7705781706e-01 101.0000000000 +0.3050392983 6.7705781706e-01 101.0000000000 +0.3055022613 6.8042815428e-01 102.0000000000 +0.3059652243 6.8042815428e-01 102.0000000000 +0.3064281872 6.8334926480e-01 103.0000000000 +0.3068911502 6.8334926480e-01 103.0000000000 +0.3073541131 6.8619807839e-01 104.0000000000 +0.3078170761 6.9654503494e-01 106.0000000000 +0.3082800391 6.9654503494e-01 106.0000000000 +0.3087430020 6.9654503494e-01 106.0000000000 +0.3092059650 6.9937534078e-01 108.0000000000 +0.3096689279 6.9937534078e-01 108.0000000000 +0.3101318909 7.0216111232e-01 109.0000000000 +0.3105948539 7.0216111232e-01 109.0000000000 +0.3110578168 7.0493772287e-01 112.0000000000 +0.3115207798 7.0493772287e-01 112.0000000000 +0.3119837428 7.0762653760e-01 114.0000000000 +0.3124467057 7.0762653760e-01 114.0000000000 +0.3129096687 7.1030594699e-01 117.0000000000 +0.3133726316 7.1030594699e-01 117.0000000000 +0.3138355946 7.1298075991e-01 118.0000000000 +0.3142985576 7.1298075991e-01 118.0000000000 +0.3147615205 7.1562815243e-01 123.0000000000 +0.3152244835 7.1562815243e-01 123.0000000000 +0.3156874465 7.1820933056e-01 124.0000000000 +0.3161504094 7.1820933056e-01 124.0000000000 +0.3166133724 7.2076549767e-01 125.0000000000 +0.3170763353 7.2076549767e-01 125.0000000000 +0.3175392983 7.2326680412e-01 126.0000000000 +0.3180022613 7.2326680412e-01 126.0000000000 +0.3184652242 7.2572687663e-01 127.0000000000 +0.3189281872 7.2572687663e-01 127.0000000000 +0.3193911502 7.2815242027e-01 128.0000000000 +0.3198541131 7.2815242027e-01 128.0000000000 +0.3203170761 7.3036897080e-01 129.0000000000 +0.3207800390 7.3036897080e-01 129.0000000000 +0.3212430020 7.3255136712e-01 132.0000000000 +0.3217059650 7.3255136712e-01 132.0000000000 +0.3221689279 7.3471328106e-01 133.0000000000 +0.3226318909 7.3471328106e-01 133.0000000000 +0.3230948538 7.3684010554e-01 134.0000000000 +0.3235578168 7.3684010554e-01 134.0000000000 +0.3240207798 7.3890667989e-01 135.0000000000 +0.3244837427 7.3890667989e-01 135.0000000000 +0.3249467057 7.4084799328e-01 136.0000000000 +0.3254096687 7.4084799328e-01 136.0000000000 +0.3258726316 7.4278118267e-01 138.0000000000 +0.3263355946 7.4278118267e-01 138.0000000000 +0.3267985575 7.4464226626e-01 140.0000000000 +0.3272615205 7.4464226626e-01 140.0000000000 +0.3277244835 7.4646953023e-01 141.0000000000 +0.3281874464 7.4646953023e-01 141.0000000000 +0.3286504094 7.4826691908e-01 142.0000000000 +0.3291133724 7.4826691908e-01 142.0000000000 +0.3295763353 7.5005782020e-01 145.0000000000 +0.3300392983 7.5005782020e-01 145.0000000000 +0.3305022612 7.5169893758e-01 147.0000000000 +0.3309652242 7.5169893758e-01 147.0000000000 +0.3314281872 7.5332970009e-01 148.0000000000 +0.3318911501 7.5332970009e-01 148.0000000000 +0.3323541131 7.5494876907e-01 149.0000000000 +0.3328170761 7.5494876907e-01 149.0000000000 +0.3332800390 7.5642644783e-01 150.0000000000 +0.3337430020 7.5642644783e-01 150.0000000000 +0.3342059649 7.5783852091e-01 151.0000000000 +0.3346689279 7.5783852091e-01 151.0000000000 +0.3351318909 7.5924889407e-01 152.0000000000 +0.3355948538 7.5924889407e-01 152.0000000000 +0.3360578168 7.6063373248e-01 154.0000000000 +0.3365207797 7.6063373248e-01 154.0000000000 +0.3369837427 7.6199820611e-01 155.0000000000 +0.3374467057 7.6199820611e-01 155.0000000000 +0.3379096686 7.6335542422e-01 156.0000000000 +0.3383726316 7.6335542422e-01 156.0000000000 +0.3388355946 7.6470006370e-01 157.0000000000 +0.3392985575 7.6470006370e-01 157.0000000000 +0.3397615205 7.6599198178e-01 158.0000000000 +0.3402244834 7.6599198178e-01 158.0000000000 +0.3406874464 7.6724247936e-01 160.0000000000 +0.3411504094 7.6724247936e-01 160.0000000000 +0.3416133723 7.6849015616e-01 161.0000000000 +0.3420763353 7.6849015616e-01 161.0000000000 +0.3425392983 7.6972715746e-01 162.0000000000 +0.3430022612 7.6972715746e-01 162.0000000000 +0.3434652242 7.7095374903e-01 163.0000000000 +0.3439281871 7.7095374903e-01 163.0000000000 +0.3443911501 7.7216249760e-01 164.0000000000 +0.3448541131 7.7216249760e-01 164.0000000000 +0.3453170760 7.7335962051e-01 166.0000000000 +0.3457800390 7.7335962051e-01 166.0000000000 +0.3462430019 7.7442807638e-01 167.0000000000 +0.3467059649 7.7442807638e-01 167.0000000000 +0.3471689279 7.7539880219e-01 168.0000000000 +0.3476318908 7.7539880219e-01 168.0000000000 +0.3480948538 7.7636778866e-01 169.0000000000 +0.3485578168 7.7636778866e-01 169.0000000000 +0.3490207797 7.7731641657e-01 171.0000000000 +0.3494837427 7.7731641657e-01 171.0000000000 +0.3499467056 7.7824877371e-01 174.0000000000 +0.3504096686 7.7914572537e-01 175.0000000000 +0.3508726316 7.7914572537e-01 175.0000000000 +0.3513355945 7.8002247097e-01 177.0000000000 +0.3517985575 7.8002247097e-01 177.0000000000 +0.3522615205 7.8086067125e-01 178.0000000000 +0.3527244834 7.8086067125e-01 178.0000000000 +0.3531874464 7.8163887630e-01 179.0000000000 +0.3536504093 7.8163887630e-01 179.0000000000 +0.3541133723 7.8240600336e-01 181.0000000000 +0.3545763353 7.8240600336e-01 181.0000000000 +0.3550392982 7.8314409989e-01 184.0000000000 +0.3555022612 7.8314409989e-01 184.0000000000 +0.3559652242 7.8383169099e-01 186.0000000000 +0.3564281871 7.8383169099e-01 186.0000000000 +0.3568911501 7.8450134577e-01 187.0000000000 +0.3573541130 7.8450134577e-01 187.0000000000 +0.3578170760 7.8512805047e-01 188.0000000000 +0.3582800390 7.8512805047e-01 188.0000000000 +0.3587430019 7.8573194540e-01 189.0000000000 +0.3592059649 7.8573194540e-01 189.0000000000 +0.3596689278 7.8633393031e-01 190.0000000000 +0.3601318908 7.8633393031e-01 190.0000000000 +0.3605948538 7.8690859179e-01 191.0000000000 +0.3610578167 7.8690859179e-01 191.0000000000 +0.3615207797 7.8747996801e-01 192.0000000000 +0.3619837427 7.8747996801e-01 192.0000000000 +0.3624467056 7.8803887419e-01 193.0000000000 +0.3629096686 7.8803887419e-01 193.0000000000 +0.3633726315 7.8859456514e-01 194.0000000000 +0.3638355945 7.8859456514e-01 194.0000000000 +0.3642985575 7.8913343056e-01 195.0000000000 +0.3647615204 7.8913343056e-01 195.0000000000 +0.3652244834 7.8966555492e-01 196.0000000000 +0.3656874464 7.8966555492e-01 196.0000000000 +0.3661504093 7.9017307661e-01 197.0000000000 +0.3666133723 7.9017307661e-01 197.0000000000 +0.3670763352 7.9065559438e-01 198.0000000000 +0.3675392982 7.9065559438e-01 198.0000000000 +0.4055022611 7.9322676687e-01 199.0000000000 +0.4059652241 7.9322676687e-01 199.0000000000 +0.4087430018 8.0045415963e-01 201.0000000000 +0.4092059648 8.0045415963e-01 201.0000000000 +0.4096689278 8.0514882511e-01 202.0000000000 +0.4101318907 8.1209600019e-01 204.0000000000 +0.4105948537 8.1209600019e-01 204.0000000000 +0.4110578166 8.1737956994e-01 206.0000000000 +0.4115207796 8.1737956994e-01 206.0000000000 +0.4119837426 8.1737956994e-01 206.0000000000 +0.4156874463 8.1858761664e-01 209.0000000000 +0.4161504092 8.1858761664e-01 209.0000000000 +0.4166133722 8.2640142278e-01 211.0000000000 +0.4170763351 8.2640142278e-01 211.0000000000 +0.4175392981 8.2754028714e-01 212.0000000000 +0.4180022611 8.2754028714e-01 212.0000000000 +0.4184652240 8.2915417858e-01 214.0000000000 +0.4189281870 8.2915417858e-01 214.0000000000 +0.4193911500 8.3059096299e-01 217.0000000000 +0.4198541129 8.3059096299e-01 217.0000000000 +0.4203170759 8.3147701529e-01 220.0000000000 +0.4207800388 8.3147701529e-01 220.0000000000 +0.4397615203 8.3272365494e-01 221.0000000000 +0.4402244832 8.3272365494e-01 221.0000000000 +0.5892985570 8.3328560006e-01 222.0000000000 +0.5897615200 8.3328560006e-01 222.0000000000 +0.5957800385 8.3457630036e-01 223.0000000000 +0.5962430015 8.3457630036e-01 223.0000000000 +0.5980948533 8.3521454775e-01 226.0000000000 +0.5985578163 8.3615220309e-01 227.0000000000 +0.5990207792 8.4486928271e-01 228.0000000000 +0.5994837422 8.4486928271e-01 228.0000000000 +0.5999467052 8.5405085425e-01 230.0000000000 +0.6004096681 8.5405085425e-01 230.0000000000 +0.6008726311 8.5604472035e-01 231.0000000000 +0.6013355940 8.5604472035e-01 231.0000000000 +0.6017985570 8.5771859238e-01 232.0000000000 +0.6022615200 8.5771859238e-01 232.0000000000 +0.6027244829 8.5771859238e-01 232.0000000000 +0.6031874459 8.5839360194e-01 235.0000000000 +0.6036504088 8.5839360194e-01 235.0000000000 +0.6041133718 8.5839360194e-01 235.0000000000 +0.6064281866 8.5899444272e-01 236.0000000000 +0.6068911496 8.5899444272e-01 236.0000000000 +0.6096689274 8.5963991514e-01 237.0000000000 +0.6101318903 8.5963991514e-01 237.0000000000 +0.6110578162 8.6025000685e-01 238.0000000000 +0.6115207792 8.6025000685e-01 238.0000000000 +0.6119837422 8.6248995137e-01 239.0000000000 +0.6124467051 8.6248995137e-01 239.0000000000 diff --git a/gwemopt/tests/data/expected_results/tiles_coverage_int_ZTF.txt b/gwemopt/tests/data/expected_results/tiles_coverage_int_ZTF.txt index 146eb6c9..0e75ae64 100644 --- a/gwemopt/tests/data/expected_results/tiles_coverage_int_ZTF.txt +++ b/gwemopt/tests/data/expected_results/tiles_coverage_int_ZTF.txt @@ -1,4 +1,4 @@ -0.4709610778 3.4300508142e-01 46.4234297507 -0.4913314481 7.4957276454e-01 89.5945966263 -0.4920258926 7.4957276454e-01 89.5945966263 -0.4924888555 7.4957276454e-01 89.5945966263 +0.4709610778 3.9959847921e-01 51.9181803342 +0.4913314481 8.3231116217e-01 99.3743346875 +0.4920258926 8.3231116217e-01 99.3743346875 +0.4924888555 8.3231116217e-01 99.3743346875 diff --git a/gwemopt/tests/test_schedule.py b/gwemopt/tests/test_schedule.py index d3de778d..ef957406 100644 --- a/gwemopt/tests/test_schedule.py +++ b/gwemopt/tests/test_schedule.py @@ -38,7 +38,6 @@ def test_scheduler(): "--doCoverage", "--coverageFiles", os.path.join(temp_dir, "coverage_ZTF.dat"), - "--doObservability", "--doMovie", ], ), @@ -53,7 +52,7 @@ def test_scheduler(): "GLADE", ], ), - ("DECam", ["--doChipGaps", "--max_nb_tiles", "5", "--doMinimalTiling"]), + ("DECam", ["--max_nb_tiles", "5", "--doMinimalTiling"]), # ('TRE', []), # ("WINTER", []), # ('TNT', ["--tilesType", "galaxy", "--powerlaw_dist_exp", "1.0", "--doCatalog", "--galaxy_grade", "Sloc"]), @@ -111,7 +110,6 @@ def test_scheduler(): # Test the extra efficiency/coverage files extra_test_files = [ - "map.dat", "summary.dat", # "efficiency.txt", # "efficiency_tophat.txt", diff --git a/gwemopt/tiles.py b/gwemopt/tiles.py index a9aaa7cb..0616fb1b 100644 --- a/gwemopt/tiles.py +++ b/gwemopt/tiles.py @@ -1,12 +1,15 @@ +import bisect import copy +import astropy_healpix as ah import healpy as hp import ligo.segments as segments import numpy as np import pandas as pd from astropy import units as u -from astropy.coordinates import SkyCoord +from astropy.coordinates import ICRS, SkyCoord from astropy.time import Time +from mocpy import MOC from scipy.stats import norm from shapely.geometry import MultiPoint from tqdm import tqdm @@ -16,6 +19,10 @@ import gwemopt.segments from gwemopt.utils.geometry import angular_distance +LEVEL = MOC.MAX_ORDER +HPX = ah.HEALPix(nside=ah.level_to_nside(LEVEL), order="nested", frame=ICRS()) +PIXEL_AREA = HPX.pixel_area.to_value(u.sr) + TILE_TYPES = ["moc", "galaxy"] @@ -554,14 +561,9 @@ def perturb_tiles(params, config_struct, telescope, map_struct, tile_struct): map_struct_hold["prob"][moc_struct[key]["ipix"]] = -1 ipix_keep = np.setdiff1d(ipix_keep, moc_struct[key]["ipix"]) - if params["timeallocationType"] == "absmag": - tile_struct = absmag_tiles_struct( - params, config_struct, telescope, map_struct, moc_struct - ) - elif params["timeallocationType"] == "powerlaw": - tile_struct = powerlaw_tiles_struct( - params, config_struct, telescope, map_struct, moc_struct - ) + tile_struct = powerlaw_tiles_struct( + params, config_struct, telescope, map_struct, moc_struct + ) tile_struct = gwemopt.segments.get_segments_tiles( params, config_struct, tile_struct ) @@ -616,19 +618,8 @@ def schedule_alternating( ) if not params["tilesType"] == "galaxy": - if params["timeallocationType"] == "absmag": - tile_struct = absmag_tiles_struct( - params, config_struct, telescope, map_struct, tile_struct - ) - else: - tile_struct = powerlaw_tiles_struct( - params, config_struct, telescope, map_struct, tile_struct - ) - - if params["treasuremap_token"] is not None and previous_coverage_struct: - # erases tiles from a previous round - tile_struct = gwemopt.coverage.update_observed_tiles( - params, tile_struct_hold, previous_coverage_struct + tile_struct = powerlaw_tiles_struct( + params, config_struct, telescope, map_struct, tile_struct ) # set unbalanced fields to 0 @@ -849,21 +840,14 @@ def galaxy(params, map_struct, catalog_struct: pd.DataFrame): moc_struct[cnt]["galaxies"] = row["galaxies"] cnt = cnt + 1 - if params["timeallocationType"] == "absmag": - tile_struct = absmag_tiles_struct( - params, config_struct, telescope, map_struct, moc_struct - ) - elif params["timeallocationType"] == "powerlaw": - tile_struct = powerlaw_tiles_struct( - params, - config_struct, - telescope, - map_struct, - moc_struct, - catalog_struct=catalog_struct, - ) - else: - raise ValueError("timeallocationType not recognized") + tile_struct = powerlaw_tiles_struct( + params, + config_struct, + telescope, + map_struct, + moc_struct, + catalog_struct=catalog_struct, + ) tile_struct = gwemopt.segments.get_segments_tiles( params, config_struct, tile_struct @@ -890,176 +874,6 @@ def galaxy(params, map_struct, catalog_struct: pd.DataFrame): return tile_structs -def absmag_tiles_struct(params, config_struct, telescope, map_struct, tile_struct): - keys = tile_struct.keys() - ntiles = len(keys) - if ntiles == 0: - return tile_struct - - if "observability" in map_struct: - prob = map_struct["observability"][telescope]["prob"] - else: - prob = map_struct["prob"] - - n, cl, dist_exp = ( - params["powerlaw_n"], - params["powerlaw_cl"], - params["powerlaw_dist_exp"], - ) - - if params["tilesType"] == "galaxy": - tile_probs = compute_tiles_map( - params, tile_struct, prob, func="center", ipix_keep=map_struct["ipix_keep"] - ) - else: - tile_probs = compute_tiles_map( - params, - tile_struct, - prob, - func="np.sum(x)", - ipix_keep=map_struct["ipix_keep"], - ) - - tile_probs[tile_probs < np.max(tile_probs) * 0.01] = 0.0 - - prob_scaled = copy.deepcopy(prob) - prob_sorted = np.sort(prob_scaled)[::-1] - prob_indexes = np.argsort(prob_scaled)[::-1] - prob_cumsum = np.cumsum(prob_sorted) - index = np.argmin(np.abs(prob_cumsum - cl)) + 1 - # prob_scaled[prob_indexes[index:]] = 1e-10 - prob_scaled[prob_indexes[index:]] = 0.0 - prob_scaled = prob_scaled**n - prob_scaled = prob_scaled / np.nansum(prob_scaled) - - if params["tilesType"] == "galaxy": - ranked_tile_probs = compute_tiles_map( - params, - tile_struct, - prob_scaled, - func="center", - ipix_keep=map_struct["ipix_keep"], - ) - else: - ranked_tile_probs = compute_tiles_map( - params, - tile_struct, - prob_scaled, - func="np.sum(x)", - ipix_keep=map_struct["ipix_keep"], - ) - - ranked_tile_probs[np.isnan(ranked_tile_probs)] = 0.0 - ranked_tile_probs_thresh = np.max(ranked_tile_probs) * 0.01 - ranked_tile_probs[ranked_tile_probs <= ranked_tile_probs_thresh] = 0.0 - ranked_tile_probs = ranked_tile_probs / np.nansum(ranked_tile_probs) - - if "distmed" in map_struct: - distmed = map_struct["distmed"] + 2 * map_struct["diststd"] - distmed[distmed <= 0] = np.nan - distmed[~np.isfinite(distmed)] = np.nan - # distmed[distmed 0)[0] - ranked_tile_times[idx, ii] = config_struct["exposuretime_orig"] * nexp[idx] - if "max_exposure" in config_struct.keys(): - idy = np.where(ranked_tile_times[idx, ii] > config_struct["max_exposure"])[ - 0 - ] - ranked_tile_times[idx[idy], ii] = config_struct["max_exposure"] - - for key, prob, exposureTime, tileprob in zip( - keys, ranked_tile_probs, ranked_tile_times, tile_probs - ): - # Try to load the minimum duration of time from telescope config file - # Otherwise set it to zero - try: - min_obs_duration = config_struct["min_observability_duration"] / 24 - except: - min_obs_duration = 0.0 - - # Check that a given tile is observable a minimum amount of time - # If not set the proba associated to the tile to zero - if ( - "segmentlist" - and "prob" in tile_struct[key] - and tile_struct[key]["segmentlist"] - and min_obs_duration > 0.0 - ): - observability_duration = 0.0 - for counter in range(len(tile_struct[key]["segmentlist"])): - observability_duration += ( - tile_struct[key]["segmentlist"][counter][1] - - tile_struct[key]["segmentlist"][counter][0] - ) - if ( - tile_struct[key]["prob"] > 0.0 - and observability_duration < min_obs_duration - ): - tileprob = 0.0 - - tile_struct[key]["prob"] = tileprob - if prob == 0.0: - tile_struct[key]["exposureTime"] = 0.0 - tile_struct[key]["nexposures"] = 0 - tile_struct[key]["filt"] = [] - else: - if params["doReferences"] and (telescope in ["ZTF", "DECam"]): - tile_struct[key]["exposureTime"] = [] - tile_struct[key]["nexposures"] = [] - tile_struct[key]["filt"] = [] - if key in config_struct["reference_images"]: - for ii in range(len(params["filters"])): - if ( - params["filters"][ii] - in config_struct["reference_images"][key] - ): - tile_struct[key]["exposureTime"].append(exposureTime[ii]) - tile_struct[key]["filt"].append(params["filters"][ii]) - tile_struct[key]["nexposures"] = len( - tile_struct[key]["exposureTime"] - ) - else: - tile_struct[key]["exposureTime"] = 0.0 - tile_struct[key]["nexposures"] = 0 - tile_struct[key]["filt"] = [] - else: - tile_struct[key]["exposureTime"] = exposureTime - tile_struct[key]["nexposures"] = len(params["exposuretimes"]) - tile_struct[key]["filt"] = params["filters"] - - return tile_struct - - def powerlaw_tiles_struct( params, config_struct, telescope, map_struct, tile_struct, catalog_struct=None ): @@ -1070,11 +884,6 @@ def powerlaw_tiles_struct( tot_obs_time = config_struct["tot_obs_time"] - if "observability" in map_struct: - prob = map_struct["observability"][telescope]["prob"] - else: - prob = map_struct["prob"] - n, cl, dist_exp = ( params["powerlaw_n"], params["powerlaw_cl"], @@ -1085,96 +894,29 @@ def powerlaw_tiles_struct( tile_probs = compute_tiles_map( params, tile_struct, - prob, + map_struct["skymap"], func="galaxy", - ipix_keep=map_struct["ipix_keep"], + moc_keep=map_struct["moc_keep"], catalog_struct=catalog_struct, ) else: tile_probs = compute_tiles_map( params, tile_struct, - prob, + map_struct["skymap"], func="np.sum(x)", - ipix_keep=map_struct["ipix_keep"], + moc_keep=map_struct["moc_keep"], ) tile_probs[tile_probs < np.max(tile_probs) * 0.01] = 0.0 - prob_scaled = copy.deepcopy(prob) - prob_sorted = np.sort(prob_scaled)[::-1] - prob_indexes = np.argsort(prob_scaled)[::-1] - prob_cumsum = np.cumsum(prob_sorted) - index = np.argmin(np.abs(prob_cumsum - cl)) + 1 - # prob_scaled[prob_indexes[index:]] = 1e-10 - prob_scaled[prob_indexes[index:]] = 0.0 - prob_scaled = prob_scaled**n - prob_scaled = prob_scaled / np.nansum(prob_scaled) - - if params["tilesType"] == "galaxy": - ranked_tile_probs = compute_tiles_map( - params, - tile_struct, - prob_scaled, - func="galaxy", - ipix_keep=map_struct["ipix_keep"], - catalog_struct=catalog_struct, - ) - else: - ranked_tile_probs = compute_tiles_map( - params, - tile_struct, - prob_scaled, - func="np.sum(x)", - ipix_keep=map_struct["ipix_keep"], - ) - - ranked_tile_probs[np.isnan(ranked_tile_probs)] = 0.0 - ranked_tile_probs_thresh = np.max(ranked_tile_probs) * 0.01 - ranked_tile_probs[ranked_tile_probs <= ranked_tile_probs_thresh] = 0.0 - ranked_tile_probs = ranked_tile_probs / np.nansum(ranked_tile_probs) - - if "distmed" in map_struct: - distmed = map_struct["distmed"] - distmed[distmed <= 0] = np.nan - distmed[~np.isfinite(distmed)] = np.nan - # distmed[distmed 0, ii] = params["exposuretimes"][ii] + ranked_tile_times[tile_probs > 0, ii] = params["exposuretimes"][ii] - for key, prob, exposureTime, tileprob in zip( - keys, ranked_tile_probs, ranked_tile_times, tile_probs - ): + for key, tileprob, exposureTime in zip(keys, tile_probs, ranked_tile_times): # Try to load the minimum duration of time from telescope config file # Otherwise set it to zero try: @@ -1203,7 +945,7 @@ def powerlaw_tiles_struct( tileprob = 0.0 tile_struct[key]["prob"] = tileprob - if prob == 0.0: + if tileprob == 0.0: tile_struct[key]["exposureTime"] = 0.0 tile_struct[key]["nexposures"] = 0 tile_struct[key]["filt"] = [] @@ -1236,16 +978,14 @@ def powerlaw_tiles_struct( else: ranked_tile_times = gwemopt.utils.integrationTime( tot_obs_time, - ranked_tile_probs, + tile_probs, func=None, T_int=config_struct["exposuretime"], ) keys = tile_struct.keys() - for key, prob, exposureTime, tileprob in zip( - keys, ranked_tile_probs, ranked_tile_times, tile_probs - ): + for key, tileprob, exposureTime in zip(keys, tile_probs, ranked_tile_times): # Try to load the minimum duration of time from telescope config file # Otherwise set it to zero try: @@ -1271,9 +1011,9 @@ def powerlaw_tiles_struct( tile_struct[key]["prob"] > 0.0 and observability_duration < min_obs_duration ): - prob = 0.0 + tileprob = 0.0 - tile_struct[key]["prob"] = prob + tile_struct[key]["prob"] = tileprob tile_struct[key]["exposureTime"] = exposureTime tile_struct[key]["nexposures"] = int( np.floor(exposureTime / config_struct["exposuretime"]) @@ -1291,14 +1031,9 @@ def moc(params, map_struct, moc_structs, doSegments=True): config_struct = params["config"][telescope] moc_struct = moc_structs[telescope] - if params["timeallocationType"] == "absmag": - tile_struct = gwemopt.tiles.absmag_tiles_struct( - params, config_struct, telescope, map_struct, moc_struct - ) - elif params["timeallocationType"] == "powerlaw": - tile_struct = powerlaw_tiles_struct( - params, config_struct, telescope, map_struct, moc_struct - ) + tile_struct = powerlaw_tiles_struct( + params, config_struct, telescope, map_struct, moc_struct + ) if doSegments: tile_struct = gwemopt.segments.get_segments_tiles( @@ -1309,90 +1044,13 @@ def moc(params, map_struct, moc_structs, doSegments=True): return tile_structs -def pem_tiles_struct(params, config_struct, telescope, map_struct, tile_struct): - tot_obs_time = config_struct["tot_obs_time"] - - if "observability" in map_struct: - prob = map_struct["observability"][telescope]["prob"] - else: - prob = map_struct["prob"] - - ranked_tile_probs = compute_tiles_map( - params, tile_struct, prob, func="np.sum(x)", ipix_keep=map_struct["ipix_keep"] - ) - - lim_mag = config_struct["magnitude"] - lim_time = config_struct["exposuretime"] - - tau = np.arange(lim_time, 3600.0, lim_time) - Loftau = None - - N_ref = 9.7847e9 - L_min = 4.9370e31 - L_max = 4.9370e33 - - if "distmu" in map_struct: - prob = map_struct["prob"] - distnorm = map_struct["distnorm"] - distmu = map_struct["distmu"] - distsigma = map_struct["distsigma"] - - D_min = 1.0e7 - D_max = 1.0e10 - R = np.linspace(D_min / 1e6, D_max / 1e6) - p_R = [ - np.sum(prob * rr**2 * distnorm * norm(distmu, distsigma).pdf(rr)) - for rr in R - ] - p_R = np.array(p_R) - p_R = p_R / len(prob) - - R = R * 1e6 - D_mu = None - D_sig = None - else: - D_mu = 200.0e6 - D_sig = 60.0e6 - R = None - p_R = None - - tau, prob = gwemopt.pem.Pem( - lim_mag, - lim_time, - N_ref=N_ref, - L_min=L_min, - L_max=L_max, - tau=tau, - Loftau=Loftau, - D_mu=D_mu, - D_sig=D_sig, - R=R, - p_R=p_R, - ) - - tprob, time_allocation = gwemopt.pem.Main( - tot_obs_time, 0, 0, ranked_tile_probs, tau, prob, "Eq" - ) - - if params["doPlots"]: - gwemopt.plotting.tauprob(params, tau, prob) - - keys = tile_struct.keys() - for key, prob, exposureTime in zip(keys, ranked_tile_probs, time_allocation): - tile_struct[key]["prob"] = prob - tile_struct[key]["exposureTime"] = exposureTime - tile_struct[key]["nexposures"] = int( - np.floor(exposureTime / config_struct["exposuretime"]) - ) - tile_struct[key]["filt"] = [config_struct["filt"]] * tile_struct[key][ - "nexposures" - ] - - return tile_struct - - def compute_tiles_map( - params, tile_struct, skymap, func=None, ipix_keep=[], catalog_struct=None + params, + tile_struct, + skymap, + func=None, + moc_keep=MOC.new_empty(max_depth=10), + catalog_struct=None, ): if func is None: f = lambda x: np.sum(x) @@ -1412,7 +1070,6 @@ def compute_tiles_map( keys = tile_struct.keys() ntiles = len(keys) vals = np.nan * np.ones((ntiles,)) - nside = hp.npix2nside(len(skymap)) for ii, key in enumerate(tile_struct.keys()): galaxies = tile_struct[key]["galaxies"] val = np.sum(catalog_struct[params["galaxy_grade"]][galaxies]) @@ -1426,31 +1083,8 @@ def compute_tiles_map( keys = tile_struct.keys() ntiles = len(keys) vals = np.nan * np.ones((ntiles,)) - for ii, key in enumerate(tile_struct.keys()): - idx = np.where(prob[tile_struct[key]["ipix"]] == -1)[0] - idx = np.setdiff1d(idx, ipix_keep) - if len(prob[tile_struct[key]["ipix"]]) == 0: - rat = 0.0 - else: - rat = float(len(idx)) / float(len(prob[tile_struct[key]["ipix"]])) - if rat > params["maximumOverlap"]: - vals[ii] = 0.0 - else: - ipix = tile_struct[key]["ipix"] - if len(ipix) == 0: - vals[ii] = 0.0 - else: - vals_to_sum = prob[ipix] - if func == "np.nanmedian(x)": - vals_to_sum[vals_to_sum < 0] = np.nan - else: - vals_to_sum[vals_to_sum < 0] = 0 - vals[ii] = f(vals_to_sum) - - ipix_keep = np.setdiff1d(ipix_keep, tile_struct[key]["ipix"]) - ipix_slice = np.setdiff1d(tile_struct[key]["ipix"], ipix_keep) - if len(ipix_slice) > 0: - prob[ipix_slice] = -1 + for ii, key in tqdm(enumerate(keys), total=len(keys)): + vals[ii] = tile_struct[key]["moc"].probability_in_multiordermap(skymap) return vals diff --git a/gwemopt/utils/pixels.py b/gwemopt/utils/pixels.py index 415a7eb9..c57e38ac 100644 --- a/gwemopt/utils/pixels.py +++ b/gwemopt/utils/pixels.py @@ -43,7 +43,6 @@ def getRegionPixels( alpha=0.4, color="k", edgecolor="k", - rotation=None, ): theta = 0.5 * np.pi - np.deg2rad(dec_pointing) phi = np.deg2rad(ra_pointing) @@ -54,12 +53,9 @@ def getRegionPixels( ra_pointing, dec_pointing, unit=u.deg ).skyoffset_frame() - proj = hp.projector.MollweideProj(rot=rotation, coord=None) + moc = None # security for the periodic limit conditions - ipix = [] - radecs = np.empty((0, 2)) - patch = [] for reg in regions: ra_tmp = reg.vertices.ra dec_tmp = reg.vertices.dec @@ -75,55 +71,22 @@ def getRegionPixels( ).transform_to(coordinates.ICRS) coords = coords_icrs.transform_to(HPX.frame) - moc = MOC.from_polygon_skycoord(coords) + if moc is None: + moc = MOC.from_polygon_skycoord(coords) + else: + moc = moc + MOC.from_polygon_skycoord(coords) - pix_id = moc.degrade_to_order(int(np.log2(nside))).flatten() - if len(pix_id) > 0: - ipix.extend(hp.nest2ring(int(nside), pix_id.tolist())) - - ras = [coord.ra.deg for coord in coords] - decs = [coord.dec.deg for coord in coords] - coords = np.vstack([ras, decs]).T - - xyz = hp.ang2vec(coords[:, 0], coords[:, 1], lonlat=True) - x, y = proj.vec2xy(xyz[:, 0], xyz[:, 1], xyz[:, 2]) - xy = np.vstack([x, y]).T - path = matplotlib.path.Path(xy) - - radecs = np.vstack([radecs, coords]) - - patch.append( - matplotlib.patches.PathPatch( - path, - alpha=alpha, - facecolor=color, - fill=True, - zorder=3, - edgecolor=edgecolor, - ) - ) - - ipix = list(set(ipix)) - area = len(ipix) * hp.nside2pixarea(nside, degrees=True) - - return ipix, radecs, patch, area + return moc def getCirclePixels( ra_pointing, dec_pointing, radius, - nside, alpha=0.4, color="k", edgecolor="k", - rotation=None, ): - theta = 0.5 * np.pi - np.deg2rad(dec_pointing) - phi = np.deg2rad(ra_pointing) - - xyz = hp.ang2vec(theta, phi) - ipix = hp.query_disc(nside, xyz, np.deg2rad(radius)) radecs = get_ellipse_coords( a=radius / np.cos(np.deg2rad(dec_pointing)), @@ -143,29 +106,10 @@ def getCirclePixels( radecs[idx, 0] = 360.0 + radecs[idx, 0] radecs = np.array(radecs) - idx1 = np.where(radecs[:, 0] >= 180.0)[0] - idx2 = np.where(radecs[:, 0] < 180.0)[0] - idx3 = np.where(radecs[:, 0] > 300.0)[0] - idx4 = np.where(radecs[:, 0] < 60.0)[0] - if (len(idx1) > 0 and len(idx2) > 0) and not (len(idx3) > 0 and len(idx4) > 0): - alpha = 0.0 + coords = coordinates.SkyCoord(radecs[:, 0] * u.deg, radecs[:, 1] * u.deg) + moc = MOC.from_polygon_skycoord(coords) - xyz = hp.ang2vec(radecs[:, 0], radecs[:, 1], lonlat=True) - - proj = hp.projector.MollweideProj(rot=rotation, coord=None) - x, y = proj.vec2xy(xyz[:, 0], xyz[:, 1], xyz[:, 2]) - xy = np.zeros(radecs.shape) - xy[:, 0] = x - xy[:, 1] = y - # path = matplotlib.path.Path(xyz[:,1:3]) - path = matplotlib.path.Path(xy) - patch = matplotlib.patches.PathPatch( - path, alpha=alpha, facecolor=color, fill=True, zorder=3, edgecolor=edgecolor - ) - - area = np.pi * radius**2 - - return ipix, radecs, patch, area + return moc def getRectanglePixels( @@ -173,11 +117,9 @@ def getRectanglePixels( dec_pointing, raSide, decSide, - nside, alpha=0.4, color="k", edgecolor="k", - rotation=None, ): area = raSide * decSide @@ -213,46 +155,30 @@ def getRectanglePixels( idx1 = np.where(np.abs(radecs[:, 1]) >= 87.0)[0] if len(idx1) == 4: - return [], [], [], [] + return MOC.new_empty(max_depth=10) idx1 = np.where((radecs[:, 1] >= 87.0) | (radecs[:, 1] <= -87.0))[0] if len(idx1) > 0: radecs = np.delete(radecs, idx1[0], 0) - xyz = [] - for r, d in radecs: - xyz.append(hp.ang2vec(r, d, lonlat=True)) - npts, junk = radecs.shape if npts == 4: - xyz = [xyz[0], xyz[1], xyz[3], xyz[2]] - ipix = hp.query_polygon(nside, np.array(xyz)) - else: - ipix = hp.query_polygon(nside, np.array(xyz)) - - xyz = np.array(xyz) - proj = hp.projector.MollweideProj(rot=rotation, coord=None) - x, y = proj.vec2xy(xyz[:, 0], xyz[:, 1], xyz[:, 2]) - xy = np.zeros(radecs.shape) - xy[:, 0] = x - xy[:, 1] = y - path = matplotlib.path.Path(xy) - patch = matplotlib.patches.PathPatch( - path, alpha=alpha, facecolor=color, fill=True, zorder=3, edgecolor=edgecolor - ) + idx = [0, 1, 3, 2] + radecs = radecs[idx, :] - return ipix, radecs, patch, area + coords = coordinates.SkyCoord(radecs[:, 0] * u.deg, radecs[:, 1] * u.deg) + moc = MOC.from_polygon_skycoord(coords) + + return moc def getSquarePixels( ra_pointing, dec_pointing, tileSide, - nside, alpha=0.4, color="k", edgecolor="k", - rotation=None, ): area = tileSide * tileSide @@ -289,32 +215,18 @@ def getSquarePixels( idx1 = np.where(np.abs(radecs[:, 1]) >= 87.0)[0] if len(idx1) == 4: - return [], [], [], [] + return MOC.new_empty(max_depth=10) idx1 = np.where((radecs[:, 1] >= 87.0) | (radecs[:, 1] <= -87.0))[0] if len(idx1) > 0: radecs = np.delete(radecs, idx1[0], 0) - xyz = [] - for r, d in radecs: - xyz.append(hp.ang2vec(r, d, lonlat=True)) - npts, junk = radecs.shape if npts == 4: - xyz = [xyz[0], xyz[1], xyz[3], xyz[2]] - ipix = hp.query_polygon(nside, np.array(xyz)) - else: - ipix = hp.query_polygon(nside, np.array(xyz)) - - xyz = np.array(xyz) - proj = hp.projector.MollweideProj(rot=rotation, coord=None) - x, y = proj.vec2xy(xyz[:, 0], xyz[:, 1], xyz[:, 2]) - xy = np.zeros(radecs.shape) - xy[:, 0] = x - xy[:, 1] = y - path = matplotlib.path.Path(xy) - patch = matplotlib.patches.PathPatch( - path, alpha=alpha, facecolor=color, fill=True, zorder=3, edgecolor=edgecolor - ) + idx = [0, 1, 3, 2] + radecs = radecs[idx, :] + + coords = coordinates.SkyCoord(radecs[:, 0] * u.deg, radecs[:, 1] * u.deg) + moc = MOC.from_polygon_skycoord(coords) - return ipix, radecs, patch, area + return moc diff --git a/gwemopt/utils/treasuremap.py b/gwemopt/utils/treasuremap.py index 4b2e0aea..928e0b53 100644 --- a/gwemopt/utils/treasuremap.py +++ b/gwemopt/utils/treasuremap.py @@ -45,7 +45,7 @@ def get_treasuremap_pointings(params): coverage_struct = {} coverage_struct["data"] = np.empty((0, 8)) coverage_struct["filters"] = [] - coverage_struct["ipix"] = [] + coverage_struct["moc"] = [] # read information into coverage struct for obs in observations: @@ -60,7 +60,7 @@ def get_treasuremap_pointings(params): ra, dec = float(pointing[0]), float(pointing[1]) filteridx = obs.find("band") + 10 # jump to starting index of filter - filter = obs[filteridx:].split('"')[0][:-1] + filt = obs[filteridx:].split('"')[0][:-1] instrumentidx = ( obs.find("instrumentid") + 16 @@ -68,13 +68,9 @@ def get_treasuremap_pointings(params): instrument_id = int(obs[instrumentidx:].split(",")[0]) if instrument_id in FOV_square: - ipix, radecs, patch, area = getSquarePixels( - ra, dec, FOV_square[instrument_id], params["nside"] - ) + moc = getSquarePixels(ra, dec, FOV_square[instrument_id]) elif instrument_id in FOV_circle: - ipix, radecs, patch, area = getCirclePixels( - ra, dec, FOV_circle[instrument_id], params["nside"] - ) + moc = getCirclePixels(ra, dec, FOV_circle[instrument_id]) else: continue @@ -83,7 +79,7 @@ def get_treasuremap_pointings(params): np.array([[ra, dec, -1, -1, -1, -1, -1, -1]]), axis=0, ) - coverage_struct["filters"].append(filter) - coverage_struct["ipix"].append(ipix) + coverage_struct["filters"].append(filt) + coverage_struct["moc"].append(moc) return coverage_struct diff --git a/pyproject.toml b/pyproject.toml index 2ded42fe..aee2e193 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -35,7 +35,7 @@ dependencies = [ 'astropy>=1.1.1', 'astropy-healpix', 'python-dateutil', - 'mocpy>=0.12.2', + 'mocpy@git+https://github.com/cds-astro/mocpy', 'astroquery', 'astroplan', 'ephem',