diff --git a/.gitignore b/.gitignore new file mode 100644 index 00000000..e4e07ec8 --- /dev/null +++ b/.gitignore @@ -0,0 +1,3 @@ +*.pyc +*checkpoint* +*egg-info* diff --git a/Roman_TDS_obseq_prism.fits b/Roman_TDS_obseq_prism.fits new file mode 100644 index 00000000..e2b64d68 Binary files /dev/null and b/Roman_TDS_obseq_prism.fits differ diff --git a/config/prism.yaml b/config/prism.yaml new file mode 100644 index 00000000..e325547f --- /dev/null +++ b/config/prism.yaml @@ -0,0 +1,173 @@ +# Default settings for roman simulation +# Includes creation of noisless oversampled images (including PSF) +# -- processing of other detector and instrument effects are still handled in the +# python postprocessing layer to enable things not currently in galsim.roman + +modules: + + # Including galsim.roman in the list of modules to import will add a number of Roman-specific + # functions and classes that we will use here. + - roman_imsim + - galsim.roman + + # We need this for one of our Eval items. GalSim does not by default import datetime into + # the globals dict it uses when evaluating Eval items, so we can tell it to import it here. + - datetime + +# Define some other information about the images +image: + + # A special Image type that knows all the Roman SCA geometry, WCS, gain, etc. + # It also by default applies a number of detector effects, but these can be turned + # off if desired by setting some parameters (given below) to False. + type: roman_sca + + wcs: + type: RomanWCS + SCA: '@image.SCA' + ra: { type: ObSeqData, field: ra } + dec: { type: ObSeqData, field: dec } + pa: { type: ObSeqData, field: pa } + # mjd: { type: ObSeqData, field: mjd } + mjd: { type: Sum_float, items: [ { type: ObSeqData_float, field: mjd }, 170 ] } + + max_sun_angle: 50 + force_cvz: True + + bandpass: + type: RomanBandpass + name: { type: ObSeqData, field: filter } + + # When you want to have multiple images generate the same random galaxies, then + # you can set up multiple random number generators with different update cadences + # by making random_seed a list. + # The default behavior is just to have the random seeds for each object go in order by + # object number across all images, but this shows how to set it up so we use two separate + # cadences. + # The first one behaves normally, which will be used for things like noise on the image. + # The second one sets the initial seed for each object to repeat to the same starting value + # at the start of each filter. If we were doing more than 3 total files, it would then + # move on to another sequence for the next 3 and so on. + random_seed: + # Used for noise and nobjects. + - { type: ObSeqData, field: visit } + + # Used for objects. Repeats sequence for each filter + # Note: Don't use $ shorthand here, since that will implicitly be evaluated once and then + # treated the same way as an integer (i.e. making a regular sequence starting from that + # value). Using an explicit dict with an Eval type means GalSim will leave it alone and + # evaluate it as is for each object. + + + # We're just doing one SCA here. + # If you wanted to do all of them in each of three filters (given below), you could use: + # + # SCA: + # type: Sequence + # first: 1 + # last: 18 + # repeat: 3 # repeat each SCA num 3 times before moving on, for the 3 filters. + # + SCA: 4 + # mjd: { type: ObSeqData, field: mjd } + mjd: { type: Sum_float, items: [ { type: ObSeqData_float, field: mjd }, 170 ] } + + filter: { type: ObSeqData, field: filter } + exptime: { type: ObSeqData, field: exptime } + + # Photon shooting is way faster for chromatic objects than fft, especially when most of them + # are fairly faint. The cross-over point for achromatic objects is generally of order + # flux=1.e6 or so (depending on the profile). Most of our objects here are much fainter than + # that. The fft rendering for chromatic is a factor of 10 or so slower still, whereas + # chromatic photon shooting is only slighly slower than achromatic, so the difference + # is even more pronounced in this case. + draw_method: 'phot' + + # These are all by default turned on, but you can turn any of them off if desired: + ignore_noise: True + stray_light: False + thermal_background: False + reciprocity_failure: False + dark_current: False + nonlinearity: False + ipc: False + read_noise: False + sky_subtract: False + +stamp: + type: Roman_stamp + world_pos: + type: SkyCatWorldPos + exptime: { type: ObSeqData, field: exptime } + min_size: 1024 + flux_cap: 1000000 + skip_failures: True + + photon_ops: + - + type: WFSSSDisperser + #- + # type: ChargeDiff + +# psf: +# type: roman_psf +# # If omitted, it would figure this out automatically, because we are using the RomanSCA image +# # type. But if we weren't, you'd have to tell it which SCA to build the PSF for. +# SCA: '@image.SCA' +# # n_waves defines how finely to sample the PSF profile over the bandpass. +# # Using 10 wavelengths usually gives decent accuracy. +# n_waves: 10 + +# Define the galaxy type and positions to use +gal: + type: SkyCatObj + +input: + obseq_data: + #file_name: /hpc/group/cosmology/repos/prism_sim/Roman_TDS_obseq_prism.fits + file_name: /hpc/group/cosmology/OpenUniverse2026/RomanTDS/Roman_TDS_obseq_11_6_23.fits + visit: 4191 + SCA: '@image.SCA' + roman_psf: + SCA: '@image.SCA' + n_waves: 5 + sky_catalog: + file_name: /hpc/group/cosmology/OpenUniverse2026/roman_rubin_cats_v1.1.2_faint/skyCatalog.yaml + edge_pix: 512 + # mjd: { type: ObSeqData, field: mjd } + mjd: { type: Sum_float, items: [ { type: ObSeqData_float, field: mjd }, 170 ] } + + exptime: { type: ObSeqData, field: exptime } + obj_types: ["snana", "diffsky_galaxy"] # , "star" + +output: + + nfiles: 1 + dir: RomanTDS_prism/images/truth + file_name: + type: FormattedStr + format: "Roman_TDS_truth_test_%s_%i_%i.fits.gz" + items: + - { type: ObSeqData, field: filter } + - { type: ObSeqData, field: visit } + - '@image.SCA' + + truth: + dir: RomanTDS_prism/truth + file_name: + type: FormattedStr + format: "Roman_TDS_index_test_%s_%i_%i.txt" + items: + - { type: ObSeqData, field: filter } + - { type: ObSeqData, field: visit } + - '@image.SCA' + columns: + object_id: "@object_id" + ra: "$sky_pos.ra.deg" + dec: "$sky_pos.dec.deg" + x: "$image_pos.x" + y: "$image_pos.y" + realized_flux: "@realized_flux" + flux: "@flux" + mag: "@mag" + obj_type: "@object_type" diff --git a/config/roman_cutouts.yaml b/config/roman_cutouts.yaml new file mode 100644 index 00000000..895d9f26 --- /dev/null +++ b/config/roman_cutouts.yaml @@ -0,0 +1,181 @@ +# Default settings for roman simulation +# Includes creation of noisless oversampled images (including PSF) +# -- processing of other detector and instrument effects are still handled in the +# python postprocessing layer to enable things not currently in galsim.roman + +modules: + + # Including galsim.roman in the list of modules to import will add a number of Roman-specific + # functions and classes that we will use here. +- roman_imsim +- galsim.roman + + # We need this for one of our Eval items. GalSim does not by default import datetime into + # the globals dict it uses when evaluating Eval items, so we can tell it to import it here. +- datetime + +# Define some other information about the images +image: + + # A special Image type that knows all the Roman SCA geometry, WCS, gain, etc. + # It also by default applies a number of detector effects, but these can be turned + # off if desired by setting some parameters (given below) to False. + type: roman_sca + + wcs: + type: RomanWCS + SCA: '@image.SCA' + ra: {type: ObSeqData, field: ra} + dec: {type: ObSeqData, field: dec} + pa: {type: ObSeqData, field: pa} + mjd: &id001 + type: Sum_float + items: + - type: ObSeqData_float + field: mjd + - 105.76 + # mjd: { type: Sum_float, items: [ { type: ObSeqData_float, field: mjd }, 170 ] } + + max_sun_angle: 50 + force_cvz: true + + bandpass: + type: RomanBandpass + name: {type: ObSeqData, field: filter} + + # When you want to have multiple images generate the same random galaxies, then + # you can set up multiple random number generators with different update cadences + # by making random_seed a list. + # The default behavior is just to have the random seeds for each object go in order by + # object number across all images, but this shows how to set it up so we use two separate + # cadences. + # The first one behaves normally, which will be used for things like noise on the image. + # The second one sets the initial seed for each object to repeat to the same starting value + # at the start of each filter. If we were doing more than 3 total files, it would then + # move on to another sequence for the next 3 and so on. + random_seed: + # Used for noise and nobjects. + - {type: ObSeqData, field: visit} + + # Used for objects. Repeats sequence for each filter + # Note: Don't use $ shorthand here, since that will implicitly be evaluated once and then + # treated the same way as an integer (i.e. making a regular sequence starting from that + # value). Using an explicit dict with an Eval type means GalSim will leave it alone and + # evaluate it as is for each object. + + + # We're just doing one SCA here. + # If you wanted to do all of them in each of three filters (given below), you could use: + # + # SCA: + # type: Sequence + # first: 1 + # last: 18 + # repeat: 3 # repeat each SCA num 3 times before moving on, for the 3 filters. + # + SCA: 5 + mjd: *id001 + # mjd: { type: Sum_float, items: [ { type: ObSeqData_float, field: mjd }, 170 ] } + + filter: {type: ObSeqData, field: filter} + exptime: {type: ObSeqData, field: exptime} + + # Photon shooting is way faster for chromatic objects than fft, especially when most of them + # are fairly faint. The cross-over point for achromatic objects is generally of order + # flux=1.e6 or so (depending on the profile). Most of our objects here are much fainter than + # that. The fft rendering for chromatic is a factor of 10 or so slower still, whereas + # chromatic photon shooting is only slighly slower than achromatic, so the difference + # is even more pronounced in this case. + draw_method: 'phot' + + # These are all by default turned on, but you can turn any of them off if desired: + ignore_noise: true + stray_light: false + thermal_background: false + reciprocity_failure: false + dark_current: false + nonlinearity: false + ipc: false + read_noise: false + sky_subtract: false + +stamp: + type: Roman_stamp + world_pos: + type: SkyCatWorldPos + exptime: {type: ObSeqData, field: exptime} + min_size: 1024 + flux_cap: 10000000000.0 + + photon_ops: + - type: SlitlessSpec + #- type: ChargeDiff + +# psf: +# type: roman_psf +# # If omitted, it would figure this out automatically, because we are using the RomanSCA image +# # type. But if we weren't, you'd have to tell it which SCA to build the PSF for. +# SCA: '@image.SCA' +# # n_waves defines how finely to sample the PSF profile over the bandpass. +# # Using 10 wavelengths usually gives decent accuracy. +# n_waves: 10 + +# Define the galaxy type and positions to use +gal: + type: SkyCatObj + +input: + obseq_data: + #file_name: /hpc/group/cosmology/repos/prism_sim/Roman_TDS_obseq_prism.fits + file_name: /hpc/group/cosmology/OpenUniverse2024/RomanTDS/Roman_TDS_obseq_11_6_23.fits + visit: 55726 + SCA: '@image.SCA' + roman_psf: + SCA: '@image.SCA' + n_waves: 5 + sky_catalog: + file_name: /hpc/group/cosmology/repos/roman_imsim/roman_rubin_cats_v1.1.2_faint/skyCatalog.yaml + edge_pix: 512 + mjd: *id001 + # mjd: { type: Sum_float, items: [ { type: ObSeqData_float, field: mjd }, 170 ] } + + exptime: {type: ObSeqData, field: exptime} + obj_types: ["snana"] # diffsky_galaxy, star + +label: + lab: "1*exp_p+15d" + + +output: + + nfiles: 1 + dir: /hpc/home/jr542/SNPIT_OU24_extract1d/exptime/truth + file_name: + type: FormattedStr + format: "Roman_TDS_truth_%s_%s_%i_%i.fits.gz" + items: + - '@label.lab' + - {type: ObSeqData, field: filter} + - {type: ObSeqData, field: visit} + - '@image.SCA' + + truth: + dir: /hpc/home/jr542/SNPIT_OU24_extract1d/exptime/index + file_name: + type: FormattedStr + format: "Roman_TDS_index_%s_%s_%i_%i.txt" + items: + - '@label.lab' + - {type: ObSeqData, field: filter} + - {type: ObSeqData, field: visit} + - '@image.SCA' + columns: + object_id: "@object_id" + ra: "$sky_pos.ra.deg" + dec: "$sky_pos.dec.deg" + x: "$image_pos.x" + y: "$image_pos.y" + realized_flux: "@realized_flux" + flux: "@flux" + mag: "@mag" + obj_type: "@object_type" diff --git a/det.py b/det.py new file mode 100644 index 00000000..e8f71ce5 --- /dev/null +++ b/det.py @@ -0,0 +1,21 @@ +import sys + +import roman_imsim + +params = sys.argv[1] +visit = int(sys.argv[2]) +sca = int(sys.argv[3]) +if len(sys.argv) > 4: + dither_from_file = sys.argv[4] +else: + import warnings + + warnings.warn("Ditherlist not provided. Persistence effect will not simulated.") + dither_from_file = None + +if len(sys.argv) > 5: + sca_path = sys.argv[5] +else: + sca_path = None + +roman_imsim.detector_physics.modify_image(params, visit, sca, dither_from_file, sca_filepath=sca_path) diff --git a/grism.yaml b/grism.yaml new file mode 100644 index 00000000..dcb93c97 --- /dev/null +++ b/grism.yaml @@ -0,0 +1,175 @@ +# Default settings for roman simulation +# Includes creation of noisless oversampled images (including PSF) +# -- processing of other detector and instrument effects are still handled in the +# python postprocessing layer to enable things not currently in galsim.roman + +modules: + + # Including galsim.roman in the list of modules to import will add a number of Roman-specific + # functions and classes that we will use here. + - roman_imsim + + # We need this for one of our Eval items. GalSim does not by default import datetime into + # the globals dict it uses when evaluating Eval items, so we can tell it to import it here. + - datetime + - galsim.roman + +# Define some other information about the images +image: + + # A special Image type that knows all the Roman SCA geometry, WCS, gain, etc. + # It also by default applies a number of detector effects, but these can be turned + # off if desired by setting some parameters (given below) to False. + type: roman_sca + + wcs: + type: RomanWCS + SCA: '@image.SCA' + ra: { type: ObSeqData, field: ra } + dec: { type: ObSeqData, field: dec } + pa: { type: ObSeqData, field: pa } + # mjd: { type: ObSeqData, field: mjd } + mjd: { type: Sum_float, items: [ { type: ObSeqData_float, field: mjd }, 170 ] } + + max_sun_angle: 50 + force_cvz: True + + bandpass: + type: RomanBandpass + name: { type: ObSeqData, field: filter } + + # When you want to have multiple images generate the same random galaxies, then + # you can set up multiple random number generators with different update cadences + # by making random_seed a list. + # The default behavior is just to have the random seeds for each object go in order by + # object number across all images, but this shows how to set it up so we use two separate + # cadences. + # The first one behaves normally, which will be used for things like noise on the image. + # The second one sets the initial seed for each object to repeat to the same starting value + # at the start of each filter. If we were doing more than 3 total files, it would then + # move on to another sequence for the next 3 and so on. + random_seed: + # Used for noise and nobjects. + - { type: ObSeqData, field: visit } + + # Used for objects. Repeats sequence for each filter + # Note: Don't use $ shorthand here, since that will implicitly be evaluated once and then + # treated the same way as an integer (i.e. making a regular sequence starting from that + # value). Using an explicit dict with an Eval type means GalSim will leave it alone and + # evaluate it as is for each object. + + + # We're just doing one SCA here. + # If you wanted to do all of them in each of three filters (given below), you could use: + # + # SCA: + # type: Sequence + # first: 1 + # last: 18 + # repeat: 3 # repeat each SCA num 3 times before moving on, for the 3 filters. + # + SCA: 16 + # mjd: { type: ObSeqData, field: mjd } + mjd: { type: Sum_float, items: [ { type: ObSeqData_float, field: mjd }, 170 ] } + + filter: { type: ObSeqData, field: filter } + exptime: { type: ObSeqData, field: exptime } + + # Photon shooting is way faster for chromatic objects than fft, especially when most of them + # are fairly faint. The cross-over point for achromatic objects is generally of order + # flux=1.e6 or so (depending on the profile). Most of our objects here are much fainter than + # that. The fft rendering for chromatic is a factor of 10 or so slower still, whereas + # chromatic photon shooting is only slighly slower than achromatic, so the difference + # is even more pronounced in this case. + draw_method: 'phot' + + # These are all by default turned on, but you can turn any of them off if desired: + ignore_noise: True + stray_light: False + thermal_background: False + reciprocity_failure: False + dark_current: False + nonlinearity: False + ipc: False + read_noise: False + sky_subtract: False + +stamp: + type: Roman_stamp + world_pos: + type: SkyCatWorldPos + exptime: { type: ObSeqData, field: exptime } + min_size: 1024 + flux_cap: 10000000 + + photon_ops: + - + type: GrismV + +# GrismNV calls the non-vector version of the code while GrismV calls the vector version + + #- + # type: ChargeDiff + +# psf: +# type: roman_psf +# # If omitted, it would figure this out automatically, because we are using the RomanSCA image +# # type. But if we weren't, you'd have to tell it which SCA to build the PSF for. +# SCA: '@image.SCA' +# # n_waves defines how finely to sample the PSF profile over the bandpass. +# # Using 10 wavelengths usually gives decent accuracy. +# n_waves: 10 + +# Define the galaxy type and positions to use +gal: + type: SkyCatObj + +input: + obseq_data: + #file_name: /hpc/group/cosmology/repos/prism_sim/Roman_TDS_obseq_prism.fits + file_name: /hpc/group/cosmology/OpenUniverse2024/RomanTDS/Roman_TDS_obseq_11_6_23.fits + visit: 14234 + SCA: '@image.SCA' + roman_psf: + SCA: '@image.SCA' + n_waves: 5 + sky_catalog: + file_name: /hpc/group/cosmology/repos/roman_imsim/roman_rubin_cats_v1.1.2_faint/skyCatalog.yaml + edge_pix: 512 + # mjd: { type: ObSeqData, field: mjd } + mjd: { type: Sum_float, items: [ { type: ObSeqData_float, field: mjd }, 170 ] } + + exptime: { type: ObSeqData, field: exptime } + obj_types: ["snana", "diffsky_galaxy"] #,snana, star + +output: + + nfiles: 1 + dir: ./output + file_name: + type: FormattedStr + format: "Roman_Grism_truth_%s_%i_%i.fits.gz" + items: + - { type: ObSeqData, field: filter } + - { type: ObSeqData, field: visit } + - '@image.SCA' + + truth: + dir: ./output + file_name: + type: FormattedStr + format: "Roman_Grism_index_test_%s_%i_%i.txt" + items: + - { type: ObSeqData, field: filter } + - { type: ObSeqData, field: visit } + - '@image.SCA' + columns: + object_id: "@object_id" + ra: "$sky_pos.ra.deg" + dec: "$sky_pos.dec.deg" + x: "$image_pos.x" + y: "$image_pos.y" + realized_flux: "@realized_flux" + flux: "@flux" + mag: "@mag" + obj_type: "@object_type" diff --git a/optical_models/README.md b/optical_models/README.md new file mode 100644 index 00000000..f0881607 --- /dev/null +++ b/optical_models/README.md @@ -0,0 +1,4 @@ +# Roman Optical Models + +This directory contains versions of various optical models from IPAC for Roman Space Telescope. +This is meant to be temporary and we should rely on centralized optical model for consistency. diff --git a/optical_models/Roman_grism_OpticalModel_v0.8.yaml b/optical_models/Roman_grism_OpticalModel_v0.8.yaml new file mode 100644 index 00000000..c84352c7 --- /dev/null +++ b/optical_models/Roman_grism_OpticalModel_v0.8.yaml @@ -0,0 +1,459 @@ +meta: + # Metadata information + product_type: wfi_calibration_optical_model + filename: Roman_grism_OpticalModel_v0.8.yaml + file_date: 2025-09-17Z00:00:00 + origin: IPAC/SSC + telescope: ROMAN + instrument: WFI + version_cdp_software: 4.0.0 + optical_element: grism + unit_wl: micron + context: null + calibration_filenames: + small_scale_flatfield: null + relative_flux_calibration: null + absolute_flux_calibration: null + spectral_psf: null + +detector_model: + # Detector specifics + naxis1: 4088 + naxis2: 4088 + crpix1: 2044.5 + crpix2: 2044.5 + pos_angle_detector: 0.0 # rotation angle for the SCAs + plate_scale: 100 # pixel/mm + pixel_scale: 0.11 # arcsec/px + + # SCA centers in mm + xy_centers: + 1: [-22.14, 12.15] + 2: [-22.29, -37.03] + 3: [-22.44, -82.06] + 4: [-66.42, 20.90] + 5: [-66.92, -28.28] + 6: [-67.42, -73.06] + 7: [-110.70, 42.20] + 8: [-111.48, -6.98] + 9: [-112.64, -51.06] + 10: [22.14, 12.15] + 11: [22.29, -37.03] + 12: [22.44, -82.06] + 13: [66.42, 20.90] + 14: [66.92, -28.28] + 15: [67.42, -73.06] + 16: [110.70, 42.20] + 17: [111.48, -6.98] + 18: [112.64, -51.06] + +optical_model: + # Wavelength specifics + wl_min: 0.9 # in micron + wl_max: 2.0 # in micron + wl_reference: 1.55 # in micron + wl_transform: linear + + # Coefficient matrix dimensions + dimen_map: + i: 6 # FPAx dimension + j: 6 # FPAy dimension + + dimen_crv: + i: 6 # Wavelength dimension [micron] + j: 6 # FPAx dimension + k: 6 # FPAy dimension + + dimen_ids: + i: 6 # Delta y distance along trace [mm] + j: 6 # FPAx dimension + k: 6 # FPAy dimension + + # List of orders defined + orders_defined: ["1", "0", "2", "m1"] + + orders: + # Positive 1st order (primary trace) + "1": + xmap_ij_coeff: + [ + [-3.5955e-01, +1.1443e-01, -2.4937e-02, -8.5102e-02, -1.8489e-01, -2.9182e-02], + [+3.2123e+02, -1.5175e+01, +1.4631e+01, -1.8579e+00, -2.7797e-01, -2.4955e+00], + [-7.6085e-02, -4.5501e-03, +4.4607e-01, +1.6437e+00, -6.3049e+00, -2.1155e+01], + [+1.3912e+01, -1.7340e+00, +1.8202e+00, -2.4279e+00, -1.7718e+00, -1.0506e+00], + [-1.7851e-02, -1.0088e-01, -2.3013e+00, -5.5019e+00, +4.4948e+01, +1.2027e+02], + [+5.0604e-01, -9.9092e-01, -6.1011e+00, -9.2447e+00, +3.7067e+01, +6.6016e+02], + ] + + ymap_ij_coeff: + [ + [-3.4090e-01, +3.2865e+02, -2.4468e+01, +1.5431e+01, -2.7334e+00, -6.6736e-01], + [-1.4446e-01, -3.8034e-02, -2.9770e-02, -6.2583e-01, -5.2444e-01, +2.5263e+00], + [-8.9654e+00, +1.4550e+01, -3.3497e+00, +1.7172e+00, +1.6131e+00, -1.2168e+01], + [-2.3530e-02, -5.8258e-01, +6.6854e-01, +2.4487e+01, -9.1861e+00, -2.2217e+02], + [-5.5495e-01, +6.6502e-01, +1.2796e+00, -4.6197e+00, -7.7578e+00, -4.4983e+01], + [+5.5202e-02, +3.5054e+00, -6.8692e+00, -1.6732e+02, +1.6477e+02, +1.7667e+03], + ] + + crv_ijk_coeff: + [ + [ + [-2.8983e-07, -4.3961e-05, -1.6987e-04, +2.3215e-03, +1.5122e-02, +2.3772e-02], + [-2.3874e-05, -5.8425e-05, +4.3961e-03, +8.2145e-03, -7.7760e-02, -1.7456e-01], + [+2.5767e-05, +1.3666e-03, -3.0734e-03, -7.9161e-02, -6.4914e-02, +3.7395e-01], + [+2.8161e-04, +5.9331e-04, -5.2914e-02, -4.1554e-02, +4.9381e-01, +7.8058e-02], + [-1.3348e-04, -6.6712e-03, +2.1911e-02, +2.7211e-01, -5.0212e-02, +6.7590e-02], + [-9.5253e-04, -6.3240e-03, +1.7150e-01, +3.0918e-01, +9.3130e-02, +2.3278e-02], + ], + [ + [-6.6192e-05, +3.0415e-03, +3.3949e-04, -1.4846e-02, -8.9554e-02, -1.5297e-01], + [-1.1129e-01, +2.1182e-01, +4.9519e-02, +2.2953e-03, -1.1773e-01, -3.0732e-01], + [-2.6219e-04, +1.5217e-03, +2.4976e-02, +1.8754e-01, +4.6758e-01, +2.1753e-02], + [-1.1368e-02, +8.4280e-03, -1.6475e-01, +1.4239e-02, +1.5778e+00, -1.7416e-01], + [+4.1521e-04, -1.0636e-02, -1.4610e-01, -3.1231e-01, +1.7437e-01, -2.0777e-03], + [-3.1875e-03, +5.0777e-02, +5.1666e-01, -9.1913e-01, +7.7557e-01, -1.5968e-01], + ], + [ + [+8.0929e-06, +5.8810e-06, -4.0069e-05, -1.0914e-03, -5.7713e-03, -8.6641e-03], + [+6.1084e-04, +1.9880e-04, -1.2675e-03, +2.0368e-03, +1.0144e-02, -2.8005e-02], + [+3.3427e-06, -2.0087e-05, +7.0851e-04, +1.8830e-02, +8.1666e-02, +9.4431e-02], + [-8.3102e-05, +2.5782e-03, +1.1838e-02, -1.0661e-01, -3.1770e-02, +1.2247e+00], + [-4.7582e-05, +4.5494e-04, +2.9519e-03, -1.0713e-01, -4.5005e-01, +9.8645e-02], + [+5.2798e-04, -9.2755e-03, -4.6692e-02, +2.5322e-01, -8.2722e-03, +1.7021e-01], + ], + [ + [+2.8632e-07, -1.0987e-06, -3.0053e-05, +1.0304e-04, +1.8449e-03, +4.1021e-03], + [+2.3910e-05, +1.6959e-05, -1.5548e-03, -2.1596e-03, +5.6122e-02, +1.5939e-01], + [-7.1627e-06, -3.4643e-05, +1.3837e-03, +2.0323e-03, -5.1040e-02, -1.3269e-01], + [-9.4097e-05, +4.1811e-04, +2.1901e-02, -2.5808e-02, -8.0046e-01, -1.3381e+00], + [+4.2874e-05, +3.1903e-05, -8.0804e-03, +8.5695e-03, +2.6375e-01, +1.6273e-01], + [+3.3832e-04, -3.8522e-03, -7.2512e-02, +3.7780e-01, +2.6601e+00, -8.3232e-01], + ], + [ + [-1.3554e-08, -4.7995e-07, -1.6579e-06, +8.1837e-05, +5.7305e-04, +9.9493e-04], + [-2.1527e-06, +3.0681e-06, +5.0829e-05, -4.2452e-05, -5.3410e-04, -1.2775e-04], + [-9.7028e-08, +1.3874e-05, +9.1668e-05, -1.5909e-03, -1.1434e-02, -1.7716e-02], + [+5.7755e-06, -1.5725e-04, -7.3418e-04, +5.0343e-03, +9.3652e-03, -2.6054e-02], + [+1.0689e-06, -8.8399e-05, -6.9269e-04, +7.6433e-03, +5.3774e-02, +6.8328e-02], + [-3.6498e-05, +5.6507e-04, +4.9745e-03, -3.5827e-03, -1.1187e-01, -4.4000e-01], + ], + [ + [+3.4015e-09, -1.4074e-08, +7.4324e-07, +9.2624e-06, +3.1904e-05, +3.1805e-05], + [+2.1924e-07, +8.1510e-07, +4.3347e-05, +5.2706e-05, -1.4538e-03, -4.0538e-03], + [+1.6859e-07, +1.8654e-06, -2.9319e-05, -2.5661e-04, -1.7782e-04, +8.0380e-04], + [+2.6942e-06, -2.6872e-05, -6.3923e-04, +1.0669e-03, +2.1559e-02, +3.4569e-02], + [-8.2686e-07, -5.3586e-06, +1.5423e-04, +6.0935e-04, -2.2742e-04, +7.6023e-03], + [-1.2058e-05, +1.5003e-04, +2.5011e-03, -8.2466e-03, -8.4447e-02, -7.5705e-02], + ], + ] + + ids_ijk_coeff: + [ + [ + [-6.3925e-06, -9.7160e-05, +8.2907e-05, +1.7802e-02, +1.2679e-01, +2.4713e-01], + [+1.0625e-05, +1.4368e-04, -1.0207e-03, -8.5457e-03, -6.4935e-03, -7.0560e-05], + [+7.9235e-05, +1.9033e-03, -1.9018e-02, -3.5247e-01, -1.4999e+00, -2.1349e+00], + [-2.0575e-04, -2.2302e-03, +2.3381e-02, +1.4254e-01, -9.9474e-02, -4.7824e-01], + [-2.6874e-04, -5.2756e-03, +8.1457e-02, +1.0593e+00, +5.7024e+00, +1.7843e+01], + [+9.5070e-04, +6.8763e-03, -9.9972e-02, -3.4871e-01, +3.3761e-01, -1.9832e+00], + ], + [ + [+9.1999e+00, -1.9901e+00, +3.3368e+00, -8.4777e-01, +8.0722e-02, -1.5785e-01], + [-9.8532e-04, -1.3339e-02, +1.4356e-01, +1.3866e+00, +1.1358e+00, -5.9559e+00], + [+1.1123e+00, -3.5145e-01, +1.2331e+00, +2.5528e-01, +1.8841e+00, +1.6802e+01], + [+1.6315e-02, +3.0771e-01, -2.7953e+00, -2.9614e+01, +7.7594e+00, +2.0308e+02], + [+6.4883e-02, -4.3905e-01, -2.5857e+00, -1.6075e+00, -2.2880e+00, +1.5401e+02], + [-6.8945e-02, -1.3836e+00, +1.3587e+01, +1.0504e+02, -1.4435e+02, -2.8115e+02], + ], + [ + [+4.0035e-02, +1.0811e-01, -4.4641e-02, -7.1667e-01, -4.3675e+00, -7.7473e+00], + [+1.2035e-04, -5.3304e-04, +6.2085e-02, +1.2209e+00, +8.3043e+00, +1.7183e+01], + [+4.2420e-04, +3.7886e-02, +1.5412e+00, +1.3653e+01, +3.8006e+01, +4.1362e+01], + [+1.5998e-03, -1.3844e-02, -1.2737e+00, -1.9112e+01, -1.2611e+02, -2.2540e+02], + [+1.7101e-02, -1.5165e-01, -9.0228e+00, -4.1354e+01, -4.6212e+01, -4.2114e+02], + [-1.3909e-02, -1.9935e-01, +4.5742e+00, +1.0648e+02, +5.2649e+02, -6.2036e+01], + ], + [ + [+3.9985e-02, +1.1103e-02, +7.1275e-02, -6.8424e-01, -8.1705e+00, -1.9454e+01], + [+7.3339e-04, +7.7768e-03, -1.5104e-02, -1.8456e-01, +1.0438e+00, +4.8762e+00], + [+8.6514e-03, -1.0075e-01, +2.2393e-01, +1.8858e+01, +9.8290e+01, +1.3191e+02], + [-4.3646e-03, -1.8028e-01, -1.1193e+00, +2.5259e+00, +2.1336e+01, -3.8180e+00], + [+1.8863e-03, +3.5264e-01, -4.6551e+00, -6.9272e+01, -2.5476e+02, -7.4827e+02], + [-1.1156e-02, +7.6916e-01, +1.1816e+01, -8.0201e+00, -3.2208e+02, -1.1884e+02], + ], + [ + [-3.6887e-02, -1.8676e-02, +4.3320e-03, +3.1190e+00, +2.2654e+01, +4.5074e+01], + [+3.5402e-04, +1.5770e-02, -7.3003e-02, -1.1827e+00, -2.4068e+00, +5.1967e-01], + [+6.9193e-03, +1.8285e-01, -3.5418e+00, -5.8855e+01, -2.5867e+02, -3.6750e+02], + [-1.3362e-02, -4.7636e-01, +1.2141e+00, +3.0254e+01, +6.9357e+01, +1.3322e+01], + [-6.6571e-02, -1.9827e-01, +1.7886e+01, +1.8154e+02, +8.9364e+02, +2.7645e+03], + [+9.5130e-02, +3.0355e+00, -9.8344e+00, -1.7895e+02, -2.0621e+02, +7.2850e+01], + ], + [ + [+7.3168e-02, -6.5848e-03, -4.5981e-02, +3.9894e+00, +3.4524e+01, +7.4117e+01], + [-1.4030e-03, -2.9034e-02, -1.3007e-01, +1.6208e+00, +7.7913e+00, +3.6682e+00], + [+1.2761e-02, +3.2559e-01, -3.3843e+00, -8.6182e+01, -4.0307e+02, -5.8132e+02], + [+1.2267e-02, +6.3636e-01, +3.1934e+00, -2.9877e+01, -1.1891e+02, -1.3013e+01], + [+1.5592e-02, -7.6510e-01, +2.3629e+01, +2.6913e+02, +1.2256e+03, +4.5907e+03], + [-7.2257e-03, -2.7877e+00, -1.7217e+01, +1.3011e+02, +4.5467e+02, -9.5989e+01], + ], + ] + + # Zeroth (0th) order + "0": + xmap_ij_coeff: + [ + [-3.5520e-01, +7.3441e-02, -2.2233e-02, -9.5321e-02, -1.8400e-01, -1.4418e-02], + [+3.2303e+02, -1.8450e+01, +1.5405e+01, -2.6627e+00, -2.3285e-01, -1.8705e+00], + [-7.1961e-02, -2.4415e-02, +4.9356e-01, +1.8981e+00, -7.1596e+00, -2.3478e+01], + [+1.4341e+01, -2.4617e+00, +2.1654e+00, +3.7846e-01, -3.9609e+00, -3.4064e+01], + [-1.6296e-02, -9.1369e-02, -2.5752e+00, -5.7962e+00, +5.0331e+01, +1.1362e+02], + [+5.3915e-01, -8.0323e-01, -7.2765e+00, -2.5256e+01, +5.5508e+01, +8.6663e+02], + ] + + ymap_ij_coeff: + [ + [-1.4090e+01, +3.3204e+02, -2.9538e+01, +1.6525e+01, -3.5604e+00, -2.9662e-01], + [-1.4310e-01, -3.2706e-02, -4.2031e-02, -5.9521e-01, -1.3797e-01, +3.5097e+00], + [-1.0607e+01, +1.5292e+01, -4.2203e+00, +2.5213e+00, +1.0550e+00, -1.9619e+01], + [-2.2612e-02, -5.6553e-01, +4.6729e-01, +2.3205e+01, -1.0729e+01, -2.2343e+02], + [-6.5972e-01, +7.8540e-01, +9.1862e-01, -8.3945e+00, -4.2815e+00, +8.2066e+00], + [+3.8978e-02, +3.4435e+00, -5.4159e+00, -1.6221e+02, +1.6114e+02, +1.7698e+03], + ] + + crv_ijk_coeff: + [ + [ + [-7.3386e-06, -1.0079e-04, +8.9281e-04, +6.2359e-03, +2.5594e-03, -2.8848e-02], + [-1.7630e-04, +2.7986e-04, +1.0925e-02, +2.3358e-02, -1.0247e-01, -3.3956e-01], + [+1.6194e-04, +1.3100e-03, -1.3742e-02, -7.9290e-02, -7.4830e-02, -3.4602e-01], + [-4.7624e-04, -1.1874e-02, -2.1137e-01, -4.5601e-01, +5.2284e+00, +1.6553e+01], + [-9.3884e-04, -3.4455e-03, +7.5132e-02, +5.9108e-03, -3.1793e-01, +6.8104e+00], + [+4.2586e-03, +6.3055e-02, +9.0280e-01, +1.1299e+00, -2.6513e+01, -5.7607e+01], + ], + [ + [+4.7465e-03, +8.7850e-03, -1.8075e-02, -2.3208e-01, -1.8851e+00, -4.1732e+00], + [+2.5570e-01, +2.5569e-01, -1.0518e-01, -3.1435e-01, +5.2373e+00, +1.9984e+01], + [-1.5207e-02, -3.9681e-02, +7.2673e-01, +4.2899e+00, +1.4575e+00, -1.9713e+01], + [+1.5670e-01, -7.3226e-01, -2.8970e+00, +3.2640e+00, -6.2603e+01, -1.5231e+02], + [+3.1569e-02, -9.6844e-02, -4.2057e+00, +1.5860e+00, +6.2597e+01, -1.4220e+02], + [-5.4343e-01, +2.6069e+00, +7.9851e+00, +5.9713e+01, +3.4599e+02, -1.5260e+03], + ], + [ + [-2.2767e-03, -3.6759e-02, -1.8103e-01, -1.8331e+00, +1.6644e+00, +1.7588e+01], + [-8.4095e-02, -5.4924e-01, -3.8454e-01, +2.3890e+00, -3.2930e+01, -1.3356e+02], + [+9.8624e-02, +4.1427e-01, +5.8715e+00, +2.5294e+01, -7.2680e+01, -2.8867e+02], + [-1.1080e+00, +6.0832e+00, +2.2163e+01, -2.8486e+01, +2.5792e+02, +1.0001e+03], + [+3.5028e-02, -1.6829e+00, -2.7029e+01, -8.9128e+01, +3.9215e+02, +1.5116e+03], + [+3.8173e+00, -2.2010e+01, -2.6147e+01, +1.0447e+02, -2.3991e+03, -3.6454e+03], + ], + [ + [-2.4411e-02, -2.8930e-01, -2.3783e-01, -8.5271e+00, +1.8695e+01, +1.3246e+02], + [-8.2459e-01, -3.7398e+00, +3.2354e+00, +2.5464e+01, -1.2083e+02, -5.4961e+02], + [+7.6081e-01, +3.2921e+00, +2.1971e+01, +3.5990e+01, -3.5930e+02, -5.0466e+02], + [-7.4251e+00, +3.5594e+01, +1.3748e+02, +5.4775e+01, +1.8509e+02, -7.6827e+02], + [-7.1156e-01, -1.2332e+01, -1.1907e+02, -1.2814e+02, +1.3659e+03, +6.5876e+03], + [+3.1842e+01, -1.1066e+02, -4.4214e+02, -1.0292e+03, -1.3765e+03, +6.9977e+02], + ], + [ + [+6.7967e-02, +8.6269e-01, +3.4211e+00, +5.4306e+01, +1.4050e+01, -3.7273e+02], + [+2.4014e+00, +1.4841e+01, +9.3472e+00, -8.6845e+01, +4.5770e+02, +2.2525e+03], + [-2.5956e+00, -7.8302e+00, -1.1137e+02, -7.6664e+02, +1.7730e+02, +3.6538e+03], + [+2.7021e+01, -1.8757e+02, -7.2716e+02, +1.0918e+03, +2.7526e+03, -9.2270e+02], + [+7.2039e-01, +3.8507e+01, +5.0634e+02, +1.9220e+03, -1.3327e+03, -2.1095e+02], + [-1.0472e+02, +6.9259e+02, +1.9500e+03, -4.0431e+03, +1.2635e+03, -7.6503e+02], + ], + [ + [+3.9330e-01, +4.8575e+00, +1.2434e+01, +2.3074e+02, -1.3276e+02, -2.2959e+03], + [+1.3769e+01, +6.7675e+01, -3.5283e+01, -5.4753e+02, +2.8233e+03, +1.3333e+04], + [-1.3518e+01, -6.5018e+01, -4.7573e+02, -1.3932e+03, +2.5362e+03, -1.0562e+03], + [+1.2687e+02, -8.1198e+02, -3.1009e+03, +2.5114e+03, +1.8043e+03, +1.4041e+03], + [+1.5642e+01, +3.1251e+02, +2.0783e+03, +8.3944e+02, +6.8818e+02, +4.1467e+02], + [-5.5724e+02, +2.8221e+03, +1.0148e+04, +1.1468e+02, +6.4796e+02, +2.1334e+02], + ], + ] + + ids_ijk_coeff: + [ + [ + [-1.3064e-05, -9.7361e-05, +3.1119e-04, +2.1007e-02, +1.4063e-01, +2.5404e-01], + [-5.0999e-05, +3.4579e-04, +3.1380e-03, +8.8606e-04, -1.8277e-02, -3.0460e-02], + [-2.2121e-04, -3.1176e-03, +9.3511e-03, -1.3926e-01, -1.6954e+00, -3.2100e+00], + [+9.8104e-04, -4.7541e-03, -5.1133e-02, +4.2443e-03, +5.4046e-01, +9.5622e-01], + [+2.6340e-03, +2.4607e-02, -2.7267e-01, +6.6475e-02, +1.2092e+01, +1.9308e+01], + [-3.9074e-03, +2.0177e-02, +1.9971e-01, -2.2799e-01, -2.8977e+00, +4.9856e-01], + ], + [ + [+3.7215e-01, +4.9370e-02, +2.1344e-01, -5.0526e-02, -1.0056e+00, -1.0113e+00], + [+3.5574e-03, -1.6246e-02, +6.4362e-02, +1.4436e+00, +3.0699e+00, -1.2795e+00], + [+1.1372e-01, +2.3184e-01, -4.7079e-01, -1.0818e+00, +6.4642e+00, +3.2093e+01], + [-1.0650e-02, +3.9109e-01, -2.2313e+00, -3.2551e+01, -1.5348e+01, +1.6778e+02], + [-2.7432e-01, -8.1782e-01, -8.2616e-01, +7.7465e+00, +1.9772e+01, -2.0529e+01], + [-1.6923e-02, -1.7836e+00, +1.2282e+01, +1.2227e+02, -6.4945e+01, -2.0414e+02], + ], + [ + [+6.5967e-02, -5.7215e-03, -3.7220e-01, -1.5237e+00, -2.9247e+00, -5.9351e+00], + [-1.2425e-02, +1.4588e-02, +1.9168e-01, +7.9929e-01, +2.7334e+00, +3.5777e+00], + [-1.7632e-01, -4.5266e-01, +4.9436e+00, +1.2927e+01, +2.6970e+01, +3.1729e+01], + [+6.5125e-02, -2.7226e-01, -1.5073e+00, -3.7272e+00, -6.7359e+01, -1.8009e+02], + [+1.0776e+00, +4.0970e-01, -4.6146e+00, -3.6041e+01, -2.7308e+02, -3.0943e+02], + [-5.6141e-02, +1.2347e+00, +2.6271e+00, +1.8489e+01, +3.5967e+02, +7.5026e+01], + ], + [ + [+3.8933e-02, -2.7621e-02, -7.9504e-01, -2.4986e+00, -3.9138e+00, -1.2591e+01], + [-3.2151e-02, +5.4935e-02, +5.4891e-01, -4.8873e-01, -1.3911e+01, -3.4471e+01], + [-4.2452e-01, -1.3282e+00, +1.0641e+01, +2.6529e+01, +7.7372e+01, +9.0012e+01], + [+2.3681e-01, -8.5138e-01, -5.9379e+00, +2.8100e+01, +2.0781e+02, +2.2495e+02], + [+2.5370e+00, +3.5355e+00, -1.1463e+01, -1.3511e+02, -8.3719e+02, +8.5461e+02], + [-5.1866e-01, +4.3436e+00, +2.3831e+01, -1.7733e+02, -9.4052e+02, -1.3900e+02], + ], + [ + [-3.8229e-02, +4.5211e-02, +1.5333e+00, +6.4962e+00, +1.6813e+01, +3.6059e+01], + [+5.2272e-02, -5.2203e-02, -6.3523e-01, +8.5475e-01, +2.2830e+01, +6.0676e+01], + [+7.5897e-01, +2.0197e+00, -1.9042e+01, -5.8467e+01, -2.0377e+02, -2.9844e+02], + [-3.0035e-01, +5.9832e-01, +2.4586e+00, -3.7611e+01, -2.0114e+02, -1.7869e+02], + [-4.3581e+00, -2.7281e+00, +4.1824e+00, +1.7505e+02, +1.9497e+03, +1.7820e+03], + [+4.0417e-01, -2.7749e+00, -3.2334e+00, +1.9761e+02, +5.9975e+02, -5.7887e+02], + ], + [ + [+7.1583e-02, +1.0031e-01, +2.5551e+00, +9.7458e+00, +2.3241e+01, +5.4486e+01], + [+9.0885e-02, -1.5895e-01, -1.4471e+00, +4.0181e+00, +5.3674e+01, +1.1915e+02], + [+1.2751e+00, +3.5406e+00, -3.2653e+01, -9.3573e+01, -3.3065e+02, -4.8962e+02], + [-5.9982e-01, +2.5574e+00, +1.1054e+01, -1.3041e+02, -6.4691e+02, -5.4567e+02], + [-7.2819e+00, -7.3897e+00, +1.7795e+01, +3.5934e+02, +3.3242e+03, +6.8969e+02], + [+1.0708e+00, -1.2955e+01, -2.5968e+01, +7.1060e+02, +2.0744e+03, -3.2075e+02], + ], + ] + + # Postive 2nd order + "2": + xmap_ij_coeff: + [ + [-3.6212e-01, +1.5504e-01, -2.8054e-02, -8.8852e-02, -2.1216e-01, -3.5587e-02], + [+3.1957e+02, -1.1993e+01, +1.3986e+01, -1.1167e+00, -8.5828e-01, -3.9841e+00], + [-8.0971e-02, +2.8209e-03, +4.6223e-01, +2.0840e+00, -6.1005e+00, -2.4241e+01], + [+1.3530e+01, -1.0477e+00, +1.3025e+00, -3.7049e+00, +9.3403e+00, +3.7599e+01], + [-1.8458e-02, -5.1678e-02, -2.5233e+00, -8.2864e+00, +4.6667e+01, +1.4850e+02], + [+4.3689e-01, -9.4761e-01, -4.8084e+00, -5.5006e+00, -1.7311e+01, +5.2262e+02], + ] + + ymap_ij_coeff: + [ + [+1.3278e+01, +3.2569e+02, -1.9579e+01, +1.4527e+01, -1.9943e+00, -7.7728e-01], + [-1.4572e-01, -4.3308e-02, -2.1206e-02, -6.5907e-01, -5.3671e-01, +2.9597e+00], + [-7.3654e+00, +1.3929e+01, -2.5413e+00, +1.8943e+00, +1.6638e+00, -1.8248e+01], + [-2.3441e-02, -5.5629e-01, +6.3978e-01, +2.3424e+01, -8.0034e+00, -2.1166e+02], + [-4.9859e-01, +5.7841e-01, +1.2921e+00, -6.0962e+00, -5.9977e+00, -3.0737e+01], + [+5.1479e-02, +3.3069e+00, -6.4168e+00, -1.5848e+02, +1.5345e+02, +1.6672e+03], + ] + + crv_ijk_coeff: + [ + [ + [+4.1435e-07, -1.3716e-05, +2.3380e-04, +1.0853e-03, -1.3325e-02, -4.9786e-02], + [+2.1217e-04, +1.9283e-03, -3.4822e-03, -5.1086e-02, -8.5619e-02, +2.3055e-02], + [-8.8412e-05, +4.5193e-04, -1.3803e-03, -1.8496e-02, -1.2792e-02, +9.5453e-04], + [-9.6932e-04, -1.1183e-02, +8.4350e-03, -2.1116e-02, -8.4758e-03, +2.4158e-03], + [-8.0702e-04, -1.9372e-03, +8.4695e-03, -2.7328e-03, -1.4324e-03, +1.2814e-04], + [-1.2719e-04, +1.1095e-02, -2.4845e-03, -5.2771e-03, -8.5054e-04, +1.2183e-04], + ], + [ + [-3.6176e-05, +2.9045e-03, +2.5177e-04, -2.7162e-03, -3.2398e-02, -7.5714e-02], + [-1.1017e-01, +2.0323e-01, +2.2782e-02, +4.4028e-02, +7.1813e-02, -2.4047e-01], + [-4.3320e-04, +1.8383e-03, +2.8793e-03, +2.6565e-02, +1.6056e-01, -4.0465e-02], + [-9.8574e-03, +7.1573e-02, +4.3671e-02, -3.8603e-01, +1.5634e-03, -3.8651e-02], + [-2.4864e-04, -5.8789e-03, -5.8811e-03, +8.5717e-02, +1.0493e-02, +1.0117e-03], + [-1.6044e-02, -9.3692e-02, -1.2313e-01, -6.4154e-02, -2.8140e-03, -5.6315e-03], + ], + [ + [+5.9873e-06, +8.7806e-06, +1.6216e-05, +2.3300e-05, +3.5950e-03, +9.6996e-03], + [+5.1618e-04, +1.5028e-03, +2.3491e-03, -7.5162e-03, +2.7294e-02, +1.0282e-01], + [+7.4011e-05, +1.6996e-04, +3.8860e-05, -1.0396e-02, -2.1593e-02, +3.3738e-02], + [-5.2302e-04, -7.9033e-03, +8.2171e-03, +4.7706e-02, -9.7592e-02, +5.9934e-02], + [+8.1867e-05, -1.2329e-03, +3.5153e-03, +5.5391e-02, -1.3196e-02, +7.9941e-03], + [+4.3481e-03, +1.0712e-02, -2.3598e-02, -4.6113e-02, -1.5580e-02, +5.5662e-03], + ], + [ + [+1.5498e-07, +2.1878e-06, -1.8296e-06, -7.7360e-05, +2.4905e-04, +1.1857e-03], + [+1.8552e-05, +2.2354e-04, +4.6513e-04, -1.7588e-03, +4.0971e-03, +2.6903e-02], + [+4.1476e-06, -1.2576e-05, +4.0971e-04, +7.0934e-04, -7.2637e-03, -9.6213e-03], + [-6.6556e-05, -1.7121e-03, -1.5522e-03, +8.3650e-03, -5.4329e-02, -1.2900e-01], + [+1.7134e-05, -4.8951e-05, -2.5510e-03, +2.0099e-03, +3.8959e-02, -8.4060e-03], + [+4.3914e-04, +3.4553e-03, +5.7293e-03, +2.5662e-02, +1.5022e-01, -5.1056e-02], + ], + [ + [-6.9061e-09, -1.1045e-07, -5.7526e-07, -8.4626e-07, -3.6872e-05, -1.0308e-04], + [-1.5936e-06, -1.8127e-05, -3.4389e-05, +1.0623e-04, -3.5169e-04, -1.6195e-03], + [-8.6857e-07, -1.3930e-06, +1.0299e-05, +1.1958e-04, -2.4408e-05, -9.5848e-04], + [+7.6799e-06, +1.1300e-04, -7.5791e-05, -7.6781e-04, +8.3225e-04, -9.0751e-05], + [-6.8142e-07, +9.7804e-06, -7.5605e-05, -5.9888e-04, +1.1936e-03, +7.2589e-04], + [-5.5628e-05, -1.9483e-04, +5.2015e-04, +2.6555e-03, -2.2860e-03, -4.9127e-02], + ], + [ + [-7.0830e-10, -1.4534e-08, -6.3238e-09, +3.1382e-07, -3.2066e-06, -1.1361e-05], + [-1.3108e-07, -1.9209e-06, -3.9144e-06, +1.3258e-05, -4.5307e-05, -2.4276e-04], + [-5.5699e-08, -3.4095e-08, -1.7706e-06, +2.4273e-06, +4.0564e-05, +4.9691e-06], + [+7.1234e-07, +1.3994e-05, +5.3220e-06, -7.1044e-05, +3.8603e-04, +8.4065e-04], + [-1.0124e-07, +8.6381e-07, +1.0985e-05, -4.6416e-05, -1.8453e-04, +1.4204e-04], + [-4.6310e-06, -2.8568e-05, +1.2709e-06, +2.3089e-05, -1.1129e-03, -4.6298e-03], + ], + ] + + ids_ijk_coeff: + [ + [ + [+3.3777e-06, -4.4216e-05, +2.5162e-04, +1.9424e-02, +1.4459e-01, +2.8167e-01], + [-1.2587e-05, +6.8573e-04, -1.2376e-03, -2.9407e-02, -1.8562e-03, +2.0711e-01], + [+1.4350e-04, +4.0259e-03, -2.4742e-02, -5.0561e-01, -1.7310e+00, -5.3369e-01], + [+1.9495e-04, -9.6793e-03, +1.4743e-02, +3.9111e-01, +1.8723e-01, -2.1722e+00], + [-1.2638e-03, -2.0410e-02, +2.0554e-01, +2.7740e+00, +1.7953e+00, -2.9400e+01], + [-1.0615e-03, +3.7074e-02, +1.4324e-02, -1.6527e+00, -3.3593e+00, +1.0141e+01], + ], + [ + [+1.7871e+01, -3.4880e+00, +6.2794e+00, -1.4525e+00, +3.5964e-01, -5.6563e-01], + [+7.0072e-05, -3.6367e-02, +8.0683e-02, +1.6620e+00, +4.5483e+00, +1.6336e+00], + [+2.1077e+00, -6.9768e-01, +1.9629e+00, +2.4875e+00, +2.4013e+00, +5.4301e+00], + [+8.2320e-03, +6.3459e-01, -2.3703e+00, -4.0576e+01, -3.3835e+01, +1.7702e+02], + [-4.1920e-02, -4.6659e-01, -4.4942e+00, -1.1092e+01, +2.9122e+01, +2.3289e+02], + [-9.3844e-02, -2.9206e+00, +1.4394e+01, +1.6492e+02, -2.5192e+01, -3.7879e+02], + ], + [ + [-2.2726e-02, +3.5735e-01, +1.1799e-01, +3.8373e-01, -4.4672e+00, -7.8557e+00], + [-7.3456e-03, +6.0776e-02, +2.6013e-01, +7.4442e-01, -9.2491e-02, -3.1090e+00], + [-3.4911e-02, +3.0922e-01, +7.7822e-01, +7.6216e+00, +3.5721e+01, +6.4321e+01], + [+2.4851e-02, -8.8354e-01, -2.5343e+00, +1.0341e+01, -3.6663e+01, -2.2480e+02], + [+6.3623e-01, +3.1057e-01, -6.7217e+00, -3.9064e+01, +2.9996e+01, +4.8872e+02], + [+1.1844e-01, +4.0264e+00, +4.2475e-01, -5.6667e+01, +3.6240e+02, +5.3695e+02], + ], + [ + [+5.0594e-02, -9.2209e-03, +6.8448e-01, +1.5131e+00, -8.9771e+00, -2.1616e+01], + [-1.8120e-02, +1.8563e-01, +5.6472e-01, -2.0816e+00, -2.4025e+01, -5.2753e+01], + [-7.7902e-03, +4.2233e-01, -2.4718e+00, +4.7861e+00, +1.1595e+02, +1.8825e+02], + [+9.0113e-02, -2.8374e+00, -6.3628e+00, +8.5439e+01, +3.3107e+02, +1.4253e+02], + [+1.2395e+00, +8.0736e-01, +8.5332e+00, -6.5058e+01, -4.6727e+02, +4.8180e+02], + [+8.2476e-02, +1.2924e+01, +1.1242e+01, -4.3936e+02, -1.0589e+03, +9.4410e+01], + ], + [ + [-4.3383e-02, +2.9064e-02, -9.8046e-01, -8.3894e-01, +2.1465e+01, +4.1659e+01], + [+3.0792e-02, -2.5246e-01, -9.7884e-01, +1.9606e+00, +4.0783e+01, +9.8649e+01], + [+5.4324e-02, -5.7935e-01, -1.3091e+00, -4.1503e+01, -2.0588e+02, -2.1079e+02], + [-1.2402e-01, +3.3729e+00, +7.3810e+00, -1.0731e+02, -4.1841e+02, -2.0031e+02], + [-2.5055e+00, -2.2999e+00, +1.7759e+01, +2.4670e+02, +2.2920e+02, -3.1858e+03], + [-3.5026e-01, -1.4859e+01, +5.3903e+00, +5.2345e+02, +8.1674e+02, -8.9154e+02], + ], + [ + [+6.1234e-02, +5.8373e-02, -1.7673e+00, -2.6168e+00, +3.4973e+01, +7.5800e+01], + [+5.3046e-02, -5.2516e-01, -1.8456e+00, +7.1715e+00, +8.6794e+01, +1.8866e+02], + [+8.4724e-02, -1.0420e+00, +2.3613e+00, -5.5260e+01, -3.8170e+02, -4.1285e+02], + [-2.3881e-01, +8.0109e+00, +1.7665e+01, -2.7635e+02, -1.0688e+03, -5.6093e+02], + [-3.8423e+00, -3.3850e+00, +8.2534e+00, +3.8949e+02, +8.5733e+02, -4.3728e+03], + [-4.3678e-01, -3.6539e+01, -6.0470e+00, +1.3814e+03, +2.5665e+03, -1.9364e+02], + ], + ] + + # Negative 1st order not defined yet + "m1": + xmap_ij_coeff: null + ymap_ij_coeff: null + crv_ijk_coeff: null + ids_ijk_coeff: null diff --git a/optical_models/Roman_prism_OpticalModel_v0.8.yaml b/optical_models/Roman_prism_OpticalModel_v0.8.yaml new file mode 100644 index 00000000..f21050c9 --- /dev/null +++ b/optical_models/Roman_prism_OpticalModel_v0.8.yaml @@ -0,0 +1,200 @@ +meta: + # Metadata information + product_type: wfi_calibration_optical_model + filename: Roman_prism_OpticalModel_v0.8.yaml + file_date: 2025-09-17Z00:00:00 + origin: IPAC/SSC + telescope: ROMAN + instrument: WFI + version_cdp_software: 4.0.0 + optical_element: prism + unit_wl: micron + context: null + calibration_filenames: + small_scale_flatfield: null + relative_flux_calibration: null + absolute_flux_calibration: null + spectral_psf: null + +detector_model: + # Detector specifics + naxis1: 4088 + naxis2: 4088 + crpix1: 2044.5 + crpix2: 2044.5 + pos_angle_detector: 0.0 # rotation angle for the SCAs + plate_scale: 100 # pixel/mm + pixel_scale: 0.11 # arcsec/px + + # SCA centers in mm + xy_centers: + 1: [-22.14, 12.15] + 2: [-22.29, -37.03] + 3: [-22.44, -82.06] + 4: [-66.42, 20.90] + 5: [-66.92, -28.28] + 6: [-67.42, -73.06] + 7: [-110.70, 42.20] + 8: [-111.48, -6.98] + 9: [-112.64, -51.06] + 10: [22.14, 12.15] + 11: [22.29, -37.03] + 12: [22.44, -82.06] + 13: [66.42, 20.90] + 14: [66.92, -28.28] + 15: [67.42, -73.06] + 16: [110.70, 42.20] + 17: [111.48, -6.98] + 18: [112.64, -51.06] + +optical_model: + # Wavelength specifics + wl_min: 0.75 # in micron + wl_max: 1.85 # in micron + wl_reference: 1.20 # in micron + wl_transform: log + + # Coefficient matrix dimensions + dimen_map: + i: 6 # FPAx dimension + j: 6 # FPAy dimension + + dimen_crv: + i: 6 # Wavelength dimension [micron] + j: 6 # FPAx dimension + k: 6 # FPAy dimension + + dimen_ids: + i: 6 # Delta y distance along trace [mm] + j: 6 # FPAx dimension + k: 6 # FPAy dimension + + # List of orders defined (only one order for prism) + orders_defined: ["1"] + + orders: + # Primary (and only) order + "1": + xmap_ij_coeff: + [ + [-3.2586e-01, +1.0527e-01, +1.3745e-02, -4.5530e-02, -8.4277e-02, -3.4747e-02], + [+3.2955e+02, -1.7272e+01, +1.9143e+01, -2.9991e+00, +4.1477e-01, -3.1791e+00], + [+4.0631e-03, -1.2502e-02, +2.2342e-01, +6.5607e-01, -3.1760e+00, -1.0301e+01], + [+1.8226e+01, -2.7794e+00, +2.7388e+00, -3.5692e+00, +1.5661e+01, +2.2579e+01], + [+1.5243e-02, -1.1852e-01, -6.6652e-01, -1.3321e+00, +1.7271e+01, +4.4172e+01], + [+1.5206e+00, -1.2689e+00, +8.8337e-01, -4.0813e+00, -6.8761e+01, +4.5620e+02], + ] + + ymap_ij_coeff: + [ + [-1.1839e-01, +3.3706e+02, -2.7704e+01, +1.9895e+01, -4.6543e+00, +2.4596e-01], + [-1.4475e-01, -4.4152e-04, -6.8843e-02, -1.8304e-01, +1.2799e+00, +4.2355e+00], + [-9.7035e+00, +1.8933e+01, -5.1829e+00, +4.1657e+00, +4.5831e+00, -2.0464e+00], + [-2.6165e-02, -2.7830e-01, +5.4473e-01, +1.1321e+01, -2.8367e+01, -1.6913e+02], + [-8.9614e-01, +1.7092e+00, +4.0315e-01, -6.0727e+00, -2.7675e+01, -7.2198e+01], + [+2.8417e-02, +1.9525e+00, -4.3897e+00, -8.7990e+01, +1.8236e+02, +1.1727e+03], + ] + + crv_ijk_coeff: + [ + [ + [-2.2492e-07, +5.4770e-05, +3.7001e-04, -1.8985e-03, -1.7185e-02, -2.8976e-02], + [-2.9921e-05, +1.9336e-05, +2.1551e-03, -1.1691e-02, -8.6373e-02, -8.0362e-02], + [-2.4265e-05, -1.1441e-03, -6.3954e-03, +5.9434e-02, +3.1879e-01, +1.5978e-01], + [+7.3312e-04, -3.4451e-03, -3.7784e-02, +5.1387e-01, +1.0776e+00, -4.0425e+00], + [+2.6878e-04, +7.0884e-03, +1.1873e-02, -5.6176e-01, -9.2037e-01, +5.6974e+00], + [-4.0457e-03, +2.9910e-02, +1.6559e-01, -3.7062e+00, -2.9120e+00, +4.5872e+01], + ], + [ + [-1.8914e-03, +1.9878e-03, -1.4992e-03, +2.4885e-03, +9.2476e-02, +2.1250e-01], + [-1.0026e-01, +2.2022e-01, +3.2341e-02, -9.3255e-03, +5.4591e-01, +1.3763e+00], + [-6.1239e-04, +6.1684e-03, +8.0762e-02, -1.6898e-01, -2.9921e+00, -5.2911e+00], + [-1.1154e-03, +3.7798e-02, +5.7427e-03, -1.6874e+00, -3.8949e+00, +1.3505e+01], + [+3.4985e-04, -5.1929e-02, -3.0121e-01, +3.1225e+00, +1.1981e+01, -1.7165e+01], + [+7.0985e-04, -9.8062e-02, +5.2167e-01, +1.4337e+01, -1.7709e+00, -2.2775e+02], + ], + [ + [+8.1260e-05, +6.4680e-05, +3.2363e-03, +1.3289e-02, -4.4695e-02, -1.9088e-01], + [-2.1307e-03, +4.8240e-04, +9.0698e-03, +8.8183e-02, +1.4765e-01, -1.2009e-01], + [+1.1360e-03, +4.1586e-03, -1.4717e-01, -6.0979e-01, +3.3706e+00, +1.2955e+01], + [-3.9476e-04, +1.1064e-02, -2.7526e-01, -1.1884e+00, +9.1051e-01, +6.2427e+00], + [-8.3833e-03, -3.7955e-02, +1.0771e+00, +4.1456e+00, -2.7471e+01, -9.4093e+01], + [-1.3756e-02, -1.0118e-01, +2.1130e+00, +9.9602e+00, -2.7738e+01, -1.6666e+02], + ], + [ + [-2.0954e-05, +5.0599e-04, +4.3650e-03, -4.6504e-02, -4.1889e-01, -7.7276e-01], + [+4.1453e-04, -1.3484e-04, +4.7085e-02, -2.3266e-02, -2.3461e+00, -5.6852e+00], + [+2.7285e-05, -1.8833e-02, -9.2473e-02, +1.5945e+00, +7.4791e+00, +1.6571e+00], + [+3.6739e-04, -4.3507e-02, +1.4299e-02, +7.0125e+00, +1.3574e+01, -6.6835e+01], + [+1.8313e-03, +1.8364e-01, -6.2895e-02, -1.7024e+01, -1.6371e+01, +1.9872e+02], + [+1.0855e-02, +4.7122e-01, -3.3592e+00, -6.6422e+01, +5.1154e+01, +1.1641e+03], + ], + [ + [+2.7725e-05, -3.1968e-04, -1.0687e-02, -1.1746e-02, +3.2515e-01, +8.9093e-01], + [+5.0923e-04, +2.1853e-03, -3.0705e-02, -1.2211e-01, +9.0752e-01, +3.1652e+00], + [-2.6743e-03, +3.0639e-03, +4.5293e-01, +7.6048e-01, -1.3386e+01, -3.7137e+01], + [-1.9431e-03, -5.4864e-03, +3.4218e-01, -1.7856e+00, -9.5008e+00, +7.7666e+00], + [+1.8381e-02, -2.5721e-02, -2.9441e+00, -2.7184e+00, +8.6535e+01, +1.8191e+02], + [+1.7578e-02, +4.5981e-02, -1.8301e+00, +6.9647e+00, +4.3186e+01, +3.5997e+01], + ], + [ + [-1.8512e-05, -7.6897e-05, +4.6727e-03, +3.2415e-02, +7.3888e-03, -1.8027e-01], + [-2.1656e-04, -1.5651e-03, -9.6412e-03, +9.2283e-02, +7.8580e-01, +1.2256e+00], + [+1.6684e-03, +7.9484e-03, -2.5499e-01, -1.3498e+00, +5.1437e+00, +2.4511e+01], + [+8.8800e-04, +2.8966e-02, -2.1186e-01, -2.9303e+00, -1.1968e+00, +3.6124e+01], + [-1.2615e-02, -8.0622e-02, +2.0514e+00, +1.1069e+01, -5.1401e+01, -2.3810e+02], + [-1.5491e-02, -3.0223e-01, +3.1156e+00, +3.4412e+01, -6.2894e+01, -7.2035e+02], + ], + ] + + ids_ijk_coeff: + [ + [ + [+1.5564e-04, -6.5745e-05, +3.2838e-05, +3.8243e-03, +2.7012e-02, +4.9279e-02], + [-7.3713e-06, -4.9475e-04, -2.5037e-03, +2.0149e-02, +1.2379e-01, +1.3806e-01], + [+4.7461e-05, +8.1502e-04, +1.8379e-03, -8.7453e-02, -5.2722e-01, -6.1619e-01], + [+2.3577e-04, +1.1544e-02, +3.5900e-02, -5.1436e-01, -1.9827e+00, -6.9295e-02], + [-1.3724e-04, -4.7110e-03, -1.3676e-02, +4.7176e-01, +2.4883e+00, +8.1067e-01], + [-1.0368e-03, -5.0694e-02, -1.6609e-01, +2.1009e+00, +8.6584e+00, +2.3423e+00], + ], + [ + [-4.1727e+00, +9.6560e-01, -1.6491e+00, +4.3861e-01, -1.0176e-01, +4.6196e-01], + [-1.8321e-03, -2.8891e-03, -7.9245e-02, -1.7016e-01, +1.1998e-01, +4.6587e-01], + [-6.1427e-01, +1.9097e-01, -7.5287e-01, -1.0053e+00, +3.5263e+00, +1.4230e+01], + [-1.1697e-04, +1.8113e-01, +9.9313e-01, +3.0841e-01, -2.0555e+00, +4.8690e+00], + [-3.6793e-02, +1.5900e-01, +2.0225e+00, +3.4465e+00, -3.3091e+01, -8.4416e+01], + [-1.1886e-02, -1.0401e+00, -2.7914e+00, +1.9203e+01, -2.3239e+00, -3.3767e+02], + ], + [ + [+3.5222e+00, -8.1786e-01, +1.4733e+00, -6.7672e-01, -2.2251e+00, -6.6187e+00], + [+9.3949e-04, -4.6744e-03, +5.2799e-02, +7.3705e-01, +3.3325e+00, +2.5989e+00], + [+5.7645e-01, -1.0087e-01, +2.3791e+00, +5.5573e+00, -3.7615e+00, -4.9342e+01], + [+2.1899e-02, -4.0898e-02, +1.3705e+00, -8.2368e+00, -7.4531e+01, -3.6726e+01], + [+2.7903e-02, -1.0041e+00, -1.1137e+01, -3.9063e+00, +1.2384e+02, +3.1338e+02], + [-2.5200e-02, +6.3288e-01, -1.6661e+01, -2.5255e+00, +4.2378e+02, +6.7050e+02], + ], + [ + [-1.7165e+01, +3.9277e+00, -7.4015e+00, -1.9901e+00, -2.0274e+01, -4.1911e+01], + [-8.4953e-03, +3.6511e-01, +1.4948e+00, -3.8581e+00, -1.0680e-01, +5.6026e+01], + [-2.6207e+00, +2.6236e+00, +1.8021e+01, +3.7588e+01, -1.2790e+01, -3.4414e+02], + [-2.4752e-01, -5.4441e+00, +1.0060e+00, +3.5936e+01, -3.1049e+02, -3.2021e+02], + [+1.7659e-01, -1.1618e+01, -1.0486e+02, +1.3287e+01, +7.4610e+02, +1.8452e+03], + [+1.9264e+00, +1.9159e+01, -9.8241e+01, +2.9307e+01, +2.2839e+03, -6.2356e+02], + ], + [ + [+1.0028e+01, -2.1651e+00, +5.2931e+00, +1.3249e+01, +1.1627e+02, +2.8119e+02], + [+6.7981e-02, -3.4666e-02, +2.9204e-01, +7.0101e+00, -5.9556e+01, -2.7070e+02], + [+1.6322e+00, -6.1493e+00, -5.7260e+01, -1.5677e+02, -6.0060e+02, +8.2778e+01], + [-1.1247e+00, -3.3831e+00, -3.7952e+01, +4.2561e+01, +5.5881e+02, -1.0173e+02], + [+4.7853e-01, +4.6261e+01, +3.0608e+02, -3.5100e+02, -9.6277e+01, -1.9128e+02], + [+3.1343e+00, +1.4667e+01, +3.4251e+02, +6.2373e+01, -1.2111e+02, +6.5926e+01], + ], + [ + [-2.5269e+01, +6.3160e+00, +4.5407e+00, +1.1059e+02, +6.6219e+02, +1.5296e+03], + [+1.9577e-01, -8.0562e+00, -5.6109e+01, +6.0679e+01, +1.7432e+02, -1.4518e+03], + [-2.0819e+00, -5.3057e+01, -5.1399e+02, -9.6561e+02, -3.0934e+03, +8.5337e+02], + [+1.0189e+00, +1.1410e+02, +4.4388e+02, -1.4942e+01, -3.8014e+02, -5.4623e+01], + [-5.9752e+00, +3.8024e+02, +2.5169e+03, -2.0590e+03, -4.6367e+00, +4.4683e+01], + [-2.5441e+01, -4.3004e+02, +2.8415e+00, -5.1220e+01, +4.4200e+01, -3.2053e+01], + ], + ] \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 93950d8b..5924f6cc 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -16,6 +16,7 @@ dependencies = [ "galsim", "fitsio", "astropy", + "matplotlib", "numpy", "skycatalogs", "packaging", diff --git a/roman_imsim/__init__.py b/roman_imsim/__init__.py index 72da6e7a..e7ec4bc0 100644 --- a/roman_imsim/__init__.py +++ b/roman_imsim/__init__.py @@ -4,7 +4,6 @@ disable_implicit_threading() except ImportError: pass - from .bandpass import * from .detector_physics import * diff --git a/roman_imsim/bandpass.py b/roman_imsim/bandpass.py index 96472722..818db0c7 100644 --- a/roman_imsim/bandpass.py +++ b/roman_imsim/bandpass.py @@ -25,7 +25,7 @@ def buildBandpass(self, config, base, logger): kwargs, safe = GetAllParams(config, base, req=req) name = kwargs["name"] - bandpass = roman.getBandpasses(red_limit=2000)[name] + bandpass = roman.getBandpasses(red_limit=2000, include_all_bands=True)[name] return bandpass, safe diff --git a/roman_imsim/detector_physics.py b/roman_imsim/detector_physics.py index 7e58f8d4..ea6a8518 100644 --- a/roman_imsim/detector_physics.py +++ b/roman_imsim/detector_physics.py @@ -68,12 +68,12 @@ def __init__(self, params, visit, SCA): file_name = params["input"]["obseq_data"]["file_name"] obseq_data = ObSeqDataLoader(file_name, visit, SCA, logger=None) - self.filter = obseq_data.ob["filter"] + self.filter_ = obseq_data.ob["filter"] self.sca = obseq_data.ob["sca"] self.visit = obseq_data.ob["visit"] self.date = obseq_data.ob["date"] self.exptime = obseq_data.ob["exptime"] - self.bpass = roman.getBandpasses()[self.filter] + self.bpass = roman.getBandpasses(include_all_bands=True)[self.filter_] self.WCS = roman.getWCS( world_pos=galsim.CelestialCoord(ra=obseq_data.ob["ra"], dec=obseq_data.ob["dec"]), PA=obseq_data.ob["pa"], @@ -98,7 +98,7 @@ class modify_image(object): Class to simulate non-idealities and noise of roman detector images. """ - def __init__(self, params, visit, sca, dither_from_file, sca_filepath=None, use_galsim=False): + def __init__(self, params, visit, sca, dither_from_file=None, sca_filepath=None, use_galsim=False): """ Set up noise properties of image @@ -127,7 +127,7 @@ def __init__(self, params, visit, sca, dither_from_file, sca_filepath=None, use_ self.df = None print("------- Using simple detector model --------") - self.params["output"]["file_name"]["items"] = [self.pointing.filter, visit, sca] + self.params["output"]["file_name"]["items"] = [self.pointing.filter_, visit, sca] imfilename = ParseValue(self.params["output"], "file_name", self.params, str)[0] old_filename = os.path.join(self.params["output"]["dir"], imfilename) @@ -768,7 +768,7 @@ def setup_sky(self, im, pointing, rng, rng_iter, force_cvz=False): # This image is in units of e-/pix. Finally we add the expected thermal backgrounds in this # band. These are provided in e-/pix/s, so we have to multiply by the exposure time. - self.sky += roman.thermal_backgrounds[pointing.filter] * roman.exptime + self.sky += roman.thermal_backgrounds[pointing.filter_] * roman.exptime # Median of dark current is used here instead of mean since hot pixels contribute significantly to # the mean. Stastistics of dark current for the current test detector file: (mean, std, median, max) @@ -972,7 +972,7 @@ def add_persistence(self, im, pointing): dt = ( pointing.date - p.date ).total_seconds() - roman.exptime / 2 # avg time since end of exposures - self.params["output"]["file_name"]["items"] = [p.filter, p.visit, p.sca] + self.params["output"]["file_name"]["items"] = [p.filter_, p.visit, p.sca] imfilename = ParseValue(self.params["output"], "file_name", self.params, str)[0] fn = os.path.join(self.params["output"]["dir"], imfilename) @@ -1007,7 +1007,7 @@ def add_persistence(self, im, pointing): fac_dt = ( roman.exptime / 2.0 ) / dt # linear time dependence (approximate until we get t1 and Delat t of the data) - self.params["output"]["file_name"]["items"] = [p.filter, p.visit, p.sca] + self.params["output"]["file_name"]["items"] = [p.filter_, p.visit, p.sca] imfilename = ParseValue(self.params["output"], "file_name", self.params, str)[0] fn = os.path.join(self.params["output"]["dir"], imfilename) @@ -1265,7 +1265,8 @@ def finalize_sky_im(self, im, pointing): im_pad = self.qe(im_pad) im_pad = self.bfe(im_pad) - im_pad = self.add_persistence(im_pad, pointing) + if self.params["dither_from_file"]: + im_pad = self.add_persistence(im_pad, pointing) im_pad.quantize() im_pad += self.im_dark im_pad = self.saturate(im_pad) diff --git a/roman_imsim/obseq.py b/roman_imsim/obseq.py index 79c05263..103dce85 100644 --- a/roman_imsim/obseq.py +++ b/roman_imsim/obseq.py @@ -42,7 +42,7 @@ def read_obseq(self): self.ob["pa"] = ob["pa"] * galsim.degrees self.ob["date"] = Time(ob["date"], format="mjd").datetime self.ob["mjd"] = ob["date"] - self.ob["filter"] = ob["filter"] + self.ob["filter"] = "SNPrism" # ob["filter"] # HACK ALERT self.ob["exptime"] = ob["exptime"] def get(self, field, default=None): diff --git a/roman_imsim/photonOps/README.md b/roman_imsim/photonOps/README.md new file mode 100644 index 00000000..b67535e9 --- /dev/null +++ b/roman_imsim/photonOps/README.md @@ -0,0 +1,5 @@ +# Photon Operators + +This directory contains custom photon operators that can be specified in galsim configuration files. +Each photon operator is defined in its own module, with the file name starting with an underscore. +This hides the implementation details, while keeping the codebase modular for maintainability. \ No newline at end of file diff --git a/roman_imsim/photonOps/__init__.py b/roman_imsim/photonOps/__init__.py new file mode 100644 index 00000000..4428afa9 --- /dev/null +++ b/roman_imsim/photonOps/__init__.py @@ -0,0 +1,4 @@ +from ._charge_diffusion import * +from ._grism import * +from ._slitless_spec import * +from ._wfss_disperser import * diff --git a/roman_imsim/photonOps.py b/roman_imsim/photonOps/_charge_diffusion.py similarity index 94% rename from roman_imsim/photonOps.py rename to roman_imsim/photonOps/_charge_diffusion.py index 76c9faea..45fa6cc1 100644 --- a/roman_imsim/photonOps.py +++ b/roman_imsim/photonOps/_charge_diffusion.py @@ -1,12 +1,8 @@ import numpy as np from galsim import GaussianDeviate, PhotonOp, UniformDeviate -from galsim.config import ( - GetAllParams, - GetRNG, - PhotonOpBuilder, - RegisterPhotonOpType, - get_cls_params, -) +from galsim.config import GetAllParams, GetRNG, PhotonOpBuilder, RegisterPhotonOpType, get_cls_params + +__all__ = ["ChargeDiff"] _w1 = 0.17519 _w2 = 0.53146 @@ -18,9 +14,7 @@ class ChargeDiff(PhotonOp): - """A photon operator that applies the effect of charge diffusion via a - probabilistic model limit. - """ + """A photon operator that applies the effect of charge diffusion via a probablistic model limit.""" def __init__(self, rng=None, **kwargs): diff --git a/roman_imsim/photonOps/_grism/__init__.py b/roman_imsim/photonOps/_grism/__init__.py new file mode 100644 index 00000000..ae377956 --- /dev/null +++ b/roman_imsim/photonOps/_grism/__init__.py @@ -0,0 +1,2 @@ +from ._grism_operators import * +from .optical_model_utils import * diff --git a/roman_imsim/photonOps/_grism/_grism_operators.py b/roman_imsim/photonOps/_grism/_grism_operators.py new file mode 100644 index 00000000..8eb8772b --- /dev/null +++ b/roman_imsim/photonOps/_grism/_grism_operators.py @@ -0,0 +1,324 @@ +import numpy as np +import yaml +from galsim import GalSimError, PhotonOp +from galsim.config import GetAllParams, PhotonOpBuilder, RegisterPhotonOpType, get_cls_params + +from .optical_model_utils import RomanDetectorCoordinates + +__all__ = ["GrismNV", "GrismV"] + + +class GrismNV(PhotonOp): + """A photon operator that applies the dispersion effects of the + Roman Grism (nonvector form). + + The photons will need to have wavelengths defined in order to work. + + Parameters + base_wavelength: Wavelength (in nm) represented by the fiducial photon positions + """ + + # what parameters are tunable + # _req_params = {"base_wavelength": float, "barycenter": list} + # _opt_params = {"resolution": list} + + def __init__(self, config=None): + if config is None: + self.config = "optical_models/Roman_grism_OpticalModel_v0.8.yaml" + else: + self.config = config + self.order = "1" + self.sca = 16 + # self.base_wavelength = base_wavelength + # self.resolution = np.array(resolution) + with open(self.config) as f: + data = yaml.safe_load(f) + for order_key, order_data in data["optical_model"]["orders"].items(): + if order_key == self.order: + order_dat = order_data + break + + self.xmap_coeff = np.array(order_dat["xmap_ij_coeff"]) + self.ymap_coeff = np.array(order_dat["ymap_ij_coeff"]) + self.crv_coeff = np.array(order_dat["crv_ijk_coeff"]) + self.ids_coeff = np.array(order_dat["ids_ijk_coeff"]) + self.wl_min = 0.9 + self.wl_max = 2.0 + self.wl_ref = 1.55 + + # initialize RomanDetectorCoordinates + # need these for some of the functions called + self.detector_coords = RomanDetectorCoordinates( + naxis1=data["detector_model"].get("naxis1", 4096), + naxis2=data["detector_model"].get("naxis2", 4096), + crpix1=data["detector_model"].get("crpix1", 2048), + crpix2=data["detector_model"].get("crpix2", 2048), + pos_angle_detector=data["detector_model"].get("pos_angle_detector", 0.0), + pixel_scale=data["detector_model"].get("pixel_scale", 0.11), + plate_scale=data["detector_model"].get("plate_scale", 100.0), + xy_centers=data["detector_model"].get("xy_centers", {}), + ) + + def _disperse(self, x0, y0, lam, sca, order="1"): + # convert to FPA degrees + xfpa_deg, yfpa_deg = self.detector_coords.convert_sca_to_fpa(x0, y0, sca=sca) + # get the number of photons and the starting stuff + nphot = len(lam) + + xref_mm = np.zeros(nphot) + yref_mm = np.zeros(nphot) + delta_y_mm = np.zeros(nphot) + delta_x_mm = np.zeros(nphot) + + # loop through the photons and apply the IDS matrices + for p in range(nphot): + + x = xfpa_deg[p] + y = yfpa_deg[p] + wl = lam[p] + + xr = 0.0 + yr = 0.0 + + for i in range(6): + xi = x**i + for j in range(6): + yj = y**j + xr += self.xmap_coeff[i, j] * xi * yj + yr += self.ymap_coeff[i, j] * xi * yj + + xref_mm[p] = xr + yref_mm[p] = yr + + # compute position-dependent coeffs C_i(x,y) + C = np.zeros(6) + + for i in range(6): + ci = 0.0 + for j in range(6): + xj = x**j + for k in range(6): + yk = y**k + ci += self.ids_coeff[i, j, k] * xj * yk + C[i] = ci + + # apply wavelength polynomial (relative to wl_ref) + dy = 0.0 + for i in range(6): + dy += C[i] * (wl**i - self.wl_ref**i) + + delta_y_mm[p] = dy + + # curvature (x displacement as function of delta_y) + + D = np.zeros(6) + + for i in range(6): + di = 0.0 + for j in range(6): + xj = x**j + for k in range(6): + yk = y**k + di += self.crv_coeff[i, j, k] * xj * yk + D[i] = di + + dx = 0.0 + for i in range(6): + dx += D[i] * (dy**i) + + delta_x_mm[p] = dx + + xmpa_mm = xref_mm + delta_x_mm + ympa_mm = yref_mm + delta_y_mm + # go back to pixels + + x_pix, y_pix = self.detector_coords.convert_mpa_to_sca(xmpa=xmpa_mm, ympa=ympa_mm, sca=self.sca) + + return x_pix, y_pix + + def applyTo(self, photon_array, local_wcs=None, rng=None): + """Apply the grism dispersion to the photos + + Parameters: + photon_array: A `PhotonArray` to apply the operator to. + local_wcs: A `LocalWCS` instance defining the local WCS for the current photon + bundle in case the operator needs this information. [default: None] + rng: A random number generator is not used. + """ + # photon array has .x, .y, .wavelength, .coord, .time, ... + if not photon_array.hasAllocatedWavelengths(): + raise GalSimError("Grism requires that wavelengths be set") + + # wavelength is in nm. Roman slitless thinks in microns. + # http://galsim-developers.github.io/GalSim/_build/html/photon_array.html#galsim.PhotonArray + wavelength = photon_array.wavelength / 1000.0 + + x_pix, y_pix = self._disperse(photon_array.x, photon_array.y, wavelength, self.sca) + + # add to the photon array's positions + photon_array.x = x_pix + photon_array.y = y_pix + + def __repr__(self): + # ) + # s += ")" + s = "galsim.GrismNV()" + return s + + +class GrismNVBuilder(PhotonOpBuilder): + """Build a GrismNV""" + + # This one needs special handling for obj_coord + def buildPhotonOp(self, config, base, logger): + req, opt, single, takes_rng = get_cls_params(GrismNV) + kwargs, safe = GetAllParams(config, base, req, opt, single) + # if 'sky_pos' in base: + # kwargs['obj_coord'] = base['sky_pos'] + return GrismNV(**kwargs) + + +class GrismV(PhotonOp): + """A photon operator that applies the dispersion effects of the + Roman Grism (vector version) + + The photons will need to have wavelengths defined in order to work. + + Parameters: + config_file: Path to the YAML config file with optical model + order: Grism order (e.g., '1') + sca: SCA number + """ + + # what parameters are tunable + _req_params = {} + _opt_params = {} + _single_params = [] + _takes_rng = False + + def __init__(self, config=None, order=None, sca=None): + if config is None: + self.config = "optical_models/Roman_grism_OpticalModel_v0.8.yaml" + else: + self.config = config + self.order = "1" + self.sca = 16 + # self.base_wavelength = base_wavelength + # self.resolution = np.array(resolution) + with open(self.config) as f: + data = yaml.safe_load(f) + for order_key, order_data in data["optical_model"]["orders"].items(): + if order_key == self.order: + order_dat = order_data + break + + self.xmap_coeff = np.array(order_dat["xmap_ij_coeff"]) + self.ymap_coeff = np.array(order_dat["ymap_ij_coeff"]) + self.crv_coeff = np.array(order_dat["crv_ijk_coeff"]) + self.ids_coeff = np.array(order_dat["ids_ijk_coeff"]) + self.wl_min = 0.9 + self.wl_max = 2.0 + self.wl_ref = 1.55 + # initialize RomanDetectorCoordinates + # need these for some of the functions called + self.detector_coords = RomanDetectorCoordinates( + naxis1=data["detector_model"].get("naxis1", 4096), # Adjust these defaults as needed + naxis2=data["detector_model"].get("naxis2", 4096), + crpix1=data["detector_model"].get("crpix1", 2048), + crpix2=data["detector_model"].get("crpix2", 2048), + pos_angle_detector=data["detector_model"].get("pos_angle_detector", 0.0), + pixel_scale=data["detector_model"].get("pixel_scale", 0.11), + plate_scale=data["detector_model"].get("plate_scale", 100.0), + xy_centers=data["detector_model"].get("xy_centers", {}), # Dictionary of SCA centers + ) + + def _disperse(self, x0, y0, lam, sca, order="1"): + # get the coords in degrees to match the matrices + xfpa_deg, yfpa_deg = self.detector_coords.convert_sca_to_fpa(x0, y0, sca=self.sca) + + # we start by calculating the intial offset + xfpa_deg = np.atleast_1d(xfpa_deg) + yfpa_deg = np.atleast_1d(yfpa_deg) + # Create power matrices + x_powers = xfpa_deg[:, np.newaxis] ** np.arange(6) + y_powers = yfpa_deg[:, np.newaxis] ** np.arange(6) + # Compute the offset in mm for x + xref_mm = np.diagonal(x_powers @ self.xmap_coeff @ y_powers.T) + # Compute the offset in mm for y + yref_mm = np.diagonal(x_powers @ self.ymap_coeff @ y_powers.T) + + # now we get the y displacement (dispersion by wavelength) + + wavelength = np.atleast_1d(lam) + # create power matrices, note that we still need x and y in fpa degrees + wl_powers = wavelength[:, np.newaxis] ** np.arange(6) + ref_powers = self.wl_ref ** np.arange(6) + x_pow = ( + xfpa_deg[:, np.newaxis, np.newaxis, np.newaxis] + ** np.arange(6)[np.newaxis, np.newaxis, :, np.newaxis] + ) + y_pow = ( + yfpa_deg[:, np.newaxis, np.newaxis, np.newaxis] + ** np.arange(6)[np.newaxis, np.newaxis, np.newaxis, :] + ) + ids_at_pos = np.sum(self.ids_coeff[np.newaxis, :, :, :] * x_pow * y_pow, axis=(2, 3)) + delta_y_mm = wl_powers @ ids_at_pos.T - (ref_powers @ ids_at_pos.T) + delta_y_mm = np.diagonal(delta_y_mm) + + # now we get the x displacement (curvature) + crv_at_pos = np.sum(self.crv_coeff[np.newaxis, :, :, :] * x_pow * y_pow, axis=(2, 3)) + dy_powers = delta_y_mm[:, np.newaxis] ** np.arange(6) + delta_x_mm = dy_powers @ crv_at_pos.T + delta_x_mm = np.diagonal(delta_x_mm) + + # now combine them and then convert to pixels and set the photon array coords to them + xmpa_mm = xref_mm + delta_x_mm + ympa_mm = yref_mm + delta_y_mm + + x_pix, y_pix = self.detector_coords.convert_mpa_to_sca(xmpa=xmpa_mm, ympa=ympa_mm, sca=self.sca) + + return x_pix, y_pix + + def applyTo(self, photon_array, local_wcs=None, rng=None): + """Apply the grism dispersion to the photos + + Parameters: + photon_array: A `PhotonArray` to apply the operator to. + local_wcs: A `LocalWCS` instance defining the local WCS for the current photon + bundle in case the operator needs this information. [default: None] + rng: A random number generator is not used. + """ + # photon array has .x, .y, .wavelength, .coord, .time, ... + if not photon_array.hasAllocatedWavelengths(): + raise GalSimError("Grism requires that wavelengths be set") + + # wavelength is in nm. Roman slitless thinks in microns. + # http://galsim-developers.github.io/GalSim/_build/html/photon_array.html#galsim.PhotonArray + w = photon_array.wavelength / 1000.0 + + x_pix, y_pix = self._disperse(photon_array.x, photon_array.y, w, self.sca, self.order) + + photon_array.x = x_pix + photon_array.y = y_pix + + def __repr__(self): + + # ) + # s += ")" + return "galsim.GrismV()" + + +class GrismVBuilder(PhotonOpBuilder): + """Build a Grism op""" + + # This one needs special handling for obj_coord + def buildPhotonOp(self, config, base, logger): + req, opt, single, takes_rng = get_cls_params(GrismV) + kwargs, safe = GetAllParams(config, base, req, opt, single) + # if 'sky_pos' in base: + # kwargs['obj_coord'] = base['sky_pos'] + return GrismV(**kwargs) + + +RegisterPhotonOpType("GrismV", GrismVBuilder()) +RegisterPhotonOpType("GrismNV", GrismNVBuilder()) diff --git a/roman_imsim/photonOps/_grism/optical_model_utils.py b/roman_imsim/photonOps/_grism/optical_model_utils.py new file mode 100644 index 00000000..66faacb3 --- /dev/null +++ b/roman_imsim/photonOps/_grism/optical_model_utils.py @@ -0,0 +1,259 @@ +""" +Utilities for Roman Optical Model +""" + +import warnings + +import numpy as np +import pandas as pd +from matplotlib.patches import Rectangle + +__all__ = ["RomanDetectorCoordinates"] + + +class RomanDetectorCoordinates: + """ + Helper class to handle all the coordinate transformations for WFI + """ + + def __init__( + self, + naxis1, + naxis2, + crpix1, + crpix2, + pos_angle_detector, + pixel_scale, + plate_scale, + xy_centers, + ): + # Basic parameters + self.naxis1 = naxis1 + self.naxis2 = naxis2 + self.crpix1 = crpix1 + self.crpix2 = crpix2 + self.pos_angle_detector = pos_angle_detector + self.pixel_scale = pixel_scale + self.plate_scale = plate_scale + self.xy_centers = xy_centers + + if isinstance(self.xy_centers, dict): + # Convert the dictionary of (x,y) centers to a table + self.sca_list = np.array( + sorted(list(self.xy_centers.keys())), + dtype=int, + ) + self.xy_centers_tbl = pd.DataFrame(self.xy_centers, index=["x", "y"]).T + elif isinstance(self.xy_centers, pd.DataFrame): + self.sca_list = self.xy_centers.index.to_numpy() + self.xy_centers_tbl = self.xy_centers + else: + raise Exception("Provide xy_centers as a dict or pd.Dataframe") + + # Define the Polygons for each SCA + self.sca_polygons = self.define_sca_polygons() + + def get_sca_center(self, sca): + """ + Returns the center(s) for given SCAs + """ + xcen = self.xy_centers_tbl.loc[sca, "x"] + ycen = self.xy_centers_tbl.loc[sca, "y"] + if len(np.atleast_1d(sca)) > 1: + return xcen.to_numpy(), ycen.to_numpy() + else: + return xcen, ycen + + def define_sca_polygons(self, lw=2, facecolor="none", edgecolor="k", alpha=0.9, zorder=10): + """ + Defines the SCA locations on the MPA (in mm) and their sizes (in pix) + """ + sca_polygons = {} + for sca in self.sca_list: + + xy_cen = self.xy_centers_tbl.loc[sca] + xy_vrt = np.array(xy_cen) - ( + self.crpix1 / self.plate_scale, + self.crpix2 / self.plate_scale, + ) + dx = self.naxis1 / self.plate_scale + dy = self.naxis2 / self.plate_scale + + sca_polygons[sca] = Rectangle( + xy_vrt, + dx, + dy, + angle=self.pos_angle_detector, + rotation_point="center", + lw=lw, + facecolor=facecolor, + edgecolor=edgecolor, + alpha=alpha, + zorder=zorder, + ) + return sca_polygons + + def match_pos_to_sca(self, xmpa, ympa): + """ + Returns matching SCA in which the given point(s) lie + Units are in 'mm' + """ + xmpa, ympa = np.atleast_1d(xmpa), np.atleast_1d(ympa) + pos = np.array([xmpa, ympa]).T + cond = np.array( + [ + self.sca_polygons[sca].contains_points( + self.sca_polygons[sca].get_data_transform().transform(pos) + ) + for sca in self.sca_list + ], + dtype=bool, + ) + j, i = np.where(cond) + + sca_match = np.zeros(len(xmpa), dtype=int) + sca_match[i] = self.sca_list[j] + + if len(sca_match) > 1: + raise Exception("Position matched to multiple SCAs!") + elif sca_match == 0: + warnings.warn("No matching SCA; returning closest one.") + xdist = np.abs(xmpa - self.xy_centers_tbl["x"].to_numpy()) + sca_column = self.sca_list[np.argsort(xdist)[:3]] + ydist = np.abs(ympa - self.xy_centers_tbl.loc[sca_column, "y"].to_numpy()) + return sca_column[np.argmin(ydist)] + else: + return sca_match[0] + + def convert_sca_to_mpa(self, xsca, ysca, sca): # unused and broken + """ + Returns MPA position [mm] for reference position in the SCA [pixel] + Input: xsca, ysca are source position, in px, in detector plane (SCA) + Output: xmpa, ympa are source position, in mm, in focal plane (MPA) + """ + dx = (xsca - self.crpix1) * -1 + dy = ysca - self.crpix2 + + # Rotation terms might be needed here + xoff = dx / self.plate_scale + yoff = dy / self.plate_scale + + xmpa = xoff + self.xy_centers_tbl.loc[sca, "x"].to_numpy() + ympa = yoff + self.xy_centers_tbl.loc[sca, "y"].to_numpy() + + return xmpa, ympa + + def convert_mpa_to_sca(self, xmpa, ympa, sca=None): + """ + Returns SCA position [pixel] for reference position in the MPA [mm] + Input: xmpa, ympa are source position, in px, in detector plane (SCA) + Output: xfpa, yfpa are source position, in deg, in focal plane (FPA) + """ + if sca is None: + sca = self.match_pos_to_sca(xmpa=xmpa, ympa=ympa) + + xcen, ycen = self.get_sca_center(sca=sca) + xoff = xmpa - xcen + yoff = ympa - ycen + + # Rotation terms might be needed here + xrot = xoff * self.plate_scale + yrot = yoff * self.plate_scale + + xdet = -xrot + self.crpix1 + ydet = yrot + self.crpix2 + + return xdet, ydet + + def convert_sca_to_fpa(self, xsca, ysca, sca): + """ + Returns FPA position [degree] for reference position in the SCA [px] + Input: xsca, ysca are source position, in px, in detector plane (SCA) + Output: xfpa, yfpa are source position, in deg, in focal plane (FPA) + """ + dx = (xsca - self.crpix1) * -1 + dy = ysca - self.crpix2 + + xcen, ycen = self.get_sca_center(sca=sca) + xfpa = (xcen * self.plate_scale + dx) * self.pixel_scale / 3600 + yfpa = (ycen * self.plate_scale + dy) * self.pixel_scale / 3600 + + return xfpa, yfpa + + def convert_fpa_to_sca(self, xfpa, yfpa, sca=None): + """ + Returns SCA position [pixel] for reference position in the FPA [degree] + Input: xfpa, yfpa are source position, in deg, in focal plane (FPA) + Output: xsca, ysca are source position, in px, in detector plane (SCA) + """ + if sca is None: + xmpa, ympa = self.convert_fpa_to_mpa(xfpa=xfpa, yfpa=yfpa) + sca = self.match_pos_to_sca(xmpa=xmpa, ympa=ympa) + + xcen, ycen = self.get_sca_center(sca=sca) + dx = (xfpa * 3600 / self.pixel_scale) - (xcen * self.plate_scale) + dy = (yfpa * 3600 / self.pixel_scale) - (ycen * self.plate_scale) + + xsca = -dx + self.crpix1 + ysca = dy + self.crpix2 + + return xsca, ysca + + def convert_fpa_to_mpa(self, xfpa, yfpa): + + xmpa = xfpa * 3600 / self.pixel_scale / self.plate_scale + ympa = yfpa * 3600 / self.pixel_scale / self.plate_scale + + return xmpa, ympa + + def convert_mpa_to_fpa(self, xmpa, ympa): + + xfpa = xmpa * self.plate_scale * self.pixel_scale / 3600 + yfpa = ympa * self.plate_scale * self.pixel_scale / 3600 + + return xfpa, yfpa + + def generate_xyfpa_grid(self, npts, custom_slice=None): + """ + Utility function to generate a grid of (x,y) in FPA coords + """ + xref_fpa, yref_fpa, ref_sca = [], [], [] + for sca in self.sca_list: + + xref, yref = np.mgrid[ + 0 : self.naxis1 : npts * 1j, + 0 : self.naxis2 : npts * 1j, + ] + + if custom_slice is not None: + xref = xref[custom_slice, custom_slice] + yref = yref[custom_slice, custom_slice] + + xyref_fpa = self.convert_sca_to_fpa( + xsca=xref.ravel(), + ysca=yref.ravel(), + sca=sca, + ) + xref_fpa = np.concatenate([xref_fpa, xyref_fpa[0]]) + yref_fpa = np.concatenate([yref_fpa, xyref_fpa[1]]) + ref_sca = np.concatenate([ref_sca, [sca] * len(xyref_fpa[0])]) + + return xref_fpa, yref_fpa, ref_sca + + +def transform_wl(wl, wl_reference, kind, inverse=False): + """ + Helper function to transform the wllengths for use in the optical model + """ + if kind == "linear": + if not inverse: + return wl - wl_reference + else: + return wl + wl_reference + elif kind == "log": + if not inverse: + return np.log10(wl / wl_reference) + else: + return 10**wl * wl_reference + else: + raise Exception("Invalid kind provided: choose from 'linear' or 'log'.") diff --git a/roman_imsim/photonOps/_slitless_spec.py b/roman_imsim/photonOps/_slitless_spec.py new file mode 100644 index 00000000..53641ace --- /dev/null +++ b/roman_imsim/photonOps/_slitless_spec.py @@ -0,0 +1,74 @@ +from galsim import GalSimError, PhotonOp +from galsim.config import GetAllParams, PhotonOpBuilder, RegisterPhotonOpType, get_cls_params + +__all__ = ["SlitlessSpec"] + + +class SlitlessSpec(PhotonOp): + r"""A photon operator that applies the dispersion effects of the + Roman Prism. + + The photons will need to have wavelengths defined in order to work. + + Parameters: + base_wavelength: Wavelength (in nm) represented by the fiducial photon positions + """ + + # what parameters are tunable + # _req_params = {"base_wavelength": float, "barycenter": list} + # _opt_params = {"resolution": list} + + def __init__(self): + # self.base_wavelength = base_wavelength + # self.resolution = np.array(resolution) + pass + + def applyTo(self, photon_array, local_wcs=None, rng=None): + """Apply the slitless-spectroscopy disspersion to the photos + + Parameters: + photon_array: A `PhotonArray` to apply the operator to. + local_wcs: A `LocalWCS` instance defining the local WCS for the current photon + bundle in case the operator needs this information. [default: None] + rng: A random number generator is not used. + """ + # photon array has .x, .y, .wavelength, .coord, .time, ... + if not photon_array.hasAllocatedWavelengths(): + raise GalSimError("SlitlessSpec requires that wavelengths be set") + + # wavelength is in nm. Roman slitless thinks in microns. + # http://galsim-developers.github.io/GalSim/_build/html/photon_array.html#galsim.PhotonArray + w = photon_array.wavelength / 1000.0 + + # dx = (-12.973976 + 213.353667*(w - 1.0) + -20.254574*(w - 1.0)**2) / + # (1.0 + 1.086448*(w - 1.0) + -0.573796*(w - 1.0)**2) + dy = (-81.993865 + 138.367237 * (w - 1.0) + 19.348549 * (w - 1.0) ** 2) / ( + 1.0 + 1.086447 * (w - 1.0) + -0.573797 * (w - 1.0) ** 2 + ) + + photon_array.y += dy + + # might need to change dxdz/dydz for the angle of travel through the detector. + + def __repr__(self): + # s = "galsim.SlitlessSpec(base_wavelength=%r, " % ( + # self.base_wavelength, + # ) + # s += ")" + s = "galsim.SlitlessSpec()" + return s + + +class SlitlessSpecBuilder(PhotonOpBuilder): + """Build a SlitlessSpec""" + + # This one needs special handling for obj_coord + def buildPhotonOp(self, config, base, logger): + req, opt, single, takes_rng = get_cls_params(SlitlessSpec) + kwargs, safe = GetAllParams(config, base, req, opt, single) + # if 'sky_pos' in base: + # kwargs['obj_coord'] = base['sky_pos'] + return SlitlessSpec(**kwargs) + + +RegisterPhotonOpType("SlitlessSpec", SlitlessSpecBuilder()) diff --git a/roman_imsim/photonOps/_wfss_disperser/__init__.py b/roman_imsim/photonOps/_wfss_disperser/__init__.py new file mode 100644 index 00000000..e9a3e527 --- /dev/null +++ b/roman_imsim/photonOps/_wfss_disperser/__init__.py @@ -0,0 +1 @@ +from ._wfss_disperser_operators import * diff --git a/roman_imsim/photonOps/_wfss_disperser/_wfss_disperser_operators.py b/roman_imsim/photonOps/_wfss_disperser/_wfss_disperser_operators.py new file mode 100644 index 00000000..14cc8eee --- /dev/null +++ b/roman_imsim/photonOps/_wfss_disperser/_wfss_disperser_operators.py @@ -0,0 +1,155 @@ +import numpy as np +from galsim import GalSimError, PhotonOp +from galsim.config import GetAllParams, PhotonOpBuilder, RegisterPhotonOpType, get_cls_params + +from .snpitdispenser import SNPITDisperser + +__all__ = ["WFSSSDisperser"] + + +class WFSSSDisperser(PhotonOp): + """A photon operator that applies the dispersion effects of the + Roman Grism (vector version) + + The photons will need to have wavelengths defined in order to work. + + Parameters: + config_file: Path to the YAML config file with optical model + order: Grism order (e.g., '1') + sca: SCA number + """ + + # what parameters are tunable + _req_params = {} + _opt_params = {} + _single_params = [] + _takes_rng = False + + def __init__(self, config=None, order=None, sca=None): + if config is None: + self.config = "optical_models/Roman_prism_OpticalModel_v0.8.yaml" + else: + self.config = config + self.order = "1" + self.sca = sca or 16 # TODO: This is no good. + # # self.base_wavelength = base_wavelength + # # self.resolution = np.array(resolution) + # with open(self.config) as f: + # data = yaml.safe_load(f) + # for order_key, order_data in data["optical_model"]["orders"].items(): + # if order_key == self.order: + # order_dat = order_data + # break + + # self.xmap_coeff = np.array(order_dat['xmap_ij_coeff']) + # self.ymap_coeff = np.array(order_dat['ymap_ij_coeff']) + # self.crv_coeff = np.array(order_dat['crv_ijk_coeff']) + # self.ids_coeff = np.array(order_dat['ids_ijk_coeff']) + # self.wl_min = 0.9 + # self.wl_max = 2.0 + # self.wl_ref = 1.55 + # # initialize RomanDetectorCoordinates + # # need these for some of the functions called + # self.detector_coords = RomanDetectorCoordinates( + # naxis1=data['detector_model'].get('naxis1', 4096), # Adjust these defaults as needed + # naxis2=data['detector_model'].get('naxis2', 4096), + # crpix1=data['detector_model'].get('crpix1', 2048), + # crpix2=data['detector_model'].get('crpix2', 2048), + # pos_angle_detector=data['detector_model'].get('pos_angle_detector', 0.0), + # pixel_scale=data['detector_model'].get('pixel_scale', 0.11), + # plate_scale=data['detector_model'].get('plate_scale', 100.0), + # xy_centers=data['detector_model'].get('xy_centers', {}) # Dictionary of SCA centers + # ) + self.snpit_disperser = SNPITDisperser(self.config) + + def _disperse(self, x0, y0, lam, sca, order="1"): + # get the coords in degrees to match the matrices + xfpa_deg, yfpa_deg = self.detector_coords.convert_sca_to_fpa(x0, y0, sca=self.sca) + + # we start by calculating the intial offset + xfpa_deg = np.atleast_1d(xfpa_deg) + yfpa_deg = np.atleast_1d(yfpa_deg) + # Create power matrices + x_powers = xfpa_deg[:, np.newaxis] ** np.arange(6) + y_powers = yfpa_deg[:, np.newaxis] ** np.arange(6) + # Compute the offset in mm for x + xref_mm = np.diagonal(x_powers @ self.xmap_coeff @ y_powers.T) + # Compute the offset in mm for y + yref_mm = np.diagonal(x_powers @ self.ymap_coeff @ y_powers.T) + + # now we get the y displacement (dispersion by wavelength) + + wavelength = np.atleast_1d(lam) + # create power matrices, note that we still need x and y in fpa degrees + wl_powers = wavelength[:, np.newaxis] ** np.arange(6) + ref_powers = self.wl_ref ** np.arange(6) + x_pow = ( + xfpa_deg[:, np.newaxis, np.newaxis, np.newaxis] + ** np.arange(6)[np.newaxis, np.newaxis, :, np.newaxis] + ) + y_pow = ( + yfpa_deg[:, np.newaxis, np.newaxis, np.newaxis] + ** np.arange(6)[np.newaxis, np.newaxis, np.newaxis, :] + ) + ids_at_pos = np.sum(self.ids_coeff[np.newaxis, :, :, :] * x_pow * y_pow, axis=(2, 3)) + delta_y_mm = wl_powers @ ids_at_pos.T - (ref_powers @ ids_at_pos.T) + delta_y_mm = np.diagonal(delta_y_mm) + + # now we get the x displacement (curvature) + crv_at_pos = np.sum(self.crv_coeff[np.newaxis, :, :, :] * x_pow * y_pow, axis=(2, 3)) + dy_powers = delta_y_mm[:, np.newaxis] ** np.arange(6) + delta_x_mm = dy_powers @ crv_at_pos.T + delta_x_mm = np.diagonal(delta_x_mm) + + # now combine them and then convert to pixels and set the photon array coords to them + xmpa_mm = xref_mm + delta_x_mm + ympa_mm = yref_mm + delta_y_mm + + x_pix, y_pix = self.detector_coords.convert_mpa_to_sca(xmpa=xmpa_mm, ympa=ympa_mm, sca=self.sca) + + return x_pix, y_pix + + def applyTo(self, photon_array, local_wcs=None, rng=None): + """Apply the grism dispersion to the photos + + Parameters: + photon_array: A `PhotonArray` to apply the operator to. + local_wcs: A `LocalWCS` instance defining the local WCS for the current photon + bundle in case the operator needs this information. [default: None] + rng: A random number generator is not used. + """ + # photon array has .x, .y, .wavelength, .coord, .time, ... + if not photon_array.hasAllocatedWavelengths(): + raise GalSimError("Grism requires that wavelengths be set") + + # wavelength is in nm. Roman slitless thinks in microns. + # http://galsim-developers.github.io/GalSim/_build/html/photon_array.html#galsim.PhotonArray + w = photon_array.wavelength / 1000.0 + + x_pix, y_pix = self.snpit_disperser.disperse( + photon_array.x, photon_array.y, w, self.sca, self.order, pairwise=True + ) + + photon_array.x = x_pix + photon_array.y = y_pix + + def __repr__(self): + + # ) + # s += ")" + return "galsim.GrismV()" + + +class WFSSSDisperserBuilder(PhotonOpBuilder): + """Build a Grism op""" + + # This one needs special handling for obj_coord + def buildPhotonOp(self, config, base, logger): + req, opt, single, takes_rng = get_cls_params(WFSSSDisperser) + kwargs, safe = GetAllParams(config, base, req, opt, single) + # if 'sky_pos' in base: + # kwargs['obj_coord'] = base['sky_pos'] + return WFSSSDisperser(**kwargs) + + +RegisterPhotonOpType("WFSSSDisperser", WFSSSDisperserBuilder()) diff --git a/roman_imsim/photonOps/_wfss_disperser/snpitdispenser.py b/roman_imsim/photonOps/_wfss_disperser/snpitdispenser.py new file mode 100644 index 00000000..9e47f8a4 --- /dev/null +++ b/roman_imsim/photonOps/_wfss_disperser/snpitdispenser.py @@ -0,0 +1,839 @@ +import numba as nb +import numpy as np +import yaml + +""" +Classes and numba-accelerated functions to implement the optical +model provided by the SSC (IPAC). + +Written by R. Ryan (STScI) + +Feb 8, 2026 +""" + + +@nb.njit +def horner1(c, z): + """ + Evaluate a 1D polynomial at z using Horner's method. + + Parameters + ---------- + c : ndarray, shape (n,) + Polynomial coefficients ordered from lowest to highest degree. + z : float or ndarray + Evaluation point(s). + + Returns + ------- + p : float or ndarray + Polynomial value at z. + """ + nz = c.shape[0] + p = c[-1] + for i in range(nz - 2, -1, -1): + p *= z + p += c[i] + return p + + +@nb.njit +def horner2(c, y, z): + """ + Evaluate a 2D polynomial at (y, z) using nested Horner's method. + + Parameters + ---------- + c : ndarray, shape (ny, nz) + Polynomial coefficients, where c[i, j] corresponds to y^i z^j. + y, z : float or ndarray + Evaluation point(s). + + Returns + ------- + p : float or ndarray + Polynomial value at (y, z). + """ + ny = c.shape[0] + p = horner1(c[-1], z) + for i in range(ny - 2, -1, -1): + p *= y + p += horner1(c[i], z) + return p + + +@nb.njit +def horner3(c, x, y, z): + """ + Evaluate a 3D polynomial at (x, y, z) using nested Horner's method. + + Parameters + ---------- + c : ndarray, shape (nx, ny, nz) + Polynomial coefficients, where c[i, j, k] corresponds to x^i y^j z^k. + x, y, z : float or ndarray + Evaluation point(s). + + Returns + ------- + p : float or ndarray + Polynomial value at (x, y, z). + """ + nx = c.shape[0] + p = horner2(c[-1], y, z) + for i in range(nx - 2, -1, -1): + p *= x + p += horner2(c[i], y, z) + return p + + +@nb.njit(parallel=True) +def numba_disperse_pairwise(Xij, Yij, Cijk, Dijk, x, y, w): + """ + Compute dispersed coordinates for pairwise inputs using Numba. + + This version assumes x, y, and w have identical shapes and computes + one output coordinate pair per input triplet. + + Parameters + ---------- + Xij, Yij : ndarray, shape (nx, ny) + Polynomial coefficients mapping undispersed coordinates to + Focal Plane Assembly (FPA) space. + Cijk, Dijk : ndarray, shape (nw, nx, ny) + Polynomial coefficients for dispersion in FPA space. + x, y : ndarray + Undispersed FPA coordinates. + w : ndarray + Transformed wavelength coordinate. + + Returns + ------- + xt, yt : ndarray + Dispersed FPA coordinates. + """ + xt = np.empty_like(x) + yt = np.empty_like(x) + + for idx in nb.prange(x.size): + xx = x.flat[idx] + yy = y.flat[idx] + ww = w.flat[idx] + + dy = horner3(Cijk, ww, xx, yy) + dx = horner3(Dijk, dy, xx, yy) + + xmpa = horner2(Xij, xx, yy) + ympa = horner2(Yij, xx, yy) + + xt.flat[idx] = xmpa + dx + yt.flat[idx] = ympa + dy + + return xt, yt + + +@nb.njit(parallel=True) +def numba_disperse(Xij, Yij, Cijk, Dijk, x, y, w): + """ + Compute dispersed coordinates on a grid using Numba. + + This version assumes x and y are 1D arrays of positions, and w is a + 1D array of wavelengths. Outputs are 2D arrays indexed by + (pixel_index, wavelength_index). + + Parameters + ---------- + Xij, Yij : ndarray, shape (nx, ny) + Polynomial coefficients mapping undispersed coordinates to + Focal Plane Assembly (FPA) space. + Cijk, Dijk : ndarray, shape (nw, nx, ny) + Polynomial coefficients for dispersion in FPA space. + x, y : ndarray, shape (npix,) + Undispersed FPA coordinates. + w : ndarray, shape (nlam,) + Transformed wavelength coordinate. + + Returns + ------- + xt, yt : ndarray, shape (npix, nlam) + Dispersed FPA coordinates. + """ + npix = x.size + nlam = w.size + + dim = (npix, nlam) + xt = np.empty(dim, dtype=float) + yt = np.empty(dim, dtype=float) + + for idx in nb.prange(npix): + xx = x.flat[idx] + yy = y.flat[idx] + + xmpa = horner2(Xij, xx, yy) + ympa = horner2(Yij, xx, yy) + + for jdx in nb.prange(nlam): + ww = w.flat[jdx] + + dy = horner3(Cijk, ww, xx, yy) + dx = horner3(Dijk, dy, xx, yy) + + xt[idx, jdx] = xmpa + dx + yt[idx, jdx] = ympa + dy + return xt, yt + + +@nb.njit(parallel=True) +def numba_deriv_pairwise(dCijk, Cijk, dDijk, x, y, w): + """ + Compute derivatives of dispersed coordinates with respect to wavelength + for pairwise inputs using Numba. + + Parameters + ---------- + dCijk : ndarray, shape (nw-1, nx, ny) + Derivative of Cijk coefficients with respect to the transformed + wavelength coordinate. + Cijk : ndarray, shape (nw, nx, ny) + Original dispersion coefficients. + dDijk : ndarray, shape (nw-1, nx, ny) + Derivative of Dijk coefficients with respect to the dispersion + coordinate. + x, y : ndarray + Undispersed FPA coordinates. + w : ndarray + Transformed wavelength coordinate. + + Returns + ------- + dxdw, dydw : ndarray + Derivatives of dispersed FPA coordinates with respect to transformed + wavelength. + """ + dxdw = np.empty_like(x) + dydw = np.empty_like(x) + for idx in nb.prange(x.size): + xx = x.flat[idx] + yy = y.flat[idx] + ww = w.flat[idx] + + dy = horner3(Cijk, ww, xx, yy) + + dy_dw = horner3(dCijk, ww, xx, yy) + dx_dy = horner3(dDijk, dy, xx, yy) + + dxdw[idx] = dy_dw * dx_dy + dydw[idx] = dy_dw + + return dxdw, dydw + + +@nb.njit(parallel=True) +def numba_deriv(dCijk, Cijk, dDijk, x, y, w): + """ + Compute derivatives of dispersed coordinates with respect to wavelength + on a grid using Numba. + + Parameters + ---------- + dCijk : ndarray, shape (nw-1, nx, ny) + Derivative of Cijk coefficients with respect to the transformed + wavelength coordinate. + Cijk : ndarray, shape (nw, nx, ny) + Original dispersion coefficients. + dDijk : ndarray, shape (nw-1, nx, ny) + Derivative of Dijk coefficients with respect to the dispersion + coordinate. + x, y : ndarray, shape (npix,) + Undispersed FPA coordinates. + w : ndarray, shape (nlam,) + Transformed wavelength coordinate. + + Returns + ------- + dxdw, dydw : ndarray, shape (npix, nlam) + Derivatives of dispersed FPA coordinates with respect to transformed + wavelength. + """ + npix = x.size + nlam = w.size + + dim = (npix, nlam) + dxdw = np.empty(dim, dtype=float) + dydw = np.empty(dim, dtype=float) + + for idx in nb.prange(x.size): + xx = x.flat[idx] + yy = y.flat[idx] + + for jdx in nb.prange(nlam): + ww = w.flat[jdx] + + dy = horner3(Cijk, ww, xx, yy) + + dy_dw = horner3(dCijk, ww, xx, yy) + dx_dy = horner3(dDijk, dy, xx, yy) + + dxdw[idx, jdx] = dy_dw * dx_dy + dydw[idx, jdx] = dy_dw + + return dxdw, dydw + + +class LogTransformer: + """ + Logarithmic wavelength transformer. + + Transforms physical wavelength values into a log-scaled coordinate + relative to a reference wavelength, and provides inverse and derivative + mappings. + """ + + ln10 = np.log(10.0) + + def __init__(self, lam0): + """ + Parameters + ---------- + lam0 : float + Reference wavelength. + """ + self.lam0 = lam0 + + def evaluate(self, lam): + """ + Transform physical wavelength values into log-scaled coordinates. + + Parameters + ---------- + lam : array_like + Physical wavelengths. + + Returns + ------- + w : ndarray + Log-scaled wavelength coordinates. + """ + return np.log10(lam / self.lam0) + + def invert(self, w): + """ + Invert the log-scaled wavelength transformation. + + Parameters + ---------- + w : array_like + Log-scaled wavelength coordinates. + + Returns + ------- + lam : ndarray + Physical wavelengths. + """ + return self.lam0 * (10.0**w) + + def deriv(self, lam): + """ + Compute dλ/dw for the log-wavelength transform. + + Parameters + ---------- + lam : array_like + Physical wavelengths. + + Returns + ------- + dldw : ndarray + Derivative of wavelength with respect to the transformed coordinate. + """ + return lam * self.ln10 + + +class LinearTransformer: + """ + Linear wavelength transformer. + + Transforms physical wavelength values into a linear coordinate + relative to a reference wavelength, and provides inverse and derivative + mappings. + """ + + def __init__(self, lam0): + """ + Parameters + ---------- + lam0 : float + Reference wavelength. + """ + self.lam0 = lam0 + + def evaluate(self, lam): + """ + Transform physical wavelength values into linear coordinates. + + Parameters + ---------- + lam : array_like + Physical wavelengths. + + Returns + ------- + w : ndarray + Linear wavelength coordinates. + """ + return lam - self.lam0 + + def invert(self, w): + """ + Invert the linear wavelength transformation. + + Parameters + ---------- + w : array_like + Linear wavelength coordinates. + + Returns + ------- + lam : ndarray + Physical wavelengths. + """ + return w + self.lam0 + + def deriv(self, lam): + """ + Compute dλ/dw for the linear wavelength transform. + + Parameters + ---------- + lam : array_like + Physical wavelengths. + + Returns + ------- + dldw : float or ndarray + Derivative of wavelength with respect to the transformed coordinate. + """ + return 1.0 + + +class SNPITDisperser: + """ + Spectral disperser model for the Roman/WFI, optimized to meet the needs + of the SN PIT. + + This class loads a configuration file describing detector geometry + and optical polynomial models, and provides methods to compute: + + - Dispersed positions on the Focal Plane Assembly (FPA), + - Derivatives with respect to wavelength, + - Local trace normals, + - Local dispersion scales. + + All intermediate optical coordinates are expressed in Focal Plane + Assembly (FPA) physical units and mapped to and from Sensor Chip + Assembly (SCA) pixel coordinates. + """ + + def __init__(self, conffile): + """ + Initialize the disperser from a configuration file. + + Parameters + ---------- + conffile : str + Path to the YAML configuration file. + """ + self.load_conffile(conffile) + + # pre-compile the numba functions + _, _ = self.disperse(1.0, 1.0, 1.0, 1, pairwise=True) + _, _ = self.disperse(1.0, 1.0, 1.0, 1, pairwise=False) + + _, _ = self.deriv(1.0, 1.0, 1.0, 1, pairwise=True) + _, _ = self.deriv(1.0, 1.0, 1.0, 1, pairwise=False) + + def load_conffile(self, conffile): + """ + Load and parse the disperser configuration file. + + This method reads detector geometry, optical polynomial coefficients, + and wavelength transformation settings, and precomputes derivative + coefficients for efficient runtime evaluation. + + Parameters + ---------- + conffile : str + Path to the YAML configuration file. + """ + self.conffile = conffile + + with open(self.conffile, "r") as fp: + cfg = yaml.safe_load(fp) + + # fetch the meta data + self.meta = cfg["meta"] + + # set the detector model + self.detector = cfg["detector_model"] + + # change some units + self.detector["pixel_scale"] /= 3600.0 + + # set the optical model + self.optical = {k: v for k, v in cfg["optical_model"].items()} + + # process each spectral order + self.optical["orders"] = {} + for order, data in cfg["optical_model"]["orders"].items(): + + # extract the parameters + Xij = np.asarray(data["xmap_ij_coeff"]) + Yij = np.asarray(data["ymap_ij_coeff"]) + Cijk = np.asarray(data["ids_ijk_coeff"]) + Dijk = np.asarray(data["crv_ijk_coeff"]) + + # check that the data are valid + if Xij.all() and Yij.all() and Cijk.all() and Dijk.all(): + + # compute derivatives + dCijk = self.deriv_coeffs(Cijk) + dDijk = self.deriv_coeffs(Dijk) + + # save the results + self.optical["orders"][order] = { + "Xij": Xij, + "Yij": Yij, + "Cijk": Cijk, + "dCijk": dCijk, + "Dijk": Dijk, + "dDijk": dDijk, + } + + # transform the wavelengths + transform = self.optical["wl_transform"].lower() + if transform == "log": + self.lam_transformer = LogTransformer(self.optical["wl_reference"]) + elif transform == "linear": + self.lam_transformer = LinearTransformer(self.optical["wl_reference"]) + else: + raise NotImplementedError(f"Invalid transform {transform}") + + @staticmethod + def deriv_coeffs(M): + """ + Compute derivative polynomial coefficients with respect to the first + dimension using Horner-compatible form. + + Parameters + ---------- + M : ndarray, shape (n, ...) + Polynomial coefficient tensor where the first axis corresponds + to the variable of differentiation. + + Returns + ------- + dM : ndarray, shape (n-1, ...) + Derivative coefficient tensor. + """ + n = M.shape[0] + ii = np.arange(1, n, dtype=float) + shape = (n - 1,) + (1,) * (M.ndim - 1) + dM = M[1:] * ii.reshape(shape) + + return dM + + def validate(self, x, y, lam, sca, order="1", pairwise=True): + """ + Validate input array dimensionality for disperser evaluation. + + Parameters + ---------- + x, y : ndarray + Undispersed FPA coordinates. + lam : ndarray + Wavelength array. + sca : int + The SCA number. + order : str + The spectral order. + pairwise : bool + If True, require x, y, and l to have identical shapes. + If False, require x and y to match and l to be one-dimensional. + + """ + if order not in self.optical["orders"]: + raise KeyError("Invalid spectral order") + + if sca not in self.detector["xy_centers"]: + raise KeyError("Invalid SCA number") + + if pairwise: + if (x.shape != y.shape) or (x.shape != lam.shape): + raise RuntimeError("Invalid x, y, lam shape for pairwise") + else: + if (x.shape != y.shape) or (lam.ndim != 1): + raise RuntimeError("Invalid x, y, lam shape") + + def sca_to_fpa(self, xsca, ysca, sca): + """ + Convert Sensor Chip Assembly (SCA) pixel coordinates to + Focal Plane Assembly (FPA) physical coordinates. + + Parameters + ---------- + xsca, ysca : array_like + Pixel coordinates on the SCA detector. + sca : int or key + Identifier for the detector segment. + + Returns + ------- + xfpa, yfpa : ndarray + Physical coordinates in the Focal Plane Assembly (FPA) frame. + """ + dx = self.detector["crpix1"] - xsca # WHY NEGATIVE? + dy = ysca - self.detector["crpix2"] + + xcen, ycen = self.detector["xy_centers"][sca] + + xfpa = (xcen * self.detector["plate_scale"] + dx) * self.detector["pixel_scale"] + yfpa = (ycen * self.detector["plate_scale"] + dy) * self.detector["pixel_scale"] + + return xfpa, yfpa + + def mpa_to_sca(self, xmpa, ympa, sca): + """ + Convert Mosaic Plate Assembly (MPA) physical coordinates back to + Sensor Chip Assembly (SCA) pixel coordinates. + + Parameters + ---------- + xmpa, ympa : array_like + Physical Mosaic Plate Assembly coordinates. + sca : int or key + Identifier for the detector segment. + + Returns + ------- + xsca, ysca : ndarray + Detector pixel coordinates. + """ + xcen, ycen = self.detector["xy_centers"][sca] + + xoff = (xmpa - xcen) * self.detector["plate_scale"] + yoff = (ympa - ycen) * self.detector["plate_scale"] + + xsca = self.detector["crpix1"] - xoff # WHY NEGATIVE? + ysca = self.detector["crpix2"] + yoff + + return xsca, ysca + + def disperse(self, x0, y0, lam, sca, order="1", pairwise=False): + """ + Compute dispersed detector coordinates for given source positions + and wavelengths. + + Parameters + ---------- + x0, y0 : array_like + Undispersed Sensor Chip Assembly (SCA) pixel coordinates. + lam : array_like + Physical wavelengths. + sca : int or key + Detector segment identifier. + order : str, optional + Spectral order identifier (default is '1'). + pairwise : bool, optional + If True, treat x0, y0, and lam as pairwise-aligned arrays. + If False, compute a full grid over x0/y0 and lam. + + Returns + ------- + xp, yp : ndarray + Dispersed SCA pixel coordinates. + """ + x0 = np.atleast_1d(x0) + y0 = np.atleast_1d(y0) + lam = np.atleast_1d(lam) + + # check inputs + self.validate(x0, y0, lam, sca, order=order, pairwise=pairwise) + + # convert SCA coordinates to FPA coordinates + xfpa, yfpa = self.sca_to_fpa(x0, y0, sca) + + # transform wavelengths + wtran = self.lam_transformer.evaluate(lam) + + if pairwise: + xt, yt = numba_disperse_pairwise( + self.optical["orders"][order]["Xij"], + self.optical["orders"][order]["Yij"], + self.optical["orders"][order]["Cijk"], + self.optical["orders"][order]["Dijk"], + xfpa, + yfpa, + wtran, + ) + else: + xt, yt = numba_disperse( + self.optical["orders"][order]["Xij"], + self.optical["orders"][order]["Yij"], + self.optical["orders"][order]["Cijk"], + self.optical["orders"][order]["Dijk"], + xfpa, + yfpa, + wtran, + ) + + # convert back to SCA coordinates (this method should deal with + # the traces that span multiple SCAs) + xp, yp = self.mpa_to_sca(xt, yt, sca) + + return xp, yp + + def deriv(self, x0, y0, lam, sca, order="1", pairwise=False): + """ + Compute derivatives of dispersed detector coordinates with respect + to wavelength. + + Parameters + ---------- + x0, y0 : array_like + Undispersed Sensor Chip Assembly (SCA) pixel coordinates. + lam : array_like + Physical wavelengths. + sca : int or key + Detector segment identifier. + order : str, optional + Spectral order identifier (default is '1'). + pairwise : bool, optional + If True, treat x0, y0, and lam as pairwise-aligned arrays. + If False, compute a full grid over x0/y0 and lam. + + Returns + ------- + dxdl, dydl : ndarray + Derivatives of dispersed SCA pixel coordinates with respect to + wavelength. + """ + x0 = np.atleast_1d(x0) + y0 = np.atleast_1d(y0) + lam = np.atleast_1d(lam) + + # check inputs + self.validate(x0, y0, lam, sca, order=order, pairwise=pairwise) + + # convert SCA coordinates to FPA coordinates + xfpa, yfpa = self.sca_to_fpa(x0, y0, sca) + + # transform wavelengths + wtran = self.lam_transformer.evaluate(lam) + dldw = self.lam_transformer.deriv(lam) + + if pairwise: + dxdw, dydw = numba_deriv_pairwise( + self.optical["orders"][order]["dCijk"], + self.optical["orders"][order]["Cijk"], + self.optical["orders"][order]["dDijk"], + xfpa, + yfpa, + wtran, + ) + + dxdl = dxdw / dldw + dydl = dydw / dldw + + else: + dxdw, dydw = numba_deriv( + self.optical["orders"][order]["dCijk"], + self.optical["orders"][order]["Cijk"], + self.optical["orders"][order]["dDijk"], + xfpa, + yfpa, + wtran, + ) + + dxdl = dxdw / dldw[np.newaxis, :] + dydl = dydw / dldw[np.newaxis, :] + + dxdl *= self.detector["plate_scale"] + dydl *= self.detector["plate_scale"] + + return dxdl, dydl + + def normal(self, *args, **kwargs): + """ + Compute unit-normal vectors to the spectral trace. + + This returns vectors perpendicular to the dispersion direction + at each sampled point. + + Returns + ------- + nx, ny : ndarray + Components of the normal vectors. + """ + dxdl, dydl = self.deriv(*args, **kwargs) + return dydl, -dxdl + + def dispersion(self, *args, **kwargs): + """ + Compute the local dispersion scale (dλ/dr) along the spectral trace. + + Returns + ------- + dldr : ndarray + Local dispersion (wavelength per pixel). + """ + dxdl, dydl = self.deriv(*args, **kwargs) + dldr = 1.0 / np.hypot(dxdl, dydl) + return dldr + + +if __name__ == "__main__": + + # instantiate the disperser + disp = SNPITDisperser("Roman_prism_OpticalModel_v0.8.yaml") + + # which SCA are we considering? + sca = 1 + + # length of test data + N = 1000 + + # create some random positions + x = np.random.uniform(low=0, high=4088, size=N) + y = np.random.uniform(low=0, high=4088, size=N) + + # assume one wavelength (in micron) for each position + lam = np.random.uniform(low=disp.optical["wl_min"], high=disp.optical["wl_max"], size=N) + + # pair-wise compute the dispersed positions: here, every (x, y) + # has exactly one corresponding wavelength (l), so that this maps the + # tuples: (x,y,l) -> (xp, yp). + xp, yp = disp.disperse(x, y, lam, sca, order="1", pairwise=True) + + # can also compute the derivatives, which is the tangent vector + # at some point: + dxdl, dydl = disp.deriv(x, y, lam, sca, order="1", pairwise=True) + + # can compute the dispersion (the rate of change of the wavelength + # on a path along the trace: + dldr = disp.dispersion(x, y, lam, sca, order="1", pairwise=True) + + # now demo the non-pairwise. Here if (x,y) are vectors of size + # N elements and lam is a vector of M elements, now the output data + # will be a 2d array of (N, M) elements + + # create a new wavelength grid + lam = np.linspace(disp.optical["wl_min"], disp.optical["wl_max"], 100) + + # new dispersed positions + xp, yp = disp.disperse(x, y, lam, sca, order="1", pairwise=False) + + # new tangent vectors + dxdl, dydl = disp.deriv(x, y, lam, sca, order="1", pairwise=False) + + # new dispersion + dldr = disp.dispersion(x, y, lam, sca, order="1", pairwise=False) diff --git a/roman_imsim/skycat.py b/roman_imsim/skycat.py index 4b337518..3bab3853 100644 --- a/roman_imsim/skycat.py +++ b/roman_imsim/skycat.py @@ -88,9 +88,9 @@ def __init__( @property def objects(self): - from skycatalogs import skyCatalogs - from skycatalogs import __version__ as skycatalogs_version from packaging.version import Version + from skycatalogs import __version__ as skycatalogs_version + from skycatalogs import skyCatalogs if Version(skycatalogs_version) < Version("2.0"): PolygonalRegion = skyCatalogs.PolygonalRegion @@ -183,9 +183,16 @@ def getObj(self, index, gsparams=None, rng=None, exptime=30): gsobjs = skycat_obj.get_gsobject_components(gsparams) # Compute the flux or get the cached value. - flux = ( - skycat_obj.get_roman_flux(self.bandpass.name, mjd=self.mjd) * self.exptime * roman.collecting_area - ) + try: + flux = ( + skycat_obj.get_roman_flux(self.bandpass.name, mjd=self.mjd) + * self.exptime + * roman.collecting_area + ) + except KeyError: + # This happens because SNANA knows nothing about SNPrism yet + self.logger.warning("Switching to H158 from %s for %s", self.bandpass.name, skycat_obj) + flux = skycat_obj.get_roman_flux("H158", mjd=self.mjd) * self.exptime * roman.collecting_area if np.isnan(flux): return None diff --git a/roman_imsim/stamp.py b/roman_imsim/stamp.py index bc50d4c3..64b062bd 100644 --- a/roman_imsim/stamp.py +++ b/roman_imsim/stamp.py @@ -60,6 +60,19 @@ def setup(self, config, base, xsize, ysize, ignore, logger): # or cached by the skyCatalogs code. gal.flux = gal.calculateFlux(bandpass) self.flux = gal.flux + # Cap (star) flux at 30M photons to avoid gross artifacts when trying + # to draw the Roman PSF in finite time and memory. + flux_cap = config.get("flux_cap", 3e7) + if self.flux > flux_cap: + if ( + hasattr(gal, "original") + and hasattr(gal.original, "original") + and isinstance(gal.original.original, galsim.DeltaFunction) + ) or (isinstance(gal, galsim.DeltaFunction)): + gal = gal.withFlux(flux_cap, bandpass) + self.flux = flux_cap + gal.flux = flux_cap + logger.info("Limiting the flux of object %d to %e", base["obj_num"], flux_cap) base["flux"] = gal.flux base["mag"] = -2.5 * np.log10(gal.flux) + bandpass.zeropoint # print('stamp setup2',process.memory_info().rss) @@ -99,8 +112,9 @@ def setup(self, config, base, xsize, ysize, ignore, logger): # # Get storead achromatic PSF # psf = galsim.config.BuildGSObject(base, 'psf', logger=logger)[0]['achromatic'] # obj = galsim.Convolve(gal_achrom, psf).withFlux(self.flux) - obj = gal_achrom.withGSParams(galsim.GSParams(stepk_minimum_hlr=20)) + obj = gal_achrom.withGSParams(galsim.GSParams(stepk_minimum_hlr=20, folding_threshold=0.001)) image_size = obj.getGoodImageSize(roman.pixel_scale) + image_size = max(image_size, config.get("min_size", 0)) # print('stamp setup3',process.memory_info().rss) base["pupil_bin"] = self.pupil_bin @@ -346,7 +360,6 @@ def draw(self, prof, image, method, offset, config, base, logger): # In case we had to make a bigger image, just copy the part we need. image += fft_image[image.bounds] # print('stamp draw3',process.memory_info().rss) - return image def add_poisson_noise(self, fft_image): diff --git a/tests/test_photonOps.py b/tests/test_photonOps.py new file mode 100644 index 00000000..eadae841 --- /dev/null +++ b/tests/test_photonOps.py @@ -0,0 +1,45 @@ +import os +import unittest + +# Import galsim and roman_imsim +import galsim +import numpy as np + +from roman_imsim.photonOps import ChargeDiff, SlitlessSpec, GrismNV, GrismV, WFSSSDisperser + + +class TestPhotonOps(unittest.TestCase): + """Unit tests for all registered photon operators.""" + + def setUp(self): + """Set up test fixtures.""" + # Get the root directory of the package + self.root_dir = os.path.dirname(os.path.dirname(__file__)) + + # Create a mock photon array + self.num_photons = 10 + self.mock_photon_array = self._create_mock_photon_array() + + def _create_mock_photon_array(self): + photon_array = galsim.PhotonArray(self.num_photons) + + photon_array.x = np.random.uniform(-5, 5, self.num_photons) + photon_array.y = np.random.uniform(-5, 5, self.num_photons) + + photon_array.flux = np.random.uniform(100, 1000, self.num_photons) + + photon_array.wavelength = np.random.uniform(500, 2000, self.num_photons) + + return photon_array + + def test_photon_operators_smoke(self): + """Test that all registered roman_imsim photon operators can be initialized and applied.""" + photon_operators = [ChargeDiff, SlitlessSpec, GrismNV, GrismV, WFSSSDisperser] + for photon_operator_cls in photon_operators: + with self.subTest(operator=photon_operator_cls.__name__): + photon_op = photon_operator_cls() + photon_op.applyTo(self.mock_photon_array) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/test_snpitdisperser.py b/tests/test_snpitdisperser.py new file mode 100644 index 00000000..eb969506 --- /dev/null +++ b/tests/test_snpitdisperser.py @@ -0,0 +1,389 @@ +import os +import tempfile +import unittest + +import numpy as np +import yaml + +import roman_imsim.photonOps._wfss_disperser.snpitdispenser + + +class TestSNPITDisperser(unittest.TestCase): + """Unit tests for the SNPITDisperser class and related components.""" + + def setUp(self): + """Set up test fixtures with a minimal valid configuration.""" + + # Initialize disperser + root_dir = os.path.dirname(roman_imsim.__path__[0]) + self.config_file_name = os.path.join(root_dir, "optical_models/Roman_prism_OpticalModel_v0.8.yaml") + self.disperser = roman_imsim.snpitdispenser.SNPITDisperser( + os.path.join(root_dir, "optical_models/Roman_prism_OpticalModel_v0.8.yaml") + ) + + def test_load_conffile_invalid_transform(self): + """Test that invalid wavelength transform raises error.""" + + with open(self.config_file_name, "r") as fp: + config = yaml.safe_load(fp) + config["optical_model"]["wl_transform"] = "invalid" + + temp_file = tempfile.NamedTemporaryFile(mode="w", suffix=".yaml", delete=False) + yaml.dump(config, temp_file) + temp_file.close() + + with self.assertRaises(NotImplementedError): + roman_imsim.snpitdispenser.SNPITDisperser(temp_file.name) + + os.unlink(temp_file.name) + + def test_deriv_coeffs(self): + """Test polynomial coefficient derivative computation.""" + # Test 1D case + M = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) + expected = np.array([[4, 5, 6], [14, 16, 18]]) # [1*4, 2*5, 3*7], [2*7, 2*8, 2*9] + result = roman_imsim.snpitdispenser.SNPITDisperser.deriv_coeffs(M) + np.testing.assert_array_equal(result, expected) + + # Test 3D case + M_3d = np.random.rand(4, 3, 2) + result_3d = roman_imsim.snpitdispenser.SNPITDisperser.deriv_coeffs(M_3d) + expected_shape = (3, 3, 2) + self.assertEqual(result_3d.shape, expected_shape) + + def test_validate_inputs_valid(self): + """Test input validation with valid inputs.""" + x = np.array([100, 200]) + y = np.array([150, 250]) + lam = np.array([1.0, 1.5]) + sca = 1 + + # Should not raise any exceptions + self.disperser.validate(x, y, lam, sca, order="1", pairwise=True) + self.disperser.validate(x, y, lam, sca, order="1", pairwise=False) + + def test_validate_inputs_invalid_order(self): + """Test input validation with invalid spectral order.""" + x = np.array([100, 180]) + y = np.array([120, 200]) + lam = np.array([1.0, 1.5]) + sca = 1 + + with self.assertRaises(KeyError): + self.disperser.validate(x, y, lam, sca, order="2", pairwise=True) + + def test_validate_inputs_invalid_sca(self): + """Test input validation with invalid SCA number.""" + x = np.array([110, 200]) + y = np.array([100, 220]) + lam = np.array([1.0, 1.5]) + sca = 99 # Invalid SCA + + with self.assertRaises(KeyError): + self.disperser.validate(x, y, lam, sca, order="1", pairwise=True) + + def test_validate_inputs_pairwise_shape_mismatch(self): + """Test input validation with shape mismatch in pairwise mode.""" + x = np.array([105, 200]) + y = np.array([100, 210, 300]) # Different shape + lam = np.array([1.0, 1.5]) + sca = 1 + + with self.assertRaises(RuntimeError): + self.disperser.validate(x, y, lam, sca, order="1", pairwise=True) + + def test_validate_inputs_non_pairwise_shape_mismatch(self): + """Test input validation with shape mismatch in non-pairwise mode.""" + x = np.array([133, 245]) + y = np.array([101, 205, 303]) # Different shape + lam = np.array([1.0, 1.2]) # Should be 1D + sca = 1 + + with self.assertRaises(RuntimeError): + self.disperser.validate(x, y, lam, sca, order="1", pairwise=False) + + def test_sca_to_fpa_conversion(self): + """Test SCA to FPA coordinate conversion.""" + xsca = np.array([2044.5]) + ysca = np.array([2044.5]) + sca = 1 + + xfpa, yfpa = self.disperser.sca_to_fpa(xsca, ysca, sca) + + # Should return arrays + self.assertIsInstance(xfpa, np.ndarray) + self.assertIsInstance(yfpa, np.ndarray) + + # Should have same length as input + self.assertEqual(len(xfpa), len(xsca)) + self.assertEqual(len(yfpa), len(ysca)) + + def test_mpa_to_sca_conversion(self): + """Test MPA to SCA coordinate conversion.""" + xmpa = np.array([0.0]) + ympa = np.array([0.0]) + sca = 1 + + xsca, ysca = self.disperser.mpa_to_sca(xmpa, ympa, sca) + + # Should return arrays + self.assertIsInstance(xsca, np.ndarray) + self.assertIsInstance(ysca, np.ndarray) + + # Should have same length as input + self.assertEqual(len(xsca), len(xmpa)) + self.assertEqual(len(ysca), len(ympa)) + + @unittest.skip("Unclear if this should pass.") + def test_coordinate_conversion_roundtrip(self): + """Test that coordinate conversion is approximately reversible.""" + xsca_orig = np.array([1234.0, 2000.0]) + ysca_orig = np.array([1000.0, 2567.0]) + sca = 1 + + # Convert SCA -> FPA -> SCA + xfpa, yfpa = self.disperser.sca_to_fpa(xsca_orig, ysca_orig, sca) + xsca_back, ysca_back = self.disperser.mpa_to_sca(xfpa, yfpa, sca) + + # Should be approximately equal (allowing for numerical precision) + np.testing.assert_allclose(xsca_orig, xsca_back, rtol=1e-10) + np.testing.assert_allclose(ysca_orig, ysca_back, rtol=1e-10) + + def test_disperse_pairwise(self): + """Test dispersion computation in pairwise mode.""" + x = np.array([1020.0, 2048.0]) + y = np.array([2040.0, 1000.0]) + lam = np.array([1.0, 1.5]) + sca = 1 + + xp, yp = self.disperser.disperse(x, y, lam, sca, order="1", pairwise=True) + + # Should return arrays + self.assertIsInstance(xp, np.ndarray) + self.assertIsInstance(yp, np.ndarray) + + # Should have same shape as input + self.assertEqual(xp.shape, x.shape) + self.assertEqual(yp.shape, y.shape) + + # Should be different from input (dispersion should change positions) + self.assertFalse(np.allclose(xp, x)) + self.assertFalse(np.allclose(yp, y)) + + def test_disperse_non_pairwise(self): + """Test dispersion computation in non-pairwise mode.""" + x = np.array([1000.0, 2000.0]) + y = np.array([1000.0, 2000.0]) + lam = np.array([1.0, 1.5]) # 1D wavelength array + sca = 1 + + xp, yp = self.disperser.disperse(x, y, lam, sca, order="1", pairwise=False) + + # Should return 2D arrays + self.assertEqual(xp.ndim, 2) + self.assertEqual(yp.ndim, 2) + + # Should have shape (len(x), len(lam)) + expected_shape = (len(x), len(lam)) + self.assertEqual(xp.shape, expected_shape) + self.assertEqual(yp.shape, expected_shape) + + def test_deriv_pairwise(self): + """Test derivative computation in pairwise mode.""" + x = np.array([1000.0, 2000.0]) + y = np.array([1000.0, 2000.0]) + lam = np.array([1.0, 1.5]) + sca = 1 + + dxdl, dydl = self.disperser.deriv(x, y, lam, sca, order="1", pairwise=True) + + # Should return arrays + self.assertIsInstance(dxdl, np.ndarray) + self.assertIsInstance(dydl, np.ndarray) + + # Should have same shape as input + self.assertEqual(dxdl.shape, x.shape) + self.assertEqual(dydl.shape, y.shape) + + # Should be non-zero (dispersion should have non-zero derivative) + self.assertFalse(np.allclose(dxdl, 0)) + self.assertFalse(np.allclose(dydl, 0)) + + def test_deriv_non_pairwise(self): + """Test derivative computation in non-pairwise mode.""" + x = np.array([1000.0, 2000.0]) + y = np.array([1000.0, 2000.0]) + lam = np.array([1.0, 1.5]) # 1D wavelength array + sca = 1 + + dxdl, dydl = self.disperser.deriv(x, y, lam, sca, order="1", pairwise=False) + + # Should return 2D arrays + self.assertEqual(dxdl.ndim, 2) + self.assertEqual(dydl.ndim, 2) + + # Should have shape (len(x), len(lam)) + expected_shape = (len(x), len(lam)) + self.assertEqual(dxdl.shape, expected_shape) + self.assertEqual(dydl.shape, expected_shape) + + def test_normal(self): + """Test normal vector computation.""" + x = np.array([1000.0, 2000.0]) + y = np.array([1000.0, 2000.0]) + lam = np.array([1.0, 1.5]) + sca = 1 + + nx, ny = self.disperser.normal(x, y, lam, sca, order="1", pairwise=True) + + # Should return arrays + self.assertIsInstance(nx, np.ndarray) + self.assertIsInstance(ny, np.ndarray) + + # Should have same shape as input + self.assertEqual(nx.shape, x.shape) + self.assertEqual(ny.shape, y.shape) + + # Normal vectors should be perpendicular to derivative vectors + dxdl, dydl = self.disperser.deriv(x, y, lam, sca, order="1", pairwise=True) + dot_product = dxdl * nx + dydl * ny + np.testing.assert_allclose(dot_product, 0, atol=1e-10) + + def test_dispersion(self): + """Test dispersion scale computation.""" + x = np.array([1000.0, 2000.0]) + y = np.array([1000.0, 2000.0]) + lam = np.array([1.0, 1.5]) + sca = 1 + + dldr = self.disperser.dispersion(x, y, lam, sca, order="1", pairwise=True) + + # Should return array + self.assertIsInstance(dldr, np.ndarray) + + # Should have same shape as input + self.assertEqual(dldr.shape, x.shape) + + # Should be positive (dispersion scale should be positive) + self.assertTrue(np.all(dldr > 0)) + + def test_scalar_inputs(self): + """Test that scalar inputs are handled correctly.""" + x = 1000.0 + y = 1000.0 + lam = 1.0 + sca = 1 + + # Should work with scalar inputs + xp, yp = self.disperser.disperse(x, y, lam, sca, order="1", pairwise=True) + dxdl, dydl = self.disperser.deriv(x, y, lam, sca, order="1", pairwise=True) + + # Should return 1D arrays + self.assertEqual(xp.ndim, 1) + self.assertEqual(yp.ndim, 1) + self.assertEqual(dxdl.ndim, 1) + self.assertEqual(dydl.ndim, 1) + + # Should have length 1 + self.assertEqual(len(xp), 1) + self.assertEqual(len(yp), 1) + self.assertEqual(len(dxdl), 1) + self.assertEqual(len(dydl), 1) + + +class TestLogTransformer(unittest.TestCase): + """Unit tests for the LogTransformer class.""" + + def setUp(self): + """Set up test fixtures.""" + self.lam0 = 1.0 + self.transformer = roman_imsim.snpitdispenser.LogTransformer(self.lam0) + + def test_initialization(self): + """Test transformer initialization.""" + self.assertEqual(self.transformer.lam0, self.lam0) + self.assertAlmostEqual(self.transformer.ln10, np.log(10.0)) + + def test_evaluate(self): + """Test wavelength transformation.""" + lam = np.array([0.5, 1.0, 2.0]) + w = self.transformer.evaluate(lam) + + expected = np.log10(lam / self.lam0) + np.testing.assert_allclose(w, expected) + + def test_invert(self): + """Test inverse transformation.""" + w = np.array([-0.3, 0.0, 0.3]) + lam = self.transformer.invert(w) + + expected = self.lam0 * (10.0**w) + np.testing.assert_allclose(lam, expected) + + def test_roundtrip(self): + """Test that evaluate and invert are inverses.""" + lam_orig = np.array([0.5, 1.0, 2.0]) + + w = self.transformer.evaluate(lam_orig) + lam_back = self.transformer.invert(w) + + np.testing.assert_allclose(lam_orig, lam_back, rtol=1e-10) + + def test_deriv(self): + """Test derivative computation.""" + lam = np.array([0.5, 1.0, 2.0]) + dldw = self.transformer.deriv(lam) + + expected = lam * np.log(10.0) + np.testing.assert_allclose(dldw, expected) + + +class TestLinearTransformer(unittest.TestCase): + """Unit tests for the LinearTransformer class.""" + + def setUp(self): + """Set up test fixtures.""" + self.lam0 = 1.0 + self.transformer = roman_imsim.snpitdispenser.LinearTransformer(self.lam0) + + def test_initialization(self): + """Test transformer initialization.""" + self.assertEqual(self.transformer.lam0, self.lam0) + + def test_evaluate(self): + """Test wavelength transformation.""" + lam = np.array([0.5, 1.0, 2.0]) + w = self.transformer.evaluate(lam) + + expected = lam - self.lam0 + np.testing.assert_allclose(w, expected) + + def test_invert(self): + """Test inverse transformation.""" + w = np.array([-0.5, 0.0, 1.0]) + lam = self.transformer.invert(w) + + expected = w + self.lam0 + np.testing.assert_allclose(lam, expected) + + def test_roundtrip(self): + """Test that evaluate and invert are inverses.""" + lam_orig = np.array([0.5, 1.0, 2.0]) + + w = self.transformer.evaluate(lam_orig) + lam_back = self.transformer.invert(w) + + np.testing.assert_allclose(lam_orig, lam_back, rtol=1e-10) + + def test_deriv(self): + """Test derivative computation.""" + lam = np.array([0.5, 1.0, 2.0]) + dldw = self.transformer.deriv(lam) + + # Derivative should always be 1 for linear transform + expected = 1.0 + np.testing.assert_allclose(dldw, expected) + + +if __name__ == "__main__": + unittest.main()