From e47010a19b9793e031bcef8e0c888aac132f8448 Mon Sep 17 00:00:00 2001 From: postoecker Date: Wed, 8 Apr 2026 14:51:59 +0200 Subject: [PATCH 1/5] Change source folder to src/curtains/ --- .gitignore | 6 +++--- pyproject.toml | 2 +- {curtains => src/curtains}/analysis.py | 0 {curtains => src/curtains}/images.py | 2 +- {curtains => src/curtains}/instrument.py | 0 {curtains => src/curtains}/models/arf.fits | Bin {curtains => src/curtains}/models/rmf.fits | Bin {curtains => src/curtains}/process.py | 0 8 files changed, 5 insertions(+), 5 deletions(-) rename {curtains => src/curtains}/analysis.py (100%) rename {curtains => src/curtains}/images.py (99%) rename {curtains => src/curtains}/instrument.py (100%) rename {curtains => src/curtains}/models/arf.fits (100%) rename {curtains => src/curtains}/models/rmf.fits (100%) rename {curtains => src/curtains}/process.py (100%) diff --git a/.gitignore b/.gitignore index 8b3c745..efd947f 100644 --- a/.gitignore +++ b/.gitignore @@ -13,9 +13,9 @@ random/ output_data/ results/ -curtains/models/* -!curtains/models/arf.fits -!curtains/models/rmf.fits +src/curtains/models/* +!src/curtains/models/arf.fits +!src/curtains/models/rmf.fits *.egg-info/ *.pdf diff --git a/pyproject.toml b/pyproject.toml index 0601c76..fb26b3e 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -14,4 +14,4 @@ dependencies = [ ] [tool.setuptools.packages.find] -where = ["curtains"] +where = ["src"] diff --git a/curtains/analysis.py b/src/curtains/analysis.py similarity index 100% rename from curtains/analysis.py rename to src/curtains/analysis.py diff --git a/curtains/images.py b/src/curtains/images.py similarity index 99% rename from curtains/images.py rename to src/curtains/images.py index 706a423..3a02572 100644 --- a/curtains/images.py +++ b/src/curtains/images.py @@ -338,7 +338,7 @@ def generate_from_image(image_path, no_photons, img_scale, energy, energy_spread source_counts = np.random.choice( np.arange(img_array.size), size=no_photons, - p=pmf.flatten() + p=pmf.flatten(), ) # generate bkg photons diff --git a/curtains/instrument.py b/src/curtains/instrument.py similarity index 100% rename from curtains/instrument.py rename to src/curtains/instrument.py diff --git a/curtains/models/arf.fits b/src/curtains/models/arf.fits similarity index 100% rename from curtains/models/arf.fits rename to src/curtains/models/arf.fits diff --git a/curtains/models/rmf.fits b/src/curtains/models/rmf.fits similarity index 100% rename from curtains/models/rmf.fits rename to src/curtains/models/rmf.fits diff --git a/curtains/process.py b/src/curtains/process.py similarity index 100% rename from curtains/process.py rename to src/curtains/process.py From cb606d672cb53d22e198f5af25510c297904545e Mon Sep 17 00:00:00 2001 From: postoecker Date: Wed, 8 Apr 2026 15:22:02 +0200 Subject: [PATCH 2/5] Add dependencies and version constraints --- pyproject.toml | 17 +++++++++++++++-- 1 file changed, 15 insertions(+), 2 deletions(-) diff --git a/pyproject.toml b/pyproject.toml index fb26b3e..c79a1b1 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -8,9 +8,22 @@ name = "curtains" version = "0.0.1" description = "Simulation of X-ray interferometric observations" readme = "README.rst" -requires-python = ">=3.8" +authors = [ + { name="Philipp Stoecker", email="p.o.stoecker@uva.nl" }, +] +requires-python = ">=3.11" dependencies = [ - + "numpy~=2.4", + "astropy~=7.2", + "matplotlib~=3.10.8", + "pandas~=3.0.1", + "scipy~=1.17.1", + "pillow~=12.1.1", +] + +[project.optional-dependencies] +dev = [ + "ruff~=0.13.2", ] [tool.setuptools.packages.find] From 1f88fd2d19a8bec4dacfefa53a270236460b63ff Mon Sep 17 00:00:00 2001 From: postoecker Date: Wed, 8 Apr 2026 15:26:16 +0200 Subject: [PATCH 3/5] Dynamically add version info to the package from pyproject.toml --- .gitignore | 3 ++- pyproject.toml | 24 ++++++++++++++++++------ 2 files changed, 20 insertions(+), 7 deletions(-) diff --git a/.gitignore b/.gitignore index efd947f..e2c317a 100644 --- a/.gitignore +++ b/.gitignore @@ -19,4 +19,5 @@ src/curtains/models/* *.egg-info/ *.pdf -simulation.py \ No newline at end of file +simulation.py +src/curtains/version.py diff --git a/pyproject.toml b/pyproject.toml index c79a1b1..5719785 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,11 +1,6 @@ -[build-system] -requires = ["setuptools>=61.0"] - -build-backend = "setuptools.build_meta" - [project] name = "curtains" -version = "0.0.1" +dynamic = ["version"] description = "Simulation of X-ray interferometric observations" readme = "README.rst" authors = [ @@ -26,5 +21,22 @@ dev = [ "ruff~=0.13.2", ] +[build-system] +requires = ["setuptools>=77.0","setuptools_scm>=9.2"] +build-backend = "setuptools.build_meta" + +[tool.setuptools_scm] +fallback_version = "0.1" +write_to = "src/curtains/version.py" +write_to_template = """ +# coding: utf-8 +# file generated by setuptools_scm +# don't change, don't track in version control +__version__ = version = "{version}" +__version_tuple__ = version_tuple = {version_tuple!r} +__git_hash__ = "{scm_version.node}" +__semantic_version__ = "{scm_version.tag}" +""" + [tool.setuptools.packages.find] where = ["src"] From 4d41e60ef08c54c6a4e63ef084610eb0f06b992c Mon Sep 17 00:00:00 2001 From: postoecker Date: Wed, 8 Apr 2026 17:18:15 +0200 Subject: [PATCH 4/5] Add continuous integration --- .github/workflows/ci.yml | 42 ++++ pyproject.toml | 13 +- src/curtains/analysis.py | 360 ++++++++++++++++++++---------- src/curtains/images.py | 237 +++++++++++--------- src/curtains/instrument.py | 290 +++++++++++++++--------- src/curtains/process.py | 439 ++++++++++++++++++++++--------------- 6 files changed, 872 insertions(+), 509 deletions(-) create mode 100644 .github/workflows/ci.yml diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..d1841b7 --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,42 @@ +# This workflow will install Python dependencies, run tests and lint +# with a variety of Python versions +# For more information see: +# https://docs.github.com/en/actions/automating-builds-and-tests/building-and-testing-python + +name: Curtains + +on: + push: + branches: + - main + pull_request: + branches: + - main + workflow_dispatch: + + +jobs: + build: + + runs-on: ubuntu-latest + strategy: + fail-fast: false + matrix: + python-version: ["3.11", "3.12", "3.13", "3.14"] + + steps: + - uses: actions/checkout@v4 + - name: Set up Python ${{ matrix.python-version }} + uses: actions/setup-python@v3 + with: + python-version: ${{ matrix.python-version }} + - name: Install dependencies + run: | + python -m pip install --upgrade pip + python -m pip install '.[dev]' + - name: Check formatting with ruff + run: | + ruff format --check + - name: Lint with ruff + run: | + ruff check diff --git a/pyproject.toml b/pyproject.toml index 5719785..1784356 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -10,10 +10,10 @@ requires-python = ">=3.11" dependencies = [ "numpy~=2.4", "astropy~=7.2", - "matplotlib~=3.10.8", - "pandas~=3.0.1", - "scipy~=1.17.1", - "pillow~=12.1.1", + "matplotlib~=3.10", + "pandas~=3.0", + "scipy~=1.17", + "pillow~=12.1", ] [project.optional-dependencies] @@ -40,3 +40,8 @@ __semantic_version__ = "{scm_version.tag}" [tool.setuptools.packages.find] where = ["src"] + +[tool.ruff] +# Ignore the auto-generated file +include = ["src/curtains/*.py"] +exclude = ["src/curtains/version.py"] diff --git a/src/curtains/analysis.py b/src/curtains/analysis.py index 7977af4..af23872 100644 --- a/src/curtains/analysis.py +++ b/src/curtains/analysis.py @@ -1,4 +1,3 @@ - import numpy as np import matplotlib.pyplot as plt import scipy.fft as ft @@ -6,7 +5,8 @@ from astropy.io import fits from matplotlib.colors import LinearSegmentedColormap -def hist_data(data, binsno, pixs = False, num = 0): + +def hist_data(data, binsno, pixs=False, num=0): """ Function that makes a histogram of direct output data from an interferometer object. @@ -17,12 +17,13 @@ def hist_data(data, binsno, pixs = False, num = 0): """ if pixs: - plt.hist(data, binsno, label=f'Baseline {num}') - plt.xlabel('Detector position (pixels)') + plt.hist(data, binsno, label=f"Baseline {num}") + plt.xlabel("Detector position (pixels)") else: - plt.hist(data, binsno, label=f'Baseline {num}') - plt.xlabel('Detector position (micrometers)') - plt.ylabel('Counts') + plt.hist(data, binsno, label=f"Baseline {num}") + plt.xlabel("Detector position (micrometers)") + plt.ylabel("Counts") + def ft_data(y_data, samples, spacing): """ @@ -38,7 +39,8 @@ def ft_data(y_data, samples, spacing): return ft_x_data, ft_y_data -def plot_ft(ft_x_data, ft_y_data, plot_obj, log=0, num= 0): + +def plot_ft(ft_x_data, ft_y_data, plot_obj, log=0, num=0): """ Function to plot fourier transformed interferometer data in a number of ways. @@ -48,16 +50,32 @@ def plot_ft(ft_x_data, ft_y_data, plot_obj, log=0, num= 0): log (int in [0,2]): indicates how many axes are to be in log scale, with 1 having only the y-axis in log. """ if log == 0: - plot_obj.plot(ft.fftshift(ft_x_data), (ft.fftshift(ft_y_data)), label=f'Baseline {num}') + plot_obj.plot( + ft.fftshift(ft_x_data), (ft.fftshift(ft_y_data)), label=f"Baseline {num}" + ) if log == 1: - plot_obj.semilogy(ft.fftshift(ft_x_data), (ft.fftshift(ft_y_data)), label=f'Baseline {num}') + plot_obj.semilogy( + ft.fftshift(ft_x_data), (ft.fftshift(ft_y_data)), label=f"Baseline {num}" + ) if log == 2: - plot_obj.loglog(ft.fftshift(ft_x_data), (ft.fftshift(ft_y_data)), label=f'Baseline {num}') - -def image_recon_smooth(data, instrument, fov, samples = np.array([512, 512]), progress = 0, error = 0.1, recon_type = "IFFT", verbose = True): + plot_obj.loglog( + ft.fftshift(ft_x_data), (ft.fftshift(ft_y_data)), label=f"Baseline {num}" + ) + + +def image_recon_smooth( + data, + instrument, + fov, + samples=np.array([512, 512]), + progress=0, + error=0.1, + recon_type="IFFT", + verbose=True, +): """ This function is to be used to reconstruct images from interferometer data. - Bins input data based on roll angle, which is important to fill out the uv-plane that will be fed into + Bins input data based on roll angle, which is important to fill out the uv-plane that will be fed into the 2d inverse fourier transform. Args: @@ -66,17 +84,17 @@ def image_recon_smooth(data, instrument, fov, samples = np.array([512, 512]), pr samples (int): N for the NxN matrix that is the uv-plane used for the 2d inverse fourier transform. TODO: outdated progress (1 or 0): whether, or not respectively, to show progess of image reconstruction calculations. TODO: add fov argument, unknown at this time; largest axis size in arcsec - + Returns: - array, array, array: Three arrays, first of which is the reconstructed image, + array, array, array: Three arrays, first of which is the reconstructed image, second and third of which are the fourier transforms and associated uv coordinates of each roll bin + energy channel combination. - """ - - def inverse_fourier(f_values, uv, fov, error, instrument, energies, recon_type, verbose): - + """ + def inverse_fourier( + f_values, uv, fov, error, instrument, energies, recon_type, verbose + ): """ - This is a helper function that calculates the inverse fourier transform of the data from all baselines, only to + This is a helper function that calculates the inverse fourier transform of the data from all baselines, only to be used at the last step of the parent function. It is sectioned off here for legibility. It first defines an image to fill in according to the provided samples size, and calculates the sum @@ -85,12 +103,29 @@ def inverse_fourier(f_values, uv, fov, error, instrument, energies, recon_type, if recon_type == "IFFT": wavelengths = spc.h * spc.c / energies - - d = (fov * np.pi / (180 * 3600)) / np.max(samples) #1 / (2 * (np.max(np.array([baseline.D for baseline in instrument.baselines])) / np.min(wavelengths))) - n = np.ceil(1 / (2 * d * (np.min(np.array([baseline.D for baseline in instrument.baselines])) / np.max(wavelengths)) * error)).astype(int) # + 1 + + d = ( + (fov * np.pi / (180 * 3600)) / np.max(samples) + ) # 1 / (2 * (np.max(np.array([baseline.D for baseline in instrument.baselines])) / np.min(wavelengths))) + n = np.ceil( + 1 + / ( + 2 + * d + * ( + np.min( + np.array([baseline.D for baseline in instrument.baselines]) + ) + / np.max(wavelengths) + ) + * error + ) + ).astype(int) # + 1 if verbose: - print(f"Calculated image parameters:\n sample spacing: {d} \n window size: {n} ") + print( + f"Calculated image parameters:\n sample spacing: {d} \n window size: {n} " + ) # print(d) # print(1 / (2 * (np.max(np.array([baseline.D for baseline in instrument.baselines])) / np.min(wavelengths)))) @@ -98,15 +133,17 @@ def inverse_fourier(f_values, uv, fov, error, instrument, energies, recon_type, # calculate the frequency grid onto which the real sampled frequencies are mapped fft_freqs = np.fft.fftfreq(n, d) if verbose: - print(f"Fourier frequencies calculated.") + print("Fourier frequencies calculated.") # find closest frequency grid point in numpy.fft.fftfreq of the real sampled frequencies, # dividing by frequency bin width: Delta f = 1 / (n * d) - u_conv = np.around(uv[:,0] * (n * d)).astype(int) - v_conv = np.around(uv[:,1] * (n * d)).astype(int) + u_conv = np.around(uv[:, 0] * (n * d)).astype(int) + v_conv = np.around(uv[:, 1] * (n * d)).astype(int) uv_conv = np.array([u_conv, v_conv]).T # create nxn matrix, fourier space image - image = np.zeros((fft_freqs.size, fft_freqs.size), dtype = np.complex64)#complex)#np.zeros((fft_freqs_u.size, fft_freqs_v.size), dtype = complex)#np.zeros(samples, dtype = complex) + image = np.zeros( + (fft_freqs.size, fft_freqs.size), dtype=np.complex64 + ) # complex)#np.zeros((fft_freqs_u.size, fft_freqs_v.size), dtype = complex)#np.zeros(samples, dtype = complex) # add all fourier values to the corresponding frequency coordinate as per numpy.fft.fftfreq np.add.at(image, (u_conv, v_conv), f_values) @@ -118,15 +155,27 @@ def inverse_fourier(f_values, uv, fov, error, instrument, energies, recon_type, fft_image = np.real(fft_image) if verbose: print("Inverse FFT complete. Shifting image...") - fft_image = np.roll(fft_image, int(np.around(fft_image.shape[0]/2)), axis=0) - fft_image = np.roll(fft_image, int(np.around(fft_image.shape[1]/2)), axis=1) + fft_image = np.roll( + fft_image, int(np.around(fft_image.shape[0] / 2)), axis=0 + ) + fft_image = np.roll( + fft_image, int(np.around(fft_image.shape[1] / 2)), axis=1 + ) x_npix_real_image = samples[0] - x_lower_bound = np.round((fft_image.shape[0] / 2) - (0.5 * x_npix_real_image)).astype(int) - x_upper_bound = np.round((fft_image.shape[0] / 2) + (0.5 * x_npix_real_image)).astype(int) + x_lower_bound = np.round( + (fft_image.shape[0] / 2) - (0.5 * x_npix_real_image) + ).astype(int) + x_upper_bound = np.round( + (fft_image.shape[0] / 2) + (0.5 * x_npix_real_image) + ).astype(int) y_npix_real_image = samples[1] - y_lower_bound = np.round((fft_image.shape[0] / 2) - (0.5 * y_npix_real_image)).astype(int) - y_upper_bound = np.round((fft_image.shape[0] / 2) + (0.5 * y_npix_real_image)).astype(int) + y_lower_bound = np.round( + (fft_image.shape[0] / 2) - (0.5 * y_npix_real_image) + ).astype(int) + y_upper_bound = np.round( + (fft_image.shape[0] / 2) + (0.5 * y_npix_real_image) + ).astype(int) # prevent indexing problems if x_lower_bound < 0: @@ -141,27 +190,38 @@ def inverse_fourier(f_values, uv, fov, error, instrument, energies, recon_type, if verbose: print("Image shifted. Returning image...") - return (fft_image[np.ix_(np.arange(x_lower_bound, x_upper_bound, 1), np.arange(y_lower_bound, y_upper_bound, 1))], (samples.max() * (1e6 * 3600 * 360 / (2 * np.pi)) * (d * n))), image, uv_conv, fft_freqs + return ( + ( + fft_image[ + np.ix_( + np.arange(x_lower_bound, x_upper_bound, 1), + np.arange(y_lower_bound, y_upper_bound, 1), + ) + ], + (samples.max() * (1e6 * 3600 * 360 / (2 * np.pi)) * (d * n)), + ), + image, + uv_conv, + fft_freqs, + ) # else: - + # def inverse_fourier_val(x, y, v, u, fourier): # # This function is the formula for an inverse fourier transform, without the integration. - # # It is included here to make clear that a discrete inverse fourier transform is what is happening, and + # # It is included here to make clear that a discrete inverse fourier transform is what is happening, and # # to make clear what argument means what and for multi-threading. - + # global re_im # # frequency shift + 1 / (d * n) # re_im += fourier * np.exp(2j * np.pi * ((u) * x + (v) * y)) - # global re_im # shape = np.array(samples, dtype = int) # re_im = np.zeros(shape, dtype = complex) - # fov_x = fov # fov_y = fov_x*(shape[0]/shape[1]) # Assumes square pixels @@ -184,28 +244,26 @@ def inverse_fourier(f_values, uv, fov, error, instrument, energies, recon_type, # # waiting for all the threads to finish # for thread in thread_list: - + # thread.join() # return np.real(re_im), f_values, uv, uv - - # These arrays are all copied locally to reduce the amount of cross-referencing to other objects required. + # These arrays are all copied locally to reduce the amount of cross-referencing to other objects required. time_data = data.toa E_data = data.energies base_ind = data.baseline_indices pointing = data.pointing positional_data = data.pos - + # Generating the arrays that will contain the uv coordinates and associated fourier values covered by the interferometer. uv = [] f_values = np.array([]) - + # Looking only at the baselines that have associated photons for k in np.unique(base_ind): - if verbose: - print(f"Processing baseline {k+1} / {len(instrument.baselines)}") + print(f"Processing baseline {k + 1} / {len(instrument.baselines)}") # Taking only relevant photons from the current baseline in_baseline = base_ind == k @@ -224,24 +282,34 @@ def inverse_fourier(f_values, uv, fov, error, instrument, energies, recon_type, data_bin_roll = positional_data[in_baseline] if verbose: - print(f"Number of photons in baseline: {data_bin_roll.size}\nCalculating uv coordinates and fourier values...") + print( + f"Number of photons in baseline: {data_bin_roll.size}\nCalculating uv coordinates and fourier values..." + ) # Calculating u and v for middle of current bin by taking a projection of the current frequency - u = (freq_baseline + 0) * np.sin(pointing[time_data, 2][in_baseline] % (2 * np.pi)) - v = (freq_baseline + 0) * np.cos(pointing[time_data, 2][in_baseline] % (2 * np.pi)) + u = (freq_baseline + 0) * np.sin( + pointing[time_data, 2][in_baseline] % (2 * np.pi) + ) + v = (freq_baseline + 0) * np.cos( + pointing[time_data, 2][in_baseline] % (2 * np.pi) + ) new_uv_pairs = np.array([u, v]).T # Doing the same with the negative frequency - uv.extend(np.column_stack((new_uv_pairs, -new_uv_pairs)).reshape(-1, new_uv_pairs.shape[1])) + uv.extend( + np.column_stack((new_uv_pairs, -new_uv_pairs)).reshape( + -1, new_uv_pairs.shape[1] + ) + ) # Calculating value of the fourier transform for the current frequency and bin f_value = np.exp(-2j * np.pi * fourier_freq * data_bin_roll) # Doing the same with the negative frequency f_values = np.append(f_values, np.ravel([f_value, np.conjugate(f_value)], "F")) - + if verbose: - print(f"Baseline {k+1} processed.\n") + print(f"Baseline {k + 1} processed.\n") # reshaping the uv coordinate matrix uv = np.array(uv).reshape(-1, 2) @@ -250,25 +318,27 @@ def inverse_fourier(f_values, uv, fov, error, instrument, energies, recon_type, print("All baselines processed. Starting IFFT...") # create image from IFT - img, ft_img, uv_conv, fft_freqs = inverse_fourier(f_values, uv, fov, error, instrument, E_data, recon_type, verbose=verbose) + img, ft_img, uv_conv, fft_freqs = inverse_fourier( + f_values, uv, fov, error, instrument, E_data, recon_type, verbose=verbose + ) if verbose: print("Image reconstruction complete.") # store results in a reconstruction_data object - data_obj = reconstruction_data(img, ft_img, f_values, uv, uv_conv, fft_freqs if recon_type == "IFFT" else None) + data_obj = reconstruction_data( + img, ft_img, f_values, uv, uv_conv, fft_freqs if recon_type == "IFFT" else None + ) return img, uv, data_obj - class reconstruction_data: """ This class is to be used to store the data that is used for the image reconstruction. It contains the image, the fourier transform of the image, and the uv coordinates of the fourier transform. """ - def __init__(self, image, ft_image, f_values, uv, uv_conv, fft_freqs = None): - + def __init__(self, image, ft_image, f_values, uv, uv_conv, fft_freqs=None): self.image = image self.ft_image = ft_image self.f_values = f_values @@ -276,134 +346,182 @@ def __init__(self, image, ft_image, f_values, uv, uv_conv, fft_freqs = None): self.uv_conv = uv_conv self.fft_freqs = fft_freqs -def plot_spec_per_baseline(event_list, exposure_time): +def plot_spec_per_baseline(event_list, exposure_time): hdul = fits.open(event_list) event_data = hdul[1].data - num_baselines = np.max(event_data['BASELINE_ID']) + 1 + num_baselines = np.max(event_data["BASELINE_ID"]) + 1 hdul.close() - plt.figure(figsize=(7,5)) + plt.figure(figsize=(7, 5)) - min_e = np.min(event_data['DETECTED_ENERGY']) - max_e = np.max(event_data['DETECTED_ENERGY']) - bins = np.linspace(min_e, max_e, 101) # 100 bins + min_e = np.min(event_data["DETECTED_ENERGY"]) + max_e = np.max(event_data["DETECTED_ENERGY"]) + bins = np.linspace(min_e, max_e, 101) # 100 bins bin_width = bins[1] - bins[0] - weights = np.ones_like(event_data['DETECTED_ENERGY']) / (exposure_time * bin_width) + weights = np.ones_like(event_data["DETECTED_ENERGY"]) / (exposure_time * bin_width) - plt.hist(event_data['DETECTED_ENERGY'], bins=bins, weights=weights, histtype='step', log=True, label='All Baselines', color='darkblue', lw=2) + plt.hist( + event_data["DETECTED_ENERGY"], + bins=bins, + weights=weights, + histtype="step", + log=True, + label="All Baselines", + color="darkblue", + lw=2, + ) for baseline_id in range(num_baselines): - mask = event_data['BASELINE_ID'] == baseline_id - subset_weights = np.ones_like(event_data['DETECTED_ENERGY'][mask]) / (exposure_time * bin_width) - plt.hist(event_data['DETECTED_ENERGY'][mask], bins=bins, weights=subset_weights, histtype='step', label=f'Baseline {baseline_id + 1}', color='darkblue', alpha=1/(baseline_id+1), lw=1) + mask = event_data["BASELINE_ID"] == baseline_id + subset_weights = np.ones_like(event_data["DETECTED_ENERGY"][mask]) / ( + exposure_time * bin_width + ) + plt.hist( + event_data["DETECTED_ENERGY"][mask], + bins=bins, + weights=subset_weights, + histtype="step", + label=f"Baseline {baseline_id + 1}", + color="darkblue", + alpha=1 / (baseline_id + 1), + lw=1, + ) plt.xlim(min_e, max_e) - plt.yscale('log') - plt.xscale('log') - plt.xlabel('Energy (keV)') - plt.ylabel('cts/s/keV') + plt.yscale("log") + plt.xscale("log") + plt.xlabel("Energy (keV)") + plt.ylabel("cts/s/keV") plt.legend(fontsize=12, loc="best") plt.grid("both") plt.show() + colors_chandra = [ "#000000", # very dark "#242037", # deep navy "#344871", # dark blue "#89A4E3", # mid-to-light blue - "#E2F6FE" # bright white-blue + "#E2F6FE", # bright white-blue ] -xri_cmap = LinearSegmentedColormap.from_list('black_to_blue', colors_chandra, N=1e4) +xri_cmap = LinearSegmentedColormap.from_list("black_to_blue", colors_chandra, N=1e4) -def plot_inf_image(image, fov_arcsec, return_fig = False): - + +def plot_inf_image(image, fov_arcsec, return_fig=False): fig = plt.figure() image = image[0] fov_uas = fov_arcsec * 1e6 - extent = [-fov_uas/2, fov_uas/2, -fov_uas/2, fov_uas/2] - plt.imshow(image/np.max(np.abs((image))), extent=extent, cmap=xri_cmap, origin='lower', vmin=0, vmax=1) - plt.xlabel('$x$ ($\\mu$as)') - plt.ylabel('$y$ ($\\mu$as)') - plt.colorbar(label='Intensity (normalized)') - plt.title('Reconstructed image') + extent = [-fov_uas / 2, fov_uas / 2, -fov_uas / 2, fov_uas / 2] + plt.imshow( + image / np.max(np.abs((image))), + extent=extent, + cmap=xri_cmap, + origin="lower", + vmin=0, + vmax=1, + ) + plt.xlabel("$x$ ($\\mu$as)") + plt.ylabel("$y$ ($\\mu$as)") + plt.colorbar(label="Intensity (normalized)") + plt.title("Reconstructed image") plt.show() if return_fig: return fig - -def plot_uv(uv, return_fig = False): + +def plot_uv(uv, return_fig=False): uv_plot = uv / 1e9 # convert to Giga-lambda # plot uv-plane sampling - fig = plt.figure(figsize=(6,6)) - plt.plot(uv_plot[:, 0], uv_plot[:, 1], 'ro', markersize=1) + fig = plt.figure(figsize=(6, 6)) + plt.plot(uv_plot[:, 0], uv_plot[:, 1], "ro", markersize=1) # plt.xlim(-np.max(uv_plot) * 1.2, np.max(uv_plot) * 1.2) # plt.ylim(-np.max(uv_plot) * 1.2, np.max(uv_plot) * 1.2) - plt.xlim(-6,6) - plt.ylim(-6,6) + plt.xlim(-6, 6) + plt.ylim(-6, 6) plt.title("$uv$-plane sampling") - plt.xlabel('$u$ (G$\\lambda$)') - plt.ylabel('$v$ (G$\\lambda$)') + plt.xlabel("$u$ (G$\\lambda$)") + plt.ylabel("$v$ (G$\\lambda$)") plt.show() if return_fig: return fig -def plot_res_hist(uv, return_fig = False): - uv_abs = np.sqrt(uv[:, 0]**2 + uv[:, 1]**2) +def plot_res_hist(uv, return_fig=False): + uv_abs = np.sqrt(uv[:, 0] ** 2 + uv[:, 1] ** 2) interf_res = 1 / (2 * uv_abs) # in microarcseconds - interf_res_uas = interf_res * 1e6 * 3600 * 180 / np.pi # convert to microarcseconds - bins = np.logspace(np.log10(np.min(interf_res_uas)), np.log10(np.max(interf_res_uas)), 200) - fig = plt.figure(figsize=(7,5)) - plt.hist(interf_res_uas, bins=bins, color='darkred', alpha=0.7, histtype='stepfilled') - plt.xlabel('Interferometric resolution ($\\mu$as)') - plt.ylabel('$\\#$Photons') - plt.title('Distribution of interferometric resolution') - plt.xscale('log') - plt.yscale('log') + interf_res_uas = interf_res * 1e6 * 3600 * 180 / np.pi # convert to microarcseconds + bins = np.logspace( + np.log10(np.min(interf_res_uas)), np.log10(np.max(interf_res_uas)), 200 + ) + fig = plt.figure(figsize=(7, 5)) + plt.hist( + interf_res_uas, bins=bins, color="darkred", alpha=0.7, histtype="stepfilled" + ) + plt.xlabel("Interferometric resolution ($\\mu$as)") + plt.ylabel("$\\#$Photons") + plt.title("Distribution of interferometric resolution") + plt.xscale("log") + plt.yscale("log") plt.grid(True) plt.show() if return_fig: return fig -def subplots(image, fov_arcsec, uv, return_fig = False): - fig, axs = plt.subplots(1, 2, figsize=(14, 6), gridspec_kw={'width_ratios': [1.2, 1]}) + +def subplots(image, fov_arcsec, uv, return_fig=False): + fig, axs = plt.subplots( + 1, 2, figsize=(14, 6), gridspec_kw={"width_ratios": [1.2, 1]} + ) # Plot the reconstructed image img_array = image[0] fov_uas = fov_arcsec * 1e6 - extent = [-fov_uas/2, fov_uas/2, -fov_uas/2, fov_uas/2] - im = axs[0].imshow(img_array/np.max(np.abs(img_array)), extent=extent, cmap=xri_cmap, origin='lower', vmin=0, vmax=1) - axs[0].set_xlabel('$x$ ($\\mu$as)') - axs[0].set_ylabel('$y$ ($\\mu$as)') - axs[0].set_title('Reconstructed image') - fig.colorbar(im, ax=axs[0], label='Intensity (normalized)', fraction=0.046, pad=0.04) + extent = [-fov_uas / 2, fov_uas / 2, -fov_uas / 2, fov_uas / 2] + im = axs[0].imshow( + img_array / np.max(np.abs(img_array)), + extent=extent, + cmap=xri_cmap, + origin="lower", + vmin=0, + vmax=1, + ) + axs[0].set_xlabel("$x$ ($\\mu$as)") + axs[0].set_ylabel("$y$ ($\\mu$as)") + axs[0].set_title("Reconstructed image") + fig.colorbar( + im, ax=axs[0], label="Intensity (normalized)", fraction=0.046, pad=0.04 + ) # Plot the uv-plane sampling uv_plot = uv / 1e9 # convert to Giga-lambda - axs[1].plot(uv_plot[:, 0], uv_plot[:, 1], 'ro', markersize=1, alpha=0.5) + axs[1].plot(uv_plot[:, 0], uv_plot[:, 1], "ro", markersize=1, alpha=0.5) axs[1].set_xlim(-6, 6) axs[1].set_ylim(-6, 6) - axs[1].set_aspect('equal') + axs[1].set_aspect("equal") axs[1].set_title("$uv$-plane sampling") - axs[1].set_xlabel('$u$ (G$\\lambda$)') - axs[1].set_ylabel('$v$ (G$\\lambda$)') + axs[1].set_xlabel("$u$ (G$\\lambda$)") + axs[1].set_ylabel("$v$ (G$\\lambda$)") plt.tight_layout() plt.show() - + if return_fig: return fig def filter_data(data, E_bounds, baseline_bounds=None): - if baseline_bounds is not None: - baseline_mask = (data.baseline_indices >= baseline_bounds[0]) & (data.baseline_indices <= baseline_bounds[1]) - else: + baseline_mask = (data.baseline_indices >= baseline_bounds[0]) & ( + data.baseline_indices <= baseline_bounds[1] + ) + else: baseline_mask = np.ones_like(data.baseline_indices, dtype=bool) - mask = ((data.energies / 1e3 / spc.e > E_bounds[0]) & (data.energies / 1e3 / spc.e < E_bounds[1])) & baseline_mask + mask = ( + (data.energies / 1e3 / spc.e > E_bounds[0]) + & (data.energies / 1e3 / spc.e < E_bounds[1]) + ) & baseline_mask data_temp = data.copy() data_temp.pos = data.pos[mask] data_temp.energies = data.energies[mask] @@ -411,4 +529,4 @@ def filter_data(data, E_bounds, baseline_bounds=None): data_temp.size = len(data_temp.pos) data_temp.toa = data.toa[mask] print(f"Events in filtered data: {len(data_temp.pos)}") - return data_temp \ No newline at end of file + return data_temp diff --git a/src/curtains/images.py b/src/curtains/images.py index 3a02572..d33fe1f 100644 --- a/src/curtains/images.py +++ b/src/curtains/images.py @@ -1,30 +1,23 @@ - from PIL import Image -import gc import numpy as np import scipy.constants as spc import pandas as pd -from astropy.io import fits -from astropy.wcs import WCS - -# for testing -import matplotlib.pyplot as plt -class image(): - """ - Class that defines a data format that all images to be processed by this package should follow. +class image: + """ + Class that defines a data format that all images to be processed by this package should follow. It consists of a number of arrrays of specified size which should contain the energy, arrival time, and ... #TODO define further Note that this class only generates an empty image class of specified size. - Generating the actual photons to fill it up should happen in a seperate function that manipulates an image class object. + Generating the actual photons to fill it up should happen in a seperate function that manipulates an image class object. """ def __init__(self, size): - """ Initiation function for the class. Generates arrays of the specified size for each parameter specified in the class docstring. + """Initiation function for the class. Generates arrays of the specified size for each parameter specified in the class docstring. - Parameters: + Parameters: - size (int) = number of photons to save in arrays.\n + size (int) = number of photons to save in arrays.\n """ # Abbreviation of 'Times Of Arrival'. @@ -47,7 +40,8 @@ def __init__(self, size): self.bkg_indices = [] self.bkg_spect = None -def point_source(size, alpha, beta, energy, spectrum = None): + +def point_source(size, alpha, beta, energy, spectrum=None): """ Function that generates an image of a monochromatic point source according to some specifications. @@ -64,27 +58,27 @@ def point_source(size, alpha, beta, energy, spectrum = None): # save spectrum when given if spectrum is not None: - # load file containing spectrum or copy the spectrum if isinstance(spectrum, str): im.spectrum = np.loadtxt(spectrum) - else: + else: im.spectrum = spectrum if energy is None: - energy = np.mean(im.spectrum[:,0]) + energy = np.mean(im.spectrum[:, 0]) elif energy is None: # define either spectrum or energy raise Exception("ERROR: define either spectrum or energy") - + im.energies[:] = energy * spc.eV * 1e3 im.loc[:] = np.array([alpha, beta]) * 2 * np.pi / (3600 * 360) im.toa = np.array([i for i in range(size)]) return im -def double_point_source(size, alpha, beta, energy, spectrum = None): + +def double_point_source(size, alpha, beta, energy, spectrum=None): """ Function that generates an image of two monochromatic point sources according to some specifications. @@ -101,29 +95,29 @@ def double_point_source(size, alpha, beta, energy, spectrum = None): # save spectrum when given if spectrum is not None: - # load file containing spectrum or copy the spectrum if isinstance(spectrum, str): im.spectrum = np.loadtxt(spectrum) - else: + else: im.spectrum = spectrum if energy is None: - energy = np.mean(im.spectrum[:,0]) + energy = np.mean(im.spectrum[:, 0]) elif energy is None: # define either spectrum or energy - raise Exception("ERROR: define either spectrum or energy") - + raise ValueError("ERROR: define either spectrum or energy") + for i in range(0, size): - source = np.random.randint(0,2) + source = np.random.randint(0, 2) im.energies[i] = energy[source] * spc.eV * 1e3 im.loc[i] = np.array([alpha[source], beta[source]]) * 2 * np.pi / (3600 * 360) im.toa[i] = i return im -def m_point_sources(size, m, alpha, beta, energy, spectrum = None): + +def m_point_sources(size, m, alpha, beta, energy, spectrum=None): """ Function that generates an image of two monochromatic point sources according to some specifications. @@ -141,29 +135,29 @@ def m_point_sources(size, m, alpha, beta, energy, spectrum = None): # save spectrum when given if spectrum is not None: - # load file containing spectrum or copy the spectrum if isinstance(spectrum, str): im.spectrum = np.loadtxt(spectrum) - else: + else: im.spectrum = spectrum if energy is None: - energy = np.mean(im.spectrum[:,0]) + energy = np.mean(im.spectrum[:, 0]) elif energy is None: # define either spectrum or energy - raise Exception("ERROR: define either spectrum or energy") - + raise ValueError("ERROR: define either spectrum or energy") + for i in range(0, size): - source = np.random.randint(0,m) + source = np.random.randint(0, m) im.energies[i] = energy[source] * spc.eV * 1e3 im.loc[i] = np.array([alpha[source], beta[source]]) * 2 * np.pi / (3600 * 360) im.toa[i] = i return im -def point_source_multichromatic_range(size, alpha, beta, energy, spectrum = None): + +def point_source_multichromatic_range(size, alpha, beta, energy, spectrum=None): """ Function that generates an image of a multichromatic point source according to some specifications. @@ -180,27 +174,29 @@ def point_source_multichromatic_range(size, alpha, beta, energy, spectrum = None # save spectrum when given if spectrum is not None: - # load file containing spectrum or copy the spectrum if isinstance(spectrum, str): im.spectrum = np.loadtxt(spectrum) - else: + else: im.spectrum = spectrum if energy is None: - energy = np.mean(im.spectrum[:,0]) + energy = np.mean(im.spectrum[:, 0]) elif energy is None: # define either spectrum or energy - raise Exception("ERROR: define either spectrum or energy") - + raise ValueError("ERROR: define either spectrum or energy") + for i in range(0, size): - im.energies[i] = (np.random.random() * (energy[1] - energy[0]) + energy[0]) * spc.eV * 1e3 + im.energies[i] = ( + (np.random.random() * (energy[1] - energy[0]) + energy[0]) * spc.eV * 1e3 + ) im.loc[i] = np.array([alpha, beta]) * 2 * np.pi / (3600 * 360) im.toa[i] = i return im + def point_source_multichromatic_gauss(size, alpha, beta, energy, energy_spread): """ Function that generates an image of a multichromatic point source according to some specifications. @@ -215,8 +211,6 @@ def point_source_multichromatic_gauss(size, alpha, beta, energy, energy_spread): """ im = image(size) - spectrum = None - im.energies = np.random.normal(energy, energy_spread, size) * spc.eV * 1e3 for i in range(0, size): im.loc[i] = np.array([alpha, beta]) * 2 * np.pi / (3600 * 360) @@ -224,7 +218,8 @@ def point_source_multichromatic_gauss(size, alpha, beta, energy, energy_spread): return im -def disc(size, alpha, beta, energy, radius, energy_spread=0., spectrum = None): + +def disc(size, alpha, beta, energy, radius, energy_spread=0.0, spectrum=None): """ A function that generates photons in the shape of a continuous disk. @@ -243,35 +238,51 @@ def disc(size, alpha, beta, energy, radius, energy_spread=0., spectrum = None): # save spectrum when given if spectrum is not None: - # load file containing spectrum or copy the spectrum if isinstance(spectrum, str): im.spectrum = np.loadtxt(spectrum) - else: + else: im.spectrum = spectrum if energy is None: - energy = np.mean(im.spectrum[:,0]) + energy = np.mean(im.spectrum[:, 0]) elif energy is None: # define either spectrum or energy - raise Exception("ERROR: define either spectrum or energy") - + raise ValueError("ERROR: define either spectrum or energy") + for i in range(0, size): im.energies[i] = energy * spc.eV * 1e3 im.toa[i] = i r = np.random.random() * radius theta = np.random.random() * 2 * np.pi - im.loc[i] = np.array([alpha + r * np.cos(theta), beta + r * np.sin(theta)]) * 2 * np.pi / (3600 * 360) - - if energy_spread > 0.: + im.loc[i] = ( + np.array([alpha + r * np.cos(theta), beta + r * np.sin(theta)]) + * 2 + * np.pi + / (3600 * 360) + ) + + if energy_spread > 0.0: im.energies += np.random.normal(0, energy_spread, size) return im -def generate_from_image(image_path, no_photons, img_scale, energy, energy_spread=0., offset=[0,0], spectrum = None, bkg_phot = None, bkg_spect = None, bkg_energy = None): + +def generate_from_image( + image_path, + no_photons, + img_scale, + energy, + energy_spread=0.0, + offset=None, + spectrum=None, + bkg_phot=None, + bkg_spect=None, + bkg_energy=None, +): """ - Function that generates an image object from any arbitrary input image. + Function that generates an image object from any arbitrary input image. useful for testing realistic astrophysical sources without having to include code to simulate them here. Just have some other simulator generate an image and use this function to read that out. This function uses relative brightness of each part of the input image to generate a pmf defined at each pixel location of the image. @@ -297,9 +308,9 @@ def generate_from_image(image_path, no_photons, img_scale, energy, energy_spread # ensuring sufficient array lengths for source and bkg counts if bkg_phot is not None: - if type(bkg_phot) is float: + if isinstance(bkg_phot, float): bkg_phot = int(np.around(bkg_phot * no_photons)) - + # create image instance photon_img = image(bkg_phot + no_photons) @@ -309,23 +320,22 @@ def generate_from_image(image_path, no_photons, img_scale, energy, energy_spread # save spectrum when given if spectrum is not None: - # load file containing spectrum or copy the spectrum if isinstance(spectrum, str): photon_img.spectrum = np.loadtxt(spectrum) - else: + else: photon_img.spectrum = spectrum if energy is None: - energy = np.mean(photon_img.spectrum[:,0]) + energy = np.mean(photon_img.spectrum[:, 0]) elif energy is None: # define either spectrum or energy - raise Exception("ERROR: define either spectrum or energy") + raise ValueError("ERROR: define either spectrum or energy") # Load the image and convert it to grayscale - img = Image.open(image_path).convert('L') - + img = Image.open(image_path).convert("L") + # Convert the image to a numpy array img_array = np.array(img) @@ -333,56 +343,57 @@ def generate_from_image(image_path, no_photons, img_scale, energy, energy_spread # Generate a probability mass function from the image pmf = img_array / np.sum(img_array) - + # Draw N samples from the probability mass function source_counts = np.random.choice( np.arange(img_array.size), size=no_photons, - p=pmf.flatten(), + p=pmf.flatten(), ) - + pixel_locations = None # generate bkg photons if bkg_phot is not None: - # random point of origin within the input image (not based on FoV) - bkg_counts = np.random.choice( - np.arange(img_array.size), - size=bkg_phot - ) - - # random bkg vs based on flux - if bkg_spect == None: + bkg_counts = np.random.choice(np.arange(img_array.size), size=bkg_phot) + # random bkg vs based on flux + if bkg_spect is None: # rondomize photon TOA, while keeping the respective orders of both source and bkg photons (not based on relative fluxes) indices = np.arange(bkg_phot + no_photons) - pixel_locations = np.concatenate((np.array(bkg_counts), np.array(source_counts))) + pixel_locations = np.concatenate( + (np.array(bkg_counts), np.array(source_counts)) + ) np.random.shuffle(indices) - + # keep track of the bkg photon indices bkg_i = np.where(np.in1d(indices, np.arange(bkg_phot)))[0] photon_img.bkg_indices.extend(bkg_i) # Generating photon energies - photon_img.energies[np.arange(bkg_phot + no_photons)] = energy * spc.eV * 1e3 + photon_img.energies[np.arange(bkg_phot + no_photons)] = ( + energy * spc.eV * 1e3 + ) if bkg_energy is not None: photon_img.energies[bkg_i] = bkg_energy * spc.eV * 1e3 # Generating times of arrival # TODO: make more realistic by making Poisson times when spectra are given - photon_img.toa[np.arange(bkg_phot + no_photons)] = np.arange(bkg_phot + no_photons) + photon_img.toa[np.arange(bkg_phot + no_photons)] = np.arange( + bkg_phot + no_photons + ) else: # save spectrum when given - """TODO: sample the spectrum in process.py""" + # TODO: sample the spectrum in process.py # load file containing spectrum or copy the spectrum if isinstance(bkg_spect, str): photon_img.bkg_spect = np.loadtxt(bkg_spect) - else: + else: photon_img.bkg_spect = bkg_spect if energy is None: - energy = np.mean(photon_img.bkg_spect[:,0]) + energy = np.mean(photon_img.bkg_spect[:, 0]) # no background else: @@ -394,23 +405,36 @@ def generate_from_image(image_path, no_photons, img_scale, energy, energy_spread # Generating times of arrival # TODO: make more realistic by making Poisson times when spectra are given photon_img.toa[np.arange(no_photons)] = np.arange(no_photons) - + # Convert the flattened indices back into (x,y) coordinates - pixel_locations = np.column_stack(np.unravel_index(pixel_locations, img_array.shape)) + pixel_locations = np.column_stack( + np.unravel_index(pixel_locations, img_array.shape) + ) # Convert the sampled pixel locations to points of origin on the sky - photon_img.loc = ((pixel_locations - (pix_scale/2)) * img_scale / pix_scale.max() + np.array(offset)) * 2 * np.pi / (3600 * 360) + photon_img.loc = ( + ( + (pixel_locations - (pix_scale / 2)) * img_scale / pix_scale.max() + + np.array(offset) + ) + * 2 + * np.pi + / (3600 * 360) + ) # Adds a spread to the energies if given. - if energy_spread > 0.: - photon_img.energies += np.random.normal(0, energy_spread * spc.eV * 1e3, no_photons) + if energy_spread > 0.0: + photon_img.energies += np.random.normal( + 0, energy_spread * spc.eV * 1e3, no_photons + ) return photon_img, pix_scale + def draw_photons_from_cube(data_cube, no_photons): """ Draw photon positions (spatial) and energies from a model data cube. - + Parameters ---------- data_cube : pd.DataFrame @@ -431,7 +455,7 @@ def draw_photons_from_cube(data_cube, no_photons): """ rng = np.random.default_rng() - grouped = data_cube.groupby('energy_bin', sort=False)["flux"].sum() + grouped = data_cube.groupby("energy_bin", sort=False)["flux"].sum() total_flux = grouped.sum() if total_flux <= 0: raise ValueError("Total flux must be positive.") @@ -440,7 +464,7 @@ def draw_photons_from_cube(data_cube, no_photons): probs = pE.to_numpy() sampled_energies = rng.choice(energies, size=no_photons, p=probs) - + slices = {E: g for E, g in data_cube.groupby("energy_bin", sort=False)} records = [] @@ -450,26 +474,27 @@ def draw_photons_from_cube(data_cube, no_photons): continue slice_df = slices[E] - weights = slice_df['flux'].to_numpy() + weights = slice_df["flux"].to_numpy() if np.sum(weights) <= 0: continue weights /= np.sum(weights) indices = rng.choice(slice_df.index, size=nE, p=weights) - subset = slice_df.loc[indices, slice_df.columns.difference(['flux'])].copy() - subset['energy_bin'] = E + subset = slice_df.loc[indices, slice_df.columns.difference(["flux"])].copy() + subset["energy_bin"] = E records.append(subset) photons = pd.concat(records, ignore_index=True) im = image(no_photons) - im.energies[:] = photons['energy_bin'].to_numpy() * spc.eV * 1e3 - im.loc[:] = np.column_stack((photons['alpha'], photons['beta'])) + im.energies[:] = photons["energy_bin"].to_numpy() * spc.eV * 1e3 + im.loc[:] = np.column_stack((photons["alpha"], photons["beta"])) im.toa = np.array([i for i in range(no_photons)]) return im + def append(image1, image2): """ Function that appends two image class objects into one. @@ -483,14 +508,14 @@ def append(image1, image2): size_new = image1.size + image2.size new_image = image(size_new) - new_image.energies[0:image1.size] = image1.energies - new_image.energies[image1.size:size_new] = image2.energies + new_image.energies[0 : image1.size] = image1.energies + new_image.energies[image1.size : size_new] = image2.energies - new_image.loc[0:image1.size,:] = image1.loc - new_image.loc[image1.size:size_new,:] = image2.loc + new_image.loc[0 : image1.size, :] = image1.loc + new_image.loc[image1.size : size_new, :] = image2.loc - new_image.toa[0:image1.size] = image1.toa - new_image.toa[image1.size:size_new] = image2.toa + image1.toa[-1] + 1 + new_image.toa[0 : image1.size] = image1.toa + new_image.toa[image1.size : size_new] = image2.toa + image1.toa[-1] + 1 return new_image @@ -505,10 +530,14 @@ def create_image_from_events(event_list): """ im = image(len(event_list)) - im.energies = event_list['DETECTED_ENERGY'] * 1e3 * spc.e - im.loc[:,0] = event_list['DEC_POS'] * np.pi / 180 - im.loc[:,1] = event_list['RA_POS'] * np.pi / 180 - im.toa = event_list['TIME'].astype(int) - im.baseline_indices = (event_list['BASELINE_ID']) if 'BASELINE_ID' in event_list.columns.names else None + im.energies = event_list["DETECTED_ENERGY"] * 1e3 * spc.e + im.loc[:, 0] = event_list["DEC_POS"] * np.pi / 180 + im.loc[:, 1] = event_list["RA_POS"] * np.pi / 180 + im.toa = event_list["TIME"].astype(int) + im.baseline_indices = ( + (event_list["BASELINE_ID"]) + if "BASELINE_ID" in event_list.columns.names + else None + ) - return im \ No newline at end of file + return im diff --git a/src/curtains/instrument.py b/src/curtains/instrument.py index d1a3d23..0268232 100644 --- a/src/curtains/instrument.py +++ b/src/curtains/instrument.py @@ -1,10 +1,7 @@ - import re import math import warnings import numpy as np -import scipy.constants as spc -import scipy.interpolate as interp import matplotlib.pyplot as plt from astropy.io import fits @@ -19,12 +16,12 @@ # pos_noise = 0., energy_noise = 0., t_noise = 0., # quant_eff = r"../XRImulator/Models/detector_qe/Si_9p5_um_transmission_data.txt", # response_matrix = None): -# """ +# """ # This function is the main function that takes the 'real' photons at the camera and -# converts them to detector output data as if the detector had just detected and measured those photons. +# converts them to detector output data as if the detector had just detected and measured those photons. # It models the detector ranges, resolutions and noise for time, posistion and energie measurements # whether they are absorbed along the way, how much noise there is. - + # TODO it can be adopted to include more realistic models of detectors; and or a seperate background class; # and or energy consumption and readout times; and or more. @@ -82,10 +79,10 @@ # if self.energy_noise > 0.: # # % is for forcing it to be impossible for photons to be measured above or below energy range, while keeping random distribution # # If you want to avoid high energies bleeding over in the case of for example an emission line you want to image, -# # simply set the energy range too big to have this contamination. -# instrument_data.energies[where] = instrument_data.energies[where] + np.random.normal(0, self.energy_noise, instrument_data.energies[where].size) +# # simply set the energy range too big to have this contamination. +# instrument_data.energies[where] = instrument_data.energies[where] + np.random.normal(0, self.energy_noise, instrument_data.energies[where].size) # # - self.E_range[0]) -# # % (self.E_range[1] - self.E_range[0]) +# # % (self.E_range[1] - self.E_range[0]) # # + self.E_range[0]) # instrument_data.energies[where][instrument_data.energies[where] < self.E_range[0]] = self.E_range[0] # instrument_data.energies[where][instrument_data.energies[where] > self.E_range[1]] = self.E_range[1] @@ -107,15 +104,15 @@ # This function processes the positions at which photons arrive and # how the instrument records them. # """ - + # # Noises up the data # if self.pos_noise > 0.: # instrument_data.pos[where] = instrument_data.pos[where] + np.random.normal(0, self.pos_noise, instrument_data.pos[where].size)# - self.pos_range[0]) -# # % (self.pos_range[1] - self.pos_range[0]) +# # % (self.pos_range[1] - self.pos_range[0]) # # + self.pos_range[0]) # instrument_data.pos[where][instrument_data.pos[where] < self.pos_range[0]] = self.pos_range[0] # instrument_data.pos[where][instrument_data.pos[where] > self.pos_range[1]] = self.pos_range[1] - + # return instrument_data # def discretize_E(self, instrument_data, where): @@ -134,8 +131,8 @@ # Method that discretizes positions of incoming photons into pixel positions. # Adds an array of these locations stored to the class under the name self.discrete_pos. # """ -# instrument_data.pos[where] = (instrument_data.pos[where] - self.pos_range[0]) // self.res_pos - +# instrument_data.pos[where] = (instrument_data.pos[where] - self.pos_range[0]) // self.res_pos + # def pixel_to_pos(self, instrument_data, where): # """ Method that turns discretized positions into the positions at the center of their respective pixels. """ # instrument_data.pos[where] = (instrument_data.pos[where] + .5) * self.res_pos + self.pos_range[0] @@ -173,29 +170,40 @@ # return instrument_data -class baseline(): +class baseline: """ This class defines a single baseline in an interferometer object, and is used as a helper for the interferometer class objects. #TODO add more relevant parameters to make this more realistic. In order to fully accurately model an observation, this class can be expanded. This would necesarily also include another conceptual shift with - consequences through the rest of the code, as at the moment the image class represents a collection of all photons that will be detected, + consequences through the rest of the code, as at the moment the image class represents a collection of all photons that will be detected, which would need to shift to being a collection of photons that could be detected, with the number of input photons likely being much greater - than the detected photons. + than the detected photons. """ - def __init__(self, num_pairs, D = None, L = None, W = None, beam_angle = None, F = None, - grazing_angle = None, bench_length = None, interferometer = None, mirr_reflec = None): - """ + def __init__( + self, + num_pairs, + D=None, + L=None, + W=None, + beam_angle=None, + F=None, + grazing_angle=None, + bench_length=None, + interferometer=None, + mirr_reflec=None, + ): + """ Function that generates a single x-ray interferometer baseline according to given specifications. - - + + Parameters:\n num_pairs (int) = Number of slit-gap pairs in the slatted mirror\n D (float) = Baseline of the interferometer (in meters)\n L (float) = Length from last mirror to CCD surface (in meters)\n W (float) = Incident photon beam width (in micrometers)\n # set by projected slat width beam_angle (float) = Angle between the two beams at the detector (in radians)\n - F (float) = Effective focal length of interferometer (in meters)\n + F (float) = Effective focal length of interferometer (in meters)\n grazing_angle (float) = Angle of the mirrors with respect to the beam (in radians)\n bench_length (float) = Length of the optical bench (in meters)\n interferometer (class interferometer) = Interferometer the baseline is a part of\n @@ -208,16 +216,20 @@ def __init__(self, num_pairs, D = None, L = None, W = None, beam_angle = None, F self.beam_angle = beam_angle else: try: - self.beam_angle = W * 1e-6 / L # W to SI units + self.beam_angle = W * 1e-6 / L # W to SI units except TypeError: - raise Exception(r"ERROR: Either define the beam angle ($\theta_b$) or both the beam width ($W$) and the combining length ($L$)!") + raise Exception( + r"ERROR: Either define the beam angle ($\theta_b$) or both the beam width ($W$) and the combining length ($L$)!" + ) if D is not None: self.D = D else: try: self.D = F * self.beam_angle except TypeError: - raise Exception(r"ERROR: Either define the baseline length ($D$) or the effective focal length ($F$)!") + raise Exception( + r"ERROR: Either define the baseline length ($D$) or the effective focal length ($F$)!" + ) if F is not None: self.F = F else: @@ -228,7 +240,9 @@ def __init__(self, num_pairs, D = None, L = None, W = None, beam_angle = None, F try: self.W = self.beam_angle * L except TypeError: - raise Exception(r"ERROR: Either define the beam width ($W$) or the combining length ($L$)!") + raise Exception( + r"ERROR: Either define the beam width ($W$) or the combining length ($L$)!" + ) if L is not None: self.L = L else: @@ -240,14 +254,16 @@ def __init__(self, num_pairs, D = None, L = None, W = None, beam_angle = None, F # self.bench_length = 0.5 * self.D / np.tan(2 * grazing_angle) # except TypeError: # warnings.warn(r"Warning: if you want to check the length contraints, either define the length of the interferometer arm projected onto the optical axis ($B$)"+ - # r" or the grazing angle ($\theta_g$)!") + # r" or the grazing angle ($\theta_g$)!") if grazing_angle is not None: self.grazing_angle = grazing_angle - else: + else: try: - self.grazing_angle = np.arctan( self.D / 2 / self.bench_length) / 2 + self.grazing_angle = np.arctan(self.D / 2 / self.bench_length) / 2 except (TypeError, AttributeError): - warnings.warn(r"Warning: grazing angle could not be calculated and is set to None.") + warnings.warn( + r"Warning: grazing angle could not be calculated and is set to None." + ) self.grazing_angle = None self.num_pairs = num_pairs @@ -256,7 +272,7 @@ def __init__(self, num_pairs, D = None, L = None, W = None, beam_angle = None, F # self.camera = detector(res_E = 0.1, res_t = 1, res_pos = 2, E_range = np.array([1, 7]), pos_range = np.array([-22000, 22000])) #np.array([-1000, 1000])) #np.array([-300, 300])) ## current CMOS # initialize the search for the sampled angles in the file names - float_finder = re.compile(f'.*([0-9]+\\.[0-9]+)') + float_finder = re.compile(r".*([0-9]+\.[0-9]+)") angles = [] reflec_data = [] @@ -268,20 +284,19 @@ def __init__(self, num_pairs, D = None, L = None, W = None, beam_angle = None, F else: # loop through all sampled angles for item in mirr_reflec: - # find and save the sampled angle angles.append(float(float_finder.search(item).group(1))) # read and save the mirror reflectivity per sampled energy, assume https://henke.lbl.gov/ file format - reflec_data.append(np.loadtxt(item, skiprows = 2).T) + reflec_data.append(np.loadtxt(item, skiprows=2).T) # save to attribute - self.mirr_reflec = np.array([angles, reflec_data], dtype = object).T + self.mirr_reflec = np.array([angles, reflec_data], dtype=object).T # shows the mirror reflectivity # for angle in self.mirr_reflec: # plt.plot(angle[1][0] / 1000, angle[1][1], label = r"$\theta_{\mathrm{g}} = $" + str(angle[0]) + r"$^{\circ}$") - + # plt.ylabel("Reflectivity") # plt.xlabel("Energy (keV)") # plt.xlim(1, 7) @@ -292,11 +307,17 @@ def __init__(self, num_pairs, D = None, L = None, W = None, beam_angle = None, F # check if input values conform to physical constraints, use math.isclose against binary fraction approximation errors if not math.isclose(self.D, self.F * self.beam_angle): - raise Exception(r"ERROR: The chosen baseline length ($D$), effective focal length ($F$) and beam angle ($\theta_b$) do NOT match!") + raise Exception( + r"ERROR: The chosen baseline length ($D$), effective focal length ($F$) and beam angle ($\theta_b$) do NOT match!" + ) if not math.isclose(self.W, self.beam_angle * self.L): - raise Exception(r"ERROR: The chosen beam width ($W$), combining length ($L$) and beam angle ($\theta_b$) do NOT match!") + raise Exception( + r"ERROR: The chosen beam width ($W$), combining length ($L$) and beam angle ($\theta_b$) do NOT match!" + ) if type(self.num_pairs) is not int or self.num_pairs <= 0: - raise Exception(r"ERROR: The number of pairs must be an integer and at least one!") + raise Exception( + r"ERROR: The number of pairs must be an integer and at least one!" + ) # check if the length is within limits, only when possible # try: @@ -305,11 +326,11 @@ def __init__(self, num_pairs, D = None, L = None, W = None, beam_angle = None, F # except (TypeError, AttributeError): # warnings.warn("Warning: missing information to check the length restrictions!") - # def add_custom_detector(self, res_E, res_t, res_pos, E_range, pos_range, pos_noise = 0., E_noise = 0., t_noise = 0., + # def add_custom_detector(self, res_E, res_t, res_pos, E_range, pos_range, pos_noise = 0., E_noise = 0., t_noise = 0., # quant_eff = r"CModels/Detector QE/Si_9p5_um_transmission_data.txt", # response_matrix = None): # """ Allows for defining a custom detector. - + # Parameters:\n # res_E (float) = Energy resolution of CCD's in instrument (in KeV)\n # res_t (float) = Time resolution of CCD's in instrument (seconds)\n @@ -324,10 +345,10 @@ def __init__(self, num_pairs, D = None, L = None, W = None, beam_angle = None, F # quant_eff (string) = A file name, including the path. The file contains the quantum efficiancy of the detector in the form (1 - QE) per energy.\n # response_matrix (string) = A file name, including the path. The file contains the response_matrix of the detector.\n - + # """ # self.camera = detector(res_E, res_t, res_pos, E_range, pos_range, pos_noise, E_noise, t_noise, quant_eff, response_matrix) - + # def path_difference_change(self, times): # """Add function calls that calculate the path length difference changing over time due to different effects""" @@ -335,7 +356,7 @@ def __init__(self, num_pairs, D = None, L = None, W = None, beam_angle = None, F # return 0 * times # def eff_area(self, photon_energies): - + # # return absolute area when no reflectivity information is given, # # asumes perfect reflectivity and camera quantum efficiency # if self.mirr_reflec == None: @@ -388,34 +409,45 @@ def __init__(self, num_pairs, D = None, L = None, W = None, beam_angle = None, F # return mirr_eff_area * cam_eff_area - -class interferometer(): - """ + +class interferometer: + """ Class defining a hypothetical x-ray interferometer. It contains the code needed to generate the interferometer and adapt some of its characteristics afterwards. """ - def __init__(self, time_step = 1, wobbler = None, wobble_I = 0., wobble_c = None, wobble_file = '', - roller = None, roll_speed = 0., roll_stop_t = 0., roll_stop_a = 0., roll_init = 0, - max_ob_length = None): - """ - Function that generates a virtual x-ray interferometer according to given specifications. - - Parameters:\n - time_step (float) = arbitrary time steps in the simulator, to be based on fluxes (s)\n - - wobbler (function) = Function to use to simmulate wobble in observation (possibly not relevant here)\n - wobble_I (float) = Intensity of wobble effect, used as sigma in normally distributed random walk steps. Default is 0, which means no wobble. (in arcsec)\n - wobble_c (function) = Function to use to correct for spacecraft wobble in observation (possibly not relevant here)\n - wobble_file (file) = File containing spacecraft wobble pointing positions.\n - - roller (function) = Function to use to simulate the spacecraft rolling. Options are 'smooth_roll' and 'discrete_roll'.\n - roll_speed (float) = Indicator for how quickly spacecraft rolls around. Default is 0, meaning no roll. (in rad/sec)\n - roll_stop_t (float) = Indicator for how long spacecraft rests at specific roll if using 'discrete_roll'. Default is 0, meaning it doesn't stop. (in seconds)\n - roll_stop_a (float) = Indicator for at what angle increments spacecraft rests at if using 'discrete_roll'. Default is 0, meaning it doesn't stop. (in rads)\n - roll_init (float) = Initial roll angle. (in radians)\n - - max_ob_length (float) = Maximum length of the optical bench of the spacecraft. (in m)\n + def __init__( + self, + time_step=1, + wobbler=None, + wobble_I=0.0, + wobble_c=None, + wobble_file="", + roller=None, + roll_speed=0.0, + roll_stop_t=0.0, + roll_stop_a=0.0, + roll_init=0, + max_ob_length=None, + ): + """ + Function that generates a virtual x-ray interferometer according to given specifications. + + Parameters:\n + time_step (float) = arbitrary time steps in the simulator, to be based on fluxes (s)\n + + wobbler (function) = Function to use to simmulate wobble in observation (possibly not relevant here)\n + wobble_I (float) = Intensity of wobble effect, used as sigma in normally distributed random walk steps. Default is 0, which means no wobble. (in arcsec)\n + wobble_c (function) = Function to use to correct for spacecraft wobble in observation (possibly not relevant here)\n + wobble_file (file) = File containing spacecraft wobble pointing positions.\n + + roller (function) = Function to use to simulate the spacecraft rolling. Options are 'smooth_roll' and 'discrete_roll'.\n + roll_speed (float) = Indicator for how quickly spacecraft rolls around. Default is 0, meaning no roll. (in rad/sec)\n + roll_stop_t (float) = Indicator for how long spacecraft rests at specific roll if using 'discrete_roll'. Default is 0, meaning it doesn't stop. (in seconds)\n + roll_stop_a (float) = Indicator for at what angle increments spacecraft rests at if using 'discrete_roll'. Default is 0, meaning it doesn't stop. (in rads)\n + roll_init (float) = Initial roll angle. (in radians)\n + + max_ob_length (float) = Maximum length of the optical bench of the spacecraft. (in m)\n """ self.baselines = [] @@ -436,7 +468,7 @@ def __init__(self, time_step = 1, wobbler = None, wobble_I = 0., wobble_c = None self.max_ob_length = max_ob_length def random_wobble(self, pointing): - """ + """ Function that adds 'wobble' to the spacecraft, slightly offsetting its pointing every timestep. It models wobble as a random walk with a given intensity that is used as the sigma for a normally distributed step size in both the pitch and yaw directions. @@ -449,16 +481,18 @@ def random_wobble(self, pointing): pointing (array): That same, but now with wobble data.\n """ - pointing[1:, :2] = pointing[:-1, :2] + np.random.normal(0, self.wobble_I, size=(len(pointing[:, 0]) - 1, 2)) * 2 * np.pi / (3600 * 360) + pointing[1:, :2] = pointing[:-1, :2] + np.random.normal( + 0, self.wobble_I, size=(len(pointing[:, 0]) - 1, 2) + ) * 2 * np.pi / (3600 * 360) return pointing - + def file_wobble(self, pointing): - """ + """ Function that adds 'wobble' to the spacecraft, slightly offsetting its pointing every timestep. - This function uses an input file in a csv format (with ',' as delimiter) to read out pointing data, + This function uses an input file in a csv format (with ',' as delimiter) to read out pointing data, probably generated with a different simulator. #TODO This function is mostly a placeholder, to be replaced later to adapt to the actual format this data - will take. This is only one way it could look, but it should how to structure an eventual replacement for + will take. This is only one way it could look, but it should how to structure an eventual replacement for whoever wants to adapt the code. Parameters: @@ -469,12 +503,12 @@ def file_wobble(self, pointing): pointing (array): That same, but now with wobble data.\n """ - pointing[:, :2] = np.genfromtxt(self.wobble_file, np.float64, delimiter=',') + pointing[:, :2] = np.genfromtxt(self.wobble_file, np.float64, delimiter=",") return pointing def smooth_roller(self, pointing): """ - Function that generates the roll portion of the pointing data for the instrument. + Function that generates the roll portion of the pointing data for the instrument. This function is used for a continuous model of rolling the instrument, with a predefined roll velocity. @@ -488,12 +522,14 @@ def smooth_roller(self, pointing): """ self.roll_speed = np.pi / pointing[:, 2].size - pointing[:, 2] = (np.arange(pointing[:, 2].size) * self.roll_speed * self.time_step) + self.roll_init + pointing[:, 2] = ( + np.arange(pointing[:, 2].size) * self.roll_speed * self.time_step + ) + self.roll_init return pointing def discrete_roller(self, pointing): """ - Function that generates the roll portion of the pointing data for the instrument. + Function that generates the roll portion of the pointing data for the instrument. This function is used for a discrete model of rolling the instrument, with starts and stops at specified roll angle intervals. @@ -506,7 +542,7 @@ def discrete_roller(self, pointing): pointing (array): That same, but now with roll data. """ - # Calculates the stopping interval in timestep units + # Calculates the stopping interval in timestep units time_to_move = self.roll_stop_t # The angle over which to move after the stopping interval angle_to_move = self.roll_stop_a @@ -517,8 +553,7 @@ def discrete_roller(self, pointing): else: pointing[i, 2] = 0 - - # # Calculates the stopping interval in timestep units + # # Calculates the stopping interval in timestep units # time_to_move = self.roll_stop_t // self.time_step # # The angle over which to move after the stopping interval # angle_to_move = self.roll_stop_a @@ -535,14 +570,14 @@ def discrete_roller(self, pointing): # if pointing[i, 2] > angle_to_move: # angle_to_move += self.roll_stop_a # time_to_move += self.roll_stop_t // self.time_step - + return pointing def gen_pointing(self, t_exp): - """ - This function generates a 3d pointing vector for each time step in an observation. It consists of - three angles, the pitch, yaw and roll. The first two are linked and generated together by the wobbler - function, while the roll is fundamentally different and thus generated differently. If no wobbler or + """ + This function generates a 3d pointing vector for each time step in an observation. It consists of + three angles, the pitch, yaw and roll. The first two are linked and generated together by the wobbler + function, while the roll is fundamentally different and thus generated differently. If no wobbler or roller are given, the corresponding pointing values will be zero, indicating stillness. """ pointing = np.zeros((t_exp + 2, 3)) @@ -555,27 +590,49 @@ def gen_pointing(self, t_exp): return pointing - def add_baseline(self, num_pairs, D = None, L = None, W = None, beam_angle = None, F = None, - grazing_angle = None, bench_length = None, interferometer = None, mirr_mater = None): + def add_baseline( + self, + num_pairs, + D=None, + L=None, + W=None, + beam_angle=None, + F=None, + grazing_angle=None, + bench_length=None, + interferometer=None, + mirr_mater=None, + ): """ Function that adds a baseline of given parameters to the interferometer object. Call this function multiple times to construct a full interferometer capable of actually observing images. Without these, no photons can be measured. - + Parameters: num_pairs (int) = Number of slit-gap pairs in the slatted mirror\n D (float) = Baseline of the interferometer (in meters)\n L (float) = Length from last mirror to CCD surface (in meters)\n W (float) = Incident photon beam width (in micrometers)\n # set by projected slat width beam_angle (float) = Angle between the two beams at the detector (in radians)\n - F (float) = Effective focal length of interferometer (in meters)\n + F (float) = Effective focal length of interferometer (in meters)\n grazing_angle (float) = Angle of the mirrors with respect to the beam (in radians)\n bench_length (float) = Length of the optical bench (in meters)\n interferometer (class interferometer) = Interferometer the baseline is a part of\n mirr_mater (float) = The refrective index of the mirror material\n """ - self.baselines.append(baseline(num_pairs, D, L, W, beam_angle, F, - grazing_angle, bench_length, interferometer, - mirr_mater)) + self.baselines.append( + baseline( + num_pairs, + D, + L, + W, + beam_angle, + F, + grazing_angle, + bench_length, + interferometer, + mirr_mater, + ) + ) def add_willingale_baseline(self, D): """ @@ -583,17 +640,16 @@ def add_willingale_baseline(self, D): Call this function multiple times to construct a full interferometer capable of actually observing images. Without these, no photons can be measured. """ - - self.baselines.append(baseline(num_pairs = 30, D = D, L = 10, W = 300, bench_length= 7)) + + self.baselines.append(baseline(num_pairs=30, D=D, L=10, W=300, bench_length=7)) def clear_baselines(self): self.baselines.clear() - -def plot_arf(arf_file, baseline_mode = False): - """ - Function that plots the effective area from a given arf file. +def plot_arf(arf_file, baseline_mode=False): + """ + Function that plots the effective area from a given arf file. Parameters: arf_file (string) = Path to the arf file to plot.\n baseline_mode (bool) = Whether to plot in baseline mode (effective area per baseline) or not (total effective area). Default is False.\n @@ -601,21 +657,39 @@ def plot_arf(arf_file, baseline_mode = False): arf_data = fits.open(arf_file) - energies = arf_data[1].data['ENERG_LO'] + (arf_data[1].data['ENERG_HI'] - arf_data[1].data['ENERG_LO']) / 2 - eff_area = arf_data[1].data['SPECRESP'] + energies = ( + arf_data[1].data["ENERG_LO"] + + (arf_data[1].data["ENERG_HI"] - arf_data[1].data["ENERG_LO"]) / 2 + ) + eff_area = arf_data[1].data["SPECRESP"] # there are additional baseline extensions to arf files, which store the effective area per baseline, named 'SPECRESP_BL{i}' where i starts from 1 and is the baseline ID, find maximum i and plot all arfs per baseline if in baseline mode - plt.figure(figsize=(7,5)) + plt.figure(figsize=(7, 5)) if baseline_mode: i = 1 - while 'SPECRESP_BL' + str(i) in arf_data: - eff_area_bl = arf_data['SPECRESP_BL' + str(i)].data['SPECRESP_BL' + str(i)] + while "SPECRESP_BL" + str(i) in arf_data: + eff_area_bl = arf_data["SPECRESP_BL" + str(i)].data["SPECRESP_BL" + str(i)] # get the baseline in meters from the baseline extension header - baseline_length = arf_data['SPECRESP_BL' + str(i)].header.get('BASELINE', None) - plt.plot(energies, eff_area_bl, label = f'Baseline {i} ({baseline_length:.2f} m)', color='mediumvioletred', linewidth=1, alpha= 1/i ) + baseline_length = arf_data["SPECRESP_BL" + str(i)].header.get( + "BASELINE", None + ) + plt.plot( + energies, + eff_area_bl, + label=f"Baseline {i} ({baseline_length:.2f} m)", + color="mediumvioletred", + linewidth=1, + alpha=1 / i, + ) i += 1 - plt.plot(energies, eff_area, label = 'Total Effective Area', color='mediumvioletred', linewidth=2) + plt.plot( + energies, + eff_area, + label="Total Effective Area", + color="mediumvioletred", + linewidth=2, + ) plt.xlabel("Energy (keV)") plt.ylabel("Effective area (cm$^2$)") plt.yscale("log") @@ -627,4 +701,4 @@ def plot_arf(arf_file, baseline_mode = False): else: plt.ylim(1e2, np.max(eff_area) * 1.1) plt.grid("both") - plt.show() \ No newline at end of file + plt.show() diff --git a/src/curtains/process.py b/src/curtains/process.py index fec8ee8..ba2f14c 100644 --- a/src/curtains/process.py +++ b/src/curtains/process.py @@ -1,4 +1,3 @@ - import warnings import numpy as np import scipy.special as sps @@ -6,36 +5,43 @@ import scipy.constants as spc import scipy.interpolate as interp from astropy.io import fits -from astropy import units as u # for testing import matplotlib.pyplot as plt -class interferometer_data(): - """ +class interferometer_data: + """ Class that serves as a container for interferometer output data. Constructed as such for ease of use by way of standardization. Does not contain manipulation methods, data inside will have to be edited via external methods. """ - def __init__(self, instrument, image, eff_area_show = True, pure_diffraction = False, pure_fringes = False): - """ + + def __init__( + self, + instrument, + image, + eff_area_show=True, + pure_diffraction=False, + pure_fringes=False, + ): + """ This function is the main function that takes an image and converts it to instrument data as - if the instrument had just observed the object the image is a representation of. - It models the individual photons coming in each timestep, at what detector they end up, - whether they are absorbed along the way, how much noise there is, and also whether the + if the instrument had just observed the object the image is a representation of. + It models the individual photons coming in each timestep, at what detector they end up, + whether they are absorbed along the way, how much noise there is, and also whether the spacecraft the instrument is on wobbles, and possible correction for this. Parameters:\n instrument (interferometer class object) = Instrument object to be used to simulate observing the image.\n image (image class object) = Image object to be observed.\n - + eff_area_show (bool) = whether or not to show the effective area of the full instrument, mirrors and detector\n pure_diffraction (bool) = sample only the diffraction pattern, without the fringe pattern\n pure_fringes (bool) = sample only the fringe pattern, without the diffraction pattern\n """ - + # Useful shorthands self.size = image.size @@ -45,35 +51,23 @@ def __init__(self, instrument, image, eff_area_show = True, pure_diffraction = F # simulates the energies of the source photons going through the optical benches until they reach the detector # self.actual_energies, self.baseline_indices = self.process_photon_e_base(instrument, image, eff_area_show) - self.actual_energies, self.baseline_indices = image.energies, image.baseline_indices + self.actual_energies, self.baseline_indices = ( + image.energies, + image.baseline_indices, + ) # detector values, which can be subjected to detector effects self.energies = self.actual_energies self.toa = self.image_toa # simulates the positions (in m) of the source photons going through the optical benches until they reach the detector - self.actual_pos = self.process_photon_dpos(instrument, image, pure_diffraction, pure_fringes) - + self.actual_pos = self.process_photon_dpos( + instrument, image, pure_diffraction, pure_fringes + ) + # detector value, which can be subjected to detector effects, digitize to pixel positions - self.pos = np.round(self.actual_pos/5e-6) * 5e-6 + self.pos = np.round(self.actual_pos / 5e-6) * 5e-6 - def copy(self): - """ - Function to create a copy of the interferometer_data object. - """ - new_copy = interferometer_data.__new__(interferometer_data) - new_copy.size = self.size - new_copy.image_energies = np.copy(self.image_energies) - new_copy.image_toa = np.copy(self.image_toa) - new_copy.actual_energies = np.copy(self.actual_energies) - new_copy.baseline_indices = np.copy(self.baseline_indices) - new_copy.energies = np.copy(self.energies) - new_copy.toa = np.copy(self.toa) - new_copy.actual_pos = np.copy(self.actual_pos) - new_copy.pos = np.copy(self.pos) - new_copy.pointing = np.copy(self.pointing) - return new_copy - - def process_photon_e_base(self, instrument, image, eff_area_show = True): + def process_photon_e_base(self, instrument, image, eff_area_show=True): """ This function is a helper function for process_image that specifically processes what energy detected photons, which impact on the detector, are expected. Not to be used outside the process_image context. @@ -82,7 +76,7 @@ def process_photon_e_base(self, instrument, image, eff_area_show = True): instrument (interferometer class object) = Instrument object to be used to simulate observing the image.\n image (image class object) = Image object to be observed.\n - + eff_area_show (bool) = whether or not to show the effective area of the full instrument, mirrors and detector\n """ @@ -92,12 +86,18 @@ def process_photon_e_base(self, instrument, image, eff_area_show = True): # Select random baseline based on effective area per basline if image.spectrum is None: - # pure effective area calculation - eff_areas = np.array([baseline.eff_area(np.unique(energies)) for baseline in instrument.baselines]) - relative_eff_areas = eff_areas / np.sum(eff_areas, axis = 0) - - baseline_indices = np.random.choice(len(instrument.baselines), array_lengths, p = relative_eff_areas.T[0]) #np.random.randint(0, len(instrument.baselines), array_lengths) # + eff_areas = np.array( + [ + baseline.eff_area(np.unique(energies)) + for baseline in instrument.baselines + ] + ) + relative_eff_areas = eff_areas / np.sum(eff_areas, axis=0) + + baseline_indices = np.random.choice( + len(instrument.baselines), array_lengths, p=relative_eff_areas.T[0] + ) # np.random.randint(0, len(instrument.baselines), array_lengths) # else: # effective area calculation, including source spectra @@ -107,23 +107,41 @@ def process_photon_e_base(self, instrument, image, eff_area_show = True): # loop through the baselines for baseline in instrument.baselines: - # pull out detector for its attributes detector = baseline.camera # define the integration energies, 100 is used to integrate finer then the energy resolution # how much better energy sampling, compared to detector energy resolution finer_binning = 1e2 - integr_energy = np.linspace(detector.E_range[0], detector.E_range[1], (finer_binning * np.ceil((detector.E_range[1] - detector.E_range[0]) / detector.res_E)).astype(int)) + integr_energy = np.linspace( + detector.E_range[0], + detector.E_range[1], + ( + finer_binning + * np.ceil( + (detector.E_range[1] - detector.E_range[0]) / detector.res_E + ) + ).astype(int), + ) # get the effective area of the baseline, per energy eff_area = baseline.eff_area(integr_energy).T # spectra assumed to use keV, turned into SI units of J - interp_spec = interp.interp1d(image.spectrum[:,0] * 1e3 * spc.eV, image.spectrum[:,1], bounds_error = False, fill_value = 0)(integr_energy) - + interp_spec = interp.interp1d( + image.spectrum[:, 0] * 1e3 * spc.eV, + image.spectrum[:, 1], + bounds_error=False, + fill_value=0, + )(integr_energy) + # integrate the spectrum multiplied with the effective areas over energies, for photons per baseline - cts_persec = np.sum(interp_spec * eff_area * (integr_energy[-1] - integr_energy[0]) / integr_energy.size) + cts_persec = np.sum( + interp_spec + * eff_area + * (integr_energy[-1] - integr_energy[0]) + / integr_energy.size + ) cts_perbaseline_persec.append(cts_persec) # otherwise difficult to find error @@ -134,11 +152,23 @@ def process_photon_e_base(self, instrument, image, eff_area_show = True): e_cumsum = np.cumsum(interp_spec * eff_area) e_eCDF = e_cumsum / e_cumsum[-1] min_val_eCDF.append(np.min(e_eCDF)) - all_eCDF.append(interp.interp1d(e_eCDF, integr_energy))#, bounds_error = False, fill_value = np.nan)) + all_eCDF.append( + interp.interp1d(e_eCDF, integr_energy) + ) # , bounds_error = False, fill_value = np.nan)) if eff_area_show: - plt.plot(integr_energy / (1e3 * spc.eV), eff_area, color = '#ff7f0e', label = "Effective area") - plt.hlines(np.power(detector.pos_range[1] - detector.pos_range[0], 2), 0, 10, label = "Real area") + plt.plot( + integr_energy / (1e3 * spc.eV), + eff_area, + color="#ff7f0e", + label="Effective area", + ) + plt.hlines( + np.power(detector.pos_range[1] - detector.pos_range[0], 2), + 0, + 10, + label="Real area", + ) # plt.title("Effective area curve of mirrors and detector") plt.ylabel("Effective area ($m^2$)") plt.xlabel("Energy ($keV$)") @@ -151,20 +181,27 @@ def process_photon_e_base(self, instrument, image, eff_area_show = True): plt.show() # Select random baseline based on effective area per basline - baseline_indices = np.random.choice(len(instrument.baselines), array_lengths, p = cts_perbaseline_persec / np.sum(cts_perbaseline_persec, axis = 0)) #np.random.randint(0, len(instrument.baselines), array_lengths) # - + baseline_indices = np.random.choice( + len(instrument.baselines), + array_lengths, + p=cts_perbaseline_persec / np.sum(cts_perbaseline_persec, axis=0), + ) # np.random.randint(0, len(instrument.baselines), array_lengths) # + # Looking only at the baselines that have associated photons for index, baseline in enumerate(instrument.baselines): - # Taking only relevant photons from the current baseline - in_baseline = np.array(instrument.baselines)[baseline_indices] == baseline + in_baseline = ( + np.array(instrument.baselines)[baseline_indices] == baseline + ) # inversion sampling of photon energies, based on the eCDFs - energies[in_baseline] = all_eCDF[index](np.random.uniform(min_val_eCDF[index], 1, np.sum(in_baseline))) - + energies[in_baseline] = all_eCDF[index]( + np.random.uniform(min_val_eCDF[index], 1, np.sum(in_baseline)) + ) + return energies, baseline_indices - def fre_dif(self, wavelength, baseline, samples, y_pos = None): + def fre_dif(self, wavelength, baseline, samples, y_pos=None): """ Helper function that calculates the Fresnell difraction pattern for a beam such as the case in the interferometer. Also aplicable for, but less efficient at, @@ -181,28 +218,31 @@ def fre_dif(self, wavelength, baseline, samples, y_pos = None): # see Willingale (2004) for the definitions of the dimensionless coordinate u, # used for the Fresnel integrals u_0 = baseline.W * np.sqrt(2 / (wavelength * baseline.L)) - u_1 = lambda u, u_0: u - u_0/2 - u_2 = lambda u, u_0: u + u_0/2 # Only sample the slit size, or set values if y_pos is None: y_pos = np.linspace(-baseline.W / 2, baseline.W / 2, int(samples)) u = y_pos * np.sqrt(2 / (wavelength * baseline.L)) + u_1 = u - u_0 / 2 + u_2 = u + u_0 / 2 + # Fresnel integrals - S_1, C_1 = sps.fresnel(u_1(u, u_0)) - S_2, C_2 = sps.fresnel(u_2(u, u_0)) + S_1, C_1 = sps.fresnel(u_1) + S_2, C_2 = sps.fresnel(u_2) # Fresnel amplitude - A = ((C_2 - C_1) + 1j*(S_2 - S_1)) + A = (C_2 - C_1) + 1j * (S_2 - S_1) # Fresnel intensity A_star = np.conjugate(A) - I = np.abs(A * A_star) + Int = np.abs(A * A_star) - return I, u + return Int, u - def process_photon_dpos(self, instrument, image, pure_diffraction = False, pure_fringes = False): + def process_photon_dpos( + self, instrument, image, pure_diffraction=False, pure_fringes=False + ): """ This function is a helper function for process_image that specifically processes the locations where photons impact on the detector (hence the d(etector)pos(ition) name). Not to be used outside the process_image context. @@ -211,10 +251,10 @@ def process_photon_dpos(self, instrument, image, pure_diffraction = False, pure_ instrument (interferometer class object) = Instrument object to be used to simulate observing the image.\n image (image class object) = Image object to be observed.\n - + pure_diffraction (bool) = sample only the diffraction pattern, without the fringe pattern\n pure_fringes (bool) = sample only the fringe pattern, without the diffraction pattern\n - + """ array_lengths = self.size @@ -225,11 +265,12 @@ def process_photon_dpos(self, instrument, image, pure_diffraction = False, pure_ # Relative position is useful for the calculation of theta, since the off-axis angle is dependent on pointing position and orientation self.pointing = instrument.gen_pointing(int(np.max(self.image_toa))) pos_rel = self.pointing[self.image_toa, :2] - image.loc - theta = np.cos(self.pointing[self.image_toa, 2] - np.arctan2(pos_rel[:, 0], pos_rel[:, 1])) * np.sqrt(pos_rel[:, 0]**2 + pos_rel[:, 1]**2) - + theta = np.cos( + self.pointing[self.image_toa, 2] - np.arctan2(pos_rel[:, 0], pos_rel[:, 1]) + ) * np.sqrt(pos_rel[:, 0] ** 2 + pos_rel[:, 1] ** 2) + # only populated baselines are selected for baseline_i in np.unique(self.baseline_indices): - # which indices apply photons_in_baseline = self.baseline_indices == baseline_i @@ -250,67 +291,118 @@ def process_photon_dpos(self, instrument, image, pure_diffraction = False, pure_ # accept reject until all photons have been sampled while len(indices) > 0: - # how many photons still need to be sampled in order to reach the specified number, # no photons are lost samples_left = len(indices) # range on the detector where photons can arrive - width_box = np.array([(-baseline.W * baseline.num_pairs / 2) - baseline.L * thetas, (baseline.W * baseline.num_pairs / 2) - baseline.L * thetas]) + width_box = np.array( + [ + (-baseline.W * baseline.num_pairs / 2) - baseline.L * thetas, + (baseline.W * baseline.num_pairs / 2) - baseline.L * thetas, + ] + ) # Tested to be just above the maximum of the Fresnel diffraction, # for the height of the box needed by the accept reject method. - height_box = [0, 2.7 * self.fre_dif(wavelengths, baseline, 1, y_pos = 0)[0]] + height_box = [ + 0, + 2.7 * self.fre_dif(wavelengths, baseline, 1, y_pos=0)[0], + ] # the fringe spacing, see Willingale (2004) Eq. 1 - fringe_spacing = (wavelengths / baseline.beam_angle) + fringe_spacing = wavelengths / baseline.beam_angle # the maximum number of visable fringes - num_fringes = np.ceil((width_box[1] - width_box[0]) * baseline.num_pairs / fringe_spacing) - + num_fringes = np.ceil( + (width_box[1] - width_box[0]) * baseline.num_pairs / fringe_spacing + ) + # location on the detector with pathlength diffrence of zero, assuming the mirrors # have been alligned to correct for the distance between the combining mirrors - off_set_fringes = -(baseline.D * np.sin(thetas)) / (2 * np.sin(baseline.beam_angle / 2)) + off_set_fringes = -(baseline.D * np.sin(thetas)) / ( + 2 * np.sin(baseline.beam_angle / 2) + ) # find the left most fringe partially or entirely able to reach the detector - left_fringe = off_set_fringes - (np.ceil((off_set_fringes - baseline.L * thetas + baseline.W * baseline.num_pairs / 2) / fringe_spacing) * fringe_spacing) + left_fringe = off_set_fringes - ( + np.ceil( + ( + off_set_fringes + - baseline.L * thetas + + baseline.W * baseline.num_pairs / 2 + ) + / fringe_spacing + ) + * fringe_spacing + ) # samples where in the fringe each photon arrives - rand_cos = 0.5 * fringe_spacing * sampler.cosine.rvs(loc = 2 * np.pi * left_fringe / fringe_spacing, size = samples_left) / np.pi - + rand_cos = ( + 0.5 + * fringe_spacing + * sampler.cosine.rvs( + loc=2 * np.pi * left_fringe / fringe_spacing, size=samples_left + ) + / np.pi + ) + if pure_diffraction: # selects a random uniform position, resulting in no fringe sampling - rand_pos = np.random.uniform(width_box[0], width_box[1], size = samples_left) + rand_pos = np.random.uniform( + width_box[0], width_box[1], size=samples_left + ) else: # selects wich of the visable fringes each photon arrives in - rand_pos = rand_cos + (np.random.randint(low = 0, high = num_fringes + 1, size = samples_left) * fringe_spacing) - + rand_pos = rand_cos + ( + np.random.randint( + low=0, high=num_fringes + 1, size=samples_left + ) + * fringe_spacing + ) + # let the user know when the fringes would in practice move outside of the FoV - if np.logical_or((off_set_fringes < width_box[0]).any(), (off_set_fringes > width_box[1]).any()): + if np.logical_or( + (off_set_fringes < width_box[0]).any(), + (off_set_fringes > width_box[1]).any(), + ): warnings.warn("There are fringe centres outside of the FOV!") # boolean for wich y posistions don't fall outside of the allowed range, # due to sampling an integer number of fringes and not fractional ammounts - valid_rand_pos = np.logical_and(rand_pos >= width_box[0], rand_pos <= width_box[1]) #np.repeat(True, samples_left) # + valid_rand_pos = np.logical_and( + rand_pos >= width_box[0], rand_pos <= width_box[1] + ) # np.repeat(True, samples_left) # # map sampled positions to the central diffraction pattern # the diffraction pattern is assumed to be multiple patterns # next to one another without interfering if (baseline.num_pairs % 2) == 0: - diffrac_pos = (abs(rand_pos + baseline.L * thetas) % baseline.W) - (baseline.W / 2) #rand_pos # + diffrac_pos = (abs(rand_pos + baseline.L * thetas) % baseline.W) - ( + baseline.W / 2 + ) # rand_pos # else: - diffrac_pos = (abs(rand_pos + baseline.L * thetas + baseline.W / 2) % baseline.W) - (baseline.W / 2) #rand_pos # + diffrac_pos = ( + abs(rand_pos + baseline.L * thetas + baseline.W / 2) + % baseline.W + ) - (baseline.W / 2) # rand_pos # # accept reject - diffraction_value, _ = self.fre_dif(wavelengths, baseline, 1, diffrac_pos) - rand_value = np.random.uniform(height_box[0], height_box[1], size = samples_left) + diffraction_value, _ = self.fre_dif( + wavelengths, baseline, 1, diffrac_pos + ) + rand_value = np.random.uniform( + height_box[0], height_box[1], size=samples_left + ) if pure_fringes: accept_reject = np.repeat(True, samples_left) else: accept_reject = rand_value <= diffraction_value # keep the accepted y posistions - y_pos[indices[accept_reject * valid_rand_pos]] = rand_pos[accept_reject * valid_rand_pos] + y_pos[indices[accept_reject * valid_rand_pos]] = rand_pos[ + accept_reject * valid_rand_pos + ] # only keep the rejected theta's and wavelengths for the next round of accept reject thetas = thetas[np.invert(accept_reject * valid_rand_pos)] @@ -322,43 +414,48 @@ def process_photon_dpos(self, instrument, image, pure_diffraction = False, pure_ actual_pos[photons_in_baseline] = y_pos return actual_pos - -def generate_event_list(src_file, arf_path, rmf_path, exposure_time, output_file='event_list.fits'): + +def generate_event_list( + src_file, arf_path, rmf_path, exposure_time, output_file="event_list.fits" +): # Load spectrum data with fits.open(src_file) as hdul: - spectrum_data = hdul['SPECTRUM'].data - src_id = hdul['SRC_CAT'].data['SRC_ID'] - coords = hdul['SRC_CAT'].data['RA'], hdul['SRC_CAT'].data['DEC'] + spectrum_data = hdul["SPECTRUM"].data + src_id = hdul["SRC_CAT"].data["SRC_ID"] + coords = hdul["SRC_CAT"].data["RA"], hdul["SRC_CAT"].data["DEC"] hdul.close() with fits.open(arf_path) as hdul_arf: num_baselines = len(hdul_arf) - 2 - multiplied_spectra = np.zeros((num_baselines, len(src_id), len(spectrum_data['ENERGY'][0]))) + multiplied_spectra = np.zeros( + (num_baselines, len(src_id), len(spectrum_data["ENERGY"][0])) + ) for i in range(num_baselines): - ext_name = f'SPECRESP_BL{i+1}' + ext_name = f"SPECRESP_BL{i + 1}" arf_data = hdul_arf[ext_name].data # Calculate ARF energy midpoints - arf_energy = (arf_data['ENERG_LO'] + arf_data['ENERG_HI']) / 2.0 + arf_energy = (arf_data["ENERG_LO"] + arf_data["ENERG_HI"]) / 2.0 arf_resp = arf_data[ext_name] - energy_bin_width = arf_data['ENERG_HI'] - arf_data['ENERG_LO'] - + energy_bin_width = arf_data["ENERG_HI"] - arf_data["ENERG_LO"] + # Iterate over sources for j in range(len(src_id)): - src_energy = spectrum_data['ENERGY'][j] - src_flux = spectrum_data['FLUXDENSITY'][j] - + src_energy = spectrum_data["ENERGY"][j] + src_flux = spectrum_data["FLUXDENSITY"][j] + # Interpolate ARF effective area to the source energy grid for multiplication - arf_resp_interp = np.interp(src_energy, arf_energy, arf_resp, left=0, right=0) - + arf_resp_interp = np.interp( + src_energy, arf_energy, arf_resp, left=0, right=0 + ) + # Multiply flux by effective area and energy bin width to get counts res = src_flux * arf_resp_interp * energy_bin_width multiplied_spectra[i, j, :] = res hdul_arf.close() - # draw photons from multiplied spectra and create event list all_energies = [] all_baselines = [] @@ -366,17 +463,16 @@ def generate_event_list(src_file, arf_path, rmf_path, exposure_time, output_file all_ra = [] all_dec = [] - for i in range(num_baselines): for j in range(len(src_id)): # Calculate expected counts per energy bin expected_counts_per_bin = multiplied_spectra[i, j, :] * exposure_time - + # Sample actual counts using Poisson statistics counts_per_bin = np.random.poisson(expected_counts_per_bin) - + # Get energy values for this source - energies = spectrum_data['ENERGY'][j] + energies = spectrum_data["ENERGY"][j] # Generate photon events: repeat energy value N times where N is the sampled count sampled_energies = np.repeat(energies, counts_per_bin) sampled_baselines = np.full(len(sampled_energies), i, dtype=int) @@ -385,7 +481,7 @@ def generate_event_list(src_file, arf_path, rmf_path, exposure_time, output_file # Add sky location to sampled photons ra_pos = coords[0][j] dec_pos = coords[1][j] - + sampled_ra = np.full(len(sampled_energies), ra_pos) sampled_dec = np.full(len(sampled_energies), dec_pos) all_times.append(sampled_times) @@ -394,89 +490,87 @@ def generate_event_list(src_file, arf_path, rmf_path, exposure_time, output_file all_ra.append(sampled_ra) all_dec.append(sampled_dec) - # Concatenate all events into single arrays if all_energies: - # Concatenate source energies first (true energies) true_energies = np.concatenate(all_energies) - + # Load RMF with fits.open(rmf_path) as hdul_rmf: - rmf_data = hdul_rmf['MATRIX'].data - ebounds_data = hdul_rmf['EBOUNDS'].data + rmf_data = hdul_rmf["MATRIX"].data + ebounds_data = hdul_rmf["EBOUNDS"].data hdul_rmf.close() # RMF energy grid - f_chan = rmf_data['F_CHAN'] - n_chan = rmf_data['N_CHAN'] - matrix = rmf_data['MATRIX'] - energ_lo = rmf_data['ENERG_LO'] - energ_hi = rmf_data['ENERG_HI'] - + f_chan = rmf_data["F_CHAN"] + n_chan = rmf_data["N_CHAN"] + matrix = rmf_data["MATRIX"] + energ_lo = rmf_data["ENERG_LO"] + energ_hi = rmf_data["ENERG_HI"] + # EBOUNDS for channel to energy conversion - channel_min = ebounds_data['E_MIN'] - channel_max = ebounds_data['E_MAX'] - + channel_min = ebounds_data["E_MIN"] + channel_max = ebounds_data["E_MAX"] + final_energies = [] - + # Redistribution Loop (simplistic approach for demonstration) # Map true energies to RMF energy bins # This can be slow for large number of photons, vectorized approaches are complex due to sparse matrix # Finding the RMF bin index for each photon bin_indices = np.searchsorted(energ_hi, true_energies) - + # Filter out photons outside RMF range valid_mask = (bin_indices < len(energ_lo)) & (true_energies >= energ_lo[0]) bin_indices = bin_indices[valid_mask] - + # Update other arrays to match valid photons final_baselines = np.concatenate(all_baselines)[valid_mask] final_times = np.concatenate(all_times)[valid_mask] - + # Sample channels sampled_channels = [] - + for idx in bin_indices: # Get matrix row for this energy bin m_row = matrix[idx] - + # In FITS RMF, sparse storage: # F_CHAN is array of start channels for groups # N_CHAN is array of number of channels in each group # MATRIX is flat array of probabilities - + if len(m_row) == 0: - sampled_channels.append(-1) # No response + sampled_channels.append(-1) # No response continue - + # Reconstruct full probability array for this energy bin is needed for choice # However, usually we can just construct the specific valid channels and probs - + current_f_chan = f_chan[idx] current_n_chan = n_chan[idx] - + # If scalar (not variable length array logic in some files), handling depends on file format nuances # Assuming standard OGIP FITS variable length arrays - + if np.isscalar(current_f_chan): current_f_chan = [current_f_chan] current_n_chan = [current_n_chan] - + possible_channels = [] probabilities = [] - + # Offset in the flattened matrix array for the current row mat_idx = 0 for k in range(len(current_f_chan)): start = current_f_chan[k] count = current_n_chan[k] chunk = m_row[mat_idx : mat_idx + count] - + possible_channels.extend(range(start, start + count)) probabilities.extend(chunk) mat_idx += count - + probabilities = np.array(probabilities) # Ensure non-negative probabilities (clip small negative values to 0) @@ -484,22 +578,22 @@ def generate_event_list(src_file, arf_path, rmf_path, exposure_time, output_file prob_sum = probabilities.sum() if prob_sum > 0: - probabilities = probabilities / prob_sum # Normalize - + probabilities = probabilities / prob_sum # Normalize + # Sample one channel sampled_chan = np.random.choice(possible_channels, p=probabilities) sampled_channels.append(sampled_chan) else: sampled_channels.append(-1) - + sampled_channels = np.array(sampled_channels) - + # Filter out failed encodings valid_chan_mask = sampled_channels != -1 sampled_channels = sampled_channels[valid_chan_mask] final_baselines = final_baselines[valid_chan_mask] final_times = final_times[valid_chan_mask] - + # Ensure channel indices are within bounds valid_bounds = sampled_channels < len(channel_min) sampled_channels = sampled_channels[valid_bounds] @@ -507,53 +601,54 @@ def generate_event_list(src_file, arf_path, rmf_path, exposure_time, output_file final_times = np.round(final_times[valid_bounds], 3) ra_pos = np.concatenate(all_ra)[valid_bounds] dec_pos = np.concatenate(all_dec)[valid_bounds] - - final_energies = (channel_min[sampled_channels] + channel_max[sampled_channels]) / 2.0 + + final_energies = ( + channel_min[sampled_channels] + channel_max[sampled_channels] + ) / 2.0 else: final_energies = np.array([]) final_baselines = np.array([]) final_times = np.array([]) - - # Create a combined event list num_events = len(final_energies) - event_list = np.column_stack((final_baselines, true_energies, final_energies, final_times)) + event_list = np.column_stack( + (final_baselines, true_energies, final_energies, final_times) + ) print(f"Total photons generated: {num_events}") - # save event list to fits file primary_hdu = fits.PrimaryHDU() cols = [ - fits.Column(name='BASELINE_ID', format='I', array=event_list[:,0]), - fits.Column(name='TRUE_ENERGY', format='E', array=event_list[:,1]), - fits.Column(name='DETECTED_ENERGY', format='E', array=event_list[:,2]), - fits.Column(name='TIME', format='E', array=event_list[:,3]), - fits.Column(name='RA_POS', format='D', array=ra_pos), - fits.Column(name='DEC_POS', format='D', array=dec_pos) + fits.Column(name="BASELINE_ID", format="I", array=event_list[:, 0]), + fits.Column(name="TRUE_ENERGY", format="E", array=event_list[:, 1]), + fits.Column(name="DETECTED_ENERGY", format="E", array=event_list[:, 2]), + fits.Column(name="TIME", format="E", array=event_list[:, 3]), + fits.Column(name="RA_POS", format="D", array=ra_pos), + fits.Column(name="DEC_POS", format="D", array=dec_pos), ] hdu = fits.BinTableHDU.from_columns(cols) hdu.header.update( { - 'EXTNAME': 'EVENTS', - 'TELESCOP': 'XRI', - 'INSTRUME': 'CIS', - 'TTYPE1': 'BASELINE_ID', - 'TTYPE2': 'TRUE_ENERGY', - 'TTYPE3': 'DETECTED_ENERGY', - 'TTYPE4': 'TIME', - 'TTYPE5': 'RA_POS', - 'TTYPE6': 'DEC_POS', - 'TUNIT1': '', - 'TUNIT2': 'keV', - 'TUNIT3': 'keV', - 'TUNIT4': 's', - 'TUNIT5': 'rad', - 'TUNIT6': 'rad', - } + "EXTNAME": "EVENTS", + "TELESCOP": "XRI", + "INSTRUME": "CIS", + "TTYPE1": "BASELINE_ID", + "TTYPE2": "TRUE_ENERGY", + "TTYPE3": "DETECTED_ENERGY", + "TTYPE4": "TIME", + "TTYPE5": "RA_POS", + "TTYPE6": "DEC_POS", + "TUNIT1": "", + "TUNIT2": "keV", + "TUNIT3": "keV", + "TUNIT4": "s", + "TUNIT5": "rad", + "TUNIT6": "rad", + } ) hdul = fits.HDUList([primary_hdu, hdu]) hdul.writeto(output_file, overwrite=True) - hdul.close() \ No newline at end of file + hdul.close() From 0916a0a30864bbf9388e49dc5f6e2fa99e221648 Mon Sep 17 00:00:00 2001 From: postoecker Date: Wed, 8 Apr 2026 17:22:03 +0200 Subject: [PATCH 5/5] Change way of git-ignoring files --- .gitignore | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/.gitignore b/.gitignore index e2c317a..26851fe 100644 --- a/.gitignore +++ b/.gitignore @@ -1,8 +1,4 @@ -__pycache__/ - -.DS_Store -/.DS_Store - +*.pyc .vscode/ ignore/