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" }, 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/run/__init__.py b/src/srcsim/run/__init__.py new file mode 100644 index 0000000..11cab66 --- /dev/null +++ b/src/srcsim/run/__init__.py @@ -0,0 +1,3 @@ +from .sky import SkyDataRun +from .fixed import FixedPointingDataRun +from .generator import generator \ No newline at end of file diff --git a/src/srcsim/run.py b/src/srcsim/run/base.py similarity index 66% rename from src/srcsim/run.py rename to src/srcsim/run/base.py index f089209..07d0745 100644 --- a/src/srcsim/run.py +++ b/src/srcsim/run/base.py @@ -1,10 +1,9 @@ -import yaml import logging 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, SkyOffsetFrame, AltAz class DataRun: @@ -19,47 +18,16 @@ def __init__(self, tel_pos, tstart, tstop, obsloc, id=0, log=None): self.log = logging.getLogger(__name__) else: self.log = log.getChild(__name__) - - def __repr__(self): - frame_start = AltAz(obstime=self.tstart, location=self.obsloc) - frame_stop = AltAz(obstime=self.tstop, location=self.obsloc) - print( -f"""{type(self).__name__} instance - {'ID':.<20s}: {self.id} - {'Tel. RA/Dec':.<20s}: {self.tel_pos} - {'Tstart':.<20s}: {self.tstart.isot} - {'Tstop':.<20s}: {self.tstop.isot} - {'Tel. azimuth':.<20s}: [{self.tel_pos.transform_to(frame_start).az.to('deg'):.2f} - {self.tel_pos.transform_to(frame_stop).az.to('deg'):.2f}] - {'Tel. alt':.<20s}: [{self.tel_pos.transform_to(frame_start).alt.to('deg'):.2f} - {self.tel_pos.transform_to(frame_stop).alt.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 - - data_run = cls( - SkyCoord( - u.Quantity(cfg['pointing']['ra']), - u.Quantity(cfg['pointing']['dec']), - frame='icrs' - ), - 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 - ) + pass - return data_run + def to_dict(self): + pass + + def tel_pos_to_altaz(self, frame): + pass @classmethod def time_sort(cls, events): @@ -89,21 +57,6 @@ def update_time_delta(cls, events): return events - def to_dict(self): - data = {'id': self.id, 'pointing': {}, 'time': {}, 'location': {}} - - data['pointing']['ra'] = str(self.tel_pos.icrs.ra.to('deg').value) + ' deg' - data['pointing']['dec'] = str(self.tel_pos.icrs.dec.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 predict(self, mccollections, source, tel_pos_tolerance=None, time_step=1*u.minute): self.log.debug(f'predicting events for {source.name}') @@ -119,14 +72,14 @@ def predict(self, mccollections, source, tel_pos_tolerance=None, time_step=1*u.m 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) - + 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): @@ -151,7 +104,7 @@ def predict(self, mccollections, source, tel_pos_tolerance=None, time_step=1*u.m obstime=Time(arrival_time, format='unix'), location=self.obsloc ) - current_tel_pos = self.tel_pos.transform_to(current_frame) + 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 @@ -172,33 +125,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, @@ -211,8 +165,8 @@ def predict(self, mccollections, source, tel_pos_tolerance=None, time_step=1*u.m 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 + ra_tel = current_tel_pos.icrs.ra.to('rad').value, + dec_tel = current_tel_pos.icrs.dec.to('rad').value ) # Reconstructed events coordinates @@ -228,6 +182,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), diff --git a/src/srcsim/run/fixed.py b/src/srcsim/run/fixed.py new file mode 100644 index 0000000..eee1c68 --- /dev/null +++ b/src/srcsim/run/fixed.py @@ -0,0 +1,99 @@ +import yaml +import numpy as np +import astropy.units as u +from astropy.time import Time +from astropy.coordinates import SkyCoord, EarthLocation, AltAz + +from .base import DataRun + + +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 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/run/sky.py b/src/srcsim/run/sky.py new file mode 100644 index 0000000..5d91a46 --- /dev/null +++ b/src/srcsim/run/sky.py @@ -0,0 +1,67 @@ +import yaml +import astropy.units as u +from astropy.time import Time +from astropy.coordinates import SkyCoord, EarthLocation, AltAz + +from .base import DataRun + + +class SkyDataRun(DataRun): + def __repr__(self): + frame_start = AltAz(obstime=self.tstart, location=self.obsloc) + frame_stop = AltAz(obstime=self.tstop, location=self.obsloc) + print( +f"""{type(self).__name__} instance + {'ID':.<20s}: {self.id} + {'Tel. RA/Dec':.<20s}: {self.tel_pos} + {'Tstart':.<20s}: {self.tstart.isot} + {'Tstop':.<20s}: {self.tstop.isot} + {'Tel. azimuth':.<20s}: [{self.tel_pos.transform_to(frame_start).az.to('deg'):.2f} - {self.tel_pos.transform_to(frame_stop).az.to('deg'):.2f}] + {'Tel. alt':.<20s}: [{self.tel_pos.transform_to(frame_start).alt.to('deg'):.2f} - {self.tel_pos.transform_to(frame_stop).alt.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 + + data_run = cls( + SkyCoord( + u.Quantity(cfg['pointing']['ra']), + u.Quantity(cfg['pointing']['dec']), + frame='icrs' + ), + 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']['ra'] = str(self.tel_pos.icrs.ra.to('deg').value) + ' deg' + data['pointing']['dec'] = str(self.tel_pos.icrs.dec.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): + return self.tel_pos.transform_to(frame) 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/rungen/fixed.py b/src/srcsim/rungen/fixed.py new file mode 100644 index 0000000..d9bb660 --- /dev/null +++ b/src/srcsim/rungen/fixed.py @@ -0,0 +1,186 @@ +import yaml +import numpy as np +import astropy.units as u +from scipy.interpolate import CubicSpline +from astropy.time import Time +from astropy.coordinates import SkyCoord, AltAz, EarthLocation + +from ..run import FixedPointingDataRun +from .helpers import get_trajectory, enforce_max_interval_length + + +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 + + 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' + ) + + 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']), + lon=u.Quantity(cfg['location']['lon']), + height=u.Quantity(cfg['location']['height']), + ) + + 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, accuracy=1*u.minute, max_run_duration=None): + 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=accuracy, + obsloc=obsloc + ) + track_lead = get_trajectory( + tel_pos_lead.icrs, + tstart, + tstop=tstart + 1 * u.d, + time_step=accuracy, + 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_total > 1: + # 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=accuracy, + obsloc=obsloc + ) + track_trail = get_trajectory( + tel_pos_trail.icrs, + tstart, + tstop=tstart + ndays_full * u.d, + time_step=accuracy, + 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=accuracy, + 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))) + + # 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( + tel_pos_lead.icrs, + tstart, + tstop=tstart + ndays_total * u.d, + time_step=accuracy, + obsloc=obsloc + ) + track_trail = get_trajectory( + tel_pos_trail.icrs, + tstart, + tstop=tstart + ndays_total * u.d, + time_step=accuracy, + 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))) + + # 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) + + 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) + for run_id, (tstart, tstop, target_altaz) in enumerate(zip(tstarts, tstops, tel_pos_altaz)) + ) + + return runs diff --git a/src/srcsim/rungen/generator.py b/src/srcsim/rungen/generator.py new file mode 100644 index 0000000..98b57d3 --- /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": + runs = FixedObsGenerator.get_runs_from_config(cfg) + elif cfg['type'] == "data_matching": + runs = 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/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.py b/src/srcsim/rungen/sky.py similarity index 74% rename from src/srcsim/rungen.py rename to src/srcsim/rungen/sky.py index 63f32f6..76b8121 100644 --- a/src/srcsim/rungen.py +++ b/src/srcsim/rungen/sky.py @@ -7,84 +7,8 @@ from astropy.time import Time from astropy.coordinates import SkyCoord, AltAz, EarthLocation -from .run import DataRun - - -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), - 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 ..run import SkyDataRun +from .helpers import get_trajectory, get_time_intervals class AltAzBoxGenerator: @@ -205,7 +129,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 +141,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 +175,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 2f4ec6b..e166d79 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 generator as rungen def main(): @@ -80,7 +80,7 @@ def main(): print(srcs) log.info('preparing the data run') - run = DataRun.from_config(cfg['run']) + run = rungen(cfg['run']) print(run) log.info('starting event sampling')