diff --git a/cosipy/response/FullDetectorResponse.py b/cosipy/response/FullDetectorResponse.py index e4a10050..a889c7a4 100644 --- a/cosipy/response/FullDetectorResponse.py +++ b/cosipy/response/FullDetectorResponse.py @@ -838,7 +838,8 @@ def get_interp_response(self, coord): def get_point_source_response(self, exposure_map = None, coord = None, - scatt_map = None): + scatt_map = None, + Earth_occ = True): """ Convolve the all-sky detector response with exposure for a source at a given sky location. @@ -854,6 +855,10 @@ def get_point_source_response(self, Source coordinate scatt_map : :py:class:`SpacecraftAttitudeMap` Spacecraft attitude map + Earth_occ : bool, optional + Option to include Earth occultation in the respeonce. + Default is True, in which case you can only pass one + coord, which must be the same as was used for the scatt map. Returns ------- @@ -863,6 +868,11 @@ def get_point_source_response(self, # TODO: deprecate exposure_map in favor of coords + scatt map for both local # and interntial coords + if Earth_occ == True: + if coord != None: + if coord.size > 1: + raise ValueError("For Earth occultation you must use the same coordinate as was used for the scatt map!") + if exposure_map is not None: if not self.conformable(exposure_map): raise ValueError( diff --git a/cosipy/spacecraftfile/SpacecraftFile.py b/cosipy/spacecraftfile/SpacecraftFile.py index ce2e3d69..4615f149 100644 --- a/cosipy/spacecraftfile/SpacecraftFile.py +++ b/cosipy/spacecraftfile/SpacecraftFile.py @@ -8,6 +8,7 @@ import matplotlib.pyplot as plt from matplotlib.colors import LogNorm from matplotlib import cm, colors +from scipy import interpolate from scoords import Attitude, SpacecraftFrame from cosipy.response import FullDetectorResponse @@ -19,24 +20,39 @@ class SpacecraftFile(): - def __init__(self, time, x_pointings = None, y_pointings = None, z_pointings = None, attitude = None, - instrument = "COSI", frame = "galactic"): + def __init__(self, time, x_pointings = None, y_pointings = None, \ + z_pointings = None, earth_zenith = None, altitude = None,\ + attitude = None, instrument = "COSI", frame = "galactic"): """ - Handles the spacecraft orientation. Calculates the dwell time map and point source response over a certain orientation period. Exports the point source response as RMF and ARF files that can be read by XSPEC. + Handles the spacecraft orientation. Calculates the dwell time + map and point source response over a certain orientation period. + Exports the point source response as RMF and ARF files that can be read by XSPEC. Parameters ---------- Time : astropy.time.Time The time stamps for each pointings. Note this is NOT the time duration. x_pointings : astropy.coordinates.SkyCoord, optional - The pointings (galactic system) of the x axis of the local coordinate system attached to the spacecraft (the default is `None`, which implies no input for the x pointings). + The pointings (galactic system) of the x axis of the local + coordinate system attached to the spacecraft (the default + is `None`, which implies no input for the x pointings). y_pointings : astropy.coordinates.SkyCoord, optional - The pointings (galactic system) of the y axis of the local coordinate system attached to the spacecraft (the default is `None`, which implies no input for the y pointings). + The pointings (galactic system) of the y axis of the local + coordinate system attached to the spacecraft (the default + is `None`, which implies no input for the y pointings). z_pointings : astropy.coordinates.SkyCoord, optional - The pointings (galactic system) of the z axis of the local coordinate system attached to the spacecraft (the default is `None`, which implies no input for the z pointings). - attitude: numpy.ndarray, optional - The attitude of the spacecraft (the default is `None`, which implies no input for the attitude of the spacecraft). + The pointings (galactic system) of the z axis of the local + coordinate system attached to the spacecraft (the default + is `None`, which implies no input for the z pointings). + earth_zenith : astropy.coordinates.SkyCoord, optional + The pointings (galactic system) of the Earth zenith (the + default is `None`, which implies no input for the earth pointings). + altitude : array, optional + Altitude of the spacecraft in km. + attitude : numpy.ndarray, optional + The attitude of the spacecraft (the default is `None`, + which implies no input for the attitude of the spacecraft). instrument : str, optional The instrument name (the default is "COSI"). frame : str, optional @@ -50,6 +66,10 @@ def __init__(self, time, x_pointings = None, y_pointings = None, z_pointings = N else: raise TypeError("The time should be a astropy.time.Time object") + # Altitude + if not isinstance(altitude, (type(None))): + self._altitude = np.array(altitude) + # x pointings if isinstance(x_pointings, (SkyCoord, type(None))): self.x_pointings = x_pointings @@ -67,6 +87,12 @@ def __init__(self, time, x_pointings = None, y_pointings = None, z_pointings = N self.z_pointings = z_pointings else: raise TypeError("The z_pointing should be a NoneType or SkyCoord object!") + + # earth pointings + if isinstance(earth_zenith, (SkyCoord, type(None))): + self.earth_zenith = earth_zenith + else: + raise TypeError("The earth_zenith should be a NoneType or SkyCoord object!") # check if the x, y and z pointings are all None (no inputs). If all None, tt will try to read from attitude parameter if self.x_pointings is None and self.y_pointings is None and self.z_pointings is None: @@ -84,9 +110,10 @@ def __init__(self, time, x_pointings = None, y_pointings = None, z_pointings = N self._load_time = self._time.to_value(format = "unix") # this is not necessary, but just to make sure evething works fine... self._x_direction = np.array([x_pointings.l.deg, x_pointings.b.deg]).T # this is not necessary, but just to make sure evething works fine... self._z_direction = np.array([z_pointings.l.deg, z_pointings.b.deg]).T # this is not necessary, but just to make sure evething works fine... + self._earth_direction = np.array([earth_zenith.l.deg, earth_zenith.b.deg]).T # this is not necessary, but just to make sure evething works fine... + self.frame = frame - - + @classmethod def parse_from_file(cls, file): @@ -104,22 +131,24 @@ def parse_from_file(cls, file): The SpacecraftFile object. """ - - orientation_file = np.loadtxt(file, usecols=(1, 2, 3, 4, 5), delimiter=' ', skiprows=1, comments=("#", "EN")) + orientation_file = np.loadtxt(file, usecols=(1, 2, 3, 4, 5, 6, 7, 8),delimiter=' ', skiprows=1, comments=("#", "EN")) time_stamps = orientation_file[:, 0] axis_1 = orientation_file[:, [2, 1]] axis_2 = orientation_file[:, [4, 3]] - + axis_3 = orientation_file[:, [7, 6]] + altitude = np.array(orientation_file[:, 5]) + time = Time(time_stamps, format = "unix") xpointings = SkyCoord(l = axis_1[:,0]*u.deg, b = axis_1[:,1]*u.deg, frame = "galactic") zpointings = SkyCoord(l = axis_2[:,0]*u.deg, b = axis_2[:,1]*u.deg, frame = "galactic") - - return cls(time, x_pointings = xpointings, z_pointings = zpointings) + earthpointings = SkyCoord(l = axis_3[:,0]*u.deg, b = axis_3[:,1]*u.deg, frame = "galactic") + + return cls(time, x_pointings = xpointings, z_pointings = zpointings, earth_zenith = earthpointings, altitude = altitude) def get_time(self, time_array = None): """ - Return the arrary pf pointing times as a astropy.Time object. + Return the array pf pointing times as a astropy.Time object. Parameters ---------- @@ -139,6 +168,21 @@ def get_time(self, time_array = None): return self._time + def get_altitude(self): + + """ + Return the array of Earth altitude. + + + + Returns + ------- + numpy array + the Earth altitude. + """ + + return self._altitude + def get_time_delta(self, time_array = None): """ @@ -201,7 +245,7 @@ def source_interval(self, start, stop): Parameters ---------- start : astropy.time.Time - The star time of the orientation period. + The start time of the orientation period. stop : astropy.time.Time The end time of the orientation period. @@ -224,11 +268,15 @@ def source_interval(self, start, stop): new_times = self._load_time[start_idx : stop_idx + 1] new_x_direction = self._x_direction[start_idx : stop_idx + 1] new_z_direction = self._z_direction[start_idx : stop_idx + 1] + new_earth_direction = self._earth_direction[start_idx : stop_idx + 1] + new_earth_altitude = self._altitude[start_idx : stop_idx + 1] + else: start_idx = self._load_time.searchsorted(start.value) - 1 x_direction_start = self.interpolate_direction(start, start_idx, self._x_direction) z_direction_start = self.interpolate_direction(start, start_idx, self._z_direction) + earth_direction_start = self.interpolate_direction(start, start_idx, self._earth_direction) new_times = self._load_time[start_idx + 1 : stop_idx + 1] new_times = np.insert(new_times, 0, start.value) @@ -238,6 +286,15 @@ def source_interval(self, start, stop): new_z_direction = self._z_direction[start_idx + 1 : stop_idx + 1] new_z_direction = np.insert(new_z_direction, 0, z_direction_start, axis = 0) + + new_earth_direction = self._earth_direction[start_idx + 1 : stop_idx + 1] + new_earth_direction = np.insert(new_earth_direction, 0, earth_direction_start, axis = 0) + + # Use linear interpolation to get starting altitude at desired time. + f = interpolate.interp1d(self._time.value, self._altitude, kind="linear") + starting_alt = f(start.value) + new_earth_altitude = self._altitude[start_idx + 1 : stop_idx + 1] + new_earth_altitude = np.insert(new_earth_altitude, 0, starting_alt) if (stop.value % 1 != 0): @@ -245,6 +302,7 @@ def source_interval(self, start, stop): x_direction_stop = self.interpolate_direction(stop, stop_idx, self._x_direction) z_direction_stop = self.interpolate_direction(stop, stop_idx, self._z_direction) + earth_direction_stop = self.interpolate_direction(stop, stop_idx, self._earth_direction) new_times = np.delete(new_times, -1) new_times = np.append(new_times, stop.value) @@ -254,12 +312,23 @@ def source_interval(self, start, stop): new_z_direction = new_z_direction[:-1] new_z_direction = np.append(new_z_direction, [z_direction_stop], axis = 0) + + new_earth_direction = new_earth_direction[:-1] + new_earth_direction = np.append(new_earth_direction, [earth_direction_stop], axis = 0) + + # Use linear interpolation to get starting altitude at desired time. + f = interpolate.interp1d(self._time.value, self._altitude, kind="linear") + stop_alt = f(stop.value) + new_earth_altitude = new_earth_altitude[:-1] + new_earth_altitude = np.append(new_earth_altitude, [stop_alt]) time = Time(new_times, format = "unix") xpointings = SkyCoord(l = new_x_direction[:,0]*u.deg, b = new_x_direction[:,1]*u.deg, frame = "galactic") zpointings = SkyCoord(l = new_z_direction[:,0]*u.deg, b = new_z_direction[:,1]*u.deg, frame = "galactic") + earthpointings = SkyCoord(l = new_earth_direction[:,0]*u.deg, b = new_earth_direction[:,1]*u.deg, frame = "galactic") + altitude = new_earth_altitude - return self.__class__(time, x_pointings = xpointings, z_pointings = zpointings) + return self.__class__(time, x_pointings = xpointings, z_pointings = zpointings, earth_zenith = earthpointings, altitude =altitude) def get_attitude(self, x_pointings = None, y_pointings = None, z_pointings = None): @@ -400,7 +469,7 @@ def get_dwell_map(self, response, src_path = None, save = False): path = src_path # check if the target source path is astropy.Skycoord object if type(path) != SkyCoord: - raise TypeError("The coordiates of the source movement in the Spacecraft frame must be a SkyCoord object") + raise TypeError("The coordinates of the source movement in the Spacecraft frame must be a SkyCoord object") if path.shape[0]-1 != self.dts.shape[0]: raise ValueError("The dimensions of the dts or source coordinates are not correct. Please check your inputs.") @@ -444,22 +513,31 @@ def get_dwell_map(self, response, src_path = None, save = False): return self.dwell_map def get_scatt_map(self, + target_coord, nside, scheme = 'ring', coordsys = 'galactic', + r_earth = 6378.0, ): + """ - Bin the spacecraft attitude history into a 4D histogram that contains the accumulated time the axes of the spacecraft where looking at a given direction. + Bin the spacecraft attitude history into a 4D histogram that + contains the accumulated time the axes of the spacecraft where + looking at a given direction. Parameters ---------- + target_coord : astropy.coordinates.SkyCoord + The coordinates of the target object. nside : int The nside of the scatt map. scheme : str, optional The scheme of the scatt map (the default is "ring") coordsys : str, optional The coordinate system used in the scatt map (the default is "galactic). - + r_earth : float, optional + Earth radius in km (default is 6378 km). + Returns ------- h_ori : cosipy.spacecraftfile.scatt_map.SpacecraftAttitudeMap @@ -470,18 +548,40 @@ def get_scatt_map(self, timestamps = self.get_time() attitudes = self.get_attitude() + # Altitude at each point in the orbit: + altitude = self._altitude + + # Earth zenith at each point in the orbit: + earth_zenith = self.earth_zenith + # Fill (only 2 axes needed to fully define the orientation) h_ori = SpacecraftAttitudeMap(nside = nside, scheme = scheme, coordsys = coordsys) x,y,z = attitudes[:-1].as_axes() + + # Get max angle based on altitude: + max_angle = np.pi - np.arcsin(r_earth/(r_earth + altitude)) + max_angle *= (180/np.pi) # angles in degree + + # Calculate angle between source direction and Earth zenith + # for each time stamp: + src_angle = target_coord.separation(earth_zenith) + + # Get pointings that are occulted by Earth: + earth_occ_index = src_angle.value >= max_angle - h_ori.fill(x, y, weight = np.diff(timestamps.gps)*u.s) + # Define weights and set to 0 if blocked by Earth: + weight = np.diff(timestamps.gps)*u.s + weight[earth_occ_index[:-1]] = 0 + # Fill histogram: + h_ori.fill(x, y, weight = weight) + return h_ori - + def get_psr_rsp(self, response = None, dwell_map = None, dts = None): """ @@ -490,7 +590,7 @@ def get_psr_rsp(self, response = None, dwell_map = None, dts = None): Parameters ---------- - response : str or pathlib.Path, optional + :response : str or pathlib.Path, optional The response for the observation (the defaul is `None`, which implies that the `response` will be read from the instance). dwell_map : str, optional The time dwell map for the source, you can load saved dwell time map using this parameter if you've saved it before (the defaul is `None`, which implies that the `dwell_map` will be read from the instance). diff --git a/cosipy/test_data/20280301_2s.ori b/cosipy/test_data/20280301_2s.ori index 1d95981f..4a9b48c4 100644 --- a/cosipy/test_data/20280301_2s.ori +++ b/cosipy/test_data/20280301_2s.ori @@ -1,5 +1,5 @@ Type OrientationsGalactic -OG 1835487300.0 68.44034002307066 44.61117227945379 -21.559659976929343 44.61117227945379 -OG 1835487301.0 68.38920658776064 44.642124027021296 -21.610793412239364 44.642124027021296 -OG 1835487302.0 68.3380787943012 44.67309722321497 -21.661921205698793 44.67309722321497 -EN \ No newline at end of file +OG 1835487300.0 68.44034002307066 44.61117227945379 -21.559659976929343 44.61117227945379 0.0 0.0 0.0 +OG 1835487301.0 68.38920658776064 44.642124027021296 -21.610793412239364 44.642124027021296 0.0 0.0 0.0 +OG 1835487302.0 68.3380787943012 44.67309722321497 -21.661921205698793 44.67309722321497 0.0 0.0 0.0 +EN diff --git a/cosipy/test_data/20280301_first_10sec.ori b/cosipy/test_data/20280301_first_10sec.ori index 1aef6fc0..89db6f54 100644 --- a/cosipy/test_data/20280301_first_10sec.ori +++ b/cosipy/test_data/20280301_first_10sec.ori @@ -1,13 +1,13 @@ Type OrientationsGalactic -OG 1835478000.0 73.14907746670937 41.85821768724895 16.85092253329064 221.85821768724895 -OG 1835478001.0 73.09517926980278 41.88225011209611 16.904820730197223 221.8822501120961 -OG 1835478002.0 73.04128380352786 41.90629597072256 16.95871619647214 221.90629597072257 -OG 1835478003.0 72.98739108131268 41.93035532675578 17.012608918687327 221.93035532675577 -OG 1835478004.0 72.9335011165853 41.954428243823145 17.066498883414702 221.95442824382317 -OG 1835478005.0 72.87961392277379 41.978514785552235 17.120386077226204 221.97851478555222 -OG 1835478006.0 72.82572951330626 42.002615015570285 17.174270486693747 222.0026150155703 -OG 1835478007.0 72.77184790161073 42.02672899750497 17.228152098389273 222.02672899750493 -OG 1835478008.0 72.7179691011153 42.05085679498347 17.282030898884702 222.05085679498347 -OG 1835478009.0 72.66409312524804 42.07499847163346 17.335906874751963 222.07499847163342 -OG 1835478010.0 72.61021998743702 42.09915409108222 17.38978001256298 222.09915409108223 -# EN +OG 1835478000.0 73.14907746670937 41.85821768724895 16.85092253329064 221.85821768724895 0.0 0.0 0.0 +OG 1835478001.0 73.09517926980278 41.88225011209611 16.904820730197223 221.8822501120961 0.0 0.0 0.0 +OG 1835478002.0 73.04128380352786 41.90629597072256 16.95871619647214 221.90629597072257 0.0 0.0 0.0 +OG 1835478003.0 72.98739108131268 41.93035532675578 17.012608918687327 221.93035532675577 0.0 0.0 0.0 +OG 1835478004.0 72.9335011165853 41.954428243823145 17.066498883414702 221.95442824382317 0.0 0.0 0.0 +OG 1835478005.0 72.87961392277379 41.978514785552235 17.120386077226204 221.97851478555222 0.0 0.0 0.0 +OG 1835478006.0 72.82572951330626 42.002615015570285 17.174270486693747 222.0026150155703 0.0 0.0 0.0 +OG 1835478007.0 72.77184790161073 42.02672899750497 17.228152098389273 222.02672899750493 0.0 0.0 0.0 +OG 1835478008.0 72.7179691011153 42.05085679498347 17.282030898884702 222.05085679498347 0.0 0.0 0.0 +OG 1835478009.0 72.66409312524804 42.07499847163346 17.335906874751963 222.07499847163342 0.0 0.0 0.0 +OG 1835478010.0 72.61021998743702 42.09915409108222 17.38978001256298 222.09915409108223 0.0 0.0 0.0 +EN diff --git a/docs/tutorials/other_examples.rst b/docs/tutorials/other_examples.rst index ac9ddc9e..e31c7c7f 100644 --- a/docs/tutorials/other_examples.rst +++ b/docs/tutorials/other_examples.rst @@ -7,3 +7,5 @@ Other examples .. toctree:: :maxdepth: 1 + Point source response with Earth occultation + diff --git a/docs/tutorials/response/PSR_with_Earth_occultation_example.ipynb b/docs/tutorials/response/PSR_with_Earth_occultation_example.ipynb new file mode 100644 index 00000000..1a5a6e69 --- /dev/null +++ b/docs/tutorials/response/PSR_with_Earth_occultation_example.ipynb @@ -0,0 +1,456 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "30f0be8e-fcf6-48dd-b1ce-471518cbd421", + "metadata": {}, + "source": [ + "# Point Source Response with Earth Occultation\n", + "\n", + "This notebook requires downloading the following files from wasabi:\n", + "\n", + "1) DC3_final_530km_3_month_with_slew_15sbins_GalacticEarth.ori (from DC3)\n", + "\n", + "2) SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.h5 (from DC2)\n", + "\n", + "Additionally, you need test_earth_occ.ori, which is included in the same directory as this notebook. " + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "1b5af751-2243-4366-8be8-14e1091cf71a", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
16:35:19 WARNING   The naima package is not available. Models that depend on it will not be         functions.py:48\n",
+       "                  available                                                                                        \n",
+       "
\n" + ], + "text/plain": [ + "\u001b[38;5;46m16:35:19\u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m The naima package is not available. Models that depend on it will not be \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=781318;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/astromodels/functions/functions_1D/functions.py\u001b\\\u001b[2mfunctions.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=239732;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/astromodels/functions/functions_1D/functions.py#48\u001b\\\u001b[2m48\u001b[0m\u001b]8;;\u001b\\\n", + "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mavailable \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
         WARNING   The GSL library or the pygsl wrapper cannot be loaded. Models that depend on it  functions.py:69\n",
+       "                  will not be available.                                                                           \n",
+       "
\n" + ], + "text/plain": [ + "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m The GSL library or the pygsl wrapper cannot be loaded. Models that depend on it \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=157611;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/astromodels/functions/functions_1D/functions.py\u001b\\\u001b[2mfunctions.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=426238;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/astromodels/functions/functions_1D/functions.py#69\u001b\\\u001b[2m69\u001b[0m\u001b]8;;\u001b\\\n", + "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mwill not be available. \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
16:35:20 WARNING   The ebltable package is not available. Models that depend on it will not be     absorption.py:33\n",
+       "                  available                                                                                        \n",
+       "
\n" + ], + "text/plain": [ + "\u001b[38;5;46m16:35:20\u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m The ebltable package is not available. Models that depend on it will not be \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=421755;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/astromodels/functions/functions_1D/absorption.py\u001b\\\u001b[2mabsorption.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=212297;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/astromodels/functions/functions_1D/absorption.py#33\u001b\\\u001b[2m33\u001b[0m\u001b]8;;\u001b\\\n", + "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mavailable \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
16:35:21 INFO      Starting 3ML!                                                                     __init__.py:39\n",
+       "
\n" + ], + "text/plain": [ + "\u001b[38;5;46m16:35:21\u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;49mINFO \u001b[0m \u001b[1;38;5;251m Starting 3ML! \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=837874;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=231355;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/__init__.py#39\u001b\\\u001b[2m39\u001b[0m\u001b]8;;\u001b\\\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
         WARNING   WARNINGs here are NOT errors                                                      __init__.py:40\n",
+       "
\n" + ], + "text/plain": [ + "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m WARNINGs here are \u001b[0m\u001b[1;31mNOT\u001b[0m\u001b[1;38;5;251m errors \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=586991;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=166378;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/__init__.py#40\u001b\\\u001b[2m40\u001b[0m\u001b]8;;\u001b\\\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
         WARNING   but are inform you about optional packages that can be installed                  __init__.py:41\n",
+       "
\n" + ], + "text/plain": [ + "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m but are inform you about optional packages that can be installed \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=993822;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=55460;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/__init__.py#41\u001b\\\u001b[2m41\u001b[0m\u001b]8;;\u001b\\\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
         WARNING    to disable these messages, turn off start_warning in your config file            __init__.py:44\n",
+       "
\n" + ], + "text/plain": [ + "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m \u001b[0m\u001b[1;31m to disable these messages, turn off start_warning in your config file\u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=306104;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=60446;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/__init__.py#44\u001b\\\u001b[2m44\u001b[0m\u001b]8;;\u001b\\\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
         WARNING   no display variable set. using backend for graphics without display (agg)         __init__.py:50\n",
+       "
\n" + ], + "text/plain": [ + "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m no display variable set. using backend for graphics without display \u001b[0m\u001b[1;38;5;251m(\u001b[0m\u001b[1;38;5;251magg\u001b[0m\u001b[1;38;5;251m)\u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=314483;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=491477;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/__init__.py#50\u001b\\\u001b[2m50\u001b[0m\u001b]8;;\u001b\\\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
16:35:22 WARNING   ROOT minimizer not available                                                minimization.py:1345\n",
+       "
\n" + ], + "text/plain": [ + "\u001b[38;5;46m16:35:22\u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m ROOT minimizer not available \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=528503;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/minimizer/minimization.py\u001b\\\u001b[2mminimization.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=660971;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/minimizer/minimization.py#1345\u001b\\\u001b[2m1345\u001b[0m\u001b]8;;\u001b\\\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
         WARNING   Multinest minimizer not available                                           minimization.py:1357\n",
+       "
\n" + ], + "text/plain": [ + "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m Multinest minimizer not available \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=555187;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/minimizer/minimization.py\u001b\\\u001b[2mminimization.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=914070;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/minimizer/minimization.py#1357\u001b\\\u001b[2m1357\u001b[0m\u001b]8;;\u001b\\\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
         WARNING   PyGMO is not available                                                      minimization.py:1369\n",
+       "
\n" + ], + "text/plain": [ + "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m PyGMO is not available \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=853175;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/minimizer/minimization.py\u001b\\\u001b[2mminimization.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=263289;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/minimizer/minimization.py#1369\u001b\\\u001b[2m1369\u001b[0m\u001b]8;;\u001b\\\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
16:35:23 WARNING   The cthreeML package is not installed. You will not be able to use plugins which  __init__.py:94\n",
+       "                  require the C/C++ interface (currently HAWC)                                                     \n",
+       "
\n" + ], + "text/plain": [ + "\u001b[38;5;46m16:35:23\u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m The cthreeML package is not installed. You will not be able to use plugins which \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=195766;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=807920;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/__init__.py#94\u001b\\\u001b[2m94\u001b[0m\u001b]8;;\u001b\\\n", + "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mrequire the C/C++ interface \u001b[0m\u001b[1;38;5;251m(\u001b[0m\u001b[1;38;5;251mcurrently HAWC\u001b[0m\u001b[1;38;5;251m)\u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
         WARNING   Could not import plugin FermiLATLike.py. Do you have the relative instrument     __init__.py:144\n",
+       "                  software installed and configured?                                                               \n",
+       "
\n" + ], + "text/plain": [ + "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m Could not import plugin FermiLATLike.py. Do you have the relative instrument \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=393727;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=458256;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/__init__.py#144\u001b\\\u001b[2m144\u001b[0m\u001b]8;;\u001b\\\n", + "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251msoftware installed and configured? \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
16:35:24 WARNING   Could not import plugin HAWCLike.py. Do you have the relative instrument         __init__.py:144\n",
+       "                  software installed and configured?                                                               \n",
+       "
\n" + ], + "text/plain": [ + "\u001b[38;5;46m16:35:24\u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m Could not import plugin HAWCLike.py. Do you have the relative instrument \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=745007;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=898798;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/__init__.py#144\u001b\\\u001b[2m144\u001b[0m\u001b]8;;\u001b\\\n", + "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251msoftware installed and configured? \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
16:35:26 WARNING   No fermitools installed                                              lat_transient_builder.py:44\n",
+       "
\n" + ], + "text/plain": [ + "\u001b[38;5;46m16:35:26\u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m No fermitools installed \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=688121;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/utils/data_builders/fermi/lat_transient_builder.py\u001b\\\u001b[2mlat_transient_builder.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=962685;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/utils/data_builders/fermi/lat_transient_builder.py#44\u001b\\\u001b[2m44\u001b[0m\u001b]8;;\u001b\\\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
16:35:26 WARNING   Env. variable OMP_NUM_THREADS is not set. Please set it to 1 for optimal         __init__.py:387\n",
+       "                  performances in 3ML                                                                              \n",
+       "
\n" + ], + "text/plain": [ + "\u001b[38;5;46m16:35:26\u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m Env. variable OMP_NUM_THREADS is not set. Please set it to \u001b[0m\u001b[1;37m1\u001b[0m\u001b[1;38;5;251m for optimal \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=971080;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=965278;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/__init__.py#387\u001b\\\u001b[2m387\u001b[0m\u001b]8;;\u001b\\\n", + "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mperformances in 3ML \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
         WARNING   Env. variable MKL_NUM_THREADS is not set. Please set it to 1 for optimal         __init__.py:387\n",
+       "                  performances in 3ML                                                                              \n",
+       "
\n" + ], + "text/plain": [ + "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m Env. variable MKL_NUM_THREADS is not set. Please set it to \u001b[0m\u001b[1;37m1\u001b[0m\u001b[1;38;5;251m for optimal \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=329486;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=890782;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/__init__.py#387\u001b\\\u001b[2m387\u001b[0m\u001b]8;;\u001b\\\n", + "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mperformances in 3ML \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "
         WARNING   Env. variable NUMEXPR_NUM_THREADS is not set. Please set it to 1 for optimal     __init__.py:387\n",
+       "                  performances in 3ML                                                                              \n",
+       "
\n" + ], + "text/plain": [ + "\u001b[38;5;46m \u001b[0m\u001b[38;5;46m \u001b[0m\u001b[38;5;134mWARNING \u001b[0m \u001b[1;38;5;251m Env. variable NUMEXPR_NUM_THREADS is not set. Please set it to \u001b[0m\u001b[1;37m1\u001b[0m\u001b[1;38;5;251m for optimal \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b]8;id=85779;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/__init__.py\u001b\\\u001b[2m__init__.py\u001b[0m\u001b]8;;\u001b\\\u001b[2m:\u001b[0m\u001b]8;id=583251;file:///discover/nobackup/ckarwin/Software/COSIPY/lib/python3.10/site-packages/threeML/__init__.py#387\u001b\\\u001b[2m387\u001b[0m\u001b]8;;\u001b\\\n", + "\u001b[38;5;46m \u001b[0m \u001b[1;38;5;251mperformances in 3ML \u001b[0m\u001b[1;38;5;251m \u001b[0m\u001b[2m \u001b[0m\n" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Imports\n", + "from cosipy.spacecraftfile import SpacecraftFile\n", + "from cosipy.response import FullDetectorResponse\n", + "from cosipy.util import fetch_wasabi_file\n", + "import astropy.units as u\n", + "from astropy.coordinates import SkyCoord\n", + "import numpy as np" + ] + }, + { + "cell_type": "markdown", + "id": "120eeb60-7f17-4e97-83df-45b2f7aea854", + "metadata": {}, + "source": [ + "Load orientation file:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "064aad90-a5b4-4fb7-9a36-a7bd43f95db7", + "metadata": {}, + "outputs": [], + "source": [ + "ori_file = \"your/path/DC3_final_530km_3_month_with_slew_15sbins_GalacticEarth.ori\"\n", + "ori = SpacecraftFile.parse_from_file(ori_file)" + ] + }, + { + "cell_type": "markdown", + "id": "55314d24-d738-45e6-9f09-4035ad18ea76", + "metadata": {}, + "source": [ + "Define coordinate of source (Crab):" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "59969bd5-5a3d-456f-a534-43a175781c3d", + "metadata": {}, + "outputs": [], + "source": [ + "coord = SkyCoord(l=184.56*u.deg,b=-5.78*u.deg,frame=\"galactic\")" + ] + }, + { + "cell_type": "markdown", + "id": "3cbaa450-2719-4abc-83bf-eedccea347be", + "metadata": {}, + "source": [ + "Calculate scatt map for given coordinate:" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "8e89f3f9-aa47-47c6-acbb-a2fd1537cc53", + "metadata": {}, + "outputs": [], + "source": [ + "scatt_map = ori.get_scatt_map(coord, nside = 16, coordsys = 'galactic')" + ] + }, + { + "cell_type": "markdown", + "id": "e2c4ccf0-7543-4980-8256-3876cc454b5c", + "metadata": {}, + "source": [ + "Calculate psr:" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "fe07ffb2-e7c0-423b-894c-b3c161d3464a", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Works!\n" + ] + } + ], + "source": [ + "response_path = \"your/path/\"\n", + "response_file = response_path + \"SMEXv12.Continuum.HEALPixO3_10bins_log_flat.binnedimaging.imagingresponse.nonsparse_nside8.area.good_chunks_unzip.h5\"\n", + "with FullDetectorResponse.open(response_file) as response:\n", + " psr = response.get_point_source_response(coord = coord, scatt_map = scatt_map)\n", + "print(\"Works!\") " + ] + }, + { + "cell_type": "markdown", + "id": "75e50ee7-8f24-47f8-8c2a-7e99fa722891", + "metadata": {}, + "source": [ + "### Now let's do a simple test:\n", + "Load test orientation file: the z-axis of the instrument points to the Galactic center for 30 seconds. For the first 15 seconds, the instrument is in front of Earth, and then for the last 15.1 seconds the instrument is behind Earth. This should be reflected in the scatt map." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "6d14166f-887b-4097-bf62-cbc84cca0502", + "metadata": {}, + "outputs": [], + "source": [ + "ori_file = \"your/path/test_earth_occ.ori\"\n", + "ori = SpacecraftFile.parse_from_file(ori_file)" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "99efa500-4018-4fdb-9793-363a192cff81", + "metadata": {}, + "outputs": [], + "source": [ + "coord = SkyCoord(l=0*u.deg,b=0*u.deg,frame=\"galactic\")" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "9303a803-bc70-46b6-a0c3-6684e12ce2c5", + "metadata": {}, + "outputs": [], + "source": [ + "scatt_map = ori.get_scatt_map(coord, nside = 16, coordsys = 'galactic')" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "60b0962f-ab6a-4155-be56-fd1e510b5ed0", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "15.099999904632568" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.sum(scatt_map.contents.todense())" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python [conda env:COSIPY]", + "language": "python", + "name": "conda-env-COSIPY-py" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.15" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/tutorials/response/test_earth_occ.ori b/docs/tutorials/response/test_earth_occ.ori new file mode 100644 index 00000000..33ef6995 --- /dev/null +++ b/docs/tutorials/response/test_earth_occ.ori @@ -0,0 +1,6 @@ +Type OrientationsGalactic +OG 1835487300.0 0 90 0 0 530.2220974832549 0 0 +OG 1835487315.0 0 90 0 0 530.2220974832549 0 0 +OG 1835487315.1 0 90 0 0 530.2220974832549 0 180 +OG 1835487330.0 0 90 0 0 530.2220974832549 0 180 +EN