From 6344ce1eb9fdfaf18e6775ca298f9e1ac08da2ee Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Tue, 24 Oct 2023 18:01:16 +0900 Subject: [PATCH 01/28] Refactoring the "run" submodule. --- src/srcsim/run/__init__.py | 1 + src/srcsim/{run.py => run/sky.py} | 0 2 files changed, 1 insertion(+) create mode 100644 src/srcsim/run/__init__.py rename src/srcsim/{run.py => run/sky.py} (100%) diff --git a/src/srcsim/run/__init__.py b/src/srcsim/run/__init__.py new file mode 100644 index 0000000..0cafdae --- /dev/null +++ b/src/srcsim/run/__init__.py @@ -0,0 +1 @@ +from .sky import DataRun \ No newline at end of file diff --git a/src/srcsim/run.py b/src/srcsim/run/sky.py similarity index 100% rename from src/srcsim/run.py rename to src/srcsim/run/sky.py From 1f5d4c65193a77b770100c987dd8ffd82a51fbe5 Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Tue, 24 Oct 2023 18:27:15 +0900 Subject: [PATCH 02/28] Fixed pointing run prototype. --- src/srcsim/run/fixed.py | 299 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 299 insertions(+) create mode 100644 src/srcsim/run/fixed.py diff --git a/src/srcsim/run/fixed.py b/src/srcsim/run/fixed.py new file mode 100644 index 0000000..5706c7a --- /dev/null +++ b/src/srcsim/run/fixed.py @@ -0,0 +1,299 @@ +import yaml +import numpy as np +import pandas as pd +import astropy.units as u +from astropy.time import Time +from astropy.coordinates import SkyCoord, SkyOffsetFrame, EarthLocation, AltAz + + +class DataRun: + def __init__(self, tel_pos, tstart, tstop, obsloc, id=0): + self.id = id + self.tel_pos = tel_pos + self.obsloc = obsloc + self.tstart = tstart + self.tstop = tstop + + @classmethod + def from_config(cls, config): + pass + + def to_dict(self): + pass + + def predict(self, mccollections, source, tel_pos_tolerance, time_step): + pass + + @classmethod + def time_sort(cls, events): + timestamp = np.array(events["trigger_time"]) + sorting = np.argsort(timestamp) + + return events.iloc[sorting] + + @classmethod + def update_time_delta(cls, events): + event_time = events['trigger_time'].to_numpy() + delta_time = np.zeros_like(event_time) + + if len(event_time) > 0: + argsorted = event_time.argsort() + dt = np.diff(event_time[argsorted]) + delta_time[argsorted[:-1]] = dt + delta_time[argsorted[-1]] = dt[-1] + + events = events.drop( + columns=['delta_t'], + errors='ignore' + ) + events = events.assign( + delta_t = delta_time, + ) + + return events + + +class FixedPointingDataRun(DataRun): + def __repr__(self): + frame_start = AltAz(obstime=self.tstart, location=self.obsloc) + frame_stop = AltAz(obstime=self.tstop, location=self.obsloc) + tel_pos_start = SkyCoord( + self.tel_pos.az, + self.tel_pos.alt, + frame=frame_start + ) + tel_pos_stop = SkyCoord( + self.tel_pos.az, + self.tel_pos.alt, + frame=frame_stop + ) + print( +f"""{type(self).__name__} instance + {'ID':.<20s}: {self.id} + {'Tel. alt/az':.<20s}: {self.tel_pos} + {'Tstart':.<20s}: {self.tstart.isot} + {'Tstop':.<20s}: {self.tstop.isot} + {'Tel. RA':.<20s}: [{tel_pos_start.icrs.ra.to('deg'):.2f} - {tel_pos_stop.icrs.ra.to('deg'):.2f}] + {'Tel. Dec':.<20s}: [{tel_pos_start.icrs.dec.to('deg'):.2f} - {tel_pos_stop.icrs.dec.to('deg'):.2f}] +""" + ) + + return super().__repr__() + + @classmethod + def from_config(cls, config): + if isinstance(config, str): + cfg = yaml.load(open(config, "r"), Loader=yaml.FullLoader) + else: + cfg = config + + dummy_obstime = Time('2000-01-01T00:00:00') + location = EarthLocation( + lat=u.Quantity(cfg['location']['lat']), + lon=u.Quantity(cfg['location']['lon']), + height=u.Quantity(cfg['location']['height']), + ) + + data_run = cls( + SkyCoord( + u.Quantity(cfg['pointing']['az']), + u.Quantity(cfg['pointing']['alt']), + frame='altaz', + location=location, + obstime=dummy_obstime + ), + Time(cfg['time']['start']), + Time(cfg['time']['stop']), + EarthLocation( + lat=u.Quantity(cfg['location']['lat']), + lon=u.Quantity(cfg['location']['lon']), + height=u.Quantity(cfg['location']['height']), + ), + cfg['id'] if 'id' in cfg else 0 + ) + + return data_run + + def to_dict(self): + data = {'id': self.id, 'pointing': {}, 'time': {}, 'location': {}} + + data['pointing']['alt'] = str(self.tel_pos.alt.to('deg').value) + ' deg' + data['pointing']['az'] = str(self.tel_pos.az.to('deg').value) + ' deg' + + data['time']['start'] = self.tstart.isot + data['time']['stop'] = self.tstop.isot + + data['location']['lon'] = str(self.obsloc.lon.to('deg').value) + ' deg' + data['location']['lat'] = str(self.obsloc.lat.to('deg').value) + ' deg' + data['location']['height'] = str(self.obsloc.height.to('m').to_string()) + + return data + + def tel_pos_to_altaz(self, frame): + if frame.obstime.size > 1: + tel_pos = SkyCoord( + np.repeat(self.tel_pos.az, frame.obstime.size), + np.repeat(self.tel_pos.alt, frame.obstime.size), + frame=frame + ) + else: + tel_pos = SkyCoord( + self.tel_pos.az, + self.tel_pos.alt, + frame=frame + ) + return tel_pos + + def predict(self, mccollections, source, tel_pos_tolerance=None, time_step=1*u.minute): + events = [] + + unix_edges = np.arange(self.tstart.unix, self.tstop.unix, step=time_step.to('s').value) + if unix_edges[-1] < self.tstop.unix: + unix_edges = np.concatenate((unix_edges, [self.tstop.unix])) + + tedges = Time(unix_edges, format='unix') + + if len(tedges) > 1: + tdelta = np.diff(tedges) + else: + tdelta = [self.tstop - self.tstart] + + for tstart, dt in zip(tedges[:-1], tdelta): + frame = AltAz( + obstime=tstart + dt/2, + location=self.obsloc + ) + tel_pos = self.tel_pos_to_altaz(frame) + + if tel_pos_tolerance is None: + mc = mccollections[source.emission_type].get_closest(tel_pos.altaz) + elif isinstance(tel_pos_tolerance, u.Quantity): + mc = mccollections[source.emission_type].get_nearby(tel_pos, tel_pos_tolerance) + elif isinstance(tel_pos_tolerance, list) or isinstance(tel_pos_tolerance, tuple): + mc = mccollections[source.emission_type].get_in_box( + tel_pos, + max_lon_offset=tel_pos_tolerance[0], + max_lat_offset=tel_pos_tolerance[1], + ) + else: + raise ValueError(f"Data type '{type(tel_pos_tolerance)}' for argument 'tel_pos_tolerance' is not supported") + + nsamples = len(mc.samples) + + for sample in mc.samples: + # Randomly distributing events within the time bin + arrival_time = np.random.uniform(tstart.unix, (tstart+dt).unix, size=len(sample.data_table)) + + # Recalculating telescope Alt/Az for the mock event arrival times + current_frame = AltAz( + obstime=Time(arrival_time, format='unix'), + location=self.obsloc + ) + current_tel_pos = self.tel_pos_to_altaz(current_frame) + + # Astropy does not pass the location / time + # to the offset frame, need to do this manually + offset_frame = SkyOffsetFrame( + origin=current_tel_pos.altaz.skyoffset_frame().origin, + location=current_tel_pos.altaz.frame.location, + obstime=current_tel_pos.altaz.frame.obstime + ) + coords = SkyCoord( + sample.evt_coord.skyoffsetaltaz.lon, + sample.evt_coord.skyoffsetaltaz.lat, + frame=offset_frame + ) + expected_flux = source.dndedo(sample.evt_energy, coords.icrs) + model_flux = sample.dndedo(sample.evt_energy, sample.evt_coord) + + weights = (1 / nsamples * dt * expected_flux / model_flux).decompose() + + n_mc_events = len(sample.evt_energy) + n_events = np.random.poisson(weights.sum()) + p = weights / weights.sum() + idx = np.random.choice( + np.arange(n_mc_events), + size=n_events, + p=p + ) + + evt = sample.data_table.iloc[idx] + offset_frame = offset_frame[idx] + arrival_time = arrival_time[idx] + current_tel_pos = current_tel_pos[idx] + + # Dropping the columns we're going to (re-)fill + evt = evt.drop( + columns=['dragon_time', 'trigger_time'], + errors='ignore' + ) + evt = evt.drop( + columns=['mc_az_tel', 'mc_alt_tel', 'az_tel', 'alt_tel', 'ra_tel', 'dec_tel'], + errors='ignore' + ) + evt = evt.drop( + columns=['reco_az', 'reco_alt', 'reco_ra', 'reco_dec'], + errors='ignore' + ) + + if n_events > 0: + # Events arrival time + evt = evt.assign( + dragon_time = arrival_time, + trigger_time = arrival_time, + ) + + # Telescope pointing + evt = evt.assign( + mc_az_tel = current_tel_pos.az.to('rad').value, + mc_alt_tel = current_tel_pos.alt.to('rad').value, + az_tel = current_tel_pos.az.to('rad').value, + alt_tel = current_tel_pos.alt.to('rad').value, + ra_tel = current_tel_pos.icrs.ra.to('rad').value, + dec_tel = current_tel_pos.icrs.dec.to('rad').value + ) + + # Reconstructed events coordinates + reco_coords = SkyCoord( + evt['reco_src_x'].to_numpy() * sample.units['distance'] * sample.cam2angle, + evt['reco_src_y'].to_numpy() * sample.units['distance'] * sample.cam2angle, + frame=offset_frame + ) + evt = evt.assign( + reco_az = reco_coords.altaz.az.to('rad').value, + reco_alt = reco_coords.altaz.alt.to('rad').value, + reco_ra = reco_coords.icrs.ra.to('rad').value, + reco_dec = reco_coords.icrs.dec.to('rad').value, + ) + else: + evt = evt.assign( + dragon_time = np.zeros(0), + trigger_time = np.zeros(0), + ) + evt = evt.assign( + mc_az_tel = np.zeros(0), + mc_alt_tel = np.zeros(0), + az_tel = np.zeros(0), + alt_tel = np.zeros(0), + ra_tel = np.zeros(0), + dec_tel = np.zeros(0) + ) + evt = evt.assign( + reco_az = np.zeros(0), + reco_alt = np.zeros(0), + reco_ra = np.zeros(0), + reco_dec = np.zeros(0) + ) + + evt = evt.drop('mc_src_name', errors='ignore') + evt = evt.assign( + mc_src_name = np.repeat(source.name, len(evt)) + ) + + events.append(evt) + + events = pd.concat(events) + + events = self.update_time_delta(events) + + return events From 8b5fd3ba2b6afe1f51a062ba8a1816b16ec19e81 Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Wed, 25 Oct 2023 12:31:55 +0900 Subject: [PATCH 03/28] FixedPointingDataRun: fixed the issue during event prediction, when the zero predicted events was leading to NaN in the sampling probability (as weights.sum() == 0) and the subsequent crash. --- src/srcsim/run/fixed.py | 50 +++++++++++++++++++++-------------------- 1 file changed, 26 insertions(+), 24 deletions(-) diff --git a/src/srcsim/run/fixed.py b/src/srcsim/run/fixed.py index 5706c7a..f085c3e 100644 --- a/src/srcsim/run/fixed.py +++ b/src/srcsim/run/fixed.py @@ -210,33 +210,34 @@ def predict(self, mccollections, source, tel_pos_tolerance=None, time_step=1*u.m n_mc_events = len(sample.evt_energy) n_events = np.random.poisson(weights.sum()) - p = weights / weights.sum() - idx = np.random.choice( - np.arange(n_mc_events), - size=n_events, - p=p - ) - evt = sample.data_table.iloc[idx] - offset_frame = offset_frame[idx] - arrival_time = arrival_time[idx] - current_tel_pos = current_tel_pos[idx] + if n_events > 0: + p = weights / weights.sum() + idx = np.random.choice( + np.arange(n_mc_events), + size=n_events, + p=p + ) - # Dropping the columns we're going to (re-)fill - evt = evt.drop( - columns=['dragon_time', 'trigger_time'], - errors='ignore' - ) - evt = evt.drop( - columns=['mc_az_tel', 'mc_alt_tel', 'az_tel', 'alt_tel', 'ra_tel', 'dec_tel'], - errors='ignore' - ) - evt = evt.drop( - columns=['reco_az', 'reco_alt', 'reco_ra', 'reco_dec'], - errors='ignore' - ) + evt = sample.data_table.iloc[idx] + offset_frame = offset_frame[idx] + arrival_time = arrival_time[idx] + current_tel_pos = current_tel_pos[idx] + + # Dropping the columns we're going to (re-)fill + evt = evt.drop( + columns=['dragon_time', 'trigger_time'], + errors='ignore' + ) + evt = evt.drop( + columns=['mc_az_tel', 'mc_alt_tel', 'az_tel', 'alt_tel', 'ra_tel', 'dec_tel'], + errors='ignore' + ) + evt = evt.drop( + columns=['reco_az', 'reco_alt', 'reco_ra', 'reco_dec'], + errors='ignore' + ) - if n_events > 0: # Events arrival time evt = evt.assign( dragon_time = arrival_time, @@ -266,6 +267,7 @@ def predict(self, mccollections, source, tel_pos_tolerance=None, time_step=1*u.m reco_dec = reco_coords.icrs.dec.to('rad').value, ) else: + evt = sample.data_table.iloc[0:0] evt = evt.assign( dragon_time = np.zeros(0), trigger_time = np.zeros(0), From eff74b578c92eebff15601ca284d0d5457c5e00b Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Wed, 25 Oct 2023 12:42:01 +0900 Subject: [PATCH 04/28] Moving the base DataRun class to its own file. --- src/srcsim/run/base.py | 48 +++++++++++++++++++++++++++++++++++++++++ src/srcsim/run/fixed.py | 47 +--------------------------------------- 2 files changed, 49 insertions(+), 46 deletions(-) create mode 100644 src/srcsim/run/base.py diff --git a/src/srcsim/run/base.py b/src/srcsim/run/base.py new file mode 100644 index 0000000..b686172 --- /dev/null +++ b/src/srcsim/run/base.py @@ -0,0 +1,48 @@ +import numpy as np + +class DataRun: + def __init__(self, tel_pos, tstart, tstop, obsloc, id=0): + self.id = id + self.tel_pos = tel_pos + self.obsloc = obsloc + self.tstart = tstart + self.tstop = tstop + + @classmethod + def from_config(cls, config): + pass + + def to_dict(self): + pass + + def predict(self, mccollections, source, tel_pos_tolerance, time_step): + pass + + @classmethod + def time_sort(cls, events): + timestamp = np.array(events["trigger_time"]) + sorting = np.argsort(timestamp) + + return events.iloc[sorting] + + @classmethod + def update_time_delta(cls, events): + event_time = events['trigger_time'].to_numpy() + delta_time = np.zeros_like(event_time) + + if len(event_time) > 0: + argsorted = event_time.argsort() + dt = np.diff(event_time[argsorted]) + delta_time[argsorted[:-1]] = dt + delta_time[argsorted[-1]] = dt[-1] + + events = events.drop( + columns=['delta_t'], + errors='ignore' + ) + events = events.assign( + delta_t = delta_time, + ) + + return events + \ No newline at end of file diff --git a/src/srcsim/run/fixed.py b/src/srcsim/run/fixed.py index f085c3e..f24fe12 100644 --- a/src/srcsim/run/fixed.py +++ b/src/srcsim/run/fixed.py @@ -5,52 +5,7 @@ from astropy.time import Time from astropy.coordinates import SkyCoord, SkyOffsetFrame, EarthLocation, AltAz - -class DataRun: - def __init__(self, tel_pos, tstart, tstop, obsloc, id=0): - self.id = id - self.tel_pos = tel_pos - self.obsloc = obsloc - self.tstart = tstart - self.tstop = tstop - - @classmethod - def from_config(cls, config): - pass - - def to_dict(self): - pass - - def predict(self, mccollections, source, tel_pos_tolerance, time_step): - pass - - @classmethod - def time_sort(cls, events): - timestamp = np.array(events["trigger_time"]) - sorting = np.argsort(timestamp) - - return events.iloc[sorting] - - @classmethod - def update_time_delta(cls, events): - event_time = events['trigger_time'].to_numpy() - delta_time = np.zeros_like(event_time) - - if len(event_time) > 0: - argsorted = event_time.argsort() - dt = np.diff(event_time[argsorted]) - delta_time[argsorted[:-1]] = dt - delta_time[argsorted[-1]] = dt[-1] - - events = events.drop( - columns=['delta_t'], - errors='ignore' - ) - events = events.assign( - delta_t = delta_time, - ) - - return events +from .base import DataRun class FixedPointingDataRun(DataRun): From 2093497c6378b5f1c64f39e646e11f14210344a4 Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Wed, 25 Oct 2023 12:45:19 +0900 Subject: [PATCH 05/28] Moving the predict() method implementation to the base DataRun class - it's children are now expected to implement the tel_pos_to_altaz() method to make predict() work. --- src/srcsim/run/base.py | 164 +++++++++++++++++++++++++++++++++++++++- src/srcsim/run/fixed.py | 159 +------------------------------------- 2 files changed, 163 insertions(+), 160 deletions(-) diff --git a/src/srcsim/run/base.py b/src/srcsim/run/base.py index b686172..2326dce 100644 --- a/src/srcsim/run/base.py +++ b/src/srcsim/run/base.py @@ -1,4 +1,9 @@ import numpy as np +import pandas as pd +import astropy.units as u +from astropy.time import Time +from astropy.coordinates import SkyCoord, SkyOffsetFrame, AltAz + class DataRun: def __init__(self, tel_pos, tstart, tstop, obsloc, id=0): @@ -15,7 +20,7 @@ def from_config(cls, config): def to_dict(self): pass - def predict(self, mccollections, source, tel_pos_tolerance, time_step): + def tel_pos_to_altaz(self, frame): pass @classmethod @@ -45,4 +50,159 @@ def update_time_delta(cls, events): ) return events - \ No newline at end of file + + def predict(self, mccollections, source, tel_pos_tolerance=None, time_step=1*u.minute): + events = [] + + unix_edges = np.arange(self.tstart.unix, self.tstop.unix, step=time_step.to('s').value) + if unix_edges[-1] < self.tstop.unix: + unix_edges = np.concatenate((unix_edges, [self.tstop.unix])) + + tedges = Time(unix_edges, format='unix') + + if len(tedges) > 1: + tdelta = np.diff(tedges) + else: + tdelta = [self.tstop - self.tstart] + + for tstart, dt in zip(tedges[:-1], tdelta): + frame = AltAz( + obstime=tstart + dt/2, + location=self.obsloc + ) + tel_pos = self.tel_pos_to_altaz(frame) + + if tel_pos_tolerance is None: + mc = mccollections[source.emission_type].get_closest(tel_pos.altaz) + elif isinstance(tel_pos_tolerance, u.Quantity): + mc = mccollections[source.emission_type].get_nearby(tel_pos, tel_pos_tolerance) + elif isinstance(tel_pos_tolerance, list) or isinstance(tel_pos_tolerance, tuple): + mc = mccollections[source.emission_type].get_in_box( + tel_pos, + max_lon_offset=tel_pos_tolerance[0], + max_lat_offset=tel_pos_tolerance[1], + ) + else: + raise ValueError(f"Data type '{type(tel_pos_tolerance)}' for argument 'tel_pos_tolerance' is not supported") + + nsamples = len(mc.samples) + + for sample in mc.samples: + # Randomly distributing events within the time bin + arrival_time = np.random.uniform(tstart.unix, (tstart+dt).unix, size=len(sample.data_table)) + + # Recalculating telescope Alt/Az for the mock event arrival times + current_frame = AltAz( + obstime=Time(arrival_time, format='unix'), + location=self.obsloc + ) + current_tel_pos = self.tel_pos_to_altaz(current_frame) + + # Astropy does not pass the location / time + # to the offset frame, need to do this manually + offset_frame = SkyOffsetFrame( + origin=current_tel_pos.altaz.skyoffset_frame().origin, + location=current_tel_pos.altaz.frame.location, + obstime=current_tel_pos.altaz.frame.obstime + ) + coords = SkyCoord( + sample.evt_coord.skyoffsetaltaz.lon, + sample.evt_coord.skyoffsetaltaz.lat, + frame=offset_frame + ) + expected_flux = source.dndedo(sample.evt_energy, coords.icrs) + model_flux = sample.dndedo(sample.evt_energy, sample.evt_coord) + + weights = (1 / nsamples * dt * expected_flux / model_flux).decompose() + + n_mc_events = len(sample.evt_energy) + n_events = np.random.poisson(weights.sum()) + + if n_events > 0: + p = weights / weights.sum() + idx = np.random.choice( + np.arange(n_mc_events), + size=n_events, + p=p + ) + + evt = sample.data_table.iloc[idx] + offset_frame = offset_frame[idx] + arrival_time = arrival_time[idx] + current_tel_pos = current_tel_pos[idx] + + # Dropping the columns we're going to (re-)fill + evt = evt.drop( + columns=['dragon_time', 'trigger_time'], + errors='ignore' + ) + evt = evt.drop( + columns=['mc_az_tel', 'mc_alt_tel', 'az_tel', 'alt_tel', 'ra_tel', 'dec_tel'], + errors='ignore' + ) + evt = evt.drop( + columns=['reco_az', 'reco_alt', 'reco_ra', 'reco_dec'], + errors='ignore' + ) + + # Events arrival time + evt = evt.assign( + dragon_time = arrival_time, + trigger_time = arrival_time, + ) + + # Telescope pointing + evt = evt.assign( + mc_az_tel = current_tel_pos.az.to('rad').value, + mc_alt_tel = current_tel_pos.alt.to('rad').value, + az_tel = current_tel_pos.az.to('rad').value, + alt_tel = current_tel_pos.alt.to('rad').value, + ra_tel = current_tel_pos.icrs.ra.to('rad').value, + dec_tel = current_tel_pos.icrs.dec.to('rad').value + ) + + # Reconstructed events coordinates + reco_coords = SkyCoord( + evt['reco_src_x'].to_numpy() * sample.units['distance'] * sample.cam2angle, + evt['reco_src_y'].to_numpy() * sample.units['distance'] * sample.cam2angle, + frame=offset_frame + ) + evt = evt.assign( + reco_az = reco_coords.altaz.az.to('rad').value, + reco_alt = reco_coords.altaz.alt.to('rad').value, + reco_ra = reco_coords.icrs.ra.to('rad').value, + reco_dec = reco_coords.icrs.dec.to('rad').value, + ) + else: + evt = sample.data_table.iloc[0:0] + evt = evt.assign( + dragon_time = np.zeros(0), + trigger_time = np.zeros(0), + ) + evt = evt.assign( + mc_az_tel = np.zeros(0), + mc_alt_tel = np.zeros(0), + az_tel = np.zeros(0), + alt_tel = np.zeros(0), + ra_tel = np.zeros(0), + dec_tel = np.zeros(0) + ) + evt = evt.assign( + reco_az = np.zeros(0), + reco_alt = np.zeros(0), + reco_ra = np.zeros(0), + reco_dec = np.zeros(0) + ) + + evt = evt.drop('mc_src_name', errors='ignore') + evt = evt.assign( + mc_src_name = np.repeat(source.name, len(evt)) + ) + + events.append(evt) + + events = pd.concat(events) + + events = self.update_time_delta(events) + + return events diff --git a/src/srcsim/run/fixed.py b/src/srcsim/run/fixed.py index f24fe12..eee1c68 100644 --- a/src/srcsim/run/fixed.py +++ b/src/srcsim/run/fixed.py @@ -1,9 +1,8 @@ import yaml import numpy as np -import pandas as pd import astropy.units as u from astropy.time import Time -from astropy.coordinates import SkyCoord, SkyOffsetFrame, EarthLocation, AltAz +from astropy.coordinates import SkyCoord, EarthLocation, AltAz from .base import DataRun @@ -98,159 +97,3 @@ def tel_pos_to_altaz(self, frame): frame=frame ) return tel_pos - - def predict(self, mccollections, source, tel_pos_tolerance=None, time_step=1*u.minute): - events = [] - - unix_edges = np.arange(self.tstart.unix, self.tstop.unix, step=time_step.to('s').value) - if unix_edges[-1] < self.tstop.unix: - unix_edges = np.concatenate((unix_edges, [self.tstop.unix])) - - tedges = Time(unix_edges, format='unix') - - if len(tedges) > 1: - tdelta = np.diff(tedges) - else: - tdelta = [self.tstop - self.tstart] - - for tstart, dt in zip(tedges[:-1], tdelta): - frame = AltAz( - obstime=tstart + dt/2, - location=self.obsloc - ) - tel_pos = self.tel_pos_to_altaz(frame) - - if tel_pos_tolerance is None: - mc = mccollections[source.emission_type].get_closest(tel_pos.altaz) - elif isinstance(tel_pos_tolerance, u.Quantity): - mc = mccollections[source.emission_type].get_nearby(tel_pos, tel_pos_tolerance) - elif isinstance(tel_pos_tolerance, list) or isinstance(tel_pos_tolerance, tuple): - mc = mccollections[source.emission_type].get_in_box( - tel_pos, - max_lon_offset=tel_pos_tolerance[0], - max_lat_offset=tel_pos_tolerance[1], - ) - else: - raise ValueError(f"Data type '{type(tel_pos_tolerance)}' for argument 'tel_pos_tolerance' is not supported") - - nsamples = len(mc.samples) - - for sample in mc.samples: - # Randomly distributing events within the time bin - arrival_time = np.random.uniform(tstart.unix, (tstart+dt).unix, size=len(sample.data_table)) - - # Recalculating telescope Alt/Az for the mock event arrival times - current_frame = AltAz( - obstime=Time(arrival_time, format='unix'), - location=self.obsloc - ) - current_tel_pos = self.tel_pos_to_altaz(current_frame) - - # Astropy does not pass the location / time - # to the offset frame, need to do this manually - offset_frame = SkyOffsetFrame( - origin=current_tel_pos.altaz.skyoffset_frame().origin, - location=current_tel_pos.altaz.frame.location, - obstime=current_tel_pos.altaz.frame.obstime - ) - coords = SkyCoord( - sample.evt_coord.skyoffsetaltaz.lon, - sample.evt_coord.skyoffsetaltaz.lat, - frame=offset_frame - ) - expected_flux = source.dndedo(sample.evt_energy, coords.icrs) - model_flux = sample.dndedo(sample.evt_energy, sample.evt_coord) - - weights = (1 / nsamples * dt * expected_flux / model_flux).decompose() - - n_mc_events = len(sample.evt_energy) - n_events = np.random.poisson(weights.sum()) - - if n_events > 0: - p = weights / weights.sum() - idx = np.random.choice( - np.arange(n_mc_events), - size=n_events, - p=p - ) - - evt = sample.data_table.iloc[idx] - offset_frame = offset_frame[idx] - arrival_time = arrival_time[idx] - current_tel_pos = current_tel_pos[idx] - - # Dropping the columns we're going to (re-)fill - evt = evt.drop( - columns=['dragon_time', 'trigger_time'], - errors='ignore' - ) - evt = evt.drop( - columns=['mc_az_tel', 'mc_alt_tel', 'az_tel', 'alt_tel', 'ra_tel', 'dec_tel'], - errors='ignore' - ) - evt = evt.drop( - columns=['reco_az', 'reco_alt', 'reco_ra', 'reco_dec'], - errors='ignore' - ) - - # Events arrival time - evt = evt.assign( - dragon_time = arrival_time, - trigger_time = arrival_time, - ) - - # Telescope pointing - evt = evt.assign( - mc_az_tel = current_tel_pos.az.to('rad').value, - mc_alt_tel = current_tel_pos.alt.to('rad').value, - az_tel = current_tel_pos.az.to('rad').value, - alt_tel = current_tel_pos.alt.to('rad').value, - ra_tel = current_tel_pos.icrs.ra.to('rad').value, - dec_tel = current_tel_pos.icrs.dec.to('rad').value - ) - - # Reconstructed events coordinates - reco_coords = SkyCoord( - evt['reco_src_x'].to_numpy() * sample.units['distance'] * sample.cam2angle, - evt['reco_src_y'].to_numpy() * sample.units['distance'] * sample.cam2angle, - frame=offset_frame - ) - evt = evt.assign( - reco_az = reco_coords.altaz.az.to('rad').value, - reco_alt = reco_coords.altaz.alt.to('rad').value, - reco_ra = reco_coords.icrs.ra.to('rad').value, - reco_dec = reco_coords.icrs.dec.to('rad').value, - ) - else: - evt = sample.data_table.iloc[0:0] - evt = evt.assign( - dragon_time = np.zeros(0), - trigger_time = np.zeros(0), - ) - evt = evt.assign( - mc_az_tel = np.zeros(0), - mc_alt_tel = np.zeros(0), - az_tel = np.zeros(0), - alt_tel = np.zeros(0), - ra_tel = np.zeros(0), - dec_tel = np.zeros(0) - ) - evt = evt.assign( - reco_az = np.zeros(0), - reco_alt = np.zeros(0), - reco_ra = np.zeros(0), - reco_dec = np.zeros(0) - ) - - evt = evt.drop('mc_src_name', errors='ignore') - evt = evt.assign( - mc_src_name = np.repeat(source.name, len(evt)) - ) - - events.append(evt) - - events = pd.concat(events) - - events = self.update_time_delta(events) - - return events From 61f1a7727b830ef9ec9466d499f0c81529630ae4 Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Wed, 25 Oct 2023 12:49:06 +0900 Subject: [PATCH 06/28] Renaming DataRun -> SkyDataRun --- src/srcsim/run/__init__.py | 2 +- src/srcsim/run/sky.py | 2 +- src/srcsim/rungen.py | 8 ++++---- src/srcsim/scripts/sim_run.py | 4 ++-- 4 files changed, 8 insertions(+), 8 deletions(-) diff --git a/src/srcsim/run/__init__.py b/src/srcsim/run/__init__.py index 0cafdae..fe4c50d 100644 --- a/src/srcsim/run/__init__.py +++ b/src/srcsim/run/__init__.py @@ -1 +1 @@ -from .sky import DataRun \ No newline at end of file +from .sky import SkyDataRun \ No newline at end of file diff --git a/src/srcsim/run/sky.py b/src/srcsim/run/sky.py index 46b0abe..8ef38cb 100644 --- a/src/srcsim/run/sky.py +++ b/src/srcsim/run/sky.py @@ -6,7 +6,7 @@ from astropy.coordinates import SkyCoord, SkyOffsetFrame, EarthLocation, AltAz -class DataRun: +class SkyDataRun: def __init__(self, tel_pos, tstart, tstop, obsloc, id=0): self.id = id self.tel_pos = tel_pos diff --git a/src/srcsim/rungen.py b/src/srcsim/rungen.py index 63f32f6..b9a9aa2 100644 --- a/src/srcsim/rungen.py +++ b/src/srcsim/rungen.py @@ -7,7 +7,7 @@ from astropy.time import Time from astropy.coordinates import SkyCoord, AltAz, EarthLocation -from .run import DataRun +from .run import SkyDataRun def generator(config): @@ -205,7 +205,7 @@ def get_runs(cls, obsloc, tel_pos, tobs, altmin, altmax, azmin, azmax, tstart=No pos_angles = wobble_start_angle + np.linspace(0, 2* np.pi, num=wobble_count+1)[:-1] * u.rad runs = tuple( - DataRun( + SkyDataRun( tel_pos.directional_offset_by(pos_angles[run_id % wobble_count], wobble_offset), tstart, tstop, @@ -217,7 +217,7 @@ def get_runs(cls, obsloc, tel_pos, tobs, altmin, altmax, azmin, azmax, tstart=No else: runs = tuple( - DataRun(tel_pos, tstart, tstop, obsloc, run_id) + SkyDataRun(tel_pos, tstart, tstop, obsloc, run_id) for run_id, (tstart, tstop) in enumerate(zip(tstarts, tstops)) ) @@ -251,7 +251,7 @@ def get_runs(cls, file_mask, obsloc): ) runs = [ - DataRun( + SkyDataRun( run_info['tel_pos'], run_info['tstart'], run_info['tstop'], diff --git a/src/srcsim/scripts/sim_run.py b/src/srcsim/scripts/sim_run.py index a528507..590f1f0 100644 --- a/src/srcsim/scripts/sim_run.py +++ b/src/srcsim/scripts/sim_run.py @@ -7,7 +7,7 @@ from srcsim.mc import MCCollection from srcsim.src import generator as srcgen -from srcsim.run import DataRun +from srcsim.run import SkyDataRun def info_message(text): @@ -80,7 +80,7 @@ def main(): print(srcs) info_message('Preparing the data run') - run = DataRun.from_config(cfg['run']) + run = SkyDataRun.from_config(cfg['run']) print(run) info_message('Starting event sampling') From 7435429ce4567643829fa543e071138836ccb428 Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Wed, 25 Oct 2023 16:16:07 +0900 Subject: [PATCH 07/28] SkyDataRun now inherits DataRun for unification with FixedPointingDataRun. --- src/srcsim/run/sky.py | 194 ++---------------------------------------- 1 file changed, 5 insertions(+), 189 deletions(-) diff --git a/src/srcsim/run/sky.py b/src/srcsim/run/sky.py index 8ef38cb..fb241d2 100644 --- a/src/srcsim/run/sky.py +++ b/src/srcsim/run/sky.py @@ -5,15 +5,10 @@ from astropy.time import Time from astropy.coordinates import SkyCoord, SkyOffsetFrame, EarthLocation, AltAz +from .base import DataRun -class SkyDataRun: - def __init__(self, tel_pos, tstart, tstop, obsloc, id=0): - self.id = id - self.tel_pos = tel_pos - self.obsloc = obsloc - self.tstart = tstart - self.tstop = tstop - + +class SkyDataRun(DataRun): def __repr__(self): frame_start = AltAz(obstime=self.tstart, location=self.obsloc) frame_stop = AltAz(obstime=self.tstop, location=self.obsloc) @@ -55,34 +50,6 @@ def from_config(cls, config): return data_run - @classmethod - def time_sort(cls, events): - timestamp = np.array(events["trigger_time"]) - sorting = np.argsort(timestamp) - - return events.iloc[sorting] - - @classmethod - def update_time_delta(cls, events): - event_time = events['trigger_time'].to_numpy() - delta_time = np.zeros_like(event_time) - - if len(event_time) > 0: - argsorted = event_time.argsort() - dt = np.diff(event_time[argsorted]) - delta_time[argsorted[:-1]] = dt - delta_time[argsorted[-1]] = dt[-1] - - events = events.drop( - columns=['delta_t'], - errors='ignore' - ) - events = events.assign( - delta_t = delta_time, - ) - - return events - def to_dict(self): data = {'id': self.id, 'pointing': {}, 'time': {}, 'location': {}} @@ -98,156 +65,5 @@ def to_dict(self): return data - def predict(self, mccollections, source, tel_pos_tolerance=None, time_step=1*u.minute): - events = [] - - unix_edges = np.arange(self.tstart.unix, self.tstop.unix, step=time_step.to('s').value) - if unix_edges[-1] < self.tstop.unix: - unix_edges = np.concatenate((unix_edges, [self.tstop.unix])) - - tedges = Time(unix_edges, format='unix') - - if len(tedges) > 1: - tdelta = np.diff(tedges) - else: - tdelta = [self.tstop - self.tstart] - - for tstart, dt in zip(tedges[:-1], tdelta): - frame = AltAz( - obstime=tstart + dt/2, - location=self.obsloc - ) - tel_pos = self.tel_pos.transform_to(frame) - - if tel_pos_tolerance is None: - mc = mccollections[source.emission_type].get_closest(tel_pos.altaz) - elif isinstance(tel_pos_tolerance, u.Quantity): - mc = mccollections[source.emission_type].get_nearby(tel_pos, tel_pos_tolerance) - elif isinstance(tel_pos_tolerance, list) or isinstance(tel_pos_tolerance, tuple): - mc = mccollections[source.emission_type].get_in_box( - tel_pos, - max_lon_offset=tel_pos_tolerance[0], - max_lat_offset=tel_pos_tolerance[1], - ) - else: - raise ValueError(f"Data type '{type(tel_pos_tolerance)}' for argument 'tel_pos_tolerance' is not supported") - - nsamples = len(mc.samples) - - for sample in mc.samples: - # Randomly distributing events within the time bin - arrival_time = np.random.uniform(tstart.unix, (tstart+dt).unix, size=len(sample.data_table)) - - # Recalculating telescope Alt/Az for the mock event arrival times - current_frame = AltAz( - obstime=Time(arrival_time, format='unix'), - location=self.obsloc - ) - current_tel_pos = self.tel_pos.transform_to(current_frame) - - # Astropy does not pass the location / time - # to the offset frame, need to do this manually - offset_frame = SkyOffsetFrame( - origin=current_tel_pos.altaz.skyoffset_frame().origin, - location=current_tel_pos.altaz.frame.location, - obstime=current_tel_pos.altaz.frame.obstime - ) - coords = SkyCoord( - sample.evt_coord.skyoffsetaltaz.lon, - sample.evt_coord.skyoffsetaltaz.lat, - frame=offset_frame - ) - expected_flux = source.dndedo(sample.evt_energy, coords.icrs) - model_flux = sample.dndedo(sample.evt_energy, sample.evt_coord) - - weights = (1 / nsamples * dt * expected_flux / model_flux).decompose() - - n_mc_events = len(sample.evt_energy) - n_events = np.random.poisson(weights.sum()) - p = weights / weights.sum() - idx = np.random.choice( - np.arange(n_mc_events), - size=n_events, - p=p - ) - - evt = sample.data_table.iloc[idx] - offset_frame = offset_frame[idx] - arrival_time = arrival_time[idx] - current_tel_pos = current_tel_pos[idx] - - # Dropping the columns we're going to (re-)fill - evt = evt.drop( - columns=['dragon_time', 'trigger_time'], - errors='ignore' - ) - evt = evt.drop( - columns=['mc_az_tel', 'mc_alt_tel', 'az_tel', 'alt_tel', 'ra_tel', 'dec_tel'], - errors='ignore' - ) - evt = evt.drop( - columns=['reco_az', 'reco_alt', 'reco_ra', 'reco_dec'], - errors='ignore' - ) - - if n_events > 0: - # Events arrival time - evt = evt.assign( - dragon_time = arrival_time, - trigger_time = arrival_time, - ) - - # Telescope pointing - evt = evt.assign( - mc_az_tel = current_tel_pos.az.to('rad').value, - mc_alt_tel = current_tel_pos.alt.to('rad').value, - az_tel = current_tel_pos.az.to('rad').value, - alt_tel = current_tel_pos.alt.to('rad').value, - ra_tel = self.tel_pos.icrs.ra.to('rad').value, - dec_tel = self.tel_pos.icrs.dec.to('rad').value - ) - - # Reconstructed events coordinates - reco_coords = SkyCoord( - evt['reco_src_x'].to_numpy() * sample.units['distance'] * sample.cam2angle, - evt['reco_src_y'].to_numpy() * sample.units['distance'] * sample.cam2angle, - frame=offset_frame - ) - evt = evt.assign( - reco_az = reco_coords.altaz.az.to('rad').value, - reco_alt = reco_coords.altaz.alt.to('rad').value, - reco_ra = reco_coords.icrs.ra.to('rad').value, - reco_dec = reco_coords.icrs.dec.to('rad').value, - ) - else: - evt = evt.assign( - dragon_time = np.zeros(0), - trigger_time = np.zeros(0), - ) - evt = evt.assign( - mc_az_tel = np.zeros(0), - mc_alt_tel = np.zeros(0), - az_tel = np.zeros(0), - alt_tel = np.zeros(0), - ra_tel = np.zeros(0), - dec_tel = np.zeros(0) - ) - evt = evt.assign( - reco_az = np.zeros(0), - reco_alt = np.zeros(0), - reco_ra = np.zeros(0), - reco_dec = np.zeros(0) - ) - - evt = evt.drop('mc_src_name', errors='ignore') - evt = evt.assign( - mc_src_name = np.repeat(source.name, len(evt)) - ) - - events.append(evt) - - events = pd.concat(events) - - events = self.update_time_delta(events) - - return events + def tel_pos_to_altaz(self, frame): + return self.tel_pos.transform_to(frame) From bb4fbf28bff609b6f2744afa76ebab6fa8828215 Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Wed, 25 Oct 2023 16:34:40 +0900 Subject: [PATCH 08/28] Dropping no longer necessary imports. --- src/srcsim/run/sky.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/srcsim/run/sky.py b/src/srcsim/run/sky.py index fb241d2..5d91a46 100644 --- a/src/srcsim/run/sky.py +++ b/src/srcsim/run/sky.py @@ -1,9 +1,7 @@ import yaml -import numpy as np -import pandas as pd import astropy.units as u from astropy.time import Time -from astropy.coordinates import SkyCoord, SkyOffsetFrame, EarthLocation, AltAz +from astropy.coordinates import SkyCoord, EarthLocation, AltAz from .base import DataRun From 945a8ed80e0272bae9bcaddae47d980bf74e3714 Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Wed, 25 Oct 2023 16:35:40 +0900 Subject: [PATCH 09/28] Incrementing package version due to newly added functionality. --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index b5440dc..99dc1c4 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -4,7 +4,7 @@ build-backend = "setuptools.build_meta" [project] name = "srcsim" -version = "1.0.0" +version = "1.1.0" authors = [ { name="Ievgen Vovk", email="vovk@icrr.u-tokyo.ac.jp" }, { name="Marcel Strzys", email="strzys@icrr.u-tokyo.ac.jp" }, From cc07f78a421830bf6c00bc6ac1fb2733840ab4bd Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Wed, 25 Oct 2023 16:40:14 +0900 Subject: [PATCH 10/28] Adding FixedPointingDataRun to the "run" submodule contents. --- src/srcsim/run/__init__.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/srcsim/run/__init__.py b/src/srcsim/run/__init__.py index fe4c50d..dda540e 100644 --- a/src/srcsim/run/__init__.py +++ b/src/srcsim/run/__init__.py @@ -1 +1,2 @@ -from .sky import SkyDataRun \ No newline at end of file +from .sky import SkyDataRun +from .fixed import FixedPointingDataRun \ No newline at end of file From 730f07996262d0bc48b9893ab4d271f66da8978d Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Sat, 28 Oct 2023 17:28:21 +0900 Subject: [PATCH 11/28] FixedObsGenerator class prototype --- src/srcsim/rungen/fixed.py | 167 +++++++++++++++++++++++++++++++++++++ 1 file changed, 167 insertions(+) create mode 100644 src/srcsim/rungen/fixed.py diff --git a/src/srcsim/rungen/fixed.py b/src/srcsim/rungen/fixed.py new file mode 100644 index 0000000..8526de7 --- /dev/null +++ b/src/srcsim/rungen/fixed.py @@ -0,0 +1,167 @@ +import yaml +import numpy as np +import pandas as pd +import astropy.units as u +from scipy.interpolate import CubicSpline +from astropy.time import Time +from astropy.coordinates import SkyCoord, SkyOffsetFrame, EarthLocation, AltAz + +from srcsim.rungen import get_trajectory, get_time_intervals +from srcsim.run import FixedPointingDataRun + + +class FixedObsGenerator: + @classmethod + def get_runs_from_config(cls, config): + if isinstance(config, str): + cfg = yaml.load(open(config, "r"), Loader=yaml.FullLoader) + else: + cfg = config + + tstart = Time(cfg['time']['start']) + ra_width = u.Quantity(cfg['pointing']['width']) + alt = u.Quantity(cfg['pointing']['center']['alt']) + tel_pos = SkyCoord( + u.Quantity(cfg['pointing']['center']['ra']), + u.Quantity(cfg['pointing']['center']['dec']), + frame='icrs' + ) + + duration = u.Quantity(cfg['time']['duration']) + + obsloc = EarthLocation( + lat=u.Quantity(cfg['location']['lat']), + lon=u.Quantity(cfg['location']['lon']), + height=u.Quantity(cfg['location']['height']), + ) + + return cls.get_runs(obsloc, tel_pos, duration, alt, ra_width, tstart) + + @classmethod + def get_runs(cls, obsloc, tel_pos, tobs, alt, ra_width, tstart=None): + time_step = 1 * u.m + + tel_pos_trail = SkyCoord( + tel_pos.icrs.ra + ra_width / 2, + tel_pos.icrs.dec, + frame='icrs' + ) + tel_pos_lead = SkyCoord( + tel_pos.icrs.ra - ra_width / 2, + tel_pos.icrs.dec, + frame='icrs' + ) + + # First pass - first day + track_trail = get_trajectory( + tel_pos_trail.icrs, + tstart, + tstop=tstart + 1 * u.d, + time_step=time_step, + obsloc=obsloc + ) + track_lead = get_trajectory( + tel_pos_lead.icrs, + tstart, + tstop=tstart + 1 * u.d, + time_step=time_step, + obsloc=obsloc + ) + cs_lead = CubicSpline( + track_lead.obstime.mjd, + track_lead.alt - alt + ) + cs_trail = CubicSpline( + track_trail.obstime.mjd, + track_trail.alt - alt + ) + + tstarts = Time(cs_lead.roots(extrapolate=False), format='mjd') + tstops = Time(cs_trail.roots(extrapolate=False), format='mjd') + + # Total number of the required observation days + run_durations = tstops - tstarts + nsequences = (tobs / np.sum(run_durations)).decompose() + ndays_full = np.floor(nsequences) + ndays_total = np.ceil(nsequences) + + if ndays_full > 0: + # Next pass - to cover the simulation time interval with fully completed runs + track_lead = get_trajectory( + tel_pos_lead.icrs, + tstart, + tstop=tstart + ndays_full * u.d, + time_step=time_step, + obsloc=obsloc + ) + track_trail = get_trajectory( + tel_pos_trail.icrs, + tstart, + tstop=tstart + ndays_full * u.d, + time_step=time_step, + obsloc=obsloc + ) + + cs_lead = CubicSpline( + track_lead.obstime.mjd, + track_lead.alt - alt + ) + cs_trail = CubicSpline( + track_trail.obstime.mjd, + track_trail.alt - alt + ) + tstarts = Time(cs_lead.roots(extrapolate=False), format='mjd') + tstops = Time(cs_trail.roots(extrapolate=False), format='mjd') + remaining_tobs = tobs - np.sum(tstops - tstarts) + + # Next pass - additional incomplete runs + track_lead = get_trajectory( + tel_pos_lead.icrs, + tstart=tstart + ndays_full * u.d, + tstop=tstart + ndays_total * u.d, + time_step=time_step, + obsloc=obsloc + ) + cs_lead = CubicSpline( + track_lead.obstime.mjd, + track_lead.alt - alt + ) + _tstarts = Time(cs_lead.roots(extrapolate=False), format='mjd') + _tstops = Time(_tstarts + remaining_tobs / (len(_tstarts))) + + tstarts = Time([tstarts, _tstarts]) + tstops = Time([tstops, _tstops]) + + else: + track_lead = get_trajectory( + tel_pos_lead.icrs, + tstart, + tstop=tstart + ndays_full * u.d, + time_step=time_step, + obsloc=obsloc + ) + track_trail = get_trajectory( + tel_pos_trail.icrs, + tstart, + tstop=tstart + ndays_full * u.d, + time_step=time_step, + obsloc=obsloc + ) + + cs_lead = CubicSpline( + track_lead.obstime.mjd, + track_lead.alt - alt + ) + cs_trail = CubicSpline( + track_trail.obstime.mjd, + track_trail.alt - alt + ) + tstarts = Time(cs_lead.roots(extrapolate=False), format='mjd') + tstops = Time(tstarts + tobs / (len(tstarts))) + + runs = tuple( + FixedPointingDataRun(tel_pos, tstart, tstop, obsloc, run_id) + for run_id, (tstart, tstop) in enumerate(zip(tstarts, tstops)) + ) + + return runs From 6ab978d41ccebc4a1a3d3a85c68e194785ae474e Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Sat, 28 Oct 2023 17:32:08 +0900 Subject: [PATCH 12/28] "rungen" submodule refactoring. --- sim_in_box.py | 2 +- src/srcsim/{rungen.py => rungen/sky.py} | 2 +- src/srcsim/scripts/get_runs.py | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) rename src/srcsim/{rungen.py => rungen/sky.py} (99%) diff --git a/sim_in_box.py b/sim_in_box.py index 68b3619..007f892 100644 --- a/sim_in_box.py +++ b/sim_in_box.py @@ -7,7 +7,7 @@ from srcsim.mc import MCCollection from srcsim.src import generator as srcgen -from srcsim.rungen import AltAzBoxGenerator +from srcsim.rungen.sky import AltAzBoxGenerator def info_message(text): diff --git a/src/srcsim/rungen.py b/src/srcsim/rungen/sky.py similarity index 99% rename from src/srcsim/rungen.py rename to src/srcsim/rungen/sky.py index b9a9aa2..126dd27 100644 --- a/src/srcsim/rungen.py +++ b/src/srcsim/rungen/sky.py @@ -7,7 +7,7 @@ from astropy.time import Time from astropy.coordinates import SkyCoord, AltAz, EarthLocation -from .run import SkyDataRun +from ..run import SkyDataRun def generator(config): diff --git a/src/srcsim/scripts/get_runs.py b/src/srcsim/scripts/get_runs.py index 47ea066..3fd14b8 100644 --- a/src/srcsim/scripts/get_runs.py +++ b/src/srcsim/scripts/get_runs.py @@ -2,7 +2,7 @@ import datetime import argparse -from srcsim.rungen import generator as rungen +from srcsim.rungen.sky import generator as rungen def info_message(text): From b0c83f821a12e8577f62ed78cf5a822a22e3a5ca Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Sat, 28 Oct 2023 17:38:27 +0900 Subject: [PATCH 13/28] Run generator: moved to a separate file and added FixedObsGenerator support. --- src/srcsim/rungen/generator.py | 24 ++++++++++++++++++++++++ src/srcsim/rungen/sky.py | 18 ------------------ src/srcsim/scripts/get_runs.py | 2 +- 3 files changed, 25 insertions(+), 19 deletions(-) create mode 100644 src/srcsim/rungen/generator.py diff --git a/src/srcsim/rungen/generator.py b/src/srcsim/rungen/generator.py new file mode 100644 index 0000000..62e299e --- /dev/null +++ b/src/srcsim/rungen/generator.py @@ -0,0 +1,24 @@ +import yaml + +from .sky import AltAzBoxGenerator, DataMatchingGenerator +from .fixed import FixedObsGenerator + + +def generator(config): + if isinstance(config, str): + cfg = yaml.load(open(config, "r"), Loader=yaml.FullLoader) + else: + cfg = config + + runs = [] + + if cfg['type'] == "altaz_box": + runs = AltAzBoxGenerator.get_runs_from_config(cfg) + elif cfg['type'] == "fixed_altaz": + FixedObsGenerator.get_runs_from_config(cfg) + elif cfg['type'] == "data_matching": + DataMatchingGenerator.get_runs_from_config(cfg) + else: + raise ValueError(f"Unknown run generator type '{cfg['type']}'") + + return runs \ No newline at end of file diff --git a/src/srcsim/rungen/sky.py b/src/srcsim/rungen/sky.py index 126dd27..7f30ce5 100644 --- a/src/srcsim/rungen/sky.py +++ b/src/srcsim/rungen/sky.py @@ -10,24 +10,6 @@ from ..run import SkyDataRun -def generator(config): - if isinstance(config, str): - cfg = yaml.load(open(config, "r"), Loader=yaml.FullLoader) - else: - cfg = config - - runs = [] - - if cfg['type'] == "altaz_box": - runs = AltAzBoxGenerator.get_runs_from_config(cfg) - elif cfg['type'] == "data_matching": - DataMatchingGenerator.get_runs_from_config(cfg) - else: - raise ValueError(f"Unknown run generator type '{cfg['type']}'") - - return runs - - def get_trajectory(tel_pos, tstart, tstop, time_step, obsloc): times = Time( np.arange(tstart.unix, tstop.unix, step=time_step.to('s').value), diff --git a/src/srcsim/scripts/get_runs.py b/src/srcsim/scripts/get_runs.py index 3fd14b8..6d197de 100644 --- a/src/srcsim/scripts/get_runs.py +++ b/src/srcsim/scripts/get_runs.py @@ -2,7 +2,7 @@ import datetime import argparse -from srcsim.rungen.sky import generator as rungen +from srcsim.rungen.generator import generator as rungen def info_message(text): From 9c10702e50d0a5fd6cfb346acd765ef5e90bb31b Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Sat, 28 Oct 2023 17:43:52 +0900 Subject: [PATCH 14/28] rungen: moving extra functions to the "helpers" submodule --- src/srcsim/rungen/fixed.py | 2 +- src/srcsim/rungen/helpers.py | 62 ++++++++++++++++++++++++++++++++++++ src/srcsim/rungen/sky.py | 60 +--------------------------------- 3 files changed, 64 insertions(+), 60 deletions(-) create mode 100644 src/srcsim/rungen/helpers.py diff --git a/src/srcsim/rungen/fixed.py b/src/srcsim/rungen/fixed.py index 8526de7..5ffa3c6 100644 --- a/src/srcsim/rungen/fixed.py +++ b/src/srcsim/rungen/fixed.py @@ -6,8 +6,8 @@ from astropy.time import Time from astropy.coordinates import SkyCoord, SkyOffsetFrame, EarthLocation, AltAz -from srcsim.rungen import get_trajectory, get_time_intervals from srcsim.run import FixedPointingDataRun +from .helpers import get_trajectory class FixedObsGenerator: diff --git a/src/srcsim/rungen/helpers.py b/src/srcsim/rungen/helpers.py new file mode 100644 index 0000000..24d2088 --- /dev/null +++ b/src/srcsim/rungen/helpers.py @@ -0,0 +1,62 @@ +import numpy as np +from astropy.time import Time +from astropy.coordinates import AltAz + + +def get_trajectory(tel_pos, tstart, tstop, time_step, obsloc): + times = Time( + np.arange(tstart.unix, tstop.unix, step=time_step.to('s').value), + format='unix' + ) + frame = AltAz(obstime=times, location=obsloc) + tel_altaz = tel_pos.transform_to(frame) + + return tel_altaz + + +def enforce_max_interval_length(tstarts, tstops, max_length): + tstarts_new = [] + tstops_new = [] + + for tstart, tstop in zip(tstarts, tstops): + interval_duration = tstop - tstart + + if interval_duration > max_length: + time_edges = Time( + np.arange(tstart.unix, tstop.unix, step=max_length.to('s').value), + format='unix' + ) + if tstop not in time_edges: + time_edges = Time([time_edges, tstop]) + + for tmin, tmax in zip(time_edges[:-1], time_edges[1:]): + tstarts_new.append(tmin) + tstops_new.append(tmax) + + else: + tstarts_new.append(tstart) + tstops_new.append(tstop) + + return Time(tstarts_new), Time(tstops_new) + + +def get_time_intervals(tel_altaz, altmin, altmax, azmin, azmax, time_step, max_duration=None): + in_box = (tel_altaz.az > azmin) & (tel_altaz.az <= azmax) & (tel_altaz.alt > altmin) & (tel_altaz.alt < altmax) + nodes = np.where(np.diff(tel_altaz.obstime[in_box].unix) > time_step.to('s').value)[0] + + starts = np.concatenate(([0], nodes+1)) + + if (len(tel_altaz.obstime[in_box]) - 1) not in nodes: + stops = np.concatenate((nodes, [len(tel_altaz.obstime[in_box]) - 1])) + else: + stops = nodes + + intervals = tuple((start, stop) for start, stop in zip(starts, stops)) + + tstarts = Time([tel_altaz.obstime[in_box][interval[0]] for interval in intervals]) + tstops = Time([tel_altaz.obstime[in_box][interval[1]] for interval in intervals]) + + if max_duration is not None: + tstarts, tstops = enforce_max_interval_length(tstarts, tstops, max_duration) + + return tstarts, tstops \ No newline at end of file diff --git a/src/srcsim/rungen/sky.py b/src/srcsim/rungen/sky.py index 7f30ce5..76b8121 100644 --- a/src/srcsim/rungen/sky.py +++ b/src/srcsim/rungen/sky.py @@ -8,65 +8,7 @@ from astropy.coordinates import SkyCoord, AltAz, EarthLocation from ..run import SkyDataRun - - -def get_trajectory(tel_pos, tstart, tstop, time_step, obsloc): - times = Time( - np.arange(tstart.unix, tstop.unix, step=time_step.to('s').value), - format='unix' - ) - frame = AltAz(obstime=times, location=obsloc) - tel_altaz = tel_pos.transform_to(frame) - - return tel_altaz - - -def enforce_max_interval_length(tstarts, tstops, max_length): - tstarts_new = [] - tstops_new = [] - - for tstart, tstop in zip(tstarts, tstops): - interval_duration = tstop - tstart - - if interval_duration > max_length: - time_edges = Time( - np.arange(tstart.unix, tstop.unix, step=max_length.to('s').value), - format='unix' - ) - if tstop not in time_edges: - time_edges = Time([time_edges, tstop]) - - for tmin, tmax in zip(time_edges[:-1], time_edges[1:]): - tstarts_new.append(tmin) - tstops_new.append(tmax) - - else: - tstarts_new.append(tstart) - tstops_new.append(tstop) - - return Time(tstarts_new), Time(tstops_new) - - -def get_time_intervals(tel_altaz, altmin, altmax, azmin, azmax, time_step, max_duration=None): - in_box = (tel_altaz.az > azmin) & (tel_altaz.az <= azmax) & (tel_altaz.alt > altmin) & (tel_altaz.alt < altmax) - nodes = np.where(np.diff(tel_altaz.obstime[in_box].unix) > time_step.to('s').value)[0] - - starts = np.concatenate(([0], nodes+1)) - - if (len(tel_altaz.obstime[in_box]) - 1) not in nodes: - stops = np.concatenate((nodes, [len(tel_altaz.obstime[in_box]) - 1])) - else: - stops = nodes - - intervals = tuple((start, stop) for start, stop in zip(starts, stops)) - - tstarts = Time([tel_altaz.obstime[in_box][interval[0]] for interval in intervals]) - tstops = Time([tel_altaz.obstime[in_box][interval[1]] for interval in intervals]) - - if max_duration is not None: - tstarts, tstops = enforce_max_interval_length(tstarts, tstops, max_duration) - - return tstarts, tstops +from .helpers import get_trajectory, get_time_intervals class AltAzBoxGenerator: From 26fa602ff8f8643e3a5a39aa8efbf161c604e011 Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Sat, 28 Oct 2023 17:44:23 +0900 Subject: [PATCH 15/28] Dropping unused imports. --- src/srcsim/rungen/fixed.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/srcsim/rungen/fixed.py b/src/srcsim/rungen/fixed.py index 5ffa3c6..473447e 100644 --- a/src/srcsim/rungen/fixed.py +++ b/src/srcsim/rungen/fixed.py @@ -1,10 +1,9 @@ import yaml import numpy as np -import pandas as pd import astropy.units as u from scipy.interpolate import CubicSpline from astropy.time import Time -from astropy.coordinates import SkyCoord, SkyOffsetFrame, EarthLocation, AltAz +from astropy.coordinates import SkyCoord, EarthLocation from srcsim.run import FixedPointingDataRun from .helpers import get_trajectory From b1a287327310cb4d92703c86e4c7769347541726 Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Sat, 28 Oct 2023 17:44:48 +0900 Subject: [PATCH 16/28] Relative import to highlight it is within the same module. --- src/srcsim/rungen/fixed.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/srcsim/rungen/fixed.py b/src/srcsim/rungen/fixed.py index 473447e..8fff1f4 100644 --- a/src/srcsim/rungen/fixed.py +++ b/src/srcsim/rungen/fixed.py @@ -5,7 +5,7 @@ from astropy.time import Time from astropy.coordinates import SkyCoord, EarthLocation -from srcsim.run import FixedPointingDataRun +from ..run import FixedPointingDataRun from .helpers import get_trajectory From fb5e33ea77e2490deb3d9da6540047d839d4d35e Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Sun, 29 Oct 2023 18:59:54 +0900 Subject: [PATCH 17/28] "rungen" submodule global imports --- src/srcsim/rungen/__init__.py | 3 +++ src/srcsim/scripts/get_runs.py | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) create mode 100644 src/srcsim/rungen/__init__.py diff --git a/src/srcsim/rungen/__init__.py b/src/srcsim/rungen/__init__.py new file mode 100644 index 0000000..2845ee4 --- /dev/null +++ b/src/srcsim/rungen/__init__.py @@ -0,0 +1,3 @@ +from .generator import generator +from .sky import AltAzBoxGenerator, DataMatchingGenerator +from .fixed import FixedObsGenerator \ No newline at end of file diff --git a/src/srcsim/scripts/get_runs.py b/src/srcsim/scripts/get_runs.py index 6d197de..47ea066 100644 --- a/src/srcsim/scripts/get_runs.py +++ b/src/srcsim/scripts/get_runs.py @@ -2,7 +2,7 @@ import datetime import argparse -from srcsim.rungen.generator import generator as rungen +from srcsim.rungen import generator as rungen def info_message(text): From 8cf6afe20252c2b27ec7b26953066b5fc3bbc611 Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Sun, 29 Oct 2023 19:02:04 +0900 Subject: [PATCH 18/28] Fixing wrong units used --- src/srcsim/rungen/fixed.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/srcsim/rungen/fixed.py b/src/srcsim/rungen/fixed.py index 8fff1f4..b500f23 100644 --- a/src/srcsim/rungen/fixed.py +++ b/src/srcsim/rungen/fixed.py @@ -38,7 +38,7 @@ def get_runs_from_config(cls, config): @classmethod def get_runs(cls, obsloc, tel_pos, tobs, alt, ra_width, tstart=None): - time_step = 1 * u.m + time_step = 1 * u.minute tel_pos_trail = SkyCoord( tel_pos.icrs.ra + ra_width / 2, From d856694f57189b79a02a115ddb8949aae9e092fb Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Sun, 29 Oct 2023 19:42:50 +0900 Subject: [PATCH 19/28] FixedObsGenerator: using alt/az instead of equatorial position when generating runs. --- src/srcsim/rungen/fixed.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/src/srcsim/rungen/fixed.py b/src/srcsim/rungen/fixed.py index b500f23..fac55d9 100644 --- a/src/srcsim/rungen/fixed.py +++ b/src/srcsim/rungen/fixed.py @@ -3,7 +3,7 @@ import astropy.units as u from scipy.interpolate import CubicSpline from astropy.time import Time -from astropy.coordinates import SkyCoord, EarthLocation +from astropy.coordinates import SkyCoord, AltAz, EarthLocation from ..run import FixedPointingDataRun from .helpers import get_trajectory @@ -84,7 +84,7 @@ def get_runs(cls, obsloc, tel_pos, tobs, alt, ra_width, tstart=None): ndays_full = np.floor(nsequences) ndays_total = np.ceil(nsequences) - if ndays_full > 0: + if ndays_full > 1: # Next pass - to cover the simulation time interval with fully completed runs track_lead = get_trajectory( tel_pos_lead.icrs, @@ -135,14 +135,14 @@ def get_runs(cls, obsloc, tel_pos, tobs, alt, ra_width, tstart=None): track_lead = get_trajectory( tel_pos_lead.icrs, tstart, - tstop=tstart + ndays_full * u.d, + tstop=tstart + ndays_total * u.d, time_step=time_step, obsloc=obsloc ) track_trail = get_trajectory( tel_pos_trail.icrs, tstart, - tstop=tstart + ndays_full * u.d, + tstop=tstart + ndays_total * u.d, time_step=time_step, obsloc=obsloc ) @@ -158,8 +158,14 @@ def get_runs(cls, obsloc, tel_pos, tobs, alt, ra_width, tstart=None): tstarts = Time(cs_lead.roots(extrapolate=False), format='mjd') tstops = Time(tstarts + tobs / (len(tstarts))) + frame = AltAz( + obstime=tstarts[0], + location=obsloc + ) + tel_pos_altaz = tel_pos_lead.transform_to(frame) + runs = tuple( - FixedPointingDataRun(tel_pos, tstart, tstop, obsloc, run_id) + FixedPointingDataRun(tel_pos_altaz, tstart, tstop, obsloc, run_id) for run_id, (tstart, tstop) in enumerate(zip(tstarts, tstops)) ) From bd780fe8390344bb4c704af96adc78fc3ea61c43 Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Sun, 29 Oct 2023 19:45:11 +0900 Subject: [PATCH 20/28] Fix: ensuring generator returns created runs. --- src/srcsim/rungen/generator.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/srcsim/rungen/generator.py b/src/srcsim/rungen/generator.py index 62e299e..98b57d3 100644 --- a/src/srcsim/rungen/generator.py +++ b/src/srcsim/rungen/generator.py @@ -15,9 +15,9 @@ def generator(config): if cfg['type'] == "altaz_box": runs = AltAzBoxGenerator.get_runs_from_config(cfg) elif cfg['type'] == "fixed_altaz": - FixedObsGenerator.get_runs_from_config(cfg) + runs = FixedObsGenerator.get_runs_from_config(cfg) elif cfg['type'] == "data_matching": - DataMatchingGenerator.get_runs_from_config(cfg) + runs = DataMatchingGenerator.get_runs_from_config(cfg) else: raise ValueError(f"Unknown run generator type '{cfg['type']}'") From a7618494eef9bf04ef89f2b783f9c0c3f3952efc Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Mon, 30 Oct 2023 14:24:09 +0900 Subject: [PATCH 21/28] Run type detection when generating from config. --- src/srcsim/run/__init__.py | 3 ++- src/srcsim/run/generator.py | 10 ++++++++++ src/srcsim/scripts/sim_run.py | 4 ++-- 3 files changed, 14 insertions(+), 3 deletions(-) create mode 100644 src/srcsim/run/generator.py diff --git a/src/srcsim/run/__init__.py b/src/srcsim/run/__init__.py index dda540e..11cab66 100644 --- a/src/srcsim/run/__init__.py +++ b/src/srcsim/run/__init__.py @@ -1,2 +1,3 @@ from .sky import SkyDataRun -from .fixed import FixedPointingDataRun \ No newline at end of file +from .fixed import FixedPointingDataRun +from .generator import generator \ No newline at end of file diff --git a/src/srcsim/run/generator.py b/src/srcsim/run/generator.py new file mode 100644 index 0000000..6af3b90 --- /dev/null +++ b/src/srcsim/run/generator.py @@ -0,0 +1,10 @@ +from srcsim.run import SkyDataRun, FixedPointingDataRun + + +def generator(cfg): + if 'alt' in cfg['pointing']: + run = FixedPointingDataRun.from_config(cfg) + else: + run = SkyDataRun.from_config(cfg) + + return run \ No newline at end of file diff --git a/src/srcsim/scripts/sim_run.py b/src/srcsim/scripts/sim_run.py index 590f1f0..14e30ab 100644 --- a/src/srcsim/scripts/sim_run.py +++ b/src/srcsim/scripts/sim_run.py @@ -7,7 +7,7 @@ from srcsim.mc import MCCollection from srcsim.src import generator as srcgen -from srcsim.run import SkyDataRun +from srcsim.run import generator as rungen def info_message(text): @@ -80,7 +80,7 @@ def main(): print(srcs) info_message('Preparing the data run') - run = SkyDataRun.from_config(cfg['run']) + run = rungen(cfg['run']) print(run) info_message('Starting event sampling') From f612b46dcdc2f45f9d9201d9f335ce317849907c Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Thu, 2 Nov 2023 12:39:45 +0900 Subject: [PATCH 22/28] FixedObsGenerator: adding accuracy and max run duration reading; grouping time settings. --- src/srcsim/rungen/fixed.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/src/srcsim/rungen/fixed.py b/src/srcsim/rungen/fixed.py index fac55d9..057dc15 100644 --- a/src/srcsim/rungen/fixed.py +++ b/src/srcsim/rungen/fixed.py @@ -17,7 +17,6 @@ def get_runs_from_config(cls, config): else: cfg = config - tstart = Time(cfg['time']['start']) ra_width = u.Quantity(cfg['pointing']['width']) alt = u.Quantity(cfg['pointing']['center']['alt']) tel_pos = SkyCoord( @@ -26,7 +25,10 @@ def get_runs_from_config(cls, config): frame='icrs' ) + tstart = Time(cfg['time']['start']) duration = u.Quantity(cfg['time']['duration']) + accuracy = u.Quantity(cfg['time']['accuracy']) + max_run_duration = u.Quantity(cfg['time']['max_run_duration']) if cfg['time']['max_run_duration'] is not None else None obsloc = EarthLocation( lat=u.Quantity(cfg['location']['lat']), @@ -34,10 +36,10 @@ def get_runs_from_config(cls, config): height=u.Quantity(cfg['location']['height']), ) - return cls.get_runs(obsloc, tel_pos, duration, alt, ra_width, tstart) + return cls.get_runs(obsloc, tel_pos, duration, alt, ra_width, tstart, accuracy, max_run_duration) @classmethod - def get_runs(cls, obsloc, tel_pos, tobs, alt, ra_width, tstart=None): + def get_runs(cls, obsloc, tel_pos, tobs, alt, ra_width, tstart=None, accuracy=1*u.minute, max_run_duration=None): time_step = 1 * u.minute tel_pos_trail = SkyCoord( From 18180aaacdf3a2cb5be15dd00205c9fe115109f0 Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Thu, 2 Nov 2023 12:42:08 +0900 Subject: [PATCH 23/28] FixedObsGenerator: configurable time step size. --- src/srcsim/rungen/fixed.py | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/src/srcsim/rungen/fixed.py b/src/srcsim/rungen/fixed.py index 057dc15..0cf2bdd 100644 --- a/src/srcsim/rungen/fixed.py +++ b/src/srcsim/rungen/fixed.py @@ -40,8 +40,6 @@ def get_runs_from_config(cls, config): @classmethod def get_runs(cls, obsloc, tel_pos, tobs, alt, ra_width, tstart=None, accuracy=1*u.minute, max_run_duration=None): - time_step = 1 * u.minute - tel_pos_trail = SkyCoord( tel_pos.icrs.ra + ra_width / 2, tel_pos.icrs.dec, @@ -58,14 +56,14 @@ def get_runs(cls, obsloc, tel_pos, tobs, alt, ra_width, tstart=None, accuracy=1* tel_pos_trail.icrs, tstart, tstop=tstart + 1 * u.d, - time_step=time_step, + time_step=accuracy, obsloc=obsloc ) track_lead = get_trajectory( tel_pos_lead.icrs, tstart, tstop=tstart + 1 * u.d, - time_step=time_step, + time_step=accuracy, obsloc=obsloc ) cs_lead = CubicSpline( @@ -92,14 +90,14 @@ def get_runs(cls, obsloc, tel_pos, tobs, alt, ra_width, tstart=None, accuracy=1* tel_pos_lead.icrs, tstart, tstop=tstart + ndays_full * u.d, - time_step=time_step, + time_step=accuracy, obsloc=obsloc ) track_trail = get_trajectory( tel_pos_trail.icrs, tstart, tstop=tstart + ndays_full * u.d, - time_step=time_step, + time_step=accuracy, obsloc=obsloc ) @@ -120,7 +118,7 @@ def get_runs(cls, obsloc, tel_pos, tobs, alt, ra_width, tstart=None, accuracy=1* tel_pos_lead.icrs, tstart=tstart + ndays_full * u.d, tstop=tstart + ndays_total * u.d, - time_step=time_step, + time_step=accuracy, obsloc=obsloc ) cs_lead = CubicSpline( @@ -138,14 +136,14 @@ def get_runs(cls, obsloc, tel_pos, tobs, alt, ra_width, tstart=None, accuracy=1* tel_pos_lead.icrs, tstart, tstop=tstart + ndays_total * u.d, - time_step=time_step, + time_step=accuracy, obsloc=obsloc ) track_trail = get_trajectory( tel_pos_trail.icrs, tstart, tstop=tstart + ndays_total * u.d, - time_step=time_step, + time_step=accuracy, obsloc=obsloc ) From c5d13dd519273b7f3b629efa7f190ca30b17b8b7 Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Thu, 2 Nov 2023 12:45:13 +0900 Subject: [PATCH 24/28] FixedObsGenerator: enforcing maximal run duration. --- src/srcsim/rungen/fixed.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/srcsim/rungen/fixed.py b/src/srcsim/rungen/fixed.py index 0cf2bdd..ce8d8fa 100644 --- a/src/srcsim/rungen/fixed.py +++ b/src/srcsim/rungen/fixed.py @@ -6,7 +6,7 @@ from astropy.coordinates import SkyCoord, AltAz, EarthLocation from ..run import FixedPointingDataRun -from .helpers import get_trajectory +from .helpers import get_trajectory, enforce_max_interval_length class FixedObsGenerator: @@ -158,6 +158,8 @@ def get_runs(cls, obsloc, tel_pos, tobs, alt, ra_width, tstart=None, accuracy=1* tstarts = Time(cs_lead.roots(extrapolate=False), format='mjd') tstops = Time(tstarts + tobs / (len(tstarts))) + tstarts, tstops = enforce_max_interval_length(tstarts, tstops, max_run_duration) + frame = AltAz( obstime=tstarts[0], location=obsloc From 1768ac1eada5933f336e3bec0443ffc2e98ecb8b Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Tue, 7 Nov 2023 18:39:56 +0900 Subject: [PATCH 25/28] Bug fix in FixedObsGenerator: using the right reference times in the final Alt/Az transformation (fixed time was leading to a wrong RA assigned to some of the runs). --- src/srcsim/rungen/fixed.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/srcsim/rungen/fixed.py b/src/srcsim/rungen/fixed.py index ce8d8fa..e19da29 100644 --- a/src/srcsim/rungen/fixed.py +++ b/src/srcsim/rungen/fixed.py @@ -161,14 +161,14 @@ def get_runs(cls, obsloc, tel_pos, tobs, alt, ra_width, tstart=None, accuracy=1* tstarts, tstops = enforce_max_interval_length(tstarts, tstops, max_run_duration) frame = AltAz( - obstime=tstarts[0], + obstime=tstarts, location=obsloc ) tel_pos_altaz = tel_pos_lead.transform_to(frame) runs = tuple( - FixedPointingDataRun(tel_pos_altaz, tstart, tstop, obsloc, run_id) - for run_id, (tstart, tstop) in enumerate(zip(tstarts, tstops)) + FixedPointingDataRun(target_altaz, tstart, tstop, obsloc, run_id) + for run_id, (tstart, tstop, target_altaz) in enumerate(zip(tstarts, tstops, tel_pos_altaz)) ) return runs From 01a48812c8753499c8476ed203f374b94c142efc Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Wed, 8 Nov 2023 17:35:30 +0900 Subject: [PATCH 26/28] FixedObsGenerator: ensuring short observations will get as much as possible scheduled on the first day --- src/srcsim/rungen/fixed.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/srcsim/rungen/fixed.py b/src/srcsim/rungen/fixed.py index e19da29..a4730e8 100644 --- a/src/srcsim/rungen/fixed.py +++ b/src/srcsim/rungen/fixed.py @@ -84,7 +84,7 @@ def get_runs(cls, obsloc, tel_pos, tobs, alt, ra_width, tstart=None, accuracy=1* ndays_full = np.floor(nsequences) ndays_total = np.ceil(nsequences) - if ndays_full > 1: + if ndays_total > 1: # Next pass - to cover the simulation time interval with fully completed runs track_lead = get_trajectory( tel_pos_lead.icrs, From 983109356705ad2adbce00bd34b1e09208d18282 Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Wed, 8 Nov 2023 17:38:41 +0900 Subject: [PATCH 27/28] FixedObsGenerator: fixing a likely astropy bug --- src/srcsim/rungen/fixed.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/srcsim/rungen/fixed.py b/src/srcsim/rungen/fixed.py index a4730e8..6f26060 100644 --- a/src/srcsim/rungen/fixed.py +++ b/src/srcsim/rungen/fixed.py @@ -128,8 +128,12 @@ def get_runs(cls, obsloc, tel_pos, tobs, alt, ra_width, tstart=None, accuracy=1* _tstarts = Time(cs_lead.roots(extrapolate=False), format='mjd') _tstops = Time(_tstarts + remaining_tobs / (len(_tstarts))) - tstarts = Time([tstarts, _tstarts]) - tstops = Time([tstops, _tstops]) + # Likely an astropy bug (tested with v5.3.4): + # the latter does not work if len(tstarts) == len(_tstarts) + # tstarts = Time([tstarts, _tstarts]) + # tstops = Time([tstops, _tstops]) + tstarts = Time(np.concatenate([tstarts.mjd, _tstarts.mjd]), format='mjd') + tstops = Time(np.concatenate([tstops.mjd, _tstops.mjd]), format='mjd') else: track_lead = get_trajectory( From 238a8c31c07db5f6bf8beb24dacde6d77d171bac Mon Sep 17 00:00:00 2001 From: Ievgen Vovk Date: Wed, 8 Nov 2023 18:33:32 +0900 Subject: [PATCH 28/28] FixedObsGenerator: fixing the Alt/Az assigned to the run as using tel_pos_lead directly leads to wrong position once intervals are sliced and tel_pos_lead no longer corresponds to the interval start. --- src/srcsim/rungen/fixed.py | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/srcsim/rungen/fixed.py b/src/srcsim/rungen/fixed.py index 6f26060..d9bb660 100644 --- a/src/srcsim/rungen/fixed.py +++ b/src/srcsim/rungen/fixed.py @@ -162,13 +162,21 @@ def get_runs(cls, obsloc, tel_pos, tobs, alt, ra_width, tstart=None, accuracy=1* tstarts = Time(cs_lead.roots(extrapolate=False), format='mjd') tstops = Time(tstarts + tobs / (len(tstarts))) - tstarts, tstops = enforce_max_interval_length(tstarts, tstops, max_run_duration) - - frame = AltAz( + # Telescope positions before the time interevals will be sliced + # These should contain only two different values - those before and after culmination + _frame = AltAz( obstime=tstarts, location=obsloc ) - tel_pos_altaz = tel_pos_lead.transform_to(frame) + _tel_pos_altaz = tel_pos_lead.transform_to(_frame) + + tstarts, tstops = enforce_max_interval_length(tstarts, tstops, max_run_duration) + + # Choosing the closest Alt/Az telescope position from those before interval slicing + tel_pos_altaz = [ + _tel_pos_altaz[abs(_tel_pos_altaz.obstime - tstart).argmin()] + for tstart in tstarts + ] runs = tuple( FixedPointingDataRun(target_altaz, tstart, tstop, obsloc, run_id)