diff --git a/examples/2-clearsky-radiative-transfer/2-zeeman/10-zeeman-linear-prop.py b/examples/2-clearsky-radiative-transfer/2-zeeman/10-zeeman-linear-prop.py index 0e6ae8612b..3d185bf671 100644 --- a/examples/2-clearsky-radiative-transfer/2-zeeman/10-zeeman-linear-prop.py +++ b/examples/2-clearsky-radiative-transfer/2-zeeman/10-zeeman-linear-prop.py @@ -57,16 +57,16 @@ assert np.allclose( ws.spectral_rad[::100], np.array( - [[2.27853388e+02, 1.95866945e-04, 4.87492271e-05, -2.36891567e-02], - [2.30930417e+02, 2.15802433e-04, 5.51980780e-05, -4.03735389e-02], - [2.34870122e+02, 3.57230883e-04, 9.38282269e-05, -5.71986370e-02], - [2.40419686e+02, 1.23779635e-03, 3.24558457e-04, -6.19280114e-02], - [2.49835914e+02, 7.74686776e-03, 2.06857040e-03, -1.54924878e-02], - [2.08998520e+02, 2.38211765e+01, 1.64703709e+00, 5.47032904e-06], - [2.49835352e+02, 7.74672684e-03, 2.06868749e-03, 1.55957453e-02], - [2.40418596e+02, 1.23700264e-03, 3.24403765e-04, 6.20718269e-02], - [2.34868546e+02, 3.56494587e-04, 9.36667274e-05, 5.73507670e-02], - [2.30928366e+02, 2.15286025e-04, 5.50826926e-05, 4.05218001e-02], - [2.27850872e+02, 1.95567388e-04, 4.86829472e-05, 2.38269651e-02]] + [[ 2.27720431e+02, 1.92860863e-04, 4.80226318e-05, -2.39944534e-02], + [ 2.30795788e+02, 2.05916096e-04, 5.28115272e-05, -4.11877682e-02], + [ 2.34732775e+02, 6.07407437e-04, 1.53649762e-04, -3.35111787e-02], + [ 2.40285632e+02, 1.20481776e-03, 3.16497792e-04, -6.37947118e-02], + [ 2.49706649e+02, 5.90915798e-03, 1.64151310e-03, -9.03484330e-02], + [ 2.09884932e+02, 2.41682141e+01, 1.73723873e+00, 5.54512226e-06], + [ 2.49706087e+02, 5.90803897e-03, 1.64139792e-03, 9.04789486e-02], + [ 2.40284543e+02, 1.20399960e-03, 3.16336925e-04, 6.39391756e-02], + [ 2.34731202e+02, 6.06980643e-04, 1.53562668e-04, 3.36445494e-02], + [ 2.30793741e+02, 2.05389194e-04, 5.26935796e-05, 4.13367589e-02], + [ 2.27717921e+02, 1.92557235e-04, 4.79553616e-05, 2.41329191e-02]] ), ), "Values have drifted from expected results in spectral radiance" diff --git a/examples/2-clearsky-radiative-transfer/2-zeeman/2-zeeman.py b/examples/2-clearsky-radiative-transfer/2-zeeman/2-zeeman.py index 80e9be8a33..ba846b0882 100644 --- a/examples/2-clearsky-radiative-transfer/2-zeeman/2-zeeman.py +++ b/examples/2-clearsky-radiative-transfer/2-zeeman/2-zeeman.py @@ -57,16 +57,16 @@ assert np.allclose( ws.spectral_rad[::100], np.array( - [[2.27784836e+02, 4.26114598e-04, 1.02751718e-04, 5.69266704e-02], - [2.30863855e+02, 6.59892211e-04, 1.59299017e-04, 7.04138026e-02], - [2.34806524e+02, 1.16821510e-03, 2.82508836e-04, 9.33835996e-02], - [2.40360807e+02, 2.61610579e-03, 6.34821344e-04, 1.40029927e-01], - [2.49782423e+02, 9.74658804e-03, 2.38967070e-03, 2.69735134e-01], - [2.09027223e+02, 2.38234847e+01, 1.64634819e+00, 5.44981208e-06], - [2.49781856e+02, 9.75012917e-03, 2.39056843e-03, -2.69812587e-01], - [2.40359714e+02, 2.61791138e-03, 6.35266868e-04, -1.40106505e-01], - [2.34804947e+02, 1.16941576e-03, 2.82802436e-04, -9.34596007e-02], - [2.30861804e+02, 6.60797190e-04, 1.59519308e-04, -7.04898138e-02], - [2.27782319e+02, 4.26846028e-04, 1.02929281e-04, -5.70026325e-02]] + [[ 2.27651281e+02, 4.26102083e-04, 1.02747359e-04, 5.68792774e-02], + [ 2.30728865e+02, 6.60200934e-04, 1.59372202e-04, 7.04074781e-02], + [ 2.34671710e+02, 1.16870319e-03, 2.82624430e-04, 9.34017504e-02], + [ 2.40226519e+02, 2.61619541e-03, 6.34833393e-04, 1.40041579e-01], + [ 2.49649966e+02, 9.75456580e-03, 2.39140821e-03, 2.69997496e-01], + [ 2.09913635e+02, 2.41699652e+01, 1.73644537e+00, 5.52423749e-06], + [ 2.49649399e+02, 9.75810766e-03, 2.39230583e-03, -2.70074955e-01], + [ 2.40225428e+02, 2.61800091e-03, 6.35278864e-04, -1.40118138e-01], + [ 2.34670135e+02, 1.16990396e-03, 2.82918050e-04, -9.34777298e-02], + [ 2.30726817e+02, 6.61105654e-04, 1.59592427e-04, -7.04834283e-02], + [ 2.27648771e+02, 4.26832634e-04, 1.02924704e-04, -5.69551120e-02]] ), ), "Values have drifted from expected results in spectral radiance" diff --git a/examples/2-clearsky-radiative-transfer/2-zeeman/3-zeeman-sensor.py b/examples/2-clearsky-radiative-transfer/2-zeeman/3-zeeman-sensor.py index 00dcc92f16..ad0d2cd685 100644 --- a/examples/2-clearsky-radiative-transfer/2-zeeman/3-zeeman-sensor.py +++ b/examples/2-clearsky-radiative-transfer/2-zeeman/3-zeeman-sensor.py @@ -60,8 +60,8 @@ assert np.allclose( ws.measurement_vec[::100], np.array( - [227.85626444, 230.93431141, 234.89998492, 240.50101182, - 250.05274234, 210.84064948, 249.51259663, 240.21977571, - 234.7115623, 230.79135556, 227.73970923] + [227.72265043, 230.79931459, 234.76518855, 240.36673613, + 249.92055137, 211.84612808, 249.3798804 , 240.08547841, + 234.5767321 , 230.65637527, 227.60619589] ), ) diff --git a/examples/2-clearsky-radiative-transfer/2-zeeman/4-zeeman-transmission.py b/examples/2-clearsky-radiative-transfer/2-zeeman/4-zeeman-transmission.py index 74821249a6..93171b41b9 100644 --- a/examples/2-clearsky-radiative-transfer/2-zeeman/4-zeeman-transmission.py +++ b/examples/2-clearsky-radiative-transfer/2-zeeman/4-zeeman-transmission.py @@ -57,16 +57,16 @@ assert np.allclose( ws.spectral_rad[::100], np.array( - [[3.48824591e-06, -1.47203439e-10, -3.54155511e-11, -4.59223182e-08], - [1.71707407e-06, -1.11587582e-10, -2.68557096e-11, -2.78282572e-08], - [6.98942589e-07, -7.92932217e-11, -1.90895244e-11, -1.48129013e-08], - [2.02345668e-07, -5.04729873e-11, -1.21507664e-11, -6.26440637e-09], - [2.59661793e-08, -2.51846307e-11, -6.03516645e-12, -1.54349340e-09], - [6.89347083e-13, 6.83290031e-13, -5.69298697e-14, 7.52268204e-18], - [2.59895049e-08, -2.52264844e-11, -6.04518733e-12, 1.54841075e-09], - [2.02661517e-07, -5.06298315e-11, -1.21887251e-11, 6.29847485e-09], - [7.00441394e-07, -7.96478027e-11, -1.91752558e-11, 1.49217675e-08], - [1.72166873e-06, -1.12232853e-10, -2.70115093e-11, 2.80800108e-08], - [3.49929916e-06, -1.48242055e-10, -3.56660238e-11, 4.64085834e-08]] + [[ 3.55416974e-06, -1.49626896e-10, -3.59945621e-11, -4.65675860e-08], + [ 1.75504614e-06, -1.14112664e-10, -2.74625108e-11, -2.83292389e-08], + [ 7.16577429e-07, -8.14996162e-11, -1.96228377e-11, -1.51461032e-08], + [ 2.07916647e-07, -5.19869800e-11, -1.25199070e-11, -6.43052931e-09], + [ 2.67056326e-08, -2.58535066e-11, -6.20359298e-12, -1.58692924e-09], + [ 7.87383592e-13, 7.80210615e-13, -6.46005354e-14, 8.84600141e-18], + [ 2.67295620e-08, -2.58964360e-11, -6.21391407e-12, 1.59197636e-09], + [ 2.08240338e-07, -5.21480127e-11, -1.25589127e-11, 6.46545336e-09], + [ 7.18110068e-07, -8.18625471e-11, -1.97105925e-11, 1.52574035e-08], + [ 1.75973088e-06, -1.14770061e-10, -2.76211969e-11, 2.85857722e-08], + [ 3.56540626e-06, -1.50680035e-10, -3.62484502e-11, 4.70615011e-08]] ), ), "Values have drifted from expected results in spectral radiance" diff --git a/examples/2-clearsky-radiative-transfer/2-zeeman/5-zeeman-sun.py b/examples/2-clearsky-radiative-transfer/2-zeeman/5-zeeman-sun.py index 9c459ff836..e3b17eef37 100644 --- a/examples/2-clearsky-radiative-transfer/2-zeeman/5-zeeman-sun.py +++ b/examples/2-clearsky-radiative-transfer/2-zeeman/5-zeeman-sun.py @@ -66,11 +66,11 @@ assert np.allclose( res[::3, ::7], np.array( - [[5333.26267148, 5333.26267148, 5333.26267148], - [ 17.53700177, 17.41842254, 17.45306884], - [ 17.60721343, 17.37046458, 17.44008522], - [ 17.67721853, 17.32271259, 17.4276346 ], - [ 17.7470159 , 17.27516828, 17.41571773], - [ 17.8166044 , 17.22783343, 17.40433539], - [ 17.88598287, 17.18070982, 17.39348841]]) + [[5356.88051075, 5356.88051075, 5356.88051075], + [ 16.77380467, 16.66091577, 16.69389844], + [ 16.84063115, 16.61524188, 16.68151802], + [ 16.90725049, 16.56975243, 16.6696318 ], + [ 16.9736614 , 16.52444918, 16.65824055], + [ 17.03986264, 16.47933394, 16.64734506], + [ 17.10585295, 16.43440858, 16.63694619]]) ) diff --git a/examples/2-clearsky-radiative-transfer/2-zeeman/6-zeeman-sun-scattering.py b/examples/2-clearsky-radiative-transfer/2-zeeman/6-zeeman-sun-scattering.py index da7811ad56..fe0228a8ea 100644 --- a/examples/2-clearsky-radiative-transfer/2-zeeman/6-zeeman-sun-scattering.py +++ b/examples/2-clearsky-radiative-transfer/2-zeeman/6-zeeman-sun-scattering.py @@ -66,11 +66,11 @@ assert np.allclose( res[1::3, 1::7], - np.array([[5771.99991010, 5771.99991010, 5771.99991010], - [643.69008277, 643.69008277, 643.69008277], - [643.69030757, 643.69030764, 643.69030764], - [643.69065708, 643.69065715, 643.69065728], - [643.69113428, 643.69113428, 643.69113441], - [643.69174257, 643.69174270, 643.69174290], - [643.69248695, 643.69248708, 643.69248747]]), + np.array([[5771.9999155 , 5771.9999155 , 5771.9999155 ], + [ 642.44289971, 642.44289971, 642.44289971], + [ 642.44286517, 642.44286517, 642.44286524], + [ 642.4428158 , 642.4428158 , 642.44281587], + [ 642.44275627, 642.44275634, 642.44275648], + [ 642.44269338, 642.44269352, 642.44269372], + [ 642.44263584, 642.44263598, 642.44263625]]), ) diff --git a/examples/2-clearsky-radiative-transfer/2-zeeman/7-zeeman-refractive.py b/examples/2-clearsky-radiative-transfer/2-zeeman/7-zeeman-refractive.py index fb46cf6dfb..8fb9126d08 100644 --- a/examples/2-clearsky-radiative-transfer/2-zeeman/7-zeeman-refractive.py +++ b/examples/2-clearsky-radiative-transfer/2-zeeman/7-zeeman-refractive.py @@ -85,12 +85,12 @@ def single_propmat_agenda(ws): # %% Tests -assert np.allclose(geometric.flatten()[::21], [ 1.25790011e+02, -5.83560523e-04, -2.57441702e-03, 6.19275327e-01, - 1.38778854e+02, 1.03394079e-01, -2.12544902e-02, -5.28842204e-01, - 1.29714406e+02, -4.56150073e-04]) +assert np.allclose(geometric.flatten()[::21], [1.25486103e+02, -5.56777144e-04, -2.63816009e-03, 6.12431001e-01, + 1.38429184e+02, 1.04834522e-01, -2.21549573e-02, -5.18205673e-01, + 1.29365186e+02, -4.48679095e-04]) -assert np.allclose(refractive.flatten()[::21], [ 1.25852437e+02, -5.89498878e-04, -2.56133152e-03, 6.19139312e-01, - 1.38846218e+02, 1.03204818e-01, -2.10977115e-02, -5.31745645e-01, - 1.29755973e+02, -4.48250955e-04]) +assert np.allclose(refractive.flatten()[::21], [1.25552512e+02, -5.61580828e-04, -2.61219720e-03, 6.11932940e-01, + 1.38488011e+02, 1.04448633e-01, -2.18946938e-02, -5.19985291e-01, + 1.29439037e+02, -4.53301859e-04]) assert not np.allclose(geometric, refractive) diff --git a/examples/2-clearsky-radiative-transfer/2-zeeman/8-zeeman-adaptive-path.py b/examples/2-clearsky-radiative-transfer/2-zeeman/8-zeeman-adaptive-path.py index 0939c76684..0365b33452 100644 --- a/examples/2-clearsky-radiative-transfer/2-zeeman/8-zeeman-adaptive-path.py +++ b/examples/2-clearsky-radiative-transfer/2-zeeman/8-zeeman-adaptive-path.py @@ -51,7 +51,7 @@ srad0 = copy(ws.spectral_rad) path0 = copy(ws.ray_path) -ws.spectral_propmat_pathAdaptiveHalfPath( +ws.spectral_propmat_pathAddAdaptiveHalfPath( max_stepsize=100., max_tau=0.05, cutoff_tau=3.0) ws.spectral_radSetToBackground() ws.spectral_radSinglePathEmissionFrequencyLoop() diff --git a/examples/2-clearsky-radiative-transfer/2-zeeman/9-zeeman-linear-src.py b/examples/2-clearsky-radiative-transfer/2-zeeman/9-zeeman-linear-src.py index 1aac19834d..f1b88abbcd 100644 --- a/examples/2-clearsky-radiative-transfer/2-zeeman/9-zeeman-linear-src.py +++ b/examples/2-clearsky-radiative-transfer/2-zeeman/9-zeeman-linear-src.py @@ -57,16 +57,16 @@ assert np.allclose( ws.spectral_rad[::100], np.array( - [[2.27827137e+02, 4.25383355e-04, 1.02572576e-04, 5.68786105e-02], - [2.30903675e+02, 6.58907710e-04, 1.59056181e-04, 7.03618286e-02], - [2.34843808e+02, 1.16673733e-03, 2.82139773e-04, 9.33333640e-02], - [2.40395689e+02, 2.61325977e-03, 6.34093204e-04, 1.39982788e-01], - [2.49814591e+02, 9.73578780e-03, 2.38675267e-03, 2.69659401e-01], - [2.08998520e+02, 2.38211765e+01, 1.64703709e+00, 5.47032781e-06], - [2.49814023e+02, 9.73932829e-03, 2.38765002e-03, -2.69736936e-01], - [2.40394595e+02, 2.61506483e-03, 6.34538536e-04, -1.40059437e-01], - [2.34842230e+02, 1.16793738e-03, 2.82433194e-04, -9.34094305e-02], - [2.30901623e+02, 6.59812051e-04, 1.59276300e-04, -7.04379152e-02], - [2.27824619e+02, 4.26114226e-04, 1.02749992e-04, -5.69546795e-02]] + [[ 2.27693775e+02, 4.25360933e-04, 1.02565822e-04, 5.68301097e-02], + [ 2.30768818e+02, 6.59206212e-04, 1.59126898e-04, 7.03542474e-02], + [ 2.34709071e+02, 1.16721699e-03, 2.82253340e-04, 9.33505189e-02], + [ 2.40261444e+02, 2.61334904e-03, 6.34105246e-04, 1.39994110e-01], + [ 2.49682208e+02, 9.74388124e-03, 2.38851950e-03, 2.69924367e-01], + [ 2.09884932e+02, 2.41682141e+01, 1.73723873e+00, 5.54512081e-06], + [ 2.49681641e+02, 9.74742248e-03, 2.38941674e-03, -2.70001907e-01], + [ 2.40260352e+02, 2.61515400e-03, 6.34550524e-04, -1.40070740e-01], + [ 2.34707495e+02, 1.16841714e-03, 2.82546778e-04, -9.34265633e-02], + [ 2.30766769e+02, 6.60110278e-04, 1.59346947e-04, -7.04302732e-02], + [ 2.27691263e+02, 4.26090913e-04, 1.02743016e-04, -5.69060525e-02]] ), ), "Values have drifted from expected results in spectral radiance" diff --git a/examples/6-satellite-sensors/amsu-a.py b/examples/6-satellite-sensors/amsu-a.py index 2542684519..2ddaf9f213 100644 --- a/examples/6-satellite-sensors/amsu-a.py +++ b/examples/6-satellite-sensors/amsu-a.py @@ -57,14 +57,14 @@ def sensor_channels(channels, n): for ch in channels: # The range is [f_ref, 0], [f_ref, inf] - x = pa.arts.SensorHeterodyneFrequencyRange(ch.reference_hz, [0, np.inf]) + x = pa.arts.sensor.HeterodyneFrequencyRange(ch.reference_hz, [0, np.inf]) # Zero or more mixes, which add more ranges follow. for lo in ch.lo_offsets_hz: x.mix(lo) # The channel response is a simple boxcar with the specified bandwidth. - m = pa.arts.SensorBoxChannel(0, ch.bandwidth_hz / 2, n) + m = pa.arts.sensor.BoxChannel(0, ch.bandwidth_hz / 2, n) v = x.channel_response(m) out.append(v) @@ -84,7 +84,7 @@ def extract_channels(spectral_rad, measurement_sensor): # %% Download data and set up workspace -ws = pa.workspace.Workspace() +ws = pa.Workspace() pa.data.download() # %% Use built-in spectroscopic data. @@ -101,7 +101,7 @@ def extract_channels(spectral_rad, measurement_sensor): ws.ray_path_observer_agendaSetGeometric() # %% Simple sensor setup: position, line-of-sight, polarization, and channels. -sensor = pa.arts.SensorBuilder(channels=sensor_channels(CHANNELS, SAMPLES_PER_CHANNEL)) +sensor = pa.arts.sensor.Builder(channels=sensor_channels(CHANNELS, SAMPLES_PER_CHANNEL)) ws.measurement_sensor, ws.measurement_sensor_meta = sensor(POS, LOS) for i in range(len(ws.measurement_sensor)): @@ -170,8 +170,8 @@ def extract_channels(spectral_rad, measurement_sensor): if "ARTS_HEADLESS" not in os.environ: plt.show() -ref = [294.813805535829, 297.18778420039365, 283.26228661223786, 262.0365042504266, 244.4836135697355, 226.2694109308727, 216.11284609241537, - 209.72414936956548, 210.3336073749979, 219.42825163097098, 230.06432799428234, 241.47556633001957, 252.6256887712745, 261.4587517221863, 292.20874260781585] +ref = [294.82020426974526, 297.1947623380882, 283.302312755353, 262.0990916297846, 244.54610199393989, 226.31688048929024, 216.13528849096693, + 209.71768967552475, 210.24931675027358, 219.3006268679223, 229.9146206010939, 241.32757369987274, 252.48941437009185, 261.41619227415606, 292.22500426558753] assert np.allclose(ws.measurement_vec, ref), \ f"Mismatch to reference simulations.\nReference: {ref}\nSimulations: {ws.measurement_vec:B,}" diff --git a/python/src/pyarts3/plots/SensorObsel.py b/python/src/pyarts3/plots/SensorObsel.py index 6103517d99..3172325755 100644 --- a/python/src/pyarts3/plots/SensorObsel.py +++ b/python/src/pyarts3/plots/SensorObsel.py @@ -332,10 +332,10 @@ def plot(data: pyarts.arts.SensorObsel, import pyarts3 as pyarts import numpy as np - ant = pyarts.arts.SensorGaussianAiryAntenna(np.linspace(0, 1, 11), - np.linspace(0, 330, 12), - 0.3) - ch = pyarts.arts.SensorDiracChannel(100e9) + ant = pyarts.arts.sensor.GaussianAiryAntenna(np.linspace(0, 1, 11), + np.linspace(0, 330, 12), + 0.3) + ch = pyarts.arts.sensor.DiracChannel(100e9) obsel = ant(ch, [1, 0, 0], [90, 0]) pyarts.plots.SensorObsel.plot(obsel, polar=True, point_spread=True) diff --git a/python/src/pyarts3/utils/time_report.py b/python/src/pyarts3/utils/time_report.py index 0534ceb94c..e82c1745a6 100644 --- a/python/src/pyarts3/utils/time_report.py +++ b/python/src/pyarts3/utils/time_report.py @@ -169,12 +169,22 @@ def time_report_table(res, dt, unit): A string containing the time report in Markdown table format. """ + if len(res) == 0: + return "No methods to report." + keys = np.array(list(res.keys())) vs = np.array([dt[key][1] for key in keys]) keys = keys[np.argsort(vs)][::-1] + maxl = max(6, max([len(key) for key in keys])) + N = 20 + + def pad(s, n = N, r=False): + s = str(s) + padding = ' ' * (n - len(s)) + return padding + s if r else s + padding - out = f"| Method | Total Time [{unit}] | Min Time [{unit}] | Max Time [{unit}] | Average Time [{unit}] | Times Called |\n" - out += "| --------- | ------------ | ------------ | -------- | ------------ | ------------ |\n" + out = f"| {pad('Method', maxl)} | {pad(f"Total Time [{unit}]")} | {pad(f"Min Time [{unit}]")} | {pad(f"Max Time [{unit}]", N)} | {pad(f"Average Time [{unit}]")} | {pad("Times Called")} |\n" + out += f"| {'-' * maxl} | {'-' * N} | {'-' * N} | {'-' * N} | {'-' * N} | {'-' * N} |\n" for key in keys: numc = len(res[key]) @@ -183,6 +193,6 @@ def time_report_table(res, dt, unit): mint = min(dts) maxt = max(dts) tott = dt[key][1] - out += f"| {key} | {round(tott,1)} | {round(mint,1)} | {round(maxt,1)} | {round(avgt,1)} | {numc} |\n" + out += f"| {key:<{maxl}} | {pad(round(tott,1), r=1)} | {pad(round(mint,1), r=1)} | {pad(round(maxt,1), r=1)} | {pad(round(avgt,1), r=1)} | {pad(numc, r=1)} |\n" return out diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 19e7260acd..de92c9b1d1 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -192,6 +192,7 @@ add_library(artsworkspace STATIC m_partfun.cc m_planets.cc m_ppvar.cc + m_profile.cc m_propagation_path.cc m_propagation_path_observer.cc m_propmat.cc diff --git a/src/core/absorption/cia.cc b/src/core/absorption/cia.cc index db3eff6aae..c535ce2da2 100644 --- a/src/core/absorption/cia.cc +++ b/src/core/absorption/cia.cc @@ -146,7 +146,7 @@ void cia_interpolation(VectorView result, default: T_order = 3; break; } - using id = lagrange_interp::identity; + using id = lagrange_interp::grid_identity; // Find frequency grid positions: const auto f_lag = lagrange_interp::make_lags( diff --git a/src/core/absorption/xsec_fit.cc b/src/core/absorption/xsec_fit.cc index 84bb3159c9..19433a84c5 100644 --- a/src/core/absorption/xsec_fit.cc +++ b/src/core/absorption/xsec_fit.cc @@ -160,7 +160,7 @@ void XsecRecord::Extract(VectorView result, { const auto f_gp = - lagrange_interp::make_lags<1, lagrange_interp::identity>( + lagrange_interp::make_lags<1, lagrange_interp::grid_identity>( f_grid_active, data_f_grid_active, 0.5, "Frequency"); const auto f_itw = reinterpweights(f_gp); diff --git a/src/core/atm/atm_field.cpp b/src/core/atm/atm_field.cpp index 6543c89ea5..031d4f5ddb 100644 --- a/src/core/atm/atm_field.cpp +++ b/src/core/atm/atm_field.cpp @@ -21,6 +21,8 @@ #include #include +#include "lagrange_interp.h" + bool AtmPoint::zero_wind() const noexcept { return stdr::all_of(wind, Cmp::eq<0>()); } @@ -592,60 +594,58 @@ std::vector Field::keys() const { return out; } -void Point::check_and_fix() try { - ARTS_USER_ERROR_IF(nonstd::isnan(pressure), "Pressure is NaN") - ARTS_USER_ERROR_IF(nonstd::isnan(temperature), "Temperature is NaN") +void Point::check() { + const auto isnan = [](Numeric v) { return std::isnan(v); }; - if (stdr::all_of(wind, [](auto v) { return nonstd::isnan(v); })) { - wind = {0., 0., 0.}; - } else { - ARTS_USER_ERROR_IF( - stdr::any_of(wind, [](auto v) { return nonstd::isnan(v); }), - "Cannot have partially missing wind field. Consider setting the missing field to zero or add it completely.\n" - "Wind field [wind_u, wind_v, wind_w] is: {:B,}", - wind) + std::vector errors; + + if (not(pressure >= 0.0)) { + errors.emplace_back(std::format(" - Bad pressure: {}", pressure)); } - if (stdr::all_of(mag, [](auto v) { return nonstd::isnan(v); })) { - mag = {0., 0., 0.}; - } else { - ARTS_USER_ERROR_IF( - stdr::any_of(mag, [](auto v) { return nonstd::isnan(v); }), - "Cannot have partially missing magnetic field. Consider setting the missing field to zero or add it completely.\n" - "Magnetic field [mag_u, mag_v, mag_w] is: {:B,}", - mag) + if (not(temperature >= 0.0)) { + errors.emplace_back(std::format(" - Bad temperature: {}", temperature)); + } + + if (stdr::any_of(wind, isnan)) { + errors.emplace_back( + std::format(" - Bad wind field [wind_u, wind_v, wind_w]: {:B,}", wind)); + } + + if (stdr::any_of(mag, isnan)) { + errors.emplace_back( + std::format(" - Bad magnetic field [mag_u, mag_v, mag_w]: {:B,}", mag)); + } + + if (stdr::any_of(specs | stdv::values, isnan) or + stdr::any_of(specs | stdv::values, Cmp::lt(0.0))) { + errors.emplace_back(std::format(" - Bad species VMR: {:NB,}", specs)); } - for (auto &spec : specs) { - ARTS_USER_ERROR_IF(nonstd::isnan(spec.second) or spec.second < 0.0, - "VMR for \"{}\" is {}", - toString<1>(spec.first), - spec.second) + if (stdr::any_of(isots | stdv::values, isnan) or + stdr::any_of(isots | stdv::values, Cmp::lt(0.0))) { + errors.emplace_back( + std::format(" - Bad isotopologue ratio: {:NB,}", isots)); } - for (auto &isot : isots) { - //! Cannot check isnan because it is a valid state for isotopologue ratios - ARTS_USER_ERROR_IF(isot.second < 0.0, - "Isotopologue ratio for \"{}\" is {}", - isot.first.FullName(), - isot.second) + if (stdr::any_of(nlte | stdv::values, isnan) or + stdr::any_of(nlte | stdv::values, Cmp::lt(0.0))) { + errors.emplace_back(std::format(" - Bad non-LTE ratio: {:NB,}", nlte)); } - for (auto &nl : nlte) { - ARTS_USER_ERROR_IF(nonstd::isnan(nl.second) or nl.second < 0.0, - "Non-LTE ratio for \"{}\" is {}", - nl.first, - nl.second) + if (stdr::any_of(ssprops | stdv::values, isnan)) { + errors.emplace_back( + std::format(" - Bad scattering species property: {:NB,}", ssprops)); } - for (auto &pp : ssprops) { - ARTS_USER_ERROR_IF(nonstd::isnan(pp.second), - "Scattering Species Property value for \"", - pp.first, - pp.second) + if (not errors.empty()) { + throw std::runtime_error( + std::format("Invalid AtmPoint:\n {}", + errors | stdv::transform([](const std::string &s) { + return s + '\n'; + }) | stdr::to>())); } } -ARTS_METHOD_ERROR_CATCH Data &Field::operator[](const AtmKey &key) { return other[key]; } @@ -881,16 +881,30 @@ struct PositionalNumeric { const Numeric alt; const Numeric lat; const Numeric lon; + const bool log; operator Numeric() const { - return std::visit([&](auto &d) { return get(d, alt, lat, lon); }, data); + if (log) { + return std::visit( + [&](auto &d) { + return get(d, alt, lat, lon); + }, + data); + } + + return std::visit( + [&](auto &d) { + return get(d, alt, lat, lon); + }, + data); } }; std::optional get_optional_limit(const Data &data, const Numeric alt, const Numeric lat, - const Numeric lon) { + const Numeric lon, + const bool log_interpolation) { const auto lim = Atm::find_limit( data, std::visit([](auto &d) { return Atm::find_limits(d); }, data.data), @@ -908,8 +922,11 @@ std::optional get_optional_limit(const Data &data, if (lim.type == InterpolationExtrapolation::Zero) return 0.0; if (lim.type == InterpolationExtrapolation::Nearest) - return PositionalNumeric{ - .data = data.data, .alt = lim.alt, .lat = lim.lat, .lon = lim.lon}; + return PositionalNumeric{.data = data.data, + .alt = lim.alt, + .lat = lim.lat, + .lon = lim.lon, + .log = log_interpolation}; return std::nullopt; } @@ -919,8 +936,12 @@ std::optional get_optional_limit(const Data &data, Numeric Data::at(const Numeric alt, const Numeric lat, const Numeric lon) const { - return interp::get_optional_limit(*this, alt, lat, lon) - .value_or(interp::PositionalNumeric{data, alt, lat, lon}); + return interp::get_optional_limit(*this, alt, lat, lon, log_interpolation) + .value_or(interp::PositionalNumeric{.data = data, + .alt = alt, + .lat = lat, + .lon = lon, + .log = log_interpolation}); } Numeric Data::at(const Vector3 pos) const { return at(pos[0], pos[1], pos[2]); } @@ -935,8 +956,22 @@ Point Field::at(const Numeric alt, const Numeric lat, const Numeric lon) const alt) Point out; - for (auto &&key : keys()) out[key] = operator[](key).at(alt, lat, lon); - out.check_and_fix(); + + static_assert( + sizeof(AtmField) == + sizeof(std::unordered_map) * 5 + sizeof(Numeric), + "The loops below must be over all keys, a size change of AtmField indicates that the number of keys have changed"); + out.specs.reserve(specs.size()); + out.isots.reserve(isots.size()); + out.nlte.reserve(nlte.size()); + out.ssprops.reserve(ssprops.size()); + + for (auto &[key, data] : other) out[key] = data.at(alt, lat, lon); + for (auto &[key, data] : specs) out[key] = data.at(alt, lat, lon); + for (auto &[key, data] : isots) out[key] = data.at(alt, lat, lon); + for (auto &[key, data] : nlte) out[key] = data.at(alt, lat, lon); + for (auto &[key, data] : ssprops) out[key] = data.at(alt, lat, lon); + return out; } ARTS_METHOD_ERROR_CATCH @@ -944,6 +979,7 @@ ARTS_METHOD_ERROR_CATCH Point Field::at(const Vector3 pos) const try { return at(pos[0], pos[1], pos[2]); } + ARTS_METHOD_ERROR_CATCH } // namespace Atm diff --git a/src/core/atm/atm_field.h b/src/core/atm/atm_field.h index 9ca15952b0..d4041d2482 100644 --- a/src/core/atm/atm_field.h +++ b/src/core/atm/atm_field.h @@ -14,6 +14,8 @@ #include +#include "lagrange_interp.h" + AtmKey to_wind(const String &); AtmKey to_mag(const String &); @@ -47,10 +49,10 @@ using KeyVal = std::variant; //! returns whether the key is a wind key -bool is_wind(const KeyVal&); +bool is_wind(const KeyVal &); //! returns whether the key is a temperature key -bool is_temperature(const KeyVal&); +bool is_temperature(const KeyVal &); template concept ListKeyType = requires(T a) { @@ -159,7 +161,7 @@ struct Point { [[nodiscard]] std::pair levels( const QuantumIdentifier &band) const; - void check_and_fix(); + void check(); }; //! All the field data; if these types grow too much we might want to @@ -185,12 +187,14 @@ concept RawDataType = //! Hold all atmospheric data struct Data { FieldData data{FunctionalData{}}; + InterpolationExtrapolation alt_upp{InterpolationExtrapolation::None}; InterpolationExtrapolation alt_low{InterpolationExtrapolation::None}; InterpolationExtrapolation lat_upp{InterpolationExtrapolation::None}; InterpolationExtrapolation lat_low{InterpolationExtrapolation::None}; InterpolationExtrapolation lon_upp{InterpolationExtrapolation::None}; InterpolationExtrapolation lon_low{InterpolationExtrapolation::None}; + bool log_interpolation{false}; // Standard Data(); @@ -213,9 +217,15 @@ struct Data { [[nodiscard]] String data_type() const; + template + [[nodiscard]] decltype(auto) get_if(this auto &&self) { + return std::get_if>(&self.data); + } + template [[nodiscard]] const T &get(this auto &&self) { - auto *out = std::get_if>(&self.data); + auto *out = self.template get_if(); + if (out == nullptr) { throw std::invalid_argument( std::format("Expected type {} not in data. Instead type {} found.", @@ -458,23 +468,28 @@ struct std::formatter { if (tags.short_str) return tags.format(ctx, v.data); const std::string_view sep = tags.sep(); - tags.add_if_bracket(ctx, '['); - tags.format(ctx, v.data, sep); - tags.add_if_bracket(ctx, '['); - tags.format(ctx, - v.alt_upp, - sep, - v.alt_low, - sep, - v.lat_upp, - sep, - v.lat_low, - sep, - v.lon_upp, - sep, - v.lon_low); - tags.add_if_bracket(ctx, ']'); - tags.add_if_bracket(ctx, ']'); + if (tags.help) tags.add_if_bracket(ctx, '['); + tags.format(ctx, v.data); + + if (tags.help) { + tags.format(ctx, sep); + tags.add_if_bracket(ctx, '['); + tags.format(ctx, + v.alt_upp, + sep, + v.alt_low, + sep, + v.lat_upp, + sep, + v.lat_low, + sep, + v.lon_upp, + sep, + v.lon_low); + tags.add_if_bracket(ctx, ']'); + tags.add_if_bracket(ctx, ']'); + } + return ctx.out(); } }; diff --git a/src/core/functional/functional_atm_field.cpp b/src/core/functional/functional_atm_field.cpp index 58604520d2..2e45f9a377 100644 --- a/src/core/functional/functional_atm_field.cpp +++ b/src/core/functional/functional_atm_field.cpp @@ -5,14 +5,18 @@ #include #include "functional_atm_field_interp.h" +#include "lagrange_interp.h" namespace Atm { Numeric MagnitudeField::operator()(Numeric alt, Numeric lat, Numeric lon) const { - const Numeric M = Atm::interp::get(magnitude, alt, lat, lon); - const Numeric T = Atm::interp::get(theta, alt, lat, lon); - const Numeric P = Atm::interp::get(phi, alt, lat, lon); + const Numeric M = Atm::interp::get( + magnitude, alt, lat, lon); + const Numeric T = + Atm::interp::get(theta, alt, lat, lon); + const Numeric P = + Atm::interp::get(phi, alt, lat, lon); const Vector3 x = geocentric2ecef({M, T, P}); diff --git a/src/core/functional/functional_atm_field_interp.cpp b/src/core/functional/functional_atm_field_interp.cpp index 67b7549431..af73a01a3c 100644 --- a/src/core/functional/functional_atm_field_interp.cpp +++ b/src/core/functional/functional_atm_field_interp.cpp @@ -50,60 +50,6 @@ Matrix interpweights(const latlags& b, const lonlags& c) { c); } -Numeric get(const GeodeticField3& f, - const Numeric alt, - const Numeric lat, - const Numeric lon) { - if (not f.ok()) throw std::runtime_error("bad field"); - return std::visit( - [&data = f.data](auto&& al, auto&& la, auto&& lo) { - return lagrange_interp::interp(data, al, la, lo); - }, - altlag(f.grid<0>(), alt), - latlag(f.grid<1>(), lat), - lonlag(f.grid<2>(), lon)); -} - -Numeric get(const GeodeticField3& f, - const Tensor3& iw, - const altlags& alt, - const latlags& lat, - const lonlags& lon) { - if (not f.ok()) throw std::runtime_error("bad field"); - return std::visit( - [&data = f.data, &iw](auto&& al, auto&& la, auto&& lo) { - return lagrange_interp::interp(data, iw, al, la, lo); - }, - alt, - lat, - lon); -} - -Numeric get(const GeodeticField3& f, - Index ialt, - const Matrix& iw, - const latlags& lat, - const lonlags& lon) { - if (not f.ok()) throw std::runtime_error("bad field"); - return std::visit( - [data = f.data[ialt], &iw](auto&& la, auto&& lo) { - return lagrange_interp::interp(data, iw, la, lo); - }, - lat, - lon); -} - -Numeric get(const Numeric num, const Numeric, const Numeric, const Numeric) { - return num; -} - -Numeric get(const std::function& fd, - const Numeric alt, - const Numeric lat, - const Numeric lon) { - return fd(alt, lat, lon); -} - std::vector> flat_weight(const GeodeticField3& gf3, const Numeric alt, const Numeric lat, @@ -114,7 +60,7 @@ std::vector> flat_weight(const GeodeticField3& gf3, const Index nlat = gf3.grid<1>().size(); const Index nlon = gf3.grid<2>().size(); - using id = lagrange_interp::identity; + using id = lagrange_interp::grid_identity; using lc = lagrange_interp::loncross; return std::visit( diff --git a/src/core/functional/functional_atm_field_interp.h b/src/core/functional/functional_atm_field_interp.h index 3bb476bcb7..29091e7ed3 100644 --- a/src/core/functional/functional_atm_field_interp.h +++ b/src/core/functional/functional_atm_field_interp.h @@ -27,33 +27,69 @@ lonlags lonlag(const LonGrid& xs, Numeric x); Tensor3 interpweights(const altlags&, const latlags&, const lonlags&); Matrix interpweights(const latlags&, const lonlags&); +template Numeric get(const GeodeticField3& f, - const Numeric, - const Numeric, - const Numeric); + const Numeric alt, + const Numeric lat, + const Numeric lon) { + if (not f.ok()) throw std::runtime_error("bad field"); + return std::visit( + [&data = f.data](auto&& al, auto&& la, auto&& lo) { + return lagrange_interp::interp(data, al, la, lo); + }, + altlag(f.grid<0>(), alt), + latlag(f.grid<1>(), lat), + lonlag(f.grid<2>(), lon)); +} + +template Numeric get(const GeodeticField3& f, - const Tensor3&, - const altlags&, - const latlags&, - const lonlags&); + const Tensor3& iw, + const altlags& alt, + const latlags& lat, + const lonlags& lon) { + if (not f.ok()) throw std::runtime_error("bad field"); + return std::visit( + [&data = f.data, &iw](auto&& al, auto&& la, auto&& lo) { + return lagrange_interp::interp(data, iw, al, la, lo); + }, + alt, + lat, + lon); +} + +template Numeric get(const GeodeticField3& f, Index ialt, - const Matrix&, - const latlags&, - const lonlags&); + const Matrix& iw, + const latlags& lat, + const lonlags& lon) { + if (not f.ok()) throw std::runtime_error("bad field"); + return std::visit( + [data = f.data[ialt], &iw](auto&& la, auto&& lo) { + return lagrange_interp::interp(data, iw, la, lo); + }, + lat, + lon); +} // Numeric interpolation -Numeric get(const Numeric, const Numeric, const Numeric, const Numeric); +template +Numeric get(const Numeric num, const Numeric, const Numeric, const Numeric) { + return num; +} // Functional interpolation -Numeric get(const std::function&, - const Numeric, - const Numeric, - const Numeric); - -std::vector> flat_weight( - const GeodeticField3& gf3, - const Numeric alt, - const Numeric lat, - const Numeric lon); +template +Numeric get(const std::function& fd, + const Numeric alt, + const Numeric lat, + const Numeric lon) { + return fd(alt, lat, lon); +} + +std::vector> flat_weight(const GeodeticField3& gf3, + const Numeric alt, + const Numeric lat, + const Numeric lon); } // namespace Atm::interp diff --git a/src/core/functional/functional_atm_pressure_field.cpp b/src/core/functional/functional_atm_pressure_field.cpp index 1ba6e2189d..8b3d7f1814 100644 --- a/src/core/functional/functional_atm_pressure_field.cpp +++ b/src/core/functional/functional_atm_pressure_field.cpp @@ -53,8 +53,10 @@ std::pair HydrostaticPressure::level(Index alt_ind, const auto latlag = interp::latlag(pre.grid<1>(), la); const auto lonlag = interp::lonlag(pre.grid<2>(), lo); const auto iw = interp::interpweights(latlag, lonlag); - const Numeric p = interp::get(pre, alt_ind, iw, latlag, lonlag); - const Numeric d = interp::get(grad_p, alt_ind, iw, latlag, lonlag); + const Numeric p = interp::get( + pre, alt_ind, iw, latlag, lonlag); + const Numeric d = interp::get( + grad_p, alt_ind, iw, latlag, lonlag); return {p, d}; } diff --git a/src/core/functional/functional_gravity.cpp b/src/core/functional/functional_gravity.cpp index 03c6bff49b..4c12dfeebd 100644 --- a/src/core/functional/functional_gravity.cpp +++ b/src/core/functional/functional_gravity.cpp @@ -2,60 +2,43 @@ #include #include +#include #include -Numeric EllipsoidGravity::operator()(Numeric h, Numeric lat, Numeric lon) const { - using Conversion::cosd; - using Conversion::sind; - using Math::pow2; - using std::sqrt; +Numeric EllipsoidGravity::operator()(Numeric h, + Numeric lat, + Numeric lon) const { + const auto [r, lat2, lon2] = geodetic2geocentric({h, lat, lon}, {a, b}); - const Numeric N = a / sqrt(1 - pow2(e * sind(lat))); - const Numeric r2 = pow2((N + h) * cosd(lon) * cosd(lat)) + - pow2((N + h) * sind(lon) * cosd(lat)) + - pow2((N * (1 - pow2(e)) + h) * sind(lat)); - - return GM / r2; + return GM / Math::pow2(r); } EllipsoidGravity EllipsoidGravity::Earth() { - return {.GM = Body::Earth::GM, - .a = Body::Earth::a, - .e = std::sqrt(1 - Math::pow2(Body::Earth::b / Body::Earth::a))}; + return {.GM = Body::Earth::GM, .a = Body::Earth::a, .b = Body::Earth::b}; } EllipsoidGravity EllipsoidGravity::Mars() { - return {.GM = Body::Mars::GM, - .a = Body::Mars::a, - .e = std::sqrt(1 - Math::pow2(Body::Mars::b / Body::Mars::a))}; + return {.GM = Body::Mars::GM, .a = Body::Mars::a, .b = Body::Mars::b}; } EllipsoidGravity EllipsoidGravity::Venus() { - return {.GM = Body::Venus::GM, - .a = Body::Venus::a, - .e = std::sqrt(1 - Math::pow2(Body::Venus::b / Body::Venus::a))}; + return {.GM = Body::Venus::GM, .a = Body::Venus::a, .b = Body::Venus::b}; } EllipsoidGravity EllipsoidGravity::Jupiter() { - return {.GM = Body::Jupiter::GM, - .a = Body::Jupiter::a, - .e = std::sqrt(1 - Math::pow2(Body::Jupiter::b / Body::Jupiter::a))}; + return { + .GM = Body::Jupiter::GM, .a = Body::Jupiter::a, .b = Body::Jupiter::b}; } EllipsoidGravity EllipsoidGravity::Moon() { - return {.GM = Body::Moon::GM, - .a = Body::Moon::a, - .e = std::sqrt(1 - Math::pow2(Body::Moon::b / Body::Moon::a))}; + return {.GM = Body::Moon::GM, .a = Body::Moon::a, .b = Body::Moon::b}; } EllipsoidGravity EllipsoidGravity::Mercury() { - return {.GM = Body::Mercury::GM, - .a = Body::Mercury::a, - .e = std::sqrt(1 - Math::pow2(Body::Mercury::b / Body::Mercury::a))}; + return { + .GM = Body::Mercury::GM, .a = Body::Mercury::a, .b = Body::Mercury::b}; } EllipsoidGravity EllipsoidGravity::Saturn() { - return {.GM = Body::Saturn::GM, - .a = Body::Saturn::a, - .e = std::sqrt(1 - Math::pow2(Body::Saturn::b / Body::Saturn::a))}; + return {.GM = Body::Saturn::GM, .a = Body::Saturn::a, .b = Body::Saturn::b}; } diff --git a/src/core/functional/functional_gravity.h b/src/core/functional/functional_gravity.h index 6e6cc0ab43..cb62cbc569 100644 --- a/src/core/functional/functional_gravity.h +++ b/src/core/functional/functional_gravity.h @@ -5,7 +5,7 @@ struct EllipsoidGravity { Numeric GM; Numeric a; - Numeric e; + Numeric b; Numeric operator()(Numeric h, Numeric lat, Numeric lon) const; diff --git a/src/core/geodesy/geodetic.cpp b/src/core/geodesy/geodetic.cpp index 100bbb8110..f0ac757616 100644 --- a/src/core/geodesy/geodetic.cpp +++ b/src/core/geodesy/geodetic.cpp @@ -748,3 +748,22 @@ std::vector visible_coordinates(Vector2 pos, } return out; } + +Numeric refell2r(const Vector2 ell, const Numeric lat) { + using Conversion::cosd; + using Conversion::sind; + + const auto [a, b] = ell; + + assert(a > 0); + assert(b > 0); + assert(a >= b); + + const Numeric e2 = 1 - Math::pow2(b / a); + const Numeric c = 1 - e2; + + const Numeric ct = cosd(lat); + const Numeric st = sind(lat); + + return b / std::sqrt(c * ct * ct + st * st); +} diff --git a/src/core/geodesy/geodetic.h b/src/core/geodesy/geodetic.h index 7a26503f08..2fab3a3a9e 100644 --- a/src/core/geodesy/geodetic.h +++ b/src/core/geodesy/geodetic.h @@ -385,4 +385,32 @@ std::vector visible_coordinates(Vector2 pos, Vector2 ellipsoid, const GeodeticField2& hfield); +/*! + Reference ellipsoid radius, directly from *refellipsoid*. + + Gives the distance from the Earth's centre and the reference ellipsoid as a + function of geoCENTRIC latitude. + + For 1D, extract r directly as refellipsoid[0] to save time. + + This is the basic function to calculate the reference ellipsoid radius. + However, inside the atmosphere this radius is just used at the positions of + the lat_grid. A linear interpolation is applied between these points. This + is handled by other functions. For 2D and 3D and the grid position is + known, use *refell2d*. The function pos2refell_r handles all this in a + general way (but not always the fastest option). + + \return Ellispoid radius + \param refellipsoid [a, b] + \param latitude A geoecentric latitude. + + \author Patrick Eriksson + \date 2012-02-07 + + Update for new definition of ellipsoid and interface + \author Richard Larsson + \date 2024-05-15 +*/ +Numeric refell2r(const Vector2 ell, const Numeric lat); + #endif // geodetic_h diff --git a/src/core/lbl/lbl_jpl.cpp b/src/core/lbl/lbl_jpl.cpp index a54aced49e..c8cd64988b 100644 --- a/src/core/lbl/lbl_jpl.cpp +++ b/src/core/lbl/lbl_jpl.cpp @@ -20,7 +20,7 @@ struct reader { reader(const std::string& s) : it(s.begin()), end(s.end()) {} template - constexpr T read_next(Size n) { + constexpr T read_next(Size n, std::string_view error_context) { std::string_view orig(it, it + n); std::string_view sv = orig; skip(n); @@ -36,18 +36,22 @@ struct reader { if constexpr (std::same_as or std::same_as) { auto res = fast_float::from_chars(sv.data(), sv.data() + sv.size(), x); ARTS_USER_ERROR_IF(res.ec != std::errc{}, - "Failed to parse value from string \"{}\"", + "Failed to parse {} from string \"{}\"", + error_context, orig) ARTS_USER_ERROR_IF(res.ptr != sv.data() + sv.size(), - "Failed to fully parse string \"{}\"", + "Failed to fully parse {} string \"{}\"", + error_context, orig) } else { auto res = std::from_chars(sv.data(), sv.data() + sv.size(), x); ARTS_USER_ERROR_IF(res.ec != std::errc{}, - "Failed to parse value from string \"{}\"", + "Failed to parse {} from string \"{}\"", + error_context, orig) ARTS_USER_ERROR_IF(res.ptr != sv.data() + sv.size(), - "Failed to fully parse string \"{}\"", + "Failed to fully parse {} string \"{}\"", + error_context, orig) } @@ -69,15 +73,15 @@ struct reader { bool read_jpl_entry(jpl_record& record, reader& data) { using namespace Conversion; - record.f0 = data.read_next(13) * 1e6; // MHz to Hz - record.df = data.read_next(8) * 1e6; // MHz to Hz (?) - record.s = std::pow(10., data.read_next(8)) / - 1e12; // log_10(nm2 MHz at T0) to ARTS units - record.dr = data.read_next(2); - record.E = kaycm2joule(data.read_next(10)); // cm^-1 to J - record.g_upp = data.read_next(3); - record.qid = Jpl::data_lookup(std::abs(data.read_next(7))); - record.qnfmt = data.read_next(4); + record.f0 = data.read_next(13, "Central Frequency"sv) * 1e6; // MHz to Hz + record.df = data.read_next(8, "Frequency Deviation"sv) * 1e6; // MHz to Hz (?) + record.s = std::pow(10., data.read_next(8, "Line Intensity"sv)) / + 1e12; // log_10(nm2 MHz at T0) to ARTS units + record.dr = data.read_next(2, "dr"sv); + record.E = kaycm2joule(data.read_next(10, "Energy"sv)); // cm^-1 to J + record.g_upp = data.read_next(3, "Upper State Degeneracy"sv); + record.jpl_id = Jpl::data_lookup(std::abs(data.read_next(7, "JPL ID"sv))); + record.qnfmt = data.read_next(4, "Quantum Number Format"sv); return true; } @@ -118,16 +122,25 @@ line jpl_record::from() const { line out; - out.f0 = f0; - out.e0 = E; - out.gu = g_upp > 0 ? static_cast(g_upp) : -1; - out.gl = -1.0; - out.a = einstein_a(s, out.gu, out.e0, out.f0, qid.T0, qid.QT0); - out.z.on = false; - out.z.mdata = {}; - out.qn = {}; - out.ls = line_shape::model{}; - out.ls.T0 = qid.T0; + out.f0 = f0; + out.e0 = E; + out.gu = g_upp > 0 ? static_cast(g_upp) : -1; + out.gl = -1.0; + // Rescale so at T0, the line intensity matches the JPL value. + // This is needed because ARTS defines the line intensities based + // on built-in partition functions, which may differ from the JPL ones. + out.a = einstein_a(s, + out.gu, + out.e0, + out.f0, + jpl_id.T0, + Math::pow2(jpl_id.QT0) / + PartitionFunctions::Q(jpl_id.T0, jpl_id.qid.isot)); + out.z.on = false; + out.z.mdata = {}; + out.qn = {}; + out.ls = line_shape::model{}; + out.ls.T0 = jpl_id.T0; out.ls.single_models["AIR"_spec].data[G0] = {T1, Vector{25e3, 0.75}}; return out; diff --git a/src/core/lbl/lbl_jpl.h b/src/core/lbl/lbl_jpl.h index 8283e33cf8..17f3b5f260 100644 --- a/src/core/lbl/lbl_jpl.h +++ b/src/core/lbl/lbl_jpl.h @@ -10,7 +10,7 @@ namespace lbl { struct jpl_record { static constexpr Numeric T0 = 300.0; - Jpl::LineDataMod qid; // ID of the species + Jpl::LineDataMod jpl_id; // ID of the species // JPL format in order [F13.4,2F8.4,I2,F10.4,I3,I7,I4,12I2] Numeric f0; // Central frequency @@ -47,7 +47,7 @@ struct std::formatter { const auto sep = tags.sep(); tags.add_if_bracket(ctx, '['); tags.format(ctx, - v.qid, + v.jpl_id, sep, v.f0, sep, diff --git a/src/core/lbl/lbl_lineshape_model.cpp b/src/core/lbl/lbl_lineshape_model.cpp index 4e7e70b46b..479caaa99f 100644 --- a/src/core/lbl/lbl_lineshape_model.cpp +++ b/src/core/lbl/lbl_lineshape_model.cpp @@ -86,7 +86,7 @@ VARIABLE(X3); \ if (not nonstd::isnan(bth)) return res + (1.0 - vmr) * bth; \ \ - return res / vmr; \ + return vmr != 0.0 ? res / vmr : 0.0; \ } \ \ [[nodiscard]] Numeric model::d##mod##_dVMR(const AtmPoint& atm, \ @@ -144,7 +144,7 @@ VARIABLE(DV); \ if (not nonstd::isnan(bth)) return res + (1.0 - vmr) * bth; \ \ - return res / vmr; \ + return vmr != 0.0 ? res / vmr : 0.0; \ } #define VARIABLE(deriv) \ diff --git a/src/core/lbl/lbl_nlte.cpp b/src/core/lbl/lbl_nlte.cpp index 2a3201c7a0..693b934028 100644 --- a/src/core/lbl/lbl_nlte.cpp +++ b/src/core/lbl/lbl_nlte.cpp @@ -217,7 +217,7 @@ QuantumIdentifierVectorMap createCij( const auto numden = atm_point.number_density(key.isot); x.push_back(interp(coll.data, - coll.grid<0>().lag<1, lagrange_interp::identity>( + coll.grid<0>().lag<1, lagrange_interp::grid_identity>( atm_point.temperature)) * numden); } diff --git a/src/core/lookup/lookup_map.cpp b/src/core/lookup/lookup_map.cpp index 5ad2b5b07a..2ce9872594 100644 --- a/src/core/lookup/lookup_map.cpp +++ b/src/core/lookup/lookup_map.cpp @@ -135,7 +135,7 @@ std::array, 1> table::pressure_lagrange( ARTS_USER_ERROR_IF(not do_p(), "No pressure grid set."); const std::array plog_local = {std::log(pressure)}; const Vector& plog_v(*log_p_grid); - return lagrange_interp::make_lags( + return lagrange_interp::make_lags( plog_v, plog_local, interpolation_order, @@ -150,7 +150,7 @@ std::vector> table::frequency_lagrange( const Numeric& extpolation_factor) const try { ARTS_USER_ERROR_IF(not do_f(), "No frequency grid set."); const Vector& f_grid_v(*f_grid); - return lagrange_interp::make_lags( + return lagrange_interp::make_lags( f_grid_v, freq_grid, interpolation_order, @@ -168,7 +168,7 @@ std::array, 1> table::water_lagrange( const std::array x{water_vmr / interp(water_atmref, pressure_lagrange[0])}; const Vector& xi(*w_pert); - return lagrange_interp::make_lags( + return lagrange_interp::make_lags( xi, x, interpolation_order, extpolation_factor, "Water VMR"); } ARTS_METHOD_ERROR_CATCH @@ -182,7 +182,7 @@ std::array, 1> table::temperature_lagrange( const std::array x{temperature - interp(t_atmref, pressure_lagrange[0])}; const Vector& xi(*t_pert); - return lagrange_interp::make_lags( + return lagrange_interp::make_lags( xi, x, interpolation_order, extpolation_factor, "Temperature"); } ARTS_METHOD_ERROR_CATCH diff --git a/src/core/matpack/lagrange_interp.cc b/src/core/matpack/lagrange_interp.cc index e631a9ceac..0cdf99024c 100644 --- a/src/core/matpack/lagrange_interp.cc +++ b/src/core/matpack/lagrange_interp.cc @@ -1,7 +1,14 @@ #include "lagrange_interp.h" namespace lagrange_interp { +static_assert(grid_transformer); +static_assert(value_transformer); +static_assert(value_transformer); static_assert(cyclic); -static_assert(!cyclic); +static_assert(!cyclic); static_assert(constant_lagrange_type_list, 5>>); + +Numeric log_transform::forward(Numeric x) noexcept { return std::log(x); } + +Numeric log_transform::inverse(Numeric x) noexcept { return std::exp(x); } } // namespace lagrange_interp diff --git a/src/core/matpack/lagrange_interp.h b/src/core/matpack/lagrange_interp.h index 225124a2aa..ad3a203099 100644 --- a/src/core/matpack/lagrange_interp.h +++ b/src/core/matpack/lagrange_interp.h @@ -19,9 +19,9 @@ namespace lagrange_interp { namespace { -/*! Defines a transformer concept +/*! Defines a grid transformer concept * - * A transformer is a type that can be used to transform a numeric value + * A grid transformer is a type that can be used to transform a numeric value * to another numeric value. It is used to transform values in the Lagrange * interpolation mechanism. * @@ -29,33 +29,69 @@ namespace { * * wj(x) = (x - x0) * (x - x1) * ... * (x - xm) / ((xj - x0) * (xj - x1) * ... * (xj - xm)) for all j != m. * - * A transformer redefines this to: + * A grid transformer redefines this to: * * wj(x) = (f(x) - f(x0)) * (f(x) - f(x1)) * ... * (f(x) - f(xm)) / ((f(xj) - f(x0)) * (f(xj) - f(x1)) * ... * (f(xj) - f(xm))) for all j != m, * * where f is the transformation function defined as: * - * [static] [constexpr] Numeric operator()(Numeric x) [noexcept]. + * static [constexpr] Numeric op(Numeric x) [noexcept]. * - * The transformer is also responsible for cyclicity of a grid. The transformation is considered - * cyclic if the transformer defines the member methods: + * The grid transformer is also responsible for cyclicity of a grid. The transformation is considered + * cyclic if the grid transformer defines the member methods: * * static [constexpr] Numeric cycle(Numeric), and * static [consteval] Numeric cycle(), and * static [consteval] Numeric midpoint(). * - * The first method cycles the value x to the range defined by the transformer. - * The second returns the full cycle of the transformer. + * The first method cycles the value x to the range defined by the grid transformer. + * The second returns the full cycle of the grid transformer. * The third returns the midpoint of the cycle. * - * Note that cyclicity is only dealt with by the cycle method and not by the operator(). + * Note that cyclicity is only dealt with by the cycle method and not by the op(). */ template -concept transformer = std::same_as; +concept grid_transformer = std::same_as; + +/*! Defines a value transformer concept + * + * A value transformer takes a type and transforms its value before interpolation and back after interpolation. + * + * This is useful for interpolating values such as exponential curves. + * + * To achieve this, the value transformer defines two static methods: + * static [constexpr] Numeric forward(Numeric x) [noexcept], and + * static [constexpr] Numeric inverse(Numeric x) [noexcept]. + * + * These are used such that the interpolation result is given by inverse(sum_j w_j * forward(y_j)) + * + * Note that the value transformer must also define a value_type. Use Any if the forward/inverse methods + * are expected to be generic or the type if they are not. + */ +template +concept value_transformer = requires( + T a, + std::conditional_t, + Numeric, + typename T::value_type> v) { + { + T::forward(v) + } + -> std::same_as, + Numeric, + typename T::value_type>>; + { + T::inverse(v) + } + -> std::same_as, + Numeric, + typename T::value_type>>; +}; template concept cyclic = - transformer and std::same_as and + grid_transformer and + std::same_as and std::same_as and T::cycle() > 0.0 and std::same_as; @@ -100,21 +136,36 @@ concept runtime_lagrange_type_list = lagrange_type_list and not constant_lagrange_type_list; } // namespace -//! A transformer that does not transform the value -struct identity { - static constexpr Numeric operator()(Numeric x) noexcept { return x; } +//! A grid_transformer that does not transform the value of the grid +struct grid_identity { + static constexpr Numeric op(Numeric x) noexcept { return x; } +}; + +//! A value_transformer that does not transform the value of the grid +struct value_identity { + using value_type = Any; + + template + static constexpr auto forward(T&& x) noexcept { + return std::forward(x); + } + + template + static constexpr auto inverse(T&& x) noexcept { + return std::forward(x); + } }; -/*! A transformer that cycles the value to the range [lower, upper) +/*! A grid_transformer that cycles the value to the range [lower, upper) * * Note that only values [2 * lower - upper, 2 * upper - lower) are * cycled. This is to ensure that the transformation is quick. If - * you need a greater range, please implement another transformer + * you need a greater range, please implement another grid_transformer * and keep this fast. */ template requires(lower < upper) -struct cycler { +struct cycler final : grid_identity { static consteval Numeric midpoint() noexcept { return std::midpoint(upper, lower); } @@ -124,8 +175,14 @@ struct cycler { static constexpr Numeric cycle(Numeric x) noexcept { return x + ((x < lower) - (x >= upper)) * cycle(); } +}; - static constexpr Numeric operator()(Numeric x) noexcept { return x; } +//! A value_transformer that interpolates the logarithm of the value, useful for exponential curves +struct log_transform { + using value_type = Numeric; + + static Numeric forward(Numeric x) noexcept; + static Numeric inverse(Numeric x) noexcept; }; //! [-180, 180) cycler @@ -134,7 +191,7 @@ using loncross = cycler<-180.0, 180.0>; /****************************************************************** * Knowledge about the order of the interpolation is very important * for the interpolation to work correctly. The order is defined - * by the transformer and the size of the input array. + * by the grid_transformer and the size of the input array. ******************************************************************/ struct ascending_grid_t {}; @@ -172,7 +229,7 @@ consteval bool ascending() { * the previous x. ******************************************************************/ -template +template constexpr void update_pos(std::span indx, std::span xi, Numeric x) { @@ -183,14 +240,11 @@ constexpr void update_pos(std::span indx, if constexpr (cyclic) x = transform::cycle(x); - const Size N = P - 1; - const Size Of = N / 2; - std::span::iterator xp = - xi.begin() + indx[Of]; - const std::span::iterator xf = - xi.begin() + Of; - const std::span::iterator xe = - xi.end() - P / 2 - 1; + const Size N = P - 1; + const Size Of = N / 2; + std::span::iterator xp = xi.begin() + indx[Of]; + const std::span::iterator xf = xi.begin() + Of; + const std::span::iterator xe = xi.end() - P / 2 - 1; if constexpr (cyclic) { if constexpr (ascending()) { @@ -254,7 +308,7 @@ constexpr void update_pos(std::span indx, * with sorted evenly spaced arrays. ******************************************************************/ -template +template constexpr void find_pos(std::span indx, std::span xi, Numeric x) { @@ -293,7 +347,7 @@ constexpr void find_pos(std::span indx, * The Lagrange interpolation weights ******************************************************************/ -template +template constexpr void set_weights(std::span data, std::span indx, std::span xi, @@ -301,7 +355,7 @@ constexpr void set_weights(std::span data, const Size N = (M == std::dynamic_extent) ? (indx.size() - 1) : (M - 1); if constexpr (cyclic) x = transform::cycle(x); - x = transform{}(x); + x = transform::op(x); // Last value must make the sum equal to 1 // And i != j in the internal loop @@ -310,7 +364,7 @@ constexpr void set_weights(std::span data, if constexpr (ascending()) { if (x < xi.front()) { for (Size j = 0; j < N; ++j) { - Numeric xj = transform{}(xi[indx[j]]); + Numeric xj = transform::op(xi[indx[j]]); if (xj > transform::midpoint()) xj -= transform::cycle(); Numeric numer = 1.0; @@ -318,7 +372,7 @@ constexpr void set_weights(std::span data, for (Size k = 0; k < N; ++k) { Size m = indx[k + (k >= j)]; // i != j - Numeric xm = transform{}(xi[m]); + Numeric xm = transform::op(xi[m]); if (xm > transform::midpoint()) xm -= transform::cycle(); numer *= x - xm; @@ -329,7 +383,7 @@ constexpr void set_weights(std::span data, } } else if (x > xi.back()) { for (Size j = 0; j < N; ++j) { - Numeric xj = transform{}(xi[indx[j]]); + Numeric xj = transform::op(xi[indx[j]]); if (xj < transform::midpoint()) xj += transform::cycle(); Numeric numer = 1.0; @@ -337,7 +391,7 @@ constexpr void set_weights(std::span data, for (Size k = 0; k < N; ++k) { Size m = indx[k + (k >= j)]; // i != j - Numeric xm = transform{}(xi[m]); + Numeric xm = transform::op(xi[m]); if (xm < transform::midpoint()) xm += transform::cycle(); numer *= x - xm; denom *= xj - xm; @@ -346,14 +400,14 @@ constexpr void set_weights(std::span data, } } else { for (Size j = 0; j < N; ++j) { - Numeric xj = transform{}(xi[indx[j]]); + Numeric xj = transform::op(xi[indx[j]]); Numeric numer = 1.0; Numeric denom = 1.0; for (Size k = 0; k < N; ++k) { Size m = indx[k + (k >= j)]; // i != j - Numeric xm = transform{}(xi[m]); + Numeric xm = transform::op(xi[m]); numer *= x - xm; denom *= xj - xm; } @@ -364,7 +418,7 @@ constexpr void set_weights(std::span data, } else { if (x < xi.back()) { for (Size j = 0; j < N; ++j) { - Numeric xj = transform{}(xi[indx[j]]); + Numeric xj = transform::op(xi[indx[j]]); if (xj > transform::midpoint()) xj -= transform::cycle(); Numeric numer = 1.0; @@ -372,7 +426,7 @@ constexpr void set_weights(std::span data, for (Size k = 0; k < N; ++k) { Size m = indx[k + (k >= j)]; // i != j - Numeric xm = transform{}(xi[m]); + Numeric xm = transform::op(xi[m]); if (xm > transform::midpoint()) xm -= transform::cycle(); numer *= x - xm; @@ -383,7 +437,7 @@ constexpr void set_weights(std::span data, } } else if (x > xi.front()) { for (Size j = 0; j < N; ++j) { - Numeric xj = transform{}(xi[indx[j]]); + Numeric xj = transform::op(xi[indx[j]]); if (xj < transform::midpoint()) xj += transform::cycle(); Numeric numer = 1.0; @@ -391,7 +445,7 @@ constexpr void set_weights(std::span data, for (Size k = 0; k < N; ++k) { Size m = indx[k + (k >= j)]; // i != j - Numeric xm = transform{}(xi[m]); + Numeric xm = transform::op(xi[m]); if (xm < transform::midpoint()) xm += transform::cycle(); numer *= x - xm; denom *= xj - xm; @@ -400,14 +454,14 @@ constexpr void set_weights(std::span data, } } else { for (Size j = 0; j < N; ++j) { - Numeric xj = transform{}(xi[indx[j]]); + Numeric xj = transform::op(xi[indx[j]]); Numeric numer = 1.0; Numeric denom = 1.0; for (Size k = 0; k < N; ++k) { Size m = indx[k + (k >= j)]; // i != j - Numeric xm = transform{}(xi[m]); + Numeric xm = transform::op(xi[m]); numer *= x - xm; denom *= xj - xm; } @@ -418,14 +472,14 @@ constexpr void set_weights(std::span data, } } else { for (Size j = 0; j < N; ++j) { - Numeric xj = transform{}(xi[indx[j]]); + Numeric xj = transform::op(xi[indx[j]]); Numeric numer = 1.0; Numeric denom = 1.0; for (Size k = 0; k < N; ++k) { Size m = indx[k + (k >= j)]; // i != j - Numeric xm = transform{}(xi[m]); + Numeric xm = transform::op(xi[m]); numer *= x - xm; denom *= xj - xm; } @@ -447,7 +501,7 @@ constexpr void set_weights(std::span data, * to N + 1, where N is the polynomial order. ******************************************************************/ -template +template struct lag_t { static constexpr bool cyclic = lagrange_interp::cyclic; static constexpr bool runtime = N == -1; @@ -474,9 +528,7 @@ struct lag_t { constexpr lag_t& operator=(lag_t&&) noexcept = default; template - constexpr lag_t(std::span xi, - Numeric x, - grid) + constexpr lag_t(std::span xi, Numeric x, grid) requires(not runtime) { find_pos(indx, xi, x); @@ -484,10 +536,7 @@ struct lag_t { } template - constexpr lag_t(std::span xi, - Numeric x, - Index M, - grid) + constexpr lag_t(std::span xi, Numeric x, Index M, grid) requires(runtime) : data(M + 1), indx(M + 1) { find_pos(indx, xi, x); @@ -495,10 +544,7 @@ struct lag_t { } template - constexpr lag_t(indx_t pos, - std::span xi, - Numeric x, - grid) + constexpr lag_t(indx_t pos, std::span xi, Numeric x, grid) : indx(std::move(pos)) { if constexpr (runtime) data.resize(indx.size()); set_weights(data, indx, xi, x); @@ -522,17 +568,15 @@ struct lag_t { ******************************************************************/ namespace { -template -lag_t poly_lag( - std::span xi, Numeric x) { +template +lag_t poly_lag(std::span xi, Numeric x) { assert(xi.size() > poly); if (xi[0] < xi[1]) return lag_t(xi, x, ascending_grid_t{}); return lag_t(xi, x, descending_grid_t{}); } -template -T variant_lag_helper(std::span xi, - Numeric x) { +template +T variant_lag_helper(std::span xi, Numeric x) { if constexpr (N == 0) { return lag_t<0, transform>(xi, x, ascending_grid_t{}); } else { @@ -541,7 +585,7 @@ T variant_lag_helper(std::span xi, } } -template +template auto variant_lag_helper(std::span xi, Numeric x, std::index_sequence) { @@ -557,7 +601,7 @@ auto variant_lag_helper(std::span xi, * neighbor or linear interpolation. You can specify other values for N to * get higher order polynomial interpolation. */ -template +template auto variant_lag(std::span xi, Numeric x) { assert(xi.size() > 0); return variant_lag_helper( @@ -568,13 +612,12 @@ auto variant_lag(std::span xi, Numeric x) { * Check limits ******************************************************************/ -template -constexpr order_t check_limit( - const std::span& xi, - const std::span& xn, - Numeric extrapolation_limit, - const Index polyorder, - const char* info) try { +template +constexpr order_t check_limit(const std::span& xi, + const std::span& xn, + Numeric extrapolation_limit, + const Index polyorder, + const char* info) try { const Index n = xi.size(); const bool ascending = n <= 1 or xi[0] < xi[1]; auto retval = @@ -647,7 +690,7 @@ The actual minimum value is {8} std::format("Error in check_limit for {}:\n{}", info, e.what())); } -template Orig, matpack::ranked_md<1> New> constexpr order_t check_limit(const Orig& xi, @@ -664,7 +707,7 @@ constexpr order_t check_limit(const Orig& xi, } } -template Orig, Size N> +template Orig, Size N> constexpr order_t check_limit(const Orig& xi, const std::array& xn, Numeric extrapolation_limit, @@ -680,7 +723,7 @@ constexpr order_t check_limit(const Orig& xi, //! Fixed version of make_lags template > constexpr auto make_lags(std::span xi, @@ -736,7 +779,7 @@ constexpr auto make_lags(std::span xi, } template Orig, matpack::ranked_md<1> New, class FlagT = lag_t> @@ -754,7 +797,7 @@ constexpr auto make_lags(const Orig& xi, } template Orig, Size M, class FlagT = lag_t> @@ -766,7 +809,7 @@ constexpr auto make_lags(const Orig& xi, } //! Dynamic version of make_lags -template > constexpr auto make_lags(std::span xi, @@ -824,7 +867,7 @@ constexpr auto make_lags(std::span xi, } } -template Orig, matpack::ranked_md<1> New, class FlagT = lag_t<-1, transform>> @@ -842,7 +885,7 @@ constexpr auto make_lags(const Orig& xi, } } -template Orig, Size M, class FlagT = lag_t<-1, transform>> @@ -902,7 +945,9 @@ consteval Size inner_size() { ******************************************************************/ //! Reuse interpolation weights. -template +template constexpr auto interp(const matpack::ranked_md auto& field, const matpack::ranked_md auto& itw, const FlagTs&... lags) { @@ -914,14 +959,19 @@ constexpr auto interp(const matpack::ranked_md auto& field, for (std::array s{}; s.front() < n.front(); inc(s, n)) { out += std::apply( - [&](auto&&... i) { return field[lags.indx[i]...] * itw[i...]; }, s); + [&](auto&&... i) { + return transform::forward(field[lags.indx[i]...]) * itw[i...]; + }, + s); } - return out; + return transform::inverse(out); } //! Direct interpolation. -template +template constexpr auto interp(const matpack::ranked_md auto& field, const FlagTs&... lags) { using T = std::remove_cvref_t; @@ -933,12 +983,13 @@ constexpr auto interp(const matpack::ranked_md auto& field, for (std::array s{}; s.front() < n.front(); inc(s, n)) { out += std::apply( [&](auto&&... i) { - return (field[lags.indx[i]...] * ... * lags.data[i]); + return (transform::forward(field[lags.indx[i]...]) * ... * + lags.data[i]); }, s); } - return out; + return transform::inverse(out); } /****************************************************************** @@ -958,7 +1009,9 @@ constexpr auto interp(const matpack::ranked_md auto& field, ******************************************************************/ //! Reuse output and interpolation weights. -template +template constexpr void reinterp(matpack::mut_ranked_md auto&& out, const matpack::ranked_md auto& field, const matpack::ranked_md<2 * N> auto& itw, @@ -969,13 +1022,17 @@ constexpr void reinterp(matpack::mut_ranked_md auto&& out, for (std::array s{}; s.front() < n.front(); inc(s, n)) { std::apply( - [&](auto&&... i) { out[i...] = interp(field, itw[i...], lags[i]...); }, + [&](auto&&... i) { + out[i...] = interp(field, itw[i...], lags[i]...); + }, s); } } //! Reuse output. -template +template constexpr void reinterp(matpack::mut_ranked_md auto&& out, const matpack::ranked_md auto& field, const FlagTs&... lags) { @@ -983,12 +1040,16 @@ constexpr void reinterp(matpack::mut_ranked_md auto&& out, const std::array n{out.shape()}; for (std::array s{}; s.front() < n.front(); inc(s, n)) { - std::apply([&](auto&&... i) { out[i...] = interp(field, lags[i]...); }, s); + std::apply( + [&](auto&&... i) { out[i...] = interp(field, lags[i]...); }, + s); } } //! Reuse interpolation weights. -template +template constexpr auto reinterp(const matpack::ranked_md auto& field, const matpack::ranked_md<2 * N> auto& itw, const FlagTs&... lags) { @@ -996,28 +1057,30 @@ constexpr auto reinterp(const matpack::ranked_md auto& field, if constexpr ((constant_lagrange_type_list and ...)) { matpack::cdata_t out{}; - reinterp(out, field, itw, lags...); + reinterp(out, field, itw, lags...); return out; } else { matpack::data_t out(lags.size()...); - reinterp(out, field, itw, lags...); + reinterp(out, field, itw, lags...); return out; } } //! Direct reinterpolation. -template +template constexpr auto reinterp(const matpack::ranked_md auto& field, const FlagTs&... lags) { using T = std::remove_cvref_t; if constexpr ((constant_lagrange_type_list and ...)) { matpack::cdata_t out{}; - reinterp(out, field, lags...); + reinterp(out, field, lags...); return out; } else { matpack::data_t out(lags.size()...); - reinterp(out, field, lags...); + reinterp(out, field, lags...); return out; } } @@ -1037,26 +1100,33 @@ constexpr auto reinterp(const matpack::ranked_md auto& field, ******************************************************************/ //! Reuse output and interpolation weights. -template +template constexpr void flat_interp(matpack::mut_ranked_md<1> auto&& out, const matpack::ranked_md auto& field, const matpack::ranked_md auto& itw, const FlagTs&... lags) { const Size n = out.size(); - for (Size i = 0; i < n; ++i) out[i] = interp(field, itw[i], lags[i]...); + for (Size i = 0; i < n; ++i) + out[i] = interp(field, itw[i], lags[i]...); } //! Reuse output. -template +template constexpr void flat_interp(matpack::mut_ranked_md<1> auto&& out, const matpack::ranked_md auto& field, const FlagTs&... lags) { const Size n = out.size(); - for (Size i = 0; i < n; ++i) out[i] = interp(field, lags[i]...); + for (Size i = 0; i < n; ++i) out[i] = interp(field, lags[i]...); } //! Reuse interpolation weights. -template +template constexpr auto flat_interp(const matpack::ranked_md auto& field, const matpack::ranked_md auto& itw, const FlagTs&... lags) { @@ -1064,28 +1134,30 @@ constexpr auto flat_interp(const matpack::ranked_md auto& field, if constexpr ((constant_lagrange_type_list and ...)) { matpack::cdata_t out{}; - flat_interp(out, field, itw, lags...); + flat_interp(out, field, itw, lags...); return out; } else { matpack::data_t out(size(lags...)); - flat_interp(out, field, itw, lags...); + flat_interp(out, field, itw, lags...); return out; } } //! Direct interpolation. -template +template constexpr auto flat_interp(const matpack::ranked_md auto& field, const FlagTs&... lags) { using T = std::remove_cvref_t; if constexpr ((constant_lagrange_type_list and ...)) { matpack::cdata_t out{}; - flat_interp(out, field, lags...); + flat_interp(out, field, lags...); return out; } else { matpack::data_t out(size(lags...)); - flat_interp(out, field, lags...); + flat_interp(out, field, lags...); return out; } } @@ -1183,7 +1255,7 @@ constexpr auto reinterpweights(const FlagTs&... lags) { } } // namespace lagrange_interp -template +template struct std::formatter> { format_tags tags; @@ -1202,12 +1274,12 @@ struct std::formatter> { } }; -template +template struct xml_io_stream_name> { static constexpr std::string_view name = "lagrange_interp"; }; -template +template struct xml_io_stream> { static constexpr std::string_view type_name = xml_io_stream_name_v>; diff --git a/src/core/matpack/matpack_mdspan_helpers_grid_t.h b/src/core/matpack/matpack_mdspan_helpers_grid_t.h index e8818407a2..c9fe531785 100644 --- a/src/core/matpack/matpack_mdspan_helpers_grid_t.h +++ b/src/core/matpack/matpack_mdspan_helpers_grid_t.h @@ -133,17 +133,17 @@ class grid_t { lagrange_interp::ascending_grid_t, lagrange_interp::descending_grid_t>; - template + template [[nodiscard]] auto lag(Numeric xi) const { return lagrange_interp::lag_t{x, xi, lagsorter{}}; } - template + template [[nodiscard]] auto lag(Numeric xi, Size N) const { return lagrange_interp::lag_t<-1, transform>{x, xi, N, lagsorter{}}; } - template + template [[nodiscard]] auto lag(std::span xi, Numeric extrapolation_limit = 0.5, const char* info = "UNNAMED") const { @@ -151,7 +151,7 @@ class grid_t { x, xi, extrapolation_limit, info); } - template + template [[nodiscard]] auto lag(std::span xi, Size N, Numeric extrapolation_limit = 0.5, diff --git a/src/core/matpack/matpack_mdspan_helpers_ranged_grid_t.h b/src/core/matpack/matpack_mdspan_helpers_ranged_grid_t.h index e26510fc47..a16a757a0c 100644 --- a/src/core/matpack/matpack_mdspan_helpers_ranged_grid_t.h +++ b/src/core/matpack/matpack_mdspan_helpers_ranged_grid_t.h @@ -85,17 +85,17 @@ template return *this; } - template + template [[nodiscard]] auto lag(Numeric x) const { return lagrange_interp::lag_t{*this, x, lag_sorting_t{}}; } - template + template [[nodiscard]] auto lag(Numeric x, Size N) const { return lagrange_interp::lag_t<-1, transform>{*this, x, N, lag_sorting_t{}}; } - template + template [[nodiscard]] auto lag(std::span x, Numeric extrapolation_limit = 0.5, const char* info = "UNNAMED") const { @@ -103,7 +103,7 @@ template *this, x, extrapolation_limit, info); } - template + template [[nodiscard]] auto lag(std::span x, Size N, Numeric extrapolation_limit = 0.5, diff --git a/src/core/path/atm_path.cpp b/src/core/path/atm_path.cpp index 7a71f86c7e..1aff417417 100644 --- a/src/core/path/atm_path.cpp +++ b/src/core/path/atm_path.cpp @@ -19,12 +19,37 @@ ArrayOfAtmPoint &atm_path_resize(ArrayOfAtmPoint &atm_path, void forward_atm_path(ArrayOfAtmPoint &atm_path, const ArrayOfPropagationPathPoint &rad_path, const AtmField &atm) { - stdr::transform(rad_path, - atm_path.begin(), - [&atm](const PropagationPathPoint &pp) -> AtmPoint { - if (pp.has(PathPositionType::atm)) return atm.at(pp.pos); - return atm.at(atm.top_of_atmosphere, pp.pos[1], pp.pos[2]); - }); + std::string error{}; + + atm_path.resize(atm_path.size()); + +#pragma omp parallel for + for (Size i = 0; i < rad_path.size(); i++) { + auto &pp = rad_path[i]; + auto &ap = atm_path[i]; + try { + if (pp.has(PathPositionType::atm)) { + ap = atm.at(pp.pos); + } else { + ap = atm.at(atm.top_of_atmosphere, pp.pos[1], pp.pos[2]); + } + } catch (const std::exception &e) { +#pragma omp critical + if (error.empty()) { + error = std::format( + R"(Error extracting atmospheric point from path at index {}: +pos: {:B,} +los: {:B,} +Error message: {})", + i, + pp.pos, + pp.los, + e.what()); + } + } + } + + ARTS_USER_ERROR_IF(not error.empty(), error); } ArrayOfAtmPoint forward_atm_path(const ArrayOfPropagationPathPoint &rad_path, diff --git a/src/core/path/path_point.cpp b/src/core/path/path_point.cpp index 772c3c5746..e06012afa2 100644 --- a/src/core/path/path_point.cpp +++ b/src/core/path/path_point.cpp @@ -19,17 +19,6 @@ #include namespace path { -Numeric PropagationPathPoint::altitude() const noexcept { return pos[0]; } -Numeric& PropagationPathPoint::altitude() noexcept { return pos[0]; } -Numeric PropagationPathPoint::latitude() const noexcept { return pos[1]; } -Numeric& PropagationPathPoint::latitude() noexcept { return pos[1]; } -Numeric PropagationPathPoint::longitude() const noexcept { return pos[2]; } -Numeric& PropagationPathPoint::longitude() noexcept { return pos[2]; } -Numeric PropagationPathPoint::zenith() const noexcept { return los[0]; } -Numeric& PropagationPathPoint::zenith() noexcept { return los[0]; } -Numeric PropagationPathPoint::azimuth() const noexcept { return los[1]; } -Numeric& PropagationPathPoint::azimuth() noexcept { return los[1]; } - Vector2 mirror(const Vector2 los) { Vector2 los_mirrored; los_mirrored[0] = 180 - los[0]; diff --git a/src/core/path/path_point.h b/src/core/path/path_point.h index 5954e6e08d..f47345c26b 100644 --- a/src/core/path/path_point.h +++ b/src/core/path/path_point.h @@ -38,16 +38,30 @@ struct PropagationPathPoint { return pos_type == x or los_type == x; } - [[nodiscard]] Numeric altitude() const noexcept; - Numeric& altitude() noexcept; - [[nodiscard]] Numeric latitude() const noexcept; - Numeric& latitude() noexcept; - [[nodiscard]] Numeric longitude() const noexcept; - Numeric& longitude() noexcept; - [[nodiscard]] Numeric zenith() const noexcept; - Numeric& zenith() noexcept; - [[nodiscard]] Numeric azimuth() const noexcept; - Numeric& azimuth() noexcept; + template + [[nodiscard]] constexpr decltype(auto) altitude(this Self&& p) noexcept { + return std::forward(p).pos[0]; + } + + template + [[nodiscard]] constexpr decltype(auto) latitude(this Self&& p) noexcept { + return std::forward(p).pos[1]; + } + + template + [[nodiscard]] constexpr decltype(auto) longitude(this Self&& p) noexcept { + return std::forward(p).pos[2]; + } + + template + [[nodiscard]] constexpr decltype(auto) zenith(this Self&& p) noexcept { + return std::forward(p).los[0]; + } + + template + [[nodiscard]] constexpr decltype(auto) azimuth(this Self&& p) noexcept { + return std::forward(p).los[1]; + } }; using ArrayOfPropagationPathPoint = std::vector; diff --git a/src/core/sensor/antenna_pattern.cpp b/src/core/sensor/antenna_pattern.cpp index d25ac164cd..c634cb1bfa 100644 --- a/src/core/sensor/antenna_pattern.cpp +++ b/src/core/sensor/antenna_pattern.cpp @@ -23,8 +23,7 @@ struct AntennaBasis { struct AntennaGeometrySample { Size sample_index{}; - Numeric ant_zen{}; - Numeric ant_azi{}; + Numeric offset_zenith{}; bool has_response{}; bool is_representative{}; }; @@ -80,8 +79,6 @@ template std::vector samples; samples.reserve(zen_grid.size() * azi_grid.size()); - using Conversion::atan2d; - for (Numeric izen : zen_grid) { Size degenerate_sample_index = std::numeric_limits::max(); @@ -89,7 +86,7 @@ template const Vector3 local = antenna_frame_los({izen, iazi}); const bool degenerate = antenna_azimuth_is_degenerate(local); - Size sample_index = degenerate_sample_index; + Size sample_index = degenerate_sample_index; bool is_representative = false; if (sample_index == std::numeric_limits::max()) { const Vector3 enu = normalized(local[0] * basis.v + local[1] * basis.h + @@ -101,14 +98,13 @@ template if (degenerate) degenerate_sample_index = sample_index; } - auto& sample = samples.emplace_back(AntennaGeometrySample{ - .sample_index = sample_index, - .has_response = local[2] > 0.0, - .is_representative = is_representative}); + auto& sample = samples.emplace_back( + AntennaGeometrySample{.sample_index = sample_index, + .has_response = local[2] > 0.0, + .is_representative = is_representative}); if (sample.has_response) { - sample.ant_zen = atan2d(local[0], local[2]); - sample.ant_azi = atan2d(local[1], local[2]); + sample.offset_zenith = std::abs(izen); } } } @@ -123,13 +119,9 @@ template [[nodiscard]] AntennaPatternField make_gaussian_field(ZenGrid zen_grid, AziGrid azi_grid, - Numeric zenith_std, - Numeric azimuth_std, + Numeric std, Stokvec weight) { - ARTS_USER_ERROR_IF(zenith_std <= 0.0, - "Gaussian antenna zenith_std must be positive") - ARTS_USER_ERROR_IF(azimuth_std <= 0.0, - "Gaussian antenna azimuth_std must be positive") + ARTS_USER_ERROR_IF(std <= 0.0, "Gaussian antenna std must be positive") AntennaPatternField out{ .data_name = "gaussian"s, @@ -138,8 +130,6 @@ template .grids = {std::move(zen_grid), std::move(azi_grid)}, }; - using Conversion::atan2d; - for (Size izen = 0; izen < out.grid<0>().size(); ++izen) { for (Size iazi = 0; iazi < out.grid<1>().size(); ++iazi) { const Vector3 local = @@ -150,11 +140,8 @@ template continue; } - const Numeric ant_zen = atan2d(local[0], local[2]); - const Numeric ant_azi = atan2d(local[1], local[2]); - const Numeric exponent = - -0.5 * ((ant_zen / zenith_std) * (ant_zen / zenith_std) + - (ant_azi / azimuth_std) * (ant_azi / azimuth_std)); + const Numeric offset_zenith = out.grid<0>()[izen]; + const Numeric exponent = -0.5 * Math::pow2(offset_zenith / std); out[izen, iazi] = std::exp(exponent) * weight; } @@ -172,10 +159,9 @@ template return Conversion::hwhm2std(hwhm_deg); } -[[nodiscard]] Numeric gaussian_airy_response(Numeric ant_zen, - Numeric ant_azi, - Numeric inv_std_sq) { - return std::exp(-0.5 * (ant_zen * ant_zen + ant_azi * ant_azi) * inv_std_sq); +[[nodiscard]] Numeric gaussian_airy_response(Numeric offset_zenith, + Numeric airy_std) { + return std::exp(-0.5 * Math::pow2(offset_zenith / airy_std)); } std::shared_ptr make_single_poslos_grid( @@ -215,7 +201,7 @@ Obsel GriddedAntennaPattern::operator()(const Channel& channel, const Vector2& bore_los) const { ARTS_USER_ERROR_IF( not data.ok(), - "SensorGriddedAntennaPattern data shape does not match its grids") + "GriddedAntennaPattern data shape does not match its grids") const auto& zen_grid = data.grid<0>(); const auto& azi_grid = data.grid<1>(); @@ -232,9 +218,15 @@ Obsel GriddedAntennaPattern::operator()(const Channel& channel, for (Size izen = 0; izen < zen_grid.size(); ++izen) { for (Size iazi = 0; iazi < azi_grid.size(); ++iazi) { + const auto& sample = geometry.samples[igrid]; + if (not sample.is_representative) { + ++igrid; + continue; + } + const auto& antenna_weight = data[izen, iazi]; if (not antenna_weight.is_zero()) { - const Size sample_index = geometry.samples[igrid].sample_index; + const Size sample_index = sample.sample_index; for (Size ifreq = 0; ifreq < channel_weights.size(); ++ifreq) { if (channel_weights[ifreq] == 0.0) continue; weight_matrix[sample_index, ifreq] += @@ -253,35 +245,33 @@ std::shared_ptr GriddedAntennaPattern::clone() const { return std::make_shared(*this); } -GaussianAntenna::GaussianAntenna(ZenGrid zen_grid, - AziGrid azi_grid, - Numeric zenith_std, - Numeric azimuth_std, - Stokvec weight) { - data = make_gaussian_field(std::move(zen_grid), - std::move(azi_grid), - zenith_std, - azimuth_std, - weight); +namespace { +AziGrid make_azi_grid(Size size) { + ARTS_USER_ERROR_IF(size < 1, + "Gaussian antenna azimuth size must be at least 1") + + return Vector{nlinspace(0.0, 360.0, size + 1)[Range(0, size)]}; } +} // namespace GaussianAntenna::GaussianAntenna(ZenGrid zen_grid, - AziGrid azi_grid, Numeric std, - Stokvec weight) - : GaussianAntenna( - std::move(zen_grid), std::move(azi_grid), std, std, weight) {} + Size azi_grid_size, + Stokvec weight) { + data = make_gaussian_field( + std::move(zen_grid), make_azi_grid(azi_grid_size), std, weight); +} std::shared_ptr GaussianAntenna::clone() const { return std::make_shared(*this); } GaussianAiryAntenna::GaussianAiryAntenna(ZenGrid zen_grid, - AziGrid azi_grid, Numeric aperture_diameter, + Size azi_grid_size, Stokvec weight) : zen_grid(std::move(zen_grid)), - azi_grid(std::move(azi_grid)), + azi_grid(make_azi_grid(azi_grid_size)), aperture_diameter(aperture_diameter), weight(weight) {} @@ -297,7 +287,7 @@ Obsel GaussianAiryAntenna::operator()(const Channel& channel, ARTS_USER_ERROR_IF( freq_grid->front() <= 0.0, - "Gaussian Airy antenna requires strictly positive channel frequencies because SensorBuilder only provides the channel frequency grid") + "Gaussian Airy antenna requires strictly positive channel frequencies because sensor builder only provides the channel frequency grid") auto geometry = make_antenna_geometry_layout(zen_grid, azi_grid, pos, bore_los); @@ -305,22 +295,21 @@ Obsel GaussianAiryAntenna::operator()(const Channel& channel, channel_weights.size()); if (not weight.is_zero()) { - Vector inv_std_sq(channel_weights.size(), 0.0); + Vector airy_std(channel_weights.size(), 0.0); Vector normalization(channel_weights.size(), 0.0); Vector retained_scale(channel_weights.size(), 0.0); for (Size ifreq = 0; ifreq < channel_weights.size(); ++ifreq) { if (channel_weights[ifreq] == 0.0) continue; - const Numeric airy_std = + airy_std[ifreq] = gaussian_airy_std((*freq_grid)[ifreq], aperture_diameter); - inv_std_sq[ifreq] = 1.0 / (airy_std * airy_std); Numeric retained_mass = 0.0; for (const auto& sample : geometry.samples) { if (not sample.has_response) continue; - normalization[ifreq] += gaussian_airy_response( - sample.ant_zen, sample.ant_azi, inv_std_sq[ifreq]); + normalization[ifreq] += + gaussian_airy_response(sample.offset_zenith, airy_std[ifreq]); } if (normalization[ifreq] == 0.0) continue; @@ -330,10 +319,8 @@ Obsel GaussianAiryAntenna::operator()(const Channel& channel, if (not sample.is_representative) continue; const Numeric single_weight = - channel_weights[ifreq] * gaussian_airy_response( - sample.ant_zen, - sample.ant_azi, - inv_std_sq[ifreq]) / + channel_weights[ifreq] * + gaussian_airy_response(sample.offset_zenith, airy_std[ifreq]) / normalization[ifreq]; retained_mass += single_weight; @@ -354,10 +341,8 @@ Obsel GaussianAiryAntenna::operator()(const Channel& channel, if (retained_scale[ifreq] == 0.0) continue; const Numeric single_weight = - channel_weights[ifreq] * gaussian_airy_response( - sample.ant_zen, - sample.ant_azi, - inv_std_sq[ifreq]) / + channel_weights[ifreq] * + gaussian_airy_response(sample.offset_zenith, airy_std[ifreq]) / normalization[ifreq]; if (single_weight == 0.0) continue; diff --git a/src/core/sensor/antenna_pattern.h b/src/core/sensor/antenna_pattern.h index 3edf62cf3e..43fe3bdd94 100644 --- a/src/core/sensor/antenna_pattern.h +++ b/src/core/sensor/antenna_pattern.h @@ -76,16 +76,10 @@ struct GaussianAntenna final : GriddedAntennaPattern { GaussianAntenna(GaussianAntenna&&) = default; GaussianAntenna& operator=(GaussianAntenna&&) = default; - GaussianAntenna(ZenGrid zen_grid, - AziGrid azi_grid, - Numeric zenith_std, - Numeric azimuth_std, - Stokvec weight = {1.0, 0.0, 0.0, 0.0}); - - GaussianAntenna(ZenGrid zen_grid, - AziGrid azi_grid, + GaussianAntenna(ZenGrid grid, Numeric std, - Stokvec weight = {1.0, 0.0, 0.0, 0.0}); + Size azi_grid_size = 2, + Stokvec weight = {1.0, 0.0, 0.0, 0.0}); [[nodiscard]] std::shared_ptr clone() const override; }; @@ -103,9 +97,9 @@ struct GaussianAiryAntenna final : AntennaPattern { GaussianAiryAntenna& operator=(GaussianAiryAntenna&&) = default; GaussianAiryAntenna(ZenGrid zen_grid, - AziGrid azi_grid, Numeric aperture_diameter, - Stokvec weight = {1.0, 0.0, 0.0, 0.0}); + Size azi_grid_size = 2, + Stokvec weight = {1.0, 0.0, 0.0, 0.0}); [[nodiscard]] Obsel operator()(const Channel& channel, const Vector3& pos, diff --git a/src/core/sensor/frequency_bandpass_filters.cpp b/src/core/sensor/frequency_bandpass_filters.cpp index 74cb591578..507562f7d8 100644 --- a/src/core/sensor/frequency_bandpass_filters.cpp +++ b/src/core/sensor/frequency_bandpass_filters.cpp @@ -96,8 +96,8 @@ void deduplicate_zero_frequency_images( } using InterpolationPoint = - std::variant, - lagrange_interp::lag_t<1, lagrange_interp::identity>>; + std::variant, + lagrange_interp::lag_t<1, lagrange_interp::grid_identity>>; struct FoldedLocalPoint { InterpolationPoint interpolation; @@ -106,7 +106,7 @@ struct FoldedLocalPoint { InterpolationPoint make_interpolation_point(const AscendingGrid& grid, Numeric f) { - return lagrange_interp::variant_lag(grid, f); + return lagrange_interp::variant_lag(grid, f); } Numeric interpolate_channel_weight(const Vector& weights, @@ -189,8 +189,8 @@ void combine_samples(std::vector>& samples) { samples = std::move(combined); } -std::vector build_channels( - const FrequencyRange& range, const std::span& channels) { +std::vector build_channels(const FrequencyRange& range, + const std::span& channels) { std::vector out; out.reserve(channels.size()); @@ -235,16 +235,15 @@ std::vector build_channels( combine_samples(samples); out.push_back( - Channel{.channel = make_filter(samples, - std::format("channel-response-{}", - ichan))}); + Channel{.channel = make_filter( + samples, std::format("channel-response-{}", ichan))}); } return out; } -std::vector build_synced_channels( - const FrequencyRange& range, const Spectrometer& spectrometer) { +std::vector build_synced_channels(const FrequencyRange& range, + const Spectrometer& spectrometer) { const auto& channels = spectrometer.channels; if (channels.empty()) return {}; @@ -286,9 +285,8 @@ std::vector build_synced_channels( combine_samples(samples); out.push_back( - Channel{.channel = make_filter(samples, - std::format("channel-response-{}", - ichan))}); + Channel{.channel = make_filter( + samples, std::format("channel-response-{}", ichan))}); } return out; @@ -300,8 +298,8 @@ std::vector make_bandpass_channels( return build_channels(range, channels); } -std::vector make_bandpass_channels( - const FrequencyRange& range, const Spectrometer& spectrometer) { +std::vector make_bandpass_channels(const FrequencyRange& range, + const Spectrometer& spectrometer) { return build_synced_channels(range, spectrometer); } } // namespace sensor diff --git a/src/core/sensor/frequency_channel_selection.cpp b/src/core/sensor/frequency_channel_selection.cpp index 5e28582a2f..4afd23c5eb 100644 --- a/src/core/sensor/frequency_channel_selection.cpp +++ b/src/core/sensor/frequency_channel_selection.cpp @@ -147,6 +147,13 @@ Spectrometer::Spectrometer(const Channel& base_channel, sync_frequency_grids(); } + +Spectrometer::Spectrometer(const AscendingGrid& freq_offsets) + : channels(std::from_range, + freq_offsets | + stdv::transform([](auto f) -> DiracChannel { return f; })) { + sync_frequency_grids(); +} } // namespace sensor void xml_io_stream::write( diff --git a/src/core/sensor/frequency_channel_selection.h b/src/core/sensor/frequency_channel_selection.h index cc58dc716e..ebf3836507 100644 --- a/src/core/sensor/frequency_channel_selection.h +++ b/src/core/sensor/frequency_channel_selection.h @@ -77,6 +77,8 @@ struct Spectrometer { */ Spectrometer(const Channel& base_channel, const AscendingGrid& freq_offsets); + Spectrometer(const AscendingGrid& freq_offsets); + Spectrometer(std::vector channels) : channels(std::move(channels)) {} operator const std::vector&() const { return channels; } diff --git a/src/core/sensor/sensor_builder.cpp b/src/core/sensor/sensor_builder.cpp index f1404b1a3a..416620435f 100644 --- a/src/core/sensor/sensor_builder.cpp +++ b/src/core/sensor/sensor_builder.cpp @@ -29,28 +29,27 @@ SensorMetaInfo make_meta_info(Size nchannels, Size geometry_index) { } } // namespace -SensorBuilder::SensorBuilder() : antenna(PencilBeamAntenna{}.clone()) {} +Builder::Builder() : antenna(PencilBeamAntenna{}.clone()) {} -SensorBuilder::SensorBuilder(std::vector channels, - std::shared_ptr antenna) +Builder::Builder(std::vector channels, + std::shared_ptr antenna) : channels(std::move(channels)), antenna(std::move(antenna)) {} -SensorBuilder::SensorBuilder(const Spectrometer& spectrometer, - std::shared_ptr antenna) - : SensorBuilder(std::vector{spectrometer.channels}, - std::move(antenna)) { +Builder::Builder(const Spectrometer& spectrometer, + std::shared_ptr antenna) + : Builder(std::vector{spectrometer.channels}, std::move(antenna)) { preserve_common_frequency_grid = spectrometer.is_synced(); } -SensorBuilder::SensorBuilder(const Spectrometer& spectrometer, - const FrequencyRange& backend, - std::shared_ptr antenna) - : SensorBuilder(make_bandpass_channels(backend, spectrometer), - std::move(antenna)) { +Builder::Builder(const Spectrometer& spectrometer, + const FrequencyRange& backend, + std::shared_ptr antenna) + : Builder(make_bandpass_channels(backend, spectrometer), + std::move(antenna)) { preserve_common_frequency_grid = spectrometer.is_synced(); } -SensorBuilder& SensorBuilder::operator=(const SensorBuilder& other) { +Builder& Builder::operator=(const Builder& other) { if (this != &other) { channels = other.channels; antenna = clone_antenna(other.antenna); @@ -60,18 +59,16 @@ SensorBuilder& SensorBuilder::operator=(const SensorBuilder& other) { return *this; } -std::pair SensorBuilder::operator()( +std::pair Builder::operator()( std::span pos, std::span los) const { - ARTS_USER_ERROR_IF(channels.empty(), - "SensorBuilder requires at least one channel") - ARTS_USER_ERROR_IF(not antenna, "SensorBuilder requires an antenna pattern") + ARTS_USER_ERROR_IF(channels.empty(), "Builder requires at least one channel") + ARTS_USER_ERROR_IF(not antenna, "Builder requires an antenna pattern") ARTS_USER_ERROR_IF(pos.empty(), - "SensorBuilder requires at least one sensor position") - ARTS_USER_ERROR_IF(los.empty(), - "SensorBuilder requires at least one bore LOS") + "Builder requires at least one sensor position") + ARTS_USER_ERROR_IF(los.empty(), "Builder requires at least one bore LOS") ARTS_USER_ERROR_IF( pos.size() != los.size(), - "SensorBuilder requires matching position and LOS counts. Got {} positions and {} LOS values.", + "Builder requires matching position and LOS counts. Got {} positions and {} LOS values.", pos.size(), los.size()) @@ -116,5 +113,5 @@ std::pair SensorBuilder::operator()( return {std::move(out), std::move(meta)}; } -static_assert(SensorBuilderSelection); +static_assert(BuilderSelection); } // namespace sensor diff --git a/src/core/sensor/sensor_builder.h b/src/core/sensor/sensor_builder.h index b376934c3b..661808df64 100644 --- a/src/core/sensor/sensor_builder.h +++ b/src/core/sensor/sensor_builder.h @@ -13,30 +13,30 @@ namespace sensor { //! Combines channels with an antenna pattern and bore geometries into obsels. -struct SensorBuilder; +struct Builder; struct FrequencyRange; //! Concept for selecting sensor builders. template -concept SensorBuilderSelection = std::derived_from; +concept BuilderSelection = std::derived_from; -struct SensorBuilder { +struct Builder { std::vector channels; std::shared_ptr antenna; bool preserve_common_frequency_grid{false}; - SensorBuilder(); - SensorBuilder(std::vector channels, - std::shared_ptr antenna); - SensorBuilder(const Spectrometer& spectrometer, - std::shared_ptr antenna); - SensorBuilder(const Spectrometer& spectrometer, - const FrequencyRange& backend, - std::shared_ptr antenna); - SensorBuilder(const SensorBuilder& other) = default; - SensorBuilder(SensorBuilder&&) noexcept = default; - SensorBuilder& operator=(const SensorBuilder& other); - SensorBuilder& operator=(SensorBuilder&&) noexcept = default; + Builder(); + Builder(std::vector channels, + std::shared_ptr antenna); + Builder(const Spectrometer& spectrometer, + std::shared_ptr antenna); + Builder(const Spectrometer& spectrometer, + const FrequencyRange& backend, + std::shared_ptr antenna); + Builder(const Builder& other) = default; + Builder(Builder&&) noexcept = default; + Builder& operator=(const Builder& other); + Builder& operator=(Builder&&) noexcept = default; [[nodiscard]] std::pair operator()( std::span pos, std::span los) const; diff --git a/src/core/spec/species_tags.cc b/src/core/spec/species_tags.cc index 2b197ea1a2..c8686fb3f8 100644 --- a/src/core/spec/species_tags.cc +++ b/src/core/spec/species_tags.cc @@ -168,7 +168,10 @@ SpeciesTagTypeStatus::SpeciesTagTypeStatus( const ArrayOfSpeciesTag& abs_species) { for (auto& tag : abs_species) { switch (tag.type) { - case SpeciesTagType::Plain: Plain = true; break; + case SpeciesTagType::Plain: + Plain = true; + FreeElectrons = FreeElectrons or tag.Spec() == "free_electrons"_spec; + break; case SpeciesTagType::Predefined: Predefined = true; break; case SpeciesTagType::Cia: Cia = true; break; case SpeciesTagType::XsecFit: XsecFit = true; break; diff --git a/src/core/spec/species_tags.h b/src/core/spec/species_tags.h index 5d5531e130..1adb81baa3 100644 --- a/src/core/spec/species_tags.h +++ b/src/core/spec/species_tags.h @@ -54,7 +54,8 @@ using ArrayOfArrayOfSpeciesTag = Array; //! Struct to test of an ArrayOfArrayOfSpeciesTag contains a Speciestagtype struct SpeciesTagTypeStatus { - bool Plain{false}, Predefined{false}, Cia{false}, XsecFit{false}; + bool Plain{false}, Predefined{false}, Cia{false}, XsecFit{false}, + FreeElectrons{false}; SpeciesTagTypeStatus(const ArrayOfSpeciesTag& abs_species); }; diff --git a/src/core/subsurface/subsurf_field.cpp b/src/core/subsurface/subsurf_field.cpp index b37172749a..b502c89370 100644 --- a/src/core/subsurface/subsurf_field.cpp +++ b/src/core/subsurface/subsurf_field.cpp @@ -29,10 +29,10 @@ Numeric &Point::operator[](SubsurfaceKey x) { } Numeric Point::operator[](const SubsurfacePropertyTag &x) const { - return prop.at(x); + return props.at(x); } -Numeric &Point::operator[](const SubsurfacePropertyTag &x) { return prop[x]; } +Numeric &Point::operator[](const SubsurfacePropertyTag &x) { return props[x]; } Numeric Point::operator[](const KeyVal &x) const { return std::visit([this](auto &key) { return this->operator[](key); }, x); @@ -60,7 +60,7 @@ std::vector Point::keys() const { std::vector out; out.reserve(size()); for (auto &a : enumtyps::SubsurfaceKeyTypes) out.emplace_back(a); - for (auto &a : prop | stdv::keys) out.emplace_back(a); + for (auto &a : props | stdv::keys) out.emplace_back(a); return out; } @@ -253,7 +253,7 @@ Numeric get(const GeodeticField3 &gf3, const Numeric lon) { if (not gf3.ok()) throw std::runtime_error("bad field"); - using id = lagrange_interp::identity; + using id = lagrange_interp::grid_identity; using lc = lagrange_interp::loncross; return std::visit( @@ -330,7 +330,7 @@ std::array, 8> flat_weight_(const GeodeticField3 &gf3, const Index nlat = gf3.grid<1>().size(); const Index nlon = gf3.grid<2>().size(); - using id = lagrange_interp::identity; + using id = lagrange_interp::grid_identity; using lc = lagrange_interp::loncross; return std::visit( @@ -464,13 +464,13 @@ Index Field::size() const { return nbasic(); } std::vector out; out.reserve(size()); for (const auto &[key, _] : other) out.emplace_back(key); - for (const auto &[key, _] : prop) out.emplace_back(key); + for (const auto &[key, _] : props) out.emplace_back(key); return out; } Data &Field::operator[](const SubsurfaceKey &key) { return other[key]; } -Data &Field::operator[](const SubsurfacePropertyTag &key) { return prop[key]; } +Data &Field::operator[](const SubsurfacePropertyTag &key) { return props[key]; } Data &Field::operator[](const KeyVal &key) { return std::visit([this](auto &k) -> Data & { return this->operator[](k); }, @@ -482,7 +482,7 @@ const Data &Field::operator[](const SubsurfaceKey &key) const { } const Data &Field::operator[](const SubsurfacePropertyTag &key) const { - return prop.at(key); + return props.at(key); } const Data &Field::operator[](const KeyVal &key) const { @@ -508,7 +508,16 @@ Point Field::at(const Numeric alt, const Numeric lat, const Numeric lon) const alt) Point out; - for (auto &&key : keys()) out[key] = operator[](key).at(alt, lat, lon); + + static_assert( + sizeof(SubsurfaceField) == + sizeof(std::unordered_map) * 2 + sizeof(Numeric), + "The loops below must be over all keys, a size change of SubsurfaceField indicates that the number of keys have changed"); + out.props.reserve(props.size()); + + for (auto &[key, data] : other) out[key] = data.at(alt, lat, lon); + for (auto &[key, data] : props) out[key] = data.at(alt, lat, lon); + return out; } ARTS_METHOD_ERROR_CATCH @@ -526,7 +535,7 @@ std::string std::formatter::to_string( R"("density": )"sv, v.density, "\nproperties:\n"sv, - v.prop); + v.props); return tags.bracket ? ("{" + out + "}") : out; } @@ -545,7 +554,7 @@ std::string std::formatter::to_string( v.other.size(), sep, R"("Prop": )"sv, - v.prop.size()); + v.props.size()); } else { const std::string_view sep = tags.sep(); @@ -555,8 +564,8 @@ std::string std::formatter::to_string( out += tags.vformat(sep, R"("Other": )"sv, v.other); } - if (not v.prop.empty()) { - out += tags.vformat(sep, R"("Prop": )"sv, v.prop); + if (not v.props.empty()) { + out += tags.vformat(sep, R"("Prop": )"sv, v.props); } } diff --git a/src/core/subsurface/subsurf_field.h b/src/core/subsurface/subsurf_field.h index aa1218323c..c2f294021d 100644 --- a/src/core/subsurface/subsurf_field.h +++ b/src/core/subsurface/subsurf_field.h @@ -39,7 +39,7 @@ struct Point { Numeric temperature{0}; Numeric density{0}; - std::unordered_map prop; + std::unordered_map props; Point(); Point(const Point &); @@ -62,7 +62,7 @@ struct Point { if constexpr (isSubsurfaceKey) { return has(std::forward(keys)...); } else if constexpr (isSubsurfacePropertyTag) { - return prop.contains(key) and has(std::forward(keys)...); + return props.contains(key) and has(std::forward(keys)...); } else { static_assert( isSubsurfacePropertyTag and not isSubsurfacePropertyTag, @@ -72,7 +72,7 @@ struct Point { if constexpr (isSubsurfaceKey) { return true; } else if constexpr (isSubsurfacePropertyTag) { - return prop.contains(key); + return props.contains(key); } else { static_assert( isSubsurfacePropertyTag and not isSubsurfacePropertyTag, @@ -173,7 +173,7 @@ struct Data { struct Field final { std::unordered_map other; - std::unordered_map prop; + std::unordered_map props; Numeric bottom_depth{std::numeric_limits::max()}; diff --git a/src/core/sun/sun.cc b/src/core/sun/sun.cc index 7c5b1c44d4..6a7367ce5f 100644 --- a/src/core/sun/sun.cc +++ b/src/core/sun/sun.cc @@ -86,7 +86,7 @@ Matrix regrid_sun_spectrum(const GriddedField2& sun_spectrum_raw, const ConstVectorView data_f_grid_active = data_f_grid[active_range]; // Check if frequency is inside the range covered by the data: - lagrange_interp::check_limit( + lagrange_interp::check_limit( data_f_grid, f_grid_active, 0.5, 1, "Frequency grid for sun spectrum"); { @@ -133,53 +133,6 @@ std::ostream& operator<<(std::ostream& os, const ArrayOfSun& a) { return os; } -namespace { -/*! - Reference ellipsoid radius, directly from *refellipsoid*. - - Gives the distance from the Earth's centre and the reference ellipsoid as a - function of geoCENTRIC latitude. - - For 1D, extract r directly as refellipsoid[0] to save time. - - This is the basic function to calculate the reference ellipsoid radius. - However, inside the atmosphere this radius is just used at the positions of - the lat_grid. A linear interpolation is applied between these points. This - is handled by other functions. For 2D and 3D and the grid position is - known, use *refell2d*. The function pos2refell_r handles all this in a - general way (but not always the fastest option). - - \return Ellispoid radius - \param refellipsoid [a, b] - \param latitude A geoecentric latitude. - - \author Patrick Eriksson - \date 2012-02-07 - - Update for new definition of ellipsoid and interface - \author Richard Larsson - \date 2024-05-15 -*/ -Numeric refell2r(const Vector2 ell, const Numeric lat) { - using Conversion::cosd; - using Conversion::sind; - - const auto [a, b] = ell; - - assert(a > 0); - assert(b > 0); - assert(a >= b); - - const Numeric e2 = 1 - Math::pow2(b / a); - const Numeric c = 1 - e2; - - const Numeric ct = cosd(lat); - const Numeric st = sind(lat); - - return b / std::sqrt(c * ct * ct + st * st); -} -} // namespace - /*! Conversion from spherical to cartesian coordinates. diff --git a/src/core/surface/surf_field.cpp b/src/core/surface/surf_field.cpp index 30577034ac..f9f1ce2906 100644 --- a/src/core/surface/surf_field.cpp +++ b/src/core/surface/surf_field.cpp @@ -357,9 +357,14 @@ std::vector Field::keys() const { Point Field::at(Numeric lat, Numeric lon) const { Point out; - for (auto &key : keys()) { - out[key] = this->operator[](key).at(lat, lon); - } + static_assert( + sizeof(SurfaceField) == sizeof(std::unordered_map) * 2 + + 2 * sizeof(Numeric), + "The loops below must be over all keys, a size change of SurfaceField indicates that the number of keys have changed"); + out.prop.reserve(props.size()); + + for (auto &[key, data] : other) out[key] = data.at(lat, lon); + for (auto &[key, data] : props) out[key] = data.at(lat, lon); out.normal = normal(lat, lon); diff --git a/src/core/tests/test_laginterp.cpp b/src/core/tests/test_laginterp.cpp index a0ff8b15ec..5b9270a90a 100644 --- a/src/core/tests/test_laginterp.cpp +++ b/src/core/tests/test_laginterp.cpp @@ -11,9 +11,9 @@ #include "time_test_util.h" namespace { -template +template consteval std::string_view transform_name() { - if constexpr (std::same_as) { + if constexpr (std::same_as) { return "identity"sv; } else if constexpr (std::same_as) { return "loncross"sv; @@ -38,15 +38,16 @@ bool all_close(const std::vector& a) { } template + lagrange_interp::grid_transformer Transformer1, + lagrange_interp::grid_transformer Transformer2> void test_fixed() { using namespace lagrange_interp; constexpr auto name1 = transform_name(); constexpr auto name2 = transform_name(); - test_timer_t test_fixed_time{std::format("test_fixed N={} {}-{}", N, name1, name2)}; + test_timer_t test_fixed_time{ + std::format("test_fixed N={} {}-{}", N, name1, name2)}; const Vector grid1{-120, -90, -60, -30, 0, 30, 60, 90, 120}; const Vector grid2{120, 90, 60, 30, 0, -30, -60, -90, -120}; @@ -259,15 +260,16 @@ flat_interp: {:B,})", } } -template +template void test_runtime(Size N) { using namespace lagrange_interp; constexpr auto name1 = transform_name(); constexpr auto name2 = transform_name(); - test_timer_t test_fixed_time{std::format("test_runtime N={} {}-{}", N, name1, name2)}; + test_timer_t test_fixed_time{ + std::format("test_runtime N={} {}-{}", N, name1, name2)}; const Vector grid1{-120, -90, -60, -30, 0, 30, 60, 90, 120}; const Vector grid2{120, 90, 60, 30, 0, -30, -60, -90, -120}; @@ -533,7 +535,8 @@ void test_variant_lag() { for (Size i = 1; i < 10; ++i) { Vector grid = nlinspace(-100, 100, i); const auto lagt = - lagrange_interp::variant_lag(grid, 1); + lagrange_interp::variant_lag(grid, + 1); Size x = std::visit([](auto v) { return v.PolyOrder; }, lagt); ARTS_USER_ERROR_IF(x > N, "PolyOrder {} exceeds max lag size {}", x, N); ARTS_USER_ERROR_IF(i > N and x != N, @@ -548,23 +551,33 @@ void test_variant_lag() { int main() { for (Index i = 0; i < 5; ++i) { - test_fixed<0, lagrange_interp::identity, lagrange_interp::identity>(); - test_fixed<1, lagrange_interp::identity, lagrange_interp::identity>(); - test_fixed<2, lagrange_interp::identity, lagrange_interp::identity>(); - test_fixed<3, lagrange_interp::identity, lagrange_interp::identity>(); - test_fixed<4, lagrange_interp::identity, lagrange_interp::identity>(); - - test_fixed<0, lagrange_interp::identity, lagrange_interp::loncross>(); - test_fixed<1, lagrange_interp::identity, lagrange_interp::loncross>(); - test_fixed<2, lagrange_interp::identity, lagrange_interp::loncross>(); - test_fixed<3, lagrange_interp::identity, lagrange_interp::loncross>(); - test_fixed<4, lagrange_interp::identity, lagrange_interp::loncross>(); - - test_fixed<0, lagrange_interp::loncross, lagrange_interp::identity>(); - test_fixed<1, lagrange_interp::loncross, lagrange_interp::identity>(); - test_fixed<2, lagrange_interp::loncross, lagrange_interp::identity>(); - test_fixed<3, lagrange_interp::loncross, lagrange_interp::identity>(); - test_fixed<4, lagrange_interp::loncross, lagrange_interp::identity>(); + test_fixed<0, + lagrange_interp::grid_identity, + lagrange_interp::grid_identity>(); + test_fixed<1, + lagrange_interp::grid_identity, + lagrange_interp::grid_identity>(); + test_fixed<2, + lagrange_interp::grid_identity, + lagrange_interp::grid_identity>(); + test_fixed<3, + lagrange_interp::grid_identity, + lagrange_interp::grid_identity>(); + test_fixed<4, + lagrange_interp::grid_identity, + lagrange_interp::grid_identity>(); + + test_fixed<0, lagrange_interp::grid_identity, lagrange_interp::loncross>(); + test_fixed<1, lagrange_interp::grid_identity, lagrange_interp::loncross>(); + test_fixed<2, lagrange_interp::grid_identity, lagrange_interp::loncross>(); + test_fixed<3, lagrange_interp::grid_identity, lagrange_interp::loncross>(); + test_fixed<4, lagrange_interp::grid_identity, lagrange_interp::loncross>(); + + test_fixed<0, lagrange_interp::loncross, lagrange_interp::grid_identity>(); + test_fixed<1, lagrange_interp::loncross, lagrange_interp::grid_identity>(); + test_fixed<2, lagrange_interp::loncross, lagrange_interp::grid_identity>(); + test_fixed<3, lagrange_interp::loncross, lagrange_interp::grid_identity>(); + test_fixed<4, lagrange_interp::loncross, lagrange_interp::grid_identity>(); test_fixed<0, lagrange_interp::loncross, lagrange_interp::loncross>(); test_fixed<1, lagrange_interp::loncross, lagrange_interp::loncross>(); @@ -573,9 +586,12 @@ int main() { test_fixed<4, lagrange_interp::loncross, lagrange_interp::loncross>(); for (Size N = 0; N < 5; ++N) { - test_runtime(N); - test_runtime(N); - test_runtime(N); + test_runtime(N); + test_runtime( + N); + test_runtime( + N); test_runtime(N); } diff --git a/src/core/util/compare.h b/src/core/util/compare.h index b0f6ece091..5ffa8929a3 100644 --- a/src/core/util/compare.h +++ b/src/core/util/compare.h @@ -1,5 +1,7 @@ #pragma once +#include + namespace Cmp { //! Returns a 'less than' lambda expression for use in, e.g., std::any_of constexpr auto lt(auto v) { @@ -70,4 +72,25 @@ template constexpr auto contains() { return [](const auto& x) static { return x.contains(V); }; } + +//! Returns a `this->in` lambda expression for use in, e.g., std::any_of +template +constexpr auto in_sorted(const R& v) { + return [&v](const auto& x) { return std::ranges::binary_search(v, x); }; +} + +template +constexpr auto in(const R& v) { + return [&v](const auto& x) { return std::ranges::contains(v, x); }; +} + +template +constexpr auto not_in_sorted(const R& v) { + return [test = in_sorted(v)](const auto& x) { return not test(x); }; +} + +template +constexpr auto not_in(const R& v) { + return [test = in(v)](const auto& x) { return not test(x); }; +} } // namespace Cmp diff --git a/src/core/util/planet_data.h b/src/core/util/planet_data.h index 6dd15132ad..a5e5eaf491 100644 --- a/src/core/util/planet_data.h +++ b/src/core/util/planet_data.h @@ -1,47 +1,55 @@ #pragma once #include +#include namespace Body { namespace Earth { -inline constexpr Numeric a = 6378137.0; -inline constexpr Numeric b = 6356752.314245; -inline constexpr Numeric GM = 398600.435507; +inline constexpr Numeric a = 6378137.0; +inline constexpr Numeric b = 6356752.314245; +inline constexpr Numeric GM = 398600.435507 * 1e9; +inline constexpr Numeric day = 0.99726968 * 24 * 60 * 60; } // namespace Earth namespace Jupiter { -inline constexpr Numeric a = 71492e3; -inline constexpr Numeric b = 66854e3; -inline constexpr Numeric GM = 126712764.1; +inline constexpr Numeric a = 71492e3; +inline constexpr Numeric b = 66854e3; +inline constexpr Numeric GM = 126712764.1 * 1e9; +inline constexpr Numeric day = 0.41354 * 24 * 60 * 60; } // namespace Jupiter namespace Mars { -inline constexpr Numeric a = 3396.19e3; -inline constexpr Numeric b = 3376.20e3; -inline constexpr Numeric GM = 42828.375816; +inline constexpr Numeric a = 3396.19e3; +inline constexpr Numeric b = 3376.20e3; +inline constexpr Numeric GM = 42828.375816 * 1e9; +inline constexpr Numeric day = 1.026 * 24 * 60 * 60; } // namespace Mars namespace Moon { inline constexpr Numeric a = 1738.1e3; inline constexpr Numeric b = 1736.0e3; -inline constexpr Numeric GM = 4902.800118; +inline constexpr Numeric GM = 4902.800118 * 1e9; +inline constexpr Numeric day = 29.5 * 24 * 60 * 60; } // namespace Moon namespace Mercury { -inline constexpr Numeric a = 2439.7e3; -inline constexpr Numeric b = 2439.7e3; -inline constexpr Numeric GM = 22031.868551; +inline constexpr Numeric a = 2439.7e3; +inline constexpr Numeric b = 2439.7e3; +inline constexpr Numeric GM = 22031.868551 * 1e9; +inline constexpr Numeric day = 58.646 * 24 * 60 * 60; } // namespace Mercury namespace Venus { -inline constexpr Numeric a = 6051.8e3; -inline constexpr Numeric b = 6051.8e3; -inline constexpr Numeric GM = 324858.592; +inline constexpr Numeric a = 6051.8e3; +inline constexpr Numeric b = 6051.8e3; +inline constexpr Numeric GM = 324858.592 * 1e9; +inline constexpr Numeric day = -243.018 * 24 * 60 * 60; } // namespace Venus namespace Saturn { -inline constexpr Numeric a = 0.5 * 120536e3; -inline constexpr Numeric b = 0.5 * 108728e3; -inline constexpr Numeric GM = 37940584.841800; +inline constexpr Numeric a = 0.5 * 120536e3; +inline constexpr Numeric b = 0.5 * 108728e3; +inline constexpr Numeric GM = 37940584.841800 * 1e9; +inline constexpr Numeric day = 0.444 * 24 * 60 * 60; } // namespace Saturn } // namespace Body diff --git a/src/core/util/unique_unordered_map.h b/src/core/util/unique_unordered_map.h index 46a6b5ea92..25468a3d2c 100644 --- a/src/core/util/unique_unordered_map.h +++ b/src/core/util/unique_unordered_map.h @@ -9,35 +9,58 @@ struct UniqueMap { std::unordered_map map; T& operator[](const Key& key) { - if (map.contains(key)) + if (map.contains(key)){ throw std::runtime_error( - std::format("Duplicate key in UniqueMap: {}", key)); + std::format("Duplicate key in UniqueMap: {}", key));} return map[key]; } const T& at(const Key& key) const { return map.at(key); } - const T& back() const { return map.back(); } - const T& front() const { return map.front(); } - T& back() { return map.back(); } - T& front() { return map.front(); } - std::unordered_map::const_iterator begin() const { - return map.begin(); + template + decltype(auto) begin(this Self&& self) { + return std::forward(self).map.begin(); } - std::unordered_map::const_iterator end() const { return map.end(); } - std::unordered_map::iterator begin() { return map.begin(); } - std::unordered_map::iterator end() { return map.end(); } - std::unordered_map::iterator find(const Key& key) { - return map.find(key); + + template + decltype(auto) end(this Self&& self) { + return std::forward(self).map.end(); + } + + template + decltype(auto) cbegin(this Self&& self) { + return std::forward(self).map.cbegin(); + } + + template + decltype(auto) cend(this Self&& self) { + return std::forward(self).map.cend(); + } + + template + decltype(auto) find(this Self&& self, const Key& key) { + return std::forward(self).map.find(key); + } + + template + decltype(auto) contains(this Self&& self, const Key& key) { + return std::forward(self).map.contains(key); } - std::unordered_map::const_iterator find(const Key& key) const { - return map.find(key); + + template + decltype(auto) erase(this Self&& self, const Key& key) { + return std::forward(self).map.erase(key); + } + + template + decltype(auto) erase(this Self&& self, std::unordered_map::const_iterator iter) { + return std::forward(self).map.erase(iter); } - bool contains(const Key& key) const { return map.contains(key); } - auto erase(std::unordered_map::const_iterator iter) { - return map.erase(iter); + template + decltype(auto) erase(this Self&& self, std::unordered_map::iterator iter) { + return std::forward(self).map.erase(iter); } operator std::unordered_map&() { return map; } diff --git a/src/m_abs.cc b/src/m_abs.cc index ef85f38dc2..fe6436a8a1 100644 --- a/src/m_abs.cc +++ b/src/m_abs.cc @@ -231,9 +231,8 @@ void spectral_propmatAddFaraday(PropmatVector& spectral_propmat, } } -void spectral_propmat_agendaAuto(Agenda& spectral_propmat_agenda, - const ArrayOfSpeciesTag& abs_species, - const AbsorptionBands& abs_bands, +void spectral_propmat_agendaAuto(const Workspace& ws, + Agenda& spectral_propmat_agenda, const Index& use_abs_lookup_data, const Numeric& T_extrapolfac, const Index& ignore_errors, @@ -249,7 +248,15 @@ void spectral_propmat_agendaAuto(Agenda& spectral_propmat_agenda, AgendaCreator agenda("spectral_propmat_agenda"); - const SpeciesTagTypeStatus any_species(abs_species); + const bool has_abs_species = ws.contains("abs_species"); + const bool has_abs_bands = ws.contains("abs_bands"); + const bool has_abs_cia_data = ws.contains("abs_cia_data"); + const bool has_abs_xfit_data = ws.contains("abs_xfit_data"); + const bool has_abs_predef_data = ws.contains("abs_predef_data"); + + const SpeciesTagTypeStatus any_species( + has_abs_species ? ws.get("abs_species") + : ArrayOfSpeciesTag{}); // spectral_propmatInit agenda.add("spectral_propmatInit"); @@ -263,33 +270,32 @@ void spectral_propmat_agendaAuto(Agenda& spectral_propmat_agenda, SetWsv{"water_interp_order", water_interp_order}, SetWsv{"f_interp_order", f_interp_order}, SetWsv{"extpolfac", extpolfac}); - } else if (abs_bands.size()) { + } else if (has_abs_bands) { agenda.add("spectral_propmatAddLines", SetWsv{"no_negative_absorption", no_negative_absorption}); } //spectral_propmatAddHitranXsec - if (any_species.XsecFit) { + if (any_species.XsecFit or has_abs_xfit_data) { agenda.add("spectral_propmatAddXsecFit", SetWsv{"force_p", force_p}, SetWsv{"force_t", force_t}); } //spectral_propmatAddCIA - if (any_species.Cia) { + if (any_species.Cia or has_abs_cia_data) { agenda.add("spectral_propmatAddCIA", SetWsv{"T_extrapolfac", T_extrapolfac}, SetWsv{"ignore_errors", ignore_errors}); } //spectral_propmatAddPredefined - if (any_species.Predefined) { + if (any_species.Predefined or has_abs_predef_data) { agenda.add("spectral_propmatAddPredefined"); } //spectral_propmatAddFaraday - if (stdr::any_of( - abs_species, Cmp::eq<"free_electrons"_spec>(), &SpeciesTag::Spec)) { + if (any_species.FreeElectrons) { agenda.add("spectral_propmatAddFaraday"); } diff --git a/src/m_atm.cc b/src/m_atm.cc index 02484f3fb8..8e03989501 100644 --- a/src/m_atm.cc +++ b/src/m_atm.cc @@ -200,7 +200,8 @@ void atm_fieldAppendBaseData(AtmField &atm_field, const String &deal_with_field_component, const Index &replace_existing, const Index &allow_missing_pressure, - const Index &allow_missing_temperature) { + const Index &allow_missing_temperature, + const Index &pressure_log_interp) { ARTS_TIME_REPORT std::unordered_map keys; @@ -220,11 +221,14 @@ void atm_fieldAppendBaseData(AtmField &atm_field, using enum AtmKey; - ARTS_USER_ERROR_IF( - not atm_field.contains(p) and - not static_cast(allow_missing_pressure), - "Pressure is missing from the read atmospheric field at \"{}\"", - basename) + if (atm_field.contains(p)) { + atm_field[p].log_interpolation = pressure_log_interp; + } else { + ARTS_USER_ERROR_IF( + not static_cast(allow_missing_pressure), + "Pressure is missing from the read atmospheric field at \"{}\"", + basename) + } ARTS_USER_ERROR_IF( not atm_field.contains(t) and @@ -819,3 +823,44 @@ void atm_fieldFromProfile( atm_field = Atm::atm_from_profile( atm_profile, alt_grid, altitude_extrapolation, alt_grid.back()); } + +void atm_fieldWindIncludePlanetRotation(AtmField &atm_field, + const SurfaceField &surf_field, + const Numeric &planet_rotation_period) { + ARTS_USER_ERROR_IF(not atm_field.contains(AtmKey::wind_u), + "AtmField must contain wind_u field") + + auto lat_corr = [k1 = 2 * Constant::pi / planet_rotation_period, + ell = surf_field.ellipsoid](Numeric alt, Numeric lat) { + const Numeric k2 = k1 * Conversion::cosd(lat); + const Numeric re = prime_vertical_radius(ell, lat); + + return k2 * (re + alt); + }; + + auto &wind = atm_field[AtmKey::wind_u]; + + if (auto gf3 = wind.get_if(); gf3) { + const auto &alt_grid = gf3->grid<0>(); + const auto &lat_grid = gf3->grid<1>(); + + for (Size ialt = 0; ialt < alt_grid.size(); ialt++) { + for (Size ilat = 0; ilat < lat_grid.size(); ilat++) { + gf3->data[ialt, ilat] += lat_corr(alt_grid[ialt], lat_grid[ilat]); + } + } + } else if (auto num = wind.get_if(); num) { + atm_field[AtmKey::wind_u].data = NumericTernaryOperator{ + [x = *num, lat_corr](Numeric alt, Numeric lat, Numeric) { + return x + lat_corr(alt, lat); + }}; + } else if (auto fun = wind.get_if(); fun) { + atm_field[AtmKey::wind_u].data = NumericTernaryOperator{ + [f = fun->f, lat_corr](Numeric alt, Numeric lat, Numeric lon) { + return f(alt, lat, lon) + lat_corr(alt, lat); + }}; + } else { + ARTS_USER_ERROR( + "Unsupported data type for wind_u field when including planet rotation"); + } +} diff --git a/src/m_disort.cc b/src/m_disort.cc index 3af48e94ae..544d76b6a0 100644 --- a/src/m_disort.cc +++ b/src/m_disort.cc @@ -122,7 +122,7 @@ void spectral_radFromDisort(StokvecVector& spectral_rad, }); }, variant_lag(azi_grid, ray_point.los[1]), - variant_lag(zen_grid, ray_point.los[0])); + variant_lag(zen_grid, ray_point.los[0])); } else { std::visit( [&](auto&& alt, auto&& aa, auto&& za) { @@ -131,9 +131,9 @@ void spectral_radFromDisort(StokvecVector& spectral_rad, return interp(d, alt, aa, za); }); }, - variant_lag(std::span{alt_grid}.subspan(1), z), + variant_lag(std::span{alt_grid}.subspan(1), z), variant_lag(azi_grid, ray_point.los[1]), - variant_lag(zen_grid, ray_point.los[0])); + variant_lag(zen_grid, ray_point.los[0])); } } diff --git a/src/m_lbl.cc b/src/m_lbl.cc index 23bb052410..ba9e844963 100644 --- a/src/m_lbl.cc +++ b/src/m_lbl.cc @@ -344,7 +344,7 @@ void abs_bandsReadJPL(AbsorptionBands& abs_bands, const String& filename) try { abs_bands = {}; for (auto& line : data) { - abs_bands[line.qid.qid].lines.emplace_back(line.from()); + abs_bands[line.jpl_id.qid].lines.emplace_back(line.from()); } for (auto& [_, band] : abs_bands) { diff --git a/src/m_lookup.cc b/src/m_lookup.cc index 6a97dfc268..18c85653ad 100644 --- a/src/m_lookup.cc +++ b/src/m_lookup.cc @@ -25,8 +25,6 @@ std::conditional_t _spectral_propmatAddLookup( const Index& water_interp_order, const Index& f_interp_order, const Numeric& extpolfac) { - ARTS_TIME_REPORT - if constexpr (calc) { Vector absorption(freq_grid.size(), 0.0); diff --git a/src/m_planets.cc b/src/m_planets.cc index e9d31c5834..05460046c9 100644 --- a/src/m_planets.cc +++ b/src/m_planets.cc @@ -250,7 +250,7 @@ void surf_fieldPlanet(SurfaceField &surf_field, case Earth: return surf_fieldEarth(surf_field, "WGS84", surf_elevation); case Io: return surf_fieldIo(surf_field, "Sphere", surf_elevation); case Jupiter: - return surf_fieldJupiter(surf_field, "Sphere", surf_elevation); + return surf_fieldJupiter(surf_field, "Ellipsoid", surf_elevation); case Mars: return surf_fieldMars(surf_field, "Sphere", surf_elevation); case Venus: return surf_fieldVenus(surf_field, "Sphere", surf_elevation); case Moon: return surf_fieldMoon(surf_field, "Ellipsoid", surf_elevation); @@ -273,6 +273,6 @@ void gravity_operatorCentralMass(NumericTernaryOperator &gravity_operator, gravity_operator = NumericTernaryOperator{EllipsoidGravity{ .GM = Constant::G * mass, .a = surf_field.ellipsoid[0], - .e = std::sqrt( - 1 - Math::pow2(surf_field.ellipsoid[1] / surf_field.ellipsoid[0]))}}; + .b = surf_field.ellipsoid[1], + }}; } diff --git a/src/m_ppvar.cc b/src/m_ppvar.cc index 7a81a01f40..1f0891cdb4 100644 --- a/src/m_ppvar.cc +++ b/src/m_ppvar.cc @@ -41,6 +41,8 @@ void atm_pathFromPath(ArrayOfAtmPoint &atm_path, ARTS_TIME_REPORT forward_atm_path(atm_path_resize(atm_path, ray_path), ray_path, atm_field); + + for (auto &p : atm_path) p.check(); } ARTS_METHOD_ERROR_CATCH diff --git a/src/m_profile.cc b/src/m_profile.cc new file mode 100644 index 0000000000..d3c759a698 --- /dev/null +++ b/src/m_profile.cc @@ -0,0 +1,135 @@ +#include +#include + +#include +#include + +void Profile2Path( + ArrayOfPropmatVector& spectral_propmat_path, + ArrayOfStokvecVector& spectral_nlte_srcvec_path, + ArrayOfPropmatMatrix& spectral_propmat_jac_path, + ArrayOfStokvecMatrix& spectral_nlte_srcvec_jac_path, + ArrayOfAtmPoint& atm_path, + ArrayOfAscendingGrid& freq_grid_path, + ArrayOfVector3& freq_wind_shift_jac_path, + const AscendingGrid& freq_grid, + const AscendingGrid& alt_grid, + const ArrayOfPropagationPathPoint& ray_path, + const ArrayOfAtmPoint& atm_profile, + const ArrayOfPropmatVector& spectral_propmat_profile, + const ArrayOfPropmatMatrix& spectral_propmat_jac_profile, + const ArrayOfStokvecVector& spectral_nlte_srcvec_profile, + const ArrayOfStokvecMatrix& spectral_nlte_srcvec_jac_profile) try { + ARTS_TIME_REPORT + + ARTS_USER_ERROR_IF(not arr::same_size(alt_grid, + atm_profile, + spectral_propmat_profile, + spectral_propmat_jac_profile, + spectral_nlte_srcvec_profile, + spectral_nlte_srcvec_jac_profile), + R"(Profile input must have the same size: + +alt_grid size {} +atm_profile size {} +spectral_propmat_profile size {} +spectral_propmat_jac_profile size {} +spectral_nlte_srcvec_profile size {} +spectral_nlte_srcvec_jac_profile size {} +)", + alt_grid.size(), + atm_profile.size(), + spectral_propmat_profile.size(), + spectral_propmat_jac_profile.size(), + spectral_nlte_srcvec_profile.size(), + spectral_nlte_srcvec_jac_profile.size()) + + const std::vector ips{ + std::from_range, + ray_path | stdv::transform([&](const PropagationPathPoint& ray_point) + -> Size { + const Size ip = + stdr::distance(alt_grid.begin(), + stdr::lower_bound(alt_grid, ray_point.altitude())); + + ARTS_USER_ERROR_IF( + ip == alt_grid.size(), + "Ray path point altitude {} is out of bounds of the altitude grid [{}, {}]", + ray_point.altitude(), + alt_grid.front(), + alt_grid.back()) + + ARTS_USER_ERROR_IF( + alt_grid[ip] != ray_point.altitude(), + R"(Ray path point altitude {} does not match any point in the altitude grid. + +Grid points: {:B,} +)", + ray_point.altitude(), + alt_grid) + + return ip; + })}; + + const Size n = ips.size(); + + spectral_propmat_path.resize(n); + spectral_propmat_jac_path.resize(n); + spectral_nlte_srcvec_path.resize(n); + spectral_nlte_srcvec_jac_path.resize(n); + atm_path.resize(n); + freq_grid_path.resize(n); + freq_wind_shift_jac_path.resize(n); + + for (Size i = 0; i < n; i++) { + const Size ip = ips[i]; + atm_path[i] = atm_profile[ip]; + spectral_propmat_path[i] = spectral_propmat_profile[ip]; + spectral_propmat_jac_path[i] = spectral_propmat_jac_profile[ip]; + spectral_nlte_srcvec_path[i] = spectral_nlte_srcvec_profile[ip]; + spectral_nlte_srcvec_jac_path[i] = spectral_nlte_srcvec_jac_profile[ip]; + freq_grid_path[i] = freq_grid; + freq_wind_shift_jac_path[i] = Vector3{0, 0, 0}; + } +} +ARTS_METHOD_ERROR_CATCH + +void ProfileFromAltitude(const Workspace& ws, + ArrayOfAtmPoint& atm_profile, + ArrayOfPropmatVector& spectral_propmat_profile, + ArrayOfPropmatMatrix& spectral_propmat_jac_profile, + ArrayOfStokvecVector& spectral_nlte_srcvec_profile, + ArrayOfStokvecMatrix& spectral_nlte_srcvec_jac_profile, + const Agenda& spectral_propmat_agenda, + const JacobianTargets& jac_targets, + const AscendingGrid& alt_grid, + const AtmField& atm_field, + const AscendingGrid& freq_grid, + const Numeric& lat, + const Numeric& lon) { + const ArrayOfPropagationPathPoint ray_path{ + std::from_range, + alt_grid | stdv::transform([lat, lon](const Numeric& alt) { + PropagationPathPoint p{.pos_type = PathPositionType::atm, + .los_type = PathPositionType::atm, + .pos = {alt, lat, lon}, + .los = {0, 0}}; + return p; + })}; + + const ArrayOfAscendingGrid freq_grid_path(ray_path.size(), freq_grid); + const ArrayOfVector3 freq_wind_shift_jac_path(ray_path.size(), + Vector3{0, 0, 0}); + atm_pathFromPath(atm_profile, ray_path, atm_field); + spectral_propmat_pathFromPath(ws, + spectral_propmat_profile, + spectral_nlte_srcvec_profile, + spectral_propmat_jac_profile, + spectral_nlte_srcvec_jac_profile, + spectral_propmat_agenda, + freq_grid_path, + freq_wind_shift_jac_path, + jac_targets, + ray_path, + atm_profile); +} diff --git a/src/m_propagation_path.cc b/src/m_propagation_path.cc index f6c06baf8f..a4234f75f7 100644 --- a/src/m_propagation_path.cc +++ b/src/m_propagation_path.cc @@ -109,6 +109,26 @@ void ray_pathAddGeometricGridCrossings(ArrayOfPropagationPathPoint& ray_path, path::fill_geometric_crossings(ray_path, surf_field, alt, lat, lon); } +void ray_pathAddGeometricAltitudeGridCrossings( + ArrayOfPropagationPathPoint& ray_path, + const AscendingGrid& alt_grid, + const SurfaceField& surf_field) { + ARTS_TIME_REPORT + + ARTS_USER_ERROR_IF( + surf_field.bad_ellipsoid(), + "Surface field not properly set up - bad reference ellipsoid: {:B,}", + surf_field.ellipsoid) + + path::fill_geometric_altitude_crossings(ray_path, surf_field, alt_grid); + + ARTS_USER_ERROR_IF( + stdr::any_of(ray_path, + Cmp::not_in_sorted(alt_grid), + [](const PropagationPathPoint& p) { return p.altitude(); }), + "Failed to add only geometric altitude grid crossings"); +} + void ray_pathFillGeometricHalfStep(ArrayOfPropagationPathPoint& ray_path, const SurfaceField& surf_field, const Numeric& max_step) { diff --git a/src/m_propmat.cc b/src/m_propmat.cc index 4ae22fc6aa..55d4925d00 100644 --- a/src/m_propmat.cc +++ b/src/m_propmat.cc @@ -64,7 +64,7 @@ freq_wind_shift_jac_path size: {} element(s) } ARTS_METHOD_ERROR_CATCH -void spectral_propmat_pathAdaptiveHalfPath( +void spectral_propmat_pathAddAdaptiveHalfPath( const Workspace &ws, ArrayOfPropmatVector &spectral_propmat_path, ArrayOfStokvecVector &spectral_nlte_srcvec_path, @@ -84,18 +84,6 @@ void spectral_propmat_pathAdaptiveHalfPath( const Numeric &cutoff_tau) try { ARTS_TIME_REPORT - spectral_propmat_pathFromPath(ws, - spectral_propmat_path, - spectral_nlte_srcvec_path, - spectral_propmat_jac_path, - spectral_nlte_srcvec_jac_path, - spectral_propmat_agenda, - freq_grid_path, - freq_wind_shift_jac_path, - jac_targets, - ray_path, - atm_path); - const Size nf = freq_grid.size(); if (ray_path.size() < 2) return; @@ -154,19 +142,19 @@ void spectral_propmat_pathAdaptiveHalfPath( jac_targets, nrp, nap); - ray_path.insert_range(ray_path.begin() + ip+1, nrp); - atm_path.insert_range(atm_path.begin() + ip+1 , nap); - freq_grid_path.insert_range(freq_grid_path.begin() + ip+1 , nfgp); + ray_path.insert_range(ray_path.begin() + ip + 1, nrp); + atm_path.insert_range(atm_path.begin() + ip + 1, nap); + freq_grid_path.insert_range(freq_grid_path.begin() + ip + 1, nfgp); freq_wind_shift_jac_path.insert_range( - freq_wind_shift_jac_path.begin() + ip+1 , nfwsjp); + freq_wind_shift_jac_path.begin() + ip + 1, nfwsjp); spectral_propmat_path.insert_range( - spectral_propmat_path.begin() + ip+1 , nspp); + spectral_propmat_path.begin() + ip + 1, nspp); spectral_nlte_srcvec_path.insert_range( - spectral_nlte_srcvec_path.begin() + ip+1 , nsnsp); + spectral_nlte_srcvec_path.begin() + ip + 1, nsnsp); spectral_propmat_jac_path.insert_range( - spectral_propmat_jac_path.begin() + ip+1 , nspjp); + spectral_propmat_jac_path.begin() + ip + 1, nspjp); spectral_nlte_srcvec_jac_path.insert_range( - spectral_nlte_srcvec_jac_path.begin() + ip+1 , nsnsjp); + spectral_nlte_srcvec_jac_path.begin() + ip + 1, nsnsjp); ip += nrp.size(); } } @@ -230,3 +218,26 @@ void spectral_propmat_path_species_splitFromPath( if (not error.empty()) throw std::runtime_error(error); } ARTS_METHOD_ERROR_CATCH + +void spectral_propmat_and_atm_path_agendaSetAdaptive( + Agenda &spectral_propmat_and_atm_path_agenda, + const Numeric &max_tau, + const Numeric &cutoff_tau) { + assert(stdr::contains(internal_workspace_methods() + .at("spectral_propmat_pathAddAdaptiveHalfPath") + .gin, + "max_tau")); + assert(stdr::contains(internal_workspace_methods() + .at("spectral_propmat_pathAddAdaptiveHalfPath") + .gin, + "cutoff_tau")); + assert(internal_workspace_methods() + .at("spectral_propmat_pathAddAdaptiveHalfPath") + .gin.size() == 2); + + spectral_propmat_and_atm_path_agenda = + get_spectral_propmat_and_atm_path_agenda("AdaptiveHalfPath"); + spectral_propmat_and_atm_path_agenda.change_default("max_tau"sv, max_tau); + spectral_propmat_and_atm_path_agenda.change_default("cutoff_tau"sv, + cutoff_tau); +} \ No newline at end of file diff --git a/src/m_spectral_radiance.cc b/src/m_spectral_radiance.cc index 930bf900ff..875386a3c1 100644 --- a/src/m_spectral_radiance.cc +++ b/src/m_spectral_radiance.cc @@ -883,4 +883,4 @@ void single_radFromVector(Stokvec& single_rad, single_rad_jac.resize(spectral_rad_jac.nrows()); single_rad_jac = spectral_rad_jac[joker, f]; } -ARTS_METHOD_ERROR_CATCH +ARTS_METHOD_ERROR_CATCH \ No newline at end of file diff --git a/src/m_spectral_radiance_field.cc b/src/m_spectral_radiance_field.cc index 57e87f8ffa..632720b33c 100644 --- a/src/m_spectral_radiance_field.cc +++ b/src/m_spectral_radiance_field.cc @@ -5,10 +5,10 @@ #include namespace { -auto spectral_propmat_pathProfile(const Workspace& ws, - const Agenda& spectral_propmat_agenda, - const AscendingGrid& freq_grid, - const ArrayOfAtmPoint& atm_path) { +auto spectral_propmat_profileFromAtm(const Workspace& ws, + const Agenda& spectral_propmat_agenda, + const AscendingGrid& freq_grid, + const ArrayOfAtmPoint& atm_profile) { struct spectral_propmat_pathFromPathOutput { ArrayOfPropmatVector k; ArrayOfStokvecVector s; @@ -17,7 +17,7 @@ auto spectral_propmat_pathProfile(const Workspace& ws, }; spectral_propmat_pathFromPathOutput out; - const Size np = atm_path.size(); + const Size np = atm_profile.size(); out.k.resize(np); out.s.resize(np); out.dk.resize(np); @@ -38,7 +38,7 @@ auto spectral_propmat_pathProfile(const Workspace& ws, {}, {}, {}, - atm_path[ip], + atm_profile[ip], spectral_propmat_agenda); } catch (const std::runtime_error& e) { #pragma omp critical @@ -55,7 +55,7 @@ StokvecVector interp(const StokvecConstMatrixView& data, const Numeric& zen) { const Size NF = data.shape().back(); - const auto za = zen_grid.lag<1, lagrange_interp::identity>(zen); + const auto za = zen_grid.lag<1, lagrange_interp::grid_identity>(zen); StokvecVector out(NF); @@ -115,8 +115,8 @@ void zen_gridProfilePseudo2D(ZenGrid& zen_grid, }; if (static_cast(consider_limb)) { - Vector zenith_limb = zenith_level_limb( - alt_grid, surf_field.ellipsoid, lat, lon, azi); + Vector zenith_limb = + zenith_level_limb(alt_grid, surf_field.ellipsoid, lat, lon, azi); stdr::sort(zenith_limb); const auto [it, _] = stdr::unique(zenith_limb); @@ -150,7 +150,7 @@ void spectral_rad_fieldProfilePseudo2D( const Workspace& ws, GriddedSpectralField6& spectral_rad_field, const Agenda& spectral_propmat_agenda, - const ArrayOfAtmPoint& atm_path, + const ArrayOfAtmPoint& atm_profile, const SurfaceField& surf_field, const AscendingGrid& freq_grid, const ZenGrid& zen_grid, @@ -168,14 +168,14 @@ void spectral_rad_fieldProfilePseudo2D( surf_field.ellipsoid) ARTS_USER_ERROR_IF( - not arr::same_size(alt_grid, atm_path), + not arr::same_size(alt_grid, atm_profile), R"(Altitude grid and atmospheric point grid must have the same size Altitude grid size: {} Atmospheric point grid size: {} )", alt_grid.size(), - atm_path.size()) + atm_profile.size()) ARTS_USER_ERROR_IF(zen_grid.empty(), "Need some zenith angles") @@ -204,15 +204,14 @@ Atmospheric point grid size: {} if (NA == 0) return; - const auto propmat_data = spectral_propmat_pathProfile( - ws, spectral_propmat_agenda, freq_grid, atm_path); + const auto propmat_data = spectral_propmat_profileFromAtm( + ws, spectral_propmat_agenda, freq_grid, atm_profile); constexpr Numeric t_spac = Constant::cosmic_microwave_background_temperature; const Numeric t_surf = surf_field[SurfaceKey::t].at(lat, lon); const Vector2 ell = surf_field.ellipsoid; - const Vector zenith_limb = - zenith_level_limb(alt_grid, ell, lat, lon, azi); + const Vector zenith_limb = zenith_level_limb(alt_grid, ell, lat, lon, azi); ARTS_USER_ERROR_IF( zen_grid.front() != 0.0 or zen_grid.back() != 180.0 or @@ -257,8 +256,8 @@ Atmospheric point grid size: {} const auto& K1 = propmat_data.k[end]; const auto& J0 = propmat_data.s[beg]; const auto& J1 = propmat_data.s[end]; - const auto& T0 = atm_path[beg].temperature; - const auto& T1 = atm_path[end].temperature; + const auto& T0 = atm_profile[beg].temperature; + const auto& T1 = atm_profile[end].temperature; rtepack::nlte_step(I, freq_grid, K0, K1, J0, J1, T0, T1, r); } diff --git a/src/m_sun.cc b/src/m_sun.cc index 97a4a886fa..168ec431be 100644 --- a/src/m_sun.cc +++ b/src/m_sun.cc @@ -321,7 +321,7 @@ void spectral_rad_scat_pathSunsFirstOrderRayleigh( const AscendingGrid& freq_grid, const AtmField& atm_field, const SurfaceField& surf_field, - const Agenda& spectral_propmat_agenda, + const Agenda& spectral_propmat_and_atm_path_agenda, const TransmittanceOption& rte_option, const Numeric& depolarization_factor, const Index& hse_derivative) try { @@ -381,25 +381,26 @@ void spectral_rad_scat_pathSunsFirstOrderRayleigh( const auto& suns_path = ray_path_suns_path[ip]; for (Size isun = 0; isun < nsuns; isun++) { - const auto& sun_path = suns_path[isun]; - const auto& sun = suns[isun]; + auto sun_path = suns_path[isun]; + const auto& sun = suns[isun]; spectral_radSunOrCosmicBackground( spectral_rad_bkg, freq_grid, sun_path, sun, surf_field); - spectral_radClearskyBackgroundTransmission(ws, - spectral_rad, - spectral_rad_jac, - atm_field, - freq_grid, - jac_targets, - sun_path, - rte_option, - spectral_propmat_agenda, - spectral_rad_bkg, - spectral_rad_bkg_jac, - surf_field, - hse_derivative); + spectral_radClearskyBackgroundTransmission( + ws, + spectral_rad, + spectral_rad_jac, + sun_path, + atm_field, + freq_grid, + jac_targets, + rte_option, + spectral_propmat_and_atm_path_agenda, + spectral_rad_bkg, + spectral_rad_bkg_jac, + surf_field, + hse_derivative); ARTS_USER_ERROR_IF(spectral_rad.size() != nf, "Bad size spectral_rad (", diff --git a/src/python_interface/gen_auto_py_methods.cpp b/src/python_interface/gen_auto_py_methods.cpp index ae529c1589..b0d03096dc 100644 --- a/src/python_interface/gen_auto_py_methods.cpp +++ b/src/python_interface/gen_auto_py_methods.cpp @@ -131,7 +131,7 @@ std::string method_gout_selection(const WorkspaceMethodInternalRecord& wsm) { } else if (uses_variadic(wsm.gout_type[i])) { std::println( os, - R"-x-( auto& {0} = _{0} ? *_{0} : throw std::runtime_error(R"(Unknown output: "{0}")");)-x-", + R"-x-( auto& {0} = select_gout(_{0}, _ws, "{0}");)-x-", wsm.gout[i], wsm.gout_type[i]); } else { @@ -156,7 +156,7 @@ std::string method_gin_selection(const std::string& name, } else if (uses_variadic(wsm.gin_type[i])) { std::println( os, - R"-x-( const auto& {0} = _{0} ? *_{0} : throw std::runtime_error(R"(Unknown input: "{0}")");)-x-", + R"-x-( const auto& {0} = select_gin(_{0}, "{0}");)-x-", wsm.gin[i]); } else { if (has_default) { @@ -569,7 +569,7 @@ std::string method_error(const std::string& name, t)); } else if (uses_variadic(tt)) { arg.push_back(std::format( - R"(_{0} ? std::format("User-provided {{}}", type(_{0})) : std::string("None"))", + R"(_{0} and has_selected_value(*_{0}) ? std::format("User-provided {{}}", type(_{0})) : std::string("None"))", t)); } else { arg.push_back( @@ -602,6 +602,11 @@ std::string method_error(const std::string& name, R"(_{0} ? std::format("User-provided {{}}", _{0} -> type_name()) : std::string(R"-WSMVAR-({1})-WSMVAR-"))", t, v ? std::format("{}", to_defval_str(*v, ""sv)) : "None")); + } else if (uses_variadic(tt)) { + arg.push_back(std::format( + R"(_{0} and has_selected_value(*_{0}) ? std::format("User-provided {{}}", type(_{0})) : std::string(R"-WSMVAR-({1})-WSMVAR-"))", + t, + v ? std::format("{}", to_defval_str(*v, ""sv)) : "None")); } else { arg.push_back(std::format( R"(_{0} ? std::format("User-provided {{}}", type(_{0})) : std::string(R"-WSMVAR-({1})-WSMVAR-"))", diff --git a/src/python_interface/hpy_opaque.h b/src/python_interface/hpy_opaque.h index 8d8831910d..9ce150c0e5 100644 --- a/src/python_interface/hpy_opaque.h +++ b/src/python_interface/hpy_opaque.h @@ -3,7 +3,7 @@ #include #include -NB_MAKE_OPAQUE(Array>); +NB_MAKE_OPAQUE(Array>); NB_MAKE_OPAQUE(Array>); NB_MAKE_OPAQUE(std::unordered_map); NB_MAKE_OPAQUE(QuantumLevel); diff --git a/src/python_interface/py_agenda.cpp b/src/python_interface/py_agenda.cpp index 3d24d21998..68029fcb9c 100644 --- a/src/python_interface/py_agenda.cpp +++ b/src/python_interface/py_agenda.cpp @@ -9,9 +9,9 @@ #include #include #include -#include #include #include +#include #include #include @@ -128,7 +128,8 @@ void py_agenda(py::module_& m) try { "The output variables of the method.\n\n.. :class:`list` of :class:`str`") .doc() = "The method class of ARTS"; - auto wsvmap = py::bind_map>(m, "WsvMap"); + auto wsvmap = py::bind_map, + py::rv_policy::reference_internal>(m, "WsvMap"); wsvmap.doc() = unwrap_stars("A map from *String* to :class:`~pyarts3.arts.Wsv`"); wsvmap.def( diff --git a/src/python_interface/py_atm.cpp b/src/python_interface/py_atm.cpp index dda348e347..b61de5574a 100644 --- a/src/python_interface/py_atm.cpp +++ b/src/python_interface/py_atm.cpp @@ -72,6 +72,9 @@ void py_atm(py::module_ &m) try { "lon_low", &Atm::Data::lon_low, "Lower longitude limit\n\n.. :class:`~pyarts3.arts.InterpolationExtrapolation`") + .def_rw("log_interpolation", + &Atm::Data::log_interpolation, + "Whether to interpolate in log-space.\n\n.. :class:`bool`") .def( "set_extrapolation", [](Atm::Data &self, InterpolationExtrapolation x) { diff --git a/src/python_interface/py_interp.cpp b/src/python_interface/py_interp.cpp index 69515f3406..0fa4ea3754 100644 --- a/src/python_interface/py_interp.cpp +++ b/src/python_interface/py_interp.cpp @@ -13,7 +13,7 @@ #include "python_interface.h" namespace Python { -using id = lagrange_interp::lag_t<-1, lagrange_interp::identity>; +using id = lagrange_interp::lag_t<-1, lagrange_interp::grid_identity>; using lc = lagrange_interp::lag_t<-1, lagrange_interp::loncross>; using PythonLag = std::variant; diff --git a/src/python_interface/py_lbl.cpp b/src/python_interface/py_lbl.cpp index 07b1ec3051..16c130a182 100644 --- a/src/python_interface/py_lbl.cpp +++ b/src/python_interface/py_lbl.cpp @@ -29,13 +29,15 @@ namespace Python { void py_lbl(py::module_& m) try { auto lbl = m.def_submodule("lbl", "Line-by-line helper functions"); - auto lssmm = py::bind_map( + auto lssmm = py::bind_map( lbl, "LineShapeSpeciesModelMap"); lssmm.doc() = "A map from model variable to line shape models"; generic_interface(lssmm); auto lsmm = - py::bind_map(lbl, "LineShapeModelMap"); + py::bind_map(lbl, "LineShapeModelMap"); lsmm.doc() = "A map from species to species line shape models"; generic_interface(lsmm); @@ -989,10 +991,12 @@ fmax : ~pyarts3.arts.Numeric R"(The Omega coefficient for the ECS model)"); auto lsed = - py::bind_map(m, "LinemixingSpeciesEcsData"); + py::bind_map( + m, "LinemixingSpeciesEcsData"); generic_interface(lsed); - auto led = py::bind_map(m, "LinemixingEcsData"); + auto led = py::bind_map( + m, "LinemixingEcsData"); generic_interface(led); lbl.def( diff --git a/src/python_interface/py_lookup.cpp b/src/python_interface/py_lookup.cpp index 2810c47c35..3ef0a55f43 100644 --- a/src/python_interface/py_lookup.cpp +++ b/src/python_interface/py_lookup.cpp @@ -35,7 +35,9 @@ void py_lookup(py::module_& m) try { &AbsorptionLookupTable::xsec, "The absorption cross section table\n\n.. :class:`Tensor4`"); - auto alts = py::bind_map(m, "AbsorptionLookupTables"); + auto alts = + py::bind_map( + m, "AbsorptionLookupTables"); generic_interface(alts); alts.def( diff --git a/src/python_interface/py_nlte.cpp b/src/python_interface/py_nlte.cpp index 7c8b431fa1..110fd4dab9 100644 --- a/src/python_interface/py_nlte.cpp +++ b/src/python_interface/py_nlte.cpp @@ -7,16 +7,19 @@ namespace Python { void py_nlte(py::module_& m) try { - auto vec = - py::bind_map(m, "QuantumIdentifierVectorMap"); + auto vec = py::bind_map( + m, "QuantumIdentifierVectorMap"); generic_interface(vec); - auto num = - py::bind_map(m, "QuantumIdentifierNumericMap"); + auto num = py::bind_map( + m, "QuantumIdentifierNumericMap"); generic_interface(num); - auto gf1 = - py::bind_map(m, "QuantumIdentifierGriddedField1Map"); + auto gf1 = py::bind_map( + m, "QuantumIdentifierGriddedField1Map"); generic_interface(gf1); } catch (std::exception& e) { throw std::runtime_error( diff --git a/src/python_interface/py_planets.cpp b/src/python_interface/py_planets.cpp index 1b4bfadfdc..2dcd80e927 100644 --- a/src/python_interface/py_planets.cpp +++ b/src/python_interface/py_planets.cpp @@ -6,10 +6,21 @@ #include "hpy_numpy.h" namespace Python { +#define ADD_PLANET(PLANET) \ + auto PLANET = planets.def_submodule(#PLANET, #PLANET " constants"); \ + PLANET.attr("ellipsoid") = Vector2{Body::PLANET::a, Body::PLANET::b}; \ + PLANET.attr("GM") = Body::PLANET::GM; \ + PLANET.attr("day") = Body::PLANET::day; \ + PLANET.def( \ + "gravity_operator", \ + &EllipsoidGravity::PLANET, \ + "Get " #PLANET \ + "'s gravity object.\n\n.. :class:`~pyarts3.arts.EllipsoidGravity`"); void py_planets(py::module_& m) try { py::class_ eg(m, "EllipsoidGravity"); - eg.doc() = "Ellipsoid gravity operator - allows calculation of gravity on ellipsoidal bodies"; + eg.doc() = + "Ellipsoid gravity operator - allows calculation of gravity on ellipsoidal bodies"; generic_interface(eg); eg.def(py::init<>()) .def_rw( @@ -19,9 +30,9 @@ void py_planets(py::module_& m) try { .def_rw("a", &EllipsoidGravity::a, "Semi-major axis in m.\n\n.. :class:`~pyarts3.arts.Numeric`") - .def_rw("e", - &EllipsoidGravity::e, - "Eccentricity.\n\n.. :class:`~pyarts3.arts.Numeric`") + .def_rw("b", + &EllipsoidGravity::b, + "Semi-minor axis in m.\n\n.. :class:`~pyarts3.arts.Numeric`") .def( "__call__", [](const EllipsoidGravity& self, @@ -59,65 +70,17 @@ void py_planets(py::module_& m) try { "Saturn ellipsoid gravity operator"); auto planets = m.def_submodule("planets", "Planetary constants"); - auto earth = planets.def_submodule("Earth", "Earth constants"); - earth.attr("ellipsoid") = Vector2{Body::Earth::a, Body::Earth::b}; - earth.attr("GM") = Body::Earth::GM; - earth.def( - "gravity", - &EllipsoidGravity::Earth, - "Get Earth's gravity object.\n\n.. :class:`~pyarts3.arts.EllipsoidGravity`"); - - auto jupiter = planets.def_submodule("Jupiter", "Jupiter constants"); - jupiter.attr("ellipsoid") = Vector2{Body::Jupiter::a, Body::Jupiter::b}; - jupiter.attr("GM") = Body::Jupiter::GM; - jupiter.def( - "gravity", - &EllipsoidGravity::Jupiter, - "Get Jupiter's gravity object.\n\n.. :class:`~pyarts3.arts.EllipsoidGravity`"); - - auto mars = planets.def_submodule("Mars", "Mars constants"); - mars.attr("ellipsoid") = Vector2{Body::Mars::a, Body::Mars::b}; - mars.attr("GM") = Body::Mars::GM; - mars.def( - "gravity", - &EllipsoidGravity::Mars, - "Get Mars' gravity object.\n\n.. :class:`~pyarts3.arts.EllipsoidGravity`"); - - auto moon = planets.def_submodule("Moon", "Moon constants"); - moon.attr("ellipsoid") = Vector2{Body::Moon::a, Body::Moon::b}; - moon.attr("GM") = Body::Moon::GM; - moon.def( - "gravity", - &EllipsoidGravity::Moon, - "Get Moon's gravity object.\n\n.. :class:`~pyarts3.arts.EllipsoidGravity`"); - - auto mercury = planets.def_submodule("Mercury", "Mercury constants"); - mercury.attr("ellipsoid") = Vector2{Body::Mercury::a, Body::Mercury::b}; - mercury.attr("GM") = Body::Mercury::GM; - mercury.def( - "gravity", - &EllipsoidGravity::Mercury, - "Get Mercury's gravity object.\n\n.. :class:`~pyarts3.arts.EllipsoidGravity`"); - - auto venus = planets.def_submodule("Venus", "Venus constants"); - venus.attr("ellipsoid") = Vector2{Body::Venus::a, Body::Venus::b}; - venus.attr("GM") = Body::Venus::GM; - venus.def( - "gravity", - &EllipsoidGravity::Venus, - "Get Venus' gravity object.\n\n.. :class:`~pyarts3.arts.EllipsoidGravity`"); - - auto saturn = planets.def_submodule("Saturn", "Saturn constants"); - saturn.attr("ellipsoid") = Vector2{Body::Saturn::a, Body::Saturn::b}; - saturn.attr("GM") = Body::Saturn::GM; - saturn.def( - "gravity", - &EllipsoidGravity::Saturn, - "Get Saturn's gravity object.\n\n.. :class:`~pyarts3.arts.EllipsoidGravity`"); - + ADD_PLANET(Earth); + ADD_PLANET(Jupiter); + ADD_PLANET(Mars); + ADD_PLANET(Moon); + ADD_PLANET(Mercury); + ADD_PLANET(Venus); + ADD_PLANET(Saturn); } catch (std::exception& e) { throw std::runtime_error( std::format("DEV ERROR:\nCannot initialize planets\n{}", e.what())); } +#undef ADD_PLANET } // namespace Python \ No newline at end of file diff --git a/src/python_interface/py_predefined.cpp b/src/python_interface/py_predefined.cpp index 414d2cfb55..0be4ea8a36 100644 --- a/src/python_interface/py_predefined.cpp +++ b/src/python_interface/py_predefined.cpp @@ -34,7 +34,8 @@ void internalCKDMT430(py::module_& m) { &Absorption::PredefinedModel::MT_CKD430::WaterData::for_absco_ref, "Foreign absorption\n\n.. :class:`Vector`") .def_rw("for_closure_absco_ref", - &Absorption::PredefinedModel::MT_CKD430::WaterData::for_closure_absco_ref, + &Absorption::PredefinedModel::MT_CKD430::WaterData:: + for_closure_absco_ref, "Foreign absorption closure\n\n.. :class:`Vector`") .def_rw("wavenumbers", &Absorption::PredefinedModel::MT_CKD430::WaterData::wavenumbers, @@ -1027,7 +1028,9 @@ void py_predefined(py::module_& m) try { generic_interface(var); //! ARTS Workspace class, must live on the main (m) namespace - auto pdmd = py::bind_map(m, "PredefinedModelData"); + auto pdmd = + py::bind_map( + m, "PredefinedModelData"); generic_interface(pdmd); pdmd.def_static( "fromcatalog", diff --git a/src/python_interface/py_sensor.cpp b/src/python_interface/py_sensor.cpp index 605e3399b2..2cba96824e 100644 --- a/src/python_interface/py_sensor.cpp +++ b/src/python_interface/py_sensor.cpp @@ -421,13 +421,16 @@ See :meth:`SensorObsel.normalize` for details. generic_interface(a2); vector_interface(a2); - py::class_ sapf(m, "SensorAntennaPatternField"); + auto sen = m.def_submodule("sensor", + "Classes and functions to help create a sensor"); + + py::class_ sapf(sen, "AntennaPatternField"); sapf.doc() = "A 2D gridded antenna response on local zenith and azimuth offsets."; gridded_data_interface(sapf); generic_interface(sapf); - auto sap = py::class_(m, "SensorAntennaPattern"); + auto sap = py::class_(sen, "AntennaPattern"); sap.doc() = "Base class for angular antenna responses defined around a bore line of sight."; sap.def("__call__", @@ -439,7 +442,7 @@ See :meth:`SensorObsel.normalize` for details. // generic_interface(sap); -- this should break things. The class is abstract auto sgp = py::class_( - m, "SensorGriddedAntennaPattern"); + sen, "GriddedAntennaPattern"); sgp.doc() = "A gridded angular antenna response on local zenith and azimuth offsets."; sgp.def(py::init<>(), "Construct an empty gridded antenna pattern.") @@ -456,11 +459,11 @@ See :meth:`SensorObsel.normalize` for details. "response", &sensor::GriddedAntennaPattern::data, py::rv_policy::reference_internal, - "Local antenna response field.\n\n.. :class:`~pyarts3.arts.SensorAntennaPatternField`"); + "Local antenna response field.\n\n.. :class:`~pyarts3.arts.sensor.AntennaPatternField`"); generic_interface(sgp); auto spp = py::class_( - m, "SensorPencilBeamAntenna"); + sen, "PencilBeamAntenna"); spp.doc() = "A 1x1 antenna response centered on the bore line of sight."; spp.def( py::init(), @@ -469,28 +472,28 @@ See :meth:`SensorObsel.normalize` for details. generic_interface(spp); auto sga = py::class_( - m, "SensorGaussianAntenna"); + sen, "GaussianAntenna"); sga.doc() = - "A Gaussian antenna response defined on a local zenith/azimuth grid."; - sga.def(py::init(), - "zen_grid"_a, - "azi_grid"_a, - "zenith_std"_a, - "azimuth_std"_a, - "weight"_a = Stokvec{1.0, 0.0, 0.0, 0.0}, - "Construct a Gaussian antenna response on the supplied local grids."); + "A Gaussian antenna response defined by angular offset from the bore LOS on a local zenith/azimuth grid."; + sga.def( + py::init(), + "grid"_a, + "std"_a, + "azi_grid_size"_a = Size{2}, + "weight"_a = Stokvec{1.0, 0.0, 0.0, 0.0}, + "Construct a Gaussian antenna response on the supplied local grid. The azimuth grid is constructed from the provided size"); generic_interface(sga); auto sgairy = py::class_( - m, "SensorGaussianAiryAntenna"); + sen, "GaussianAiryAntenna"); sgairy.doc() = "A Gaussianized Airy antenna whose width scales with wavelength and aperture diameter."; sgairy.def( - py::init(), - "dzen_grid"_a, - "azi_grid"_a, + py::init(), + "grid"_a, "aperture_diameter"_a, - "weight"_a = Stokvec{1.0, 0.0, 0.0, 0.0}, + "azi_grid_size"_a = Size{2}, + "weight"_a = Stokvec{1.0, 0.0, 0.0, 0.0}, "Construct a Gaussian Airy antenna on the supplied local grids." " The channel frequency grid must be strictly positive because the" " current builder path does not carry a separate reference frequency."); @@ -513,7 +516,7 @@ See :meth:`SensorObsel.normalize` for details. generic_interface(sgairy); static_assert(arts_xml_ioable); - auto sch = py::class_(m, "SensorChannel"); + auto sch = py::class_(sen, "Channel"); sch.doc() = "Base class for relative spectrometer channel responses."; generic_interface(sch); sch.def(py::init_implicit()); @@ -536,7 +539,7 @@ See :meth:`SensorObsel.normalize` for details. "\n\n.. :class:`Vector`"); auto sbox = - py::class_(m, "SensorBoxChannel"); + py::class_(sen, "BoxChannel"); sbox.doc() = "A channel with uniform weights across a finite relative-frequency interval."; sbox.def(py::init(), @@ -554,8 +557,8 @@ See :meth:`SensorObsel.normalize` for details. "Construct a box channel directly from a relative frequency grid."); generic_interface(sbox); - auto sdirac = py::class_( - m, "SensorDiracChannel"); + auto sdirac = + py::class_(sen, "DiracChannel"); sdirac.doc() = "A single-frequency relative channel."; sdirac .def(py::init(), @@ -565,7 +568,7 @@ See :meth:`SensorObsel.normalize` for details. generic_interface(sdirac); auto sgauss = py::class_( - m, "SensorGaussianChannel"); + sen, "GaussianChannel"); sgauss.doc() = "A Gaussian relative channel response."; sgauss .def(py::init(), @@ -591,7 +594,7 @@ See :meth:`SensorObsel.normalize` for details. "Construct a zero-centered Gaussian channel on ``+/- m * std``."); generic_interface(sgauss); - auto sspec = py::class_(m, "SensorSpectrometer"); + auto sspec = py::class_(sen, "Spectrometer"); sspec.doc() = "A spectrometer made from channels arranged on one shared relative-frequency grid."; sspec @@ -600,6 +603,10 @@ See :meth:`SensorObsel.normalize` for details. "base_channel"_a, "freq_offsets"_a, "Construct a spectrometer by shifting one base channel across the supplied relative-frequency offsets.") + .def( + py::init_implicit(), + "freq_offsets"_a, + "Construct a spectrometer with Dirac channels across the supplied relative-frequency offsets.") .def(py::init>(), "channels"_a, "Construct a spectrometer from explicit channels.") @@ -608,7 +615,7 @@ See :meth:`SensorObsel.normalize` for details. [](const sensor::Spectrometer& self) -> const std::vector& { return self.channels; }, py::rv_policy::reference_internal, - "Spectrometer channels.\n\n.. :class:`list[~pyarts3.arts.SensorChannel]`") + "Spectrometer channels.\n\n.. :class:`list[~pyarts3.arts.sensor.Channel]`") .def( "is_synced", &sensor::Spectrometer::is_synced, @@ -633,38 +640,37 @@ See :meth:`SensorObsel.normalize` for details. "Iterate over the spectrometer channels."); generic_interface(sspec); - auto ssb = py::class_(m, "SensorBuilder"); + auto ssb = py::class_(sen, "Builder"); ssb.doc() = "Combines channels and an antenna pattern into sensor obsels for paired position/LOS samples."; ssb.def(py::init<>(), "Construct an empty sensor builder.") .def( "__init__", - [](sensor::SensorBuilder* self, + [](sensor::Builder* self, const std::vector& channels, const sensor::AntennaPattern& antenna) { - new (self) sensor::SensorBuilder{channels, antenna.clone()}; + new (self) sensor::Builder{channels, antenna.clone()}; }, "channels"_a, "antenna"_a = sensor::PencilBeamAntenna{}, "Construct a sensor builder from channels and one antenna pattern (defaults to pencil-beam).") .def( "__init__", - [](sensor::SensorBuilder* self, + [](sensor::Builder* self, const sensor::Spectrometer& spectrometer, const sensor::AntennaPattern& antenna) { - new (self) sensor::SensorBuilder{spectrometer, antenna.clone()}; + new (self) sensor::Builder{spectrometer, antenna.clone()}; }, "spectrometer"_a, "antenna"_a = sensor::PencilBeamAntenna{}, "Construct a sensor builder from a spectrometer and one antenna pattern (defaults to pencil-beam).") .def( "__init__", - [](sensor::SensorBuilder* self, + [](sensor::Builder* self, const sensor::AntennaPattern& antenna, const sensor::Spectrometer& spectrometer, const sensor::HeterodyneFrequencyRange& backend) { - new (self) - sensor::SensorBuilder{spectrometer, backend, antenna.clone()}; + new (self) sensor::Builder{spectrometer, backend, antenna.clone()}; }, "antenna"_a, "spectrometer"_a, @@ -672,30 +678,29 @@ See :meth:`SensorObsel.normalize` for details. "Construct a sensor builder from an antenna, a spectrometer, and a backend frequency response.") .def_rw( "channels", - &sensor::SensorBuilder::channels, - "Spectrometer channels.\n\n.. :class:`list[~pyarts3.arts.SensorChannel]`") + &sensor::Builder::channels, + "Spectrometer channels.\n\n.. :class:`list[~pyarts3.arts.sensor.Channel]`") .def_prop_rw( "antenna", - [](const sensor::SensorBuilder& self) + [](const sensor::Builder& self) -> std::shared_ptr { if (not self.antenna) - throw std::runtime_error("SensorBuilder has no antenna pattern"); + throw std::runtime_error("Sensor builder has no antenna pattern"); return self.antenna->clone(); }, - [](sensor::SensorBuilder& self, - const sensor::AntennaPattern& antenna) { + [](sensor::Builder& self, const sensor::AntennaPattern& antenna) { self.antenna = antenna.clone(); }, py::rv_policy::reference_internal, - "Angular antenna response.\n\n.. :class:`~pyarts3.arts.SensorAntennaPattern`") + "Angular antenna response.\n\n.. :class:`~pyarts3.arts.sensor.AntennaPattern`") .def( "__call__", - [](const sensor::SensorBuilder& self, + [](const sensor::Builder& self, const std::vector& pos, const std::vector& los) { if (pos.size() != los.size()) { throw std::invalid_argument(std::format( - "SensorBuilder expects matching position and LOS counts. Got {} positions and {} LOS values.", + "Sensor builder expects matching position and LOS counts. Got {} positions and {} LOS values.", pos.size(), los.size())); } @@ -711,7 +716,7 @@ geometry first and channel second, and the returned value is ``(measurement_sensor, measurement_sensor_meta)``.)") .def( "__call__", - [](const sensor::SensorBuilder& self, + [](const sensor::Builder& self, const Vector3 pos, const Vector2 los) { return self(std::span{&pos, 1}, std::span{&los, 1}); @@ -721,8 +726,7 @@ geometry first and channel second, and the returned value is R"(Build sensor obsels and metadata from paired position and bore-LOS.)") .def( "__call__", - [](const sensor::SensorBuilder& self, - const SensorPosLosVector& poslos) { + [](const sensor::Builder& self, const SensorPosLosVector& poslos) { std::vector pos(poslos.size()); std::vector los(poslos.size()); @@ -738,7 +742,7 @@ geometry first and channel second, and the returned value is generic_interface(ssb); auto shdfr = py::class_( - m, "SensorHeterodyneFrequencyRange"); + sen, "HeterodyneFrequencyRange"); shdfr.def(py::init<>(), R"(Construct an empty staged heterodyne response. @@ -822,10 +826,10 @@ geometry first and channel second, and the returned value is "channel_response", [](const sensor::HeterodyneFrequencyRange& self, const sensor::Channel& channel) { - return sensor::make_bandpass_channels( - self, std::span{&channel, 1}) - .front() - .channel; + return sensor::make_bandpass_channels( + self, std::span{&channel, 1}) + .front() + .channel; }, "channel"_a, R"(Compute the real-frequency response for one spectrometer channel. diff --git a/src/python_interface/py_species.cpp b/src/python_interface/py_species.cpp index 5bc89c7bd8..39e2ca0268 100644 --- a/src/python_interface/py_species.cpp +++ b/src/python_interface/py_species.cpp @@ -227,7 +227,9 @@ Returns ////////////////////////////////////////////////////////////////////// - auto sev = py::bind_map(m, "SpeciesEnumVectors"); + auto sev = + py::bind_map( + m, "SpeciesEnumVectors"); generic_interface(sev); py::class_ isinfo(m, "SpeciesIsotopologueInfo"); diff --git a/src/python_interface/py_surf.cpp b/src/python_interface/py_surf.cpp index e7fa2f99de..a8716d0531 100644 --- a/src/python_interface/py_surf.cpp +++ b/src/python_interface/py_surf.cpp @@ -248,8 +248,8 @@ void py_surf(py::module_ &m) try { &SubsurfaceField::other, "Other data in the subsurface field\n\n.. :class:`dict[SubsurfaceKey, SubsurfaceData]`"); ssf.def_rw( - "prop", - &SubsurfaceField::prop, + "props", + &SubsurfaceField::props, "Properties of the subsurface field\n\n.. :class:`dict[SubsurfacePropertyTag, SubsurfaceData]`"); ssf.def_rw( "bottom_depth", @@ -323,8 +323,8 @@ void py_surf(py::module_ &m) try { &SubsurfacePoint::density, "Density [kg/m^3]\n\n.. :class:`Numeric`"); ssp.def_rw( - "prop", - &SubsurfacePoint::prop, + "props", + &SubsurfacePoint::props, "Properties of the subsurface point\n\n.. :class:`dict[SubsurfacePropertyTag, Numeric]`"); ssp.def("__getitem__", [](SubsurfacePoint &self, const SubsurfaceKeyVal &key) { diff --git a/src/python_interface/python_interface.h b/src/python_interface/python_interface.h index da92387a9a..c3ad991af7 100644 --- a/src/python_interface/python_interface.h +++ b/src/python_interface/python_interface.h @@ -6,6 +6,9 @@ #include #include +#include +#include + #include "hpy_opaque.h" using ssize_t = Py_ssize_t; @@ -15,6 +18,11 @@ namespace Python { namespace py = nanobind; using namespace py::literals; +template +bool has_selected_value(const std::variant...>& x) { + return std::visit([](const auto& ptr) { return static_cast(ptr); }, x); +} + //////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////// @@ -41,6 +49,18 @@ T& select_gout(ValueHolder* const x, Workspace& ws, const char* const name) { return x ? static_cast(*x) : ws.get_or(name); } +template +std::variant...>& select_gout( + std::variant...>* const x, + Workspace&, + const char* const name) { + if (not x or not has_selected_value(*x)) { + throw std::runtime_error(std::format("Unknown output: \"{}\"", name)); + } + + return *x; +} + //////////////////////////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////////////////////////// @@ -90,6 +110,17 @@ const T& select_gin(const ValueHolder* const x, const char* const name) { std::format("Unknown input: \"{}\"", name)); } +template +const std::variant...>& select_gin( + const std::variant...>* const x, + const char* const name) { + if (not x or not has_selected_value(*x)) { + throw std::runtime_error(std::format("Unknown input: \"{}\"", name)); + } + + return *x; +} + template const T& select_gin(const T* const x, const T& defval) { return x ? *x : defval; diff --git a/src/tests/test_antenna_pattern.cc b/src/tests/test_antenna_pattern.cc index 9ea5f46cae..d2fe2a116c 100644 --- a/src/tests/test_antenna_pattern.cc +++ b/src/tests/test_antenna_pattern.cc @@ -1,5 +1,4 @@ #include - #include #include @@ -29,10 +28,9 @@ Stokvec sum_column_weights(const sensor::Obsel& obsel, Size ifreq) { return sum; } -void test_gaussian_initialization_uses_antenna_frame_offsets() { +void test_gaussian_initialization_uses_los_offset() { const Stokvec peak_weight{2.0, 1.0, 0.0, 0.0}; - const sensor::GaussianAntenna ant{ - ZenGrid{{0.0, 1.0}}, AziGrid{{0.0, 180.0}}, 1.0, 2.0, peak_weight}; + const sensor::GaussianAntenna ant{ZenGrid{{0.0, 1.0}}, 1.0, 2, peak_weight}; assert_stokvec(ant.data[0, 0], peak_weight, 1e-12, "gaussian peak weight"); @@ -41,7 +39,7 @@ void test_gaussian_initialization_uses_antenna_frame_offsets() { assert_stokvec(ant.data[1, 0], expected, 1e-12, "gaussian local [1, 0]"); assert_stokvec(ant.data[1, 1], expected, 1e-12, "gaussian local [1, 180]"); - const sensor::GaussianAntenna def{ZenGrid{{0.0}}, AziGrid{{0.0}}, 1.0, 1.0}; + const sensor::GaussianAntenna def{ZenGrid{{0.0}}, 1.0, 1}; assert_stokvec(def.data[0, 0], {1.0, 0.0, 0.0, 0.0}, 1e-12, @@ -54,28 +52,17 @@ Numeric gaussian_airy_expected_gain(Numeric zenith_deg, constexpr Numeric gaussian_airy_hwhm_factor = 3.8317059702075123156 / Constant::pi; const Numeric wavelength = Constant::speed_of_light / frequency; - const Numeric hwhm_deg = Conversion::rad2deg( - gaussian_airy_hwhm_factor * wavelength / aperture_diameter); - const Numeric ratio = zenith_deg / hwhm_deg; + const Numeric hwhm_deg = Conversion::rad2deg(gaussian_airy_hwhm_factor * + wavelength / aperture_diameter); + const Numeric ratio = zenith_deg / hwhm_deg; return std::exp(-std::log(2.0) * ratio * ratio); } -Numeric gaussian_airy_expected_std(Numeric frequency, - Numeric aperture_diameter) { - constexpr Numeric gaussian_airy_hwhm_factor = - 3.8317059702075123156 / Constant::pi; - const Numeric wavelength = Constant::speed_of_light / frequency; - const Numeric hwhm_deg = Conversion::rad2deg( - gaussian_airy_hwhm_factor * wavelength / aperture_diameter); - - return hwhm_deg / std::sqrt(2.0 * std::log(2.0)); -} - void test_gaussian_airy_is_frequency_dependent() { const Stokvec peak_weight{2.0, 0.0, 0.0, 0.0}; const sensor::GaussianAiryAntenna ant{ - ZenGrid{{0.0, 0.2}}, AziGrid{{0.0}}, 1.0, peak_weight}; + ZenGrid{{0.0, 0.2}}, 1.0, 1, peak_weight}; const sensor::BoxChannel channel{AscendingGrid{100.0e9, 200.0e9}}; const auto obsel = ant(channel, {600e3, 10.0, 20.0}, {45.0, 30.0}); @@ -83,31 +70,35 @@ void test_gaussian_airy_is_frequency_dependent() { const Numeric low_gain = gaussian_airy_expected_gain(0.2, 100.0e9, 1.0); const Numeric high_gain = gaussian_airy_expected_gain(0.2, 200.0e9, 1.0); - const Numeric low_norm = 1.0 + low_gain; - const Numeric high_norm = 1.0 + high_gain; + const Numeric low_norm = 1.0 + low_gain; + const Numeric high_norm = 1.0 + high_gain; - assert_stokvec(obsel.weight_matrix()[0, 0], - (0.5 / low_norm) * peak_weight, + assert_stokvec(obsel.weight_matrix()[0, 0], + (0.5 / low_norm) * peak_weight, 1e-12, "gaussian airy bore low frequency"); assert_stokvec(obsel.weight_matrix()[0, 1], - (0.5 / high_norm) * peak_weight, + (0.5 / high_norm) * peak_weight, 1e-12, "gaussian airy bore high frequency"); assert_stokvec(obsel.weight_matrix()[1, 0], - (0.5 * low_gain / low_norm) * peak_weight, + (0.5 * low_gain / low_norm) * peak_weight, 1e-12, "gaussian airy off-axis low frequency"); assert_stokvec(obsel.weight_matrix()[1, 1], - (0.5 * high_gain / high_norm) * peak_weight, + (0.5 * high_gain / high_norm) * peak_weight, 1e-12, "gaussian airy off-axis high frequency"); - assert_stokvec( - sum_column_weights(obsel, 0), 0.5 * peak_weight, 1e-12, "gaussian airy normalized low frequency"); - assert_stokvec( - sum_column_weights(obsel, 1), 0.5 * peak_weight, 1e-12, "gaussian airy normalized high frequency"); + assert_stokvec(sum_column_weights(obsel, 0), + 0.5 * peak_weight, + 1e-12, + "gaussian airy normalized low frequency"); + assert_stokvec(sum_column_weights(obsel, 1), + 0.5 * peak_weight, + 1e-12, + "gaussian airy normalized high frequency"); const auto low_weight = obsel.weight_matrix()[1, 0]; const auto high_weight = obsel.weight_matrix()[1, 1]; @@ -115,87 +106,73 @@ void test_gaussian_airy_is_frequency_dependent() { "Higher frequency Gaussian Airy sample must be narrower") } -void test_gaussian_airy_matches_frequency_specific_gaussian_pattern() { - const ZenGrid zen_grid{{0.1, 0.2}}; - const AziGrid azi_grid{{0.0, 0.2}}; - const Stokvec peak_weight{2.0, 0.0, 0.0, 0.0}; - const sensor::GaussianAiryAntenna ant{zen_grid, azi_grid, 1.0, peak_weight}; - const sensor::BoxChannel channel{AscendingGrid{100.0e9, 200.0e9}}; - - const auto obsel = ant(channel, {600e3, 10.0, 20.0}, {45.0, 30.0}); +void test_degenerate_azimuth_ring_is_collapsed() { + const sensor::GaussianAntenna gaussian{ZenGrid{{0.0, 1.0}}, 1.0, 3}; + const auto gaussian_obsel = gaussian( + sensor::DiracChannel{100.0e9}, {600e3, 10.0, 20.0}, {45.0, 30.0}); - for (Size ifreq = 0; ifreq < channel.freq_grid().size(); ++ifreq) { - const Numeric airy_std = - gaussian_airy_expected_std(channel.freq_grid()[ifreq], 1.0); - const sensor::GaussianAntenna gaussian{ - zen_grid, azi_grid, airy_std, airy_std, peak_weight}; - - Numeric normalization = 0.0; - for (Size izen = 0; izen < zen_grid.size(); ++izen) { - for (Size iazi = 0; iazi < azi_grid.size(); ++iazi) { - normalization += gaussian.data[izen, iazi][0] / peak_weight[0]; - } - } - - Size isample = 0; - for (Size izen = 0; izen < zen_grid.size(); ++izen) { - for (Size iazi = 0; iazi < azi_grid.size(); ++iazi) { - assert_stokvec(obsel.weight_matrix()[isample, ifreq], - (channel.weights()[ifreq] / normalization) * - gaussian.data[izen, iazi], - 1e-12, - "gaussian airy per-frequency gaussian match"); - ++isample; - } - } - - assert_stokvec(sum_column_weights(obsel, ifreq), - channel.weights()[ifreq] * peak_weight, - 1e-12, - "gaussian airy per-frequency column normalization"); - } + ARTS_USER_ERROR_IF( + gaussian_obsel.poslos_grid().size() != 4, + "Expected zero-zenith azimuth collapse to leave 4 LOS samples, got {}", + gaussian_obsel.poslos_grid().size()) + assert_stokvec(gaussian_obsel.weight_matrix()[0, 0], + {1.0, 0.0, 0.0, 0.0}, + 1e-12, + "collapsed gaussian zero-zenith ring weight"); + + const sensor::GaussianAiryAntenna airy{ZenGrid{{0.0, 0.2}}, 1.0, 3}; + const auto airy_obsel = + airy(sensor::DiracChannel{100.0e9}, {600e3, 10.0, 20.0}, {45.0, 30.0}); + const Numeric gain = gaussian_airy_expected_gain(0.2, 100.0e9, 1.0); + const Numeric retained_scale = 3.0 * (1.0 + gain) / (1.0 + 3.0 * gain); + + ARTS_USER_ERROR_IF( + airy_obsel.poslos_grid().size() != 4, + "Expected Gaussian Airy zero-zenith azimuth collapse to leave 4 LOS samples, got {}", + airy_obsel.poslos_grid().size()) + assert_stokvec( + airy_obsel.weight_matrix()[0, 0], + {retained_scale / (3.0 * (1.0 + gain)), 0.0, 0.0, 0.0}, + 1e-6, + "collapsed gaussian airy center is renormalized with retained points"); + assert_stokvec( + airy_obsel.weight_matrix()[1, 0], + {retained_scale * gain / (3.0 * (1.0 + gain)), 0.0, 0.0, 0.0}, + 1e-6, + "collapsed gaussian airy off-axis weight is renormalized uniformly"); + assert_stokvec(sum_column_weights(airy_obsel, 0), + {1.0, 0.0, 0.0, 0.0}, + 1e-6, + "collapsed gaussian airy column renormalizes to unity"); } - void test_degenerate_azimuth_ring_is_collapsed() { - const sensor::GaussianAntenna gaussian{ - ZenGrid{{0.0, 1.0}}, AziGrid{{0.0, 120.0, 240.0}}, 1.0, 1.0}; - const auto gaussian_obsel = - gaussian(sensor::DiracChannel{100.0e9}, {600e3, 10.0, 20.0}, {45.0, 30.0}); - - ARTS_USER_ERROR_IF(gaussian_obsel.poslos_grid().size() != 4, - "Expected zero-zenith azimuth collapse to leave 4 LOS samples, got {}", - gaussian_obsel.poslos_grid().size()) - assert_stokvec(gaussian_obsel.weight_matrix()[0, 0], - {3.0, 0.0, 0.0, 0.0}, - 1e-12, - "collapsed gaussian zero-zenith ring weight"); - - const sensor::GaussianAiryAntenna airy{ - ZenGrid{{0.0, 0.2}}, AziGrid{{0.0, 120.0, 240.0}}, 1.0}; - const auto airy_obsel = - airy(sensor::DiracChannel{100.0e9}, {600e3, 10.0, 20.0}, {45.0, 30.0}); - const Numeric gain = gaussian_airy_expected_gain(0.2, 100.0e9, 1.0); - const Numeric retained_scale = 3.0 * (1.0 + gain) / (1.0 + 3.0 * gain); - - ARTS_USER_ERROR_IF(airy_obsel.poslos_grid().size() != 4, - "Expected Gaussian Airy zero-zenith azimuth collapse to leave 4 LOS samples, got {}", - airy_obsel.poslos_grid().size()) - assert_stokvec(airy_obsel.weight_matrix()[0, 0], - {retained_scale / (3.0 * (1.0 + gain)), 0.0, 0.0, 0.0}, - 1e-6, - "collapsed gaussian airy center is renormalized with retained points"); - assert_stokvec(airy_obsel.weight_matrix()[1, 0], - {retained_scale * gain / (3.0 * (1.0 + gain)), 0.0, 0.0, 0.0}, - 1e-6, - "collapsed gaussian airy off-axis weight is renormalized uniformly"); - assert_stokvec(sum_column_weights(airy_obsel, 0), - {1.0, 0.0, 0.0, 0.0}, - 1e-6, - "collapsed gaussian airy column renormalizes to unity"); - } +void test_gaussian_does_not_double_count_degenerate_azimuth() { + const sensor::GaussianAntenna gaussian{ZenGrid{{0.0, 1.0}}, 1.0}; + const auto obsel = + gaussian(sensor::DiracChannel{100.0e9}, {600e3, 10.0, 20.0}, {45.0, 0.0}); + + ARTS_USER_ERROR_IF( + obsel.poslos_grid().size() != 3, + "Expected degenerate azimuth collapse to leave 3 LOS samples, got {}", + obsel.poslos_grid().size()) + + const Numeric wing = std::exp(-0.5); + assert_stokvec(obsel.weight_matrix()[0, 0], + {1.0, 0.0, 0.0, 0.0}, + 1e-12, + "gaussian center weight should be counted once"); + assert_stokvec(obsel.weight_matrix()[1, 0], + {wing, 0.0, 0.0, 0.0}, + 1e-12, + "gaussian wing at azimuth 0"); + assert_stokvec(obsel.weight_matrix()[2, 0], + {wing, 0.0, 0.0, 0.0}, + 1e-12, + "gaussian wing at azimuth 180"); +} void test_gaussian_airy_rejects_nonpositive_frequencies() { - const sensor::GaussianAiryAntenna ant{ZenGrid{{0.0}}, AziGrid{{0.0}}, 1.0}; + const sensor::GaussianAiryAntenna ant{ZenGrid{{0.0}}, 1.0, 1}; bool threw = false; try { @@ -205,16 +182,17 @@ void test_gaussian_airy_rejects_nonpositive_frequencies() { threw = true; } - ARTS_USER_ERROR_IF(not threw, - "Gaussian Airy antenna must reject nonpositive channel frequencies") + ARTS_USER_ERROR_IF( + not threw, + "Gaussian Airy antenna must reject nonpositive channel frequencies") } } // namespace int main() { - test_gaussian_initialization_uses_antenna_frame_offsets(); + test_gaussian_initialization_uses_los_offset(); test_gaussian_airy_is_frequency_dependent(); - test_gaussian_airy_matches_frequency_specific_gaussian_pattern(); test_degenerate_azimuth_ring_is_collapsed(); + test_gaussian_does_not_double_count_degenerate_azimuth(); test_gaussian_airy_rejects_nonpositive_frequencies(); return 0; -} \ No newline at end of file +} diff --git a/src/tests/test_sensor_builder.cc b/src/tests/test_sensor_builder.cc index 34e5965ba4..d2f28b449f 100644 --- a/src/tests/test_sensor_builder.cc +++ b/src/tests/test_sensor_builder.cc @@ -49,10 +49,9 @@ Numeric gaussian_airy_expected_gain(Numeric zenith_deg, } void test_sensor_builder_returns_meta_per_geometry() { - sensor::SensorBuilder builder( - {sensor::BoxChannel{AscendingGrid{100.0, 101.0}}, - sensor::DiracChannel{200.0}}, - std::make_shared()); + sensor::Builder builder({sensor::BoxChannel{AscendingGrid{100.0, 101.0}}, + sensor::DiracChannel{200.0}}, + std::make_shared()); const std::array pos{{{600e3, 10.0, 20.0}, {601e3, 11.0, 21.0}}}; const std::array los{{{20.0, 30.0}, {40.0, 50.0}}}; @@ -71,7 +70,7 @@ void test_sensor_builder_returns_meta_per_geometry() { ARTS_USER_ERROR_IF( gf0.gridname<0>() != "channel" or gf1.gridname<0>() != "channel", - "SensorBuilder meta grid must be the channel axis") + "Builder meta grid must be the channel axis") assert_close(gf0.grid<0>()[0], 0.0, 0.0, "meta[0] channel index 0"); assert_close(gf0.grid<0>()[1], 1.0, 0.0, "meta[0] channel index 1"); @@ -90,8 +89,8 @@ void test_sensor_builder_returns_meta_per_geometry() { } void test_sensor_builder_rejects_mismatched_geometry_counts() { - sensor::SensorBuilder builder({sensor::DiracChannel{}}, - std::make_shared()); + sensor::Builder builder({sensor::DiracChannel{}}, + std::make_shared()); const std::array pos{{{600e3, 10.0, 20.0}}}; const std::array los{{{20.0, 30.0}, {40.0, 50.0}}}; @@ -103,17 +102,15 @@ void test_sensor_builder_rejects_mismatched_geometry_counts() { threw = true; } - ARTS_USER_ERROR_IF( - not threw, - "SensorBuilder must reject mismatching position and LOS counts") + ARTS_USER_ERROR_IF(not threw, + "Builder must reject mismatching position and LOS counts") } void test_sensor_builder_uses_gaussian_airy_frequency_dependence() { const Stokvec peak_weight{2.0, 0.0, 0.0, 0.0}; - sensor::SensorBuilder builder( - {sensor::BoxChannel{AscendingGrid{100.0e9, 200.0e9}}}, - std::make_shared( - ZenGrid{{0.0, 0.2}}, AziGrid{{0.0}}, 1.0, peak_weight)); + sensor::Builder builder({sensor::BoxChannel{AscendingGrid{100.0e9, 200.0e9}}}, + std::make_shared( + ZenGrid{{0.0, 0.2}}, 1.0, 1, peak_weight)); const std::array pos{{{600e3, 10.0, 20.0}}}; const std::array los{{{45.0, 30.0}}}; @@ -162,8 +159,9 @@ void test_unflatten_updates_shared_poslos_grids() { ARTS_USER_ERROR_IF(obsels[0].poslos_grid_ptr() == original_poslos, "unflatten must replace the shared poslos grid") - ARTS_USER_ERROR_IF(not obsels[0].same_poslos(obsels[1]), - "obsels sharing a geometry must still share it after unflatten") + ARTS_USER_ERROR_IF( + not obsels[0].same_poslos(obsels[1]), + "obsels sharing a geometry must still share it after unflatten") assert_close(obsels[0].poslos_grid()[0].los[0], 15.0, 0.0, @@ -184,10 +182,13 @@ void test_unflatten_updates_shared_poslos_grids() { ARTS_USER_ERROR_IF(obsels[0].f_grid_ptr() == original_freq, "unflatten must replace the shared frequency grid") - ARTS_USER_ERROR_IF(not obsels[0].same_freqs(obsels[1]), - "obsels sharing frequencies must still share them after unflatten") - assert_close(obsels[0].f_grid()[0], 101.0, 0.0, "updated frequency for first obsel"); - assert_close(obsels[1].f_grid()[0], 101.0, 0.0, "updated frequency for second obsel"); + ARTS_USER_ERROR_IF( + not obsels[0].same_freqs(obsels[1]), + "obsels sharing frequencies must still share them after unflatten") + assert_close( + obsels[0].f_grid()[0], 101.0, 0.0, "updated frequency for first obsel"); + assert_close( + obsels[1].f_grid()[0], 101.0, 0.0, "updated frequency for second obsel"); } } // namespace diff --git a/src/workspace_agenda_class.cpp b/src/workspace_agenda_class.cpp index 25f9cf0ede..65c6c43e6c 100644 --- a/src/workspace_agenda_class.cpp +++ b/src/workspace_agenda_class.cpp @@ -440,6 +440,22 @@ std::string Agenda::sphinx_list(const std::string_view prep) const { return out; } +void Agenda::change_default(const std::string_view name, Wsv value) { + assert(not name.empty()); + + for (Method& method : methods) { + const auto& method_name = method.get_name(); + if (method_name.starts_with(internal_prefix) and + method_name.ends_with(name)) { + method.change_default(std::move(value)); + return; + } + } + + ARTS_USER_ERROR(std::format( + R"(Agenda "{}" does not have a default input "{}")", this->name, name)); +} + void xml_io_stream::write(std::ostream& os, const Agenda& x, bofstream* pbofs, diff --git a/src/workspace_agenda_class.h b/src/workspace_agenda_class.h index 524cf0daa0..cae05d333d 100644 --- a/src/workspace_agenda_class.h +++ b/src/workspace_agenda_class.h @@ -68,6 +68,8 @@ class Agenda { [[nodiscard]] std::string sphinx_list( const std::string_view prep = "- ") const; + + void change_default(const std::string_view name, Wsv value); }; template <> diff --git a/src/workspace_agenda_creator.cpp b/src/workspace_agenda_creator.cpp index 1d0a8fa855..07b3362747 100644 --- a/src/workspace_agenda_creator.cpp +++ b/src/workspace_agenda_creator.cpp @@ -2,6 +2,7 @@ #include #include +#include #include #include @@ -99,7 +100,8 @@ Agenda get_spectral_rad_observer_agenda(const std::string_view option) { agenda.add("spectral_rad_bkgAgendasAtEndOfPath"); agenda.add("atm_pathFromPath"); agenda.add("freq_grid_pathFromPath"); - agenda.add("spectral_propmat_pathAdaptiveHalfPath"); + agenda.add("spectral_propmat_pathFromPath"); + agenda.add("spectral_propmat_pathAddAdaptiveHalfPath"); agenda.add("spectral_radSetToBackground"); agenda.add("spectral_radSinglePathEmissionFrequencyLoop"); agenda.add("spectral_rad_jacAddSensorJacobianPerturbations"); @@ -166,10 +168,19 @@ Agenda get_measurement_inversion_agenda(const std::string_view option) { using enum measurement_inversion_agendaPredefined; switch (to(option)) { - case Standard: + case LowMemory: + agenda.add("measurement_vec_errorFromModelState"); + agenda.add("jac_targetsConditionalClear"); + agenda.add("measurement_vecFromSensor", SetWsv("kernel", "Low Memory"s)); + agenda.add("measurement_jacTransformations"); + agenda.add("measurement_vecConditionalAddError"); + agenda.add("measurement_vec_fitFromMeasurement"); + break; + case HighPerformance: agenda.add("measurement_vec_errorFromModelState"); agenda.add("jac_targetsConditionalClear"); - agenda.add("measurement_vecFromSensor"); + agenda.add("measurement_vecFromSensor", + SetWsv("kernel", "High Performance"s)); agenda.add("measurement_jacTransformations"); agenda.add("measurement_vecConditionalAddError"); agenda.add("measurement_vec_fitFromMeasurement"); @@ -249,3 +260,68 @@ Agenda get_ray_point_back_propagation_agenda(const std::string_view option) { return std::move(agenda).finalize(true); } + +Agenda get_spectral_propmat_and_atm_path_agenda(const std::string_view option) { + AgendaCreator agenda("spectral_propmat_and_atm_path_agenda"); + + using enum spectral_propmat_and_atm_path_agendaPredefined; + switch (to(option)) { + case Default: + agenda.add("atm_pathFromPath"); + agenda.add("freq_grid_pathFromPath"); + agenda.add("spectral_propmat_pathFromPath"); + agenda.add("Touch", "ray_path"); + break; + case AdaptiveHalfPath: + agenda.add("atm_pathFromPath"); + agenda.add("freq_grid_pathFromPath"); + agenda.add("spectral_propmat_pathFromPath"); + agenda.add("spectral_propmat_pathAddAdaptiveHalfPath"); + break; + case Profile2Path: + agenda.add("Profile2Path"); + agenda.add("Touch", "ray_path"); + break; + } + + return std::move(agenda).finalize(true); +} + +Agenda get_ray_path_observer_agenda(const std::string_view option) { + AgendaCreator agenda("ray_path_observer_agenda"); + + using enum ray_path_observer_agendaPredefined; + switch (to(option)) { + case GeometricProfile: + agenda.add("ray_pathInit", + SetWsv{"pos", "obs_pos"}, + SetWsv{"los", "obs_los"}, + SetWsv("as_sensor", Index{1})); + agenda.add("ray_pathSetGeometricExtremes", + SetWsv("surf_search_accuracy", Numeric{0.1}), + SetWsv("surf_safe_search", Index{1})); + agenda.add("ray_pathAddGeometricAltitudeGridCrossings"); + break; + case GeometricDefault: { + Agenda out; + auto& wsm_def = internal_workspace_methods() + .at("ray_path_observer_agendaSetGeometric") + .gin_value; + ray_path_observer_agendaSetGeometric(out, + wsm_def[0]->get(), + wsm_def[1]->get(), + wsm_def[2]->get(), + wsm_def[3]->get(), + wsm_def[4]->get(), + wsm_def[5]->get(), + wsm_def[6]->get(), + wsm_def[7]->get(), + wsm_def[8]->get(), + wsm_def[9]->get(), + wsm_def[10]->get()); + return out; + } break; + } + + return std::move(agenda).finalize(true); +} diff --git a/src/workspace_agendas.cpp b/src/workspace_agendas.cpp index a251e3e478..ca3947ec60 100644 --- a/src/workspace_agendas.cpp +++ b/src/workspace_agendas.cpp @@ -54,6 +54,40 @@ calculations that are happening deep in your ARTS method calls. }, }; + wsa_data["spectral_propmat_and_atm_path_agenda"] = { + .desc = + R"--(Computes several path parameters along the path. + +The main use of this agenda is to allow adapting the path points +based on spectral parameters. +)--", + .output = {"spectral_propmat_path", + "spectral_nlte_srcvec_path", + "spectral_propmat_jac_path", + "spectral_nlte_srcvec_jac_path", + "freq_grid_path", + "freq_wind_shift_jac_path", + "atm_path", + "ray_path"}, + .input = + {"ray_path", "jac_targets", "freq_grid", "atm_field", "surf_field"}, + .enum_options = {"Default", "AdaptiveHalfPath", "Profile2Path"}, + .enum_default = "Default", + .output_constraints = + { + {"spectral_propmat_path.size() == spectral_nlte_srcvec_path.size() and spectral_propmat_path.size() == spectral_propmat_jac_path.size() and spectral_propmat_path.size() == spectral_nlte_srcvec_jac_path.size() and spectral_propmat_path.size() == freq_grid_path.size() and spectral_propmat_path.size() == freq_wind_shift_jac_path.size() and spectral_propmat_path.size() == atm_path.size() and spectral_propmat_path.size() == ray_path.size()", + "On output, all path arrays have the same size.", + "spectral_propmat_path.size()", + "spectral_nlte_srcvec_path.size()", + "spectral_propmat_jac_path.size()", + "spectral_nlte_srcvec_jac_path.size()", + "freq_grid_path.size()", + "freq_wind_shift_jac_path.size()", + "atm_path.size()", + "ray_path.size()"}, + }, + }; + wsa_data["single_propmat_agenda"] = { .desc = R"--(Computes the propagation matrix, the non-LTE source vector, the dispersion, and their derivatives. @@ -125,7 +159,7 @@ If you do not need single-frequency-point calculations, consider using }; wsa_data["ray_path_observer_agenda"] = { - .desc = R"--(Gets the propagation path as it is observed. + .desc = R"--(Gets the propagation path as it is observed. The intent of this agenda is to provide a propagation path as seen from the observer position and line of sight. @@ -133,8 +167,10 @@ position and line of sight. .. tip:: The perhaps easiest way to set this agenda up is to use the *ray_path_observer_agendaSetGeometric* method. )--", - .output = {"ray_path"}, - .input = {"obs_pos", "obs_los"}, + .output = {"ray_path"}, + .input = {"obs_pos", "obs_los"}, + .enum_options = {"GeometricDefault", "GeometricProfile"}, + .enum_default = "GeometricDefault", }; wsa_data["ray_point_back_propagation_agenda"] = { @@ -341,8 +377,8 @@ it does a lot of unnecessary checks and operations that are not always needed. )--", .output = {"measurement_vec_fit", "measurement_jac"}, .input = {"jac_targets", "do_jac"}, - .enum_options = {"Standard"}, - .enum_default = "Standard", + .enum_options = {"LowMemory", "HighPerformance"}, + .enum_default = "LowMemory", .output_constraints = { {"do_jac != static_cast(measurement_jac.size() == 0)", diff --git a/src/workspace_group_friends.cpp b/src/workspace_group_friends.cpp index 9c9e88724e..7a754e8b98 100644 --- a/src/workspace_group_friends.cpp +++ b/src/workspace_group_friends.cpp @@ -679,8 +679,7 @@ The types are *AscendingGrid* x *LatGrid* x *LonGrid*. The grids are all sorted } for (const auto& [name, g] : internal_workspace_groups()) { - if (auto ptr = wsg_data.find(name); ptr != wsg_data.end()) - wsg_data.erase(ptr); + wsg_data.erase(name); } return wsg_data; } diff --git a/src/workspace_meta_methods.cpp b/src/workspace_meta_methods.cpp index e5df9a12d2..bb10d9ee03 100644 --- a/src/workspace_meta_methods.cpp +++ b/src/workspace_meta_methods.cpp @@ -169,15 +169,13 @@ This method simply is a convenience wrapper for that use case. .author = {"Richard Larsson"}, .methods = {"ray_pointBackground", "spectral_rad_bkgAgendasAtEndOfPath", - "atm_pathFromPath", - "freq_grid_pathFromPath", - "spectral_propmat_pathFromPath", + "spectral_propmat_and_atm_path_agendaExecute", "spectral_tramat_pathFromPath", "spectral_rad_srcvec_pathFromPropmat", "spectral_radStepByStepEmission", "spectral_rad_jacFromBackground", "spectral_rad_jacAddPathPropagation"}, - .out = {"spectral_rad", "spectral_rad_jac"}, + .out = {"spectral_rad", "spectral_rad_jac", "ray_path"}, }); wsm_meta.push_back(WorkspaceMethodInternalMetaRecord{ @@ -186,12 +184,10 @@ This method simply is a convenience wrapper for that use case. .author = {"Richard Larsson"}, .methods = {"ray_pointBackground", "spectral_rad_bkgAgendasAtEndOfPath", - "atm_pathFromPath", - "freq_grid_pathFromPath", - "spectral_propmat_pathFromPath", + "spectral_propmat_and_atm_path_agendaExecute", "spectral_radSetToBackground", "spectral_radSinglePathEmissionFrequencyLoop"}, - .out = {"spectral_rad", "spectral_rad_jac"}, + .out = {"spectral_rad", "spectral_rad_jac", "ray_path"}, }); wsm_meta.push_back(WorkspaceMethodInternalMetaRecord{ @@ -201,9 +197,7 @@ This method simply is a convenience wrapper for that use case. .author = {"Richard Larsson"}, .methods = {"ray_pointBackground", "spectral_rad_bkgAgendasAtEndOfPath", - "atm_pathFromPath", - "freq_grid_pathFromPath", - "spectral_propmat_pathFromPath", + "spectral_propmat_and_atm_path_agendaExecute", "spectral_propmat_scat_pathFromPath", "spectral_propmat_pathAddScattering", "spectral_tramat_pathFromPath", @@ -213,7 +207,7 @@ This method simply is a convenience wrapper for that use case. "spectral_radStepByStepEmission", "spectral_rad_jacFromBackground", "spectral_rad_jacAddPathPropagation"}, - .out = {"spectral_rad", "spectral_rad_jac"}, + .out = {"spectral_rad", "spectral_rad_jac", "ray_path"}, }); wsm_meta.push_back(WorkspaceMethodInternalMetaRecord{ @@ -222,14 +216,12 @@ This method simply is a convenience wrapper for that use case. .author = {"Richard Larsson"}, .methods = {"ray_pointBackground", "spectral_rad_bkgAgendasAtEndOfPath", - "atm_pathFromPath", - "freq_grid_pathFromPath", - "spectral_propmat_pathFromPath", + "spectral_propmat_and_atm_path_agendaExecute", "spectral_tramat_pathFromPath", "spectral_radCumulativeTransmission", "spectral_rad_jacFromBackground", "spectral_rad_jacAddPathPropagation"}, - .out = {"spectral_rad", "spectral_rad_jac"}, + .out = {"spectral_rad", "spectral_rad_jac", "ray_path"}, }); wsm_meta.push_back(WorkspaceMethodInternalMetaRecord{ @@ -247,14 +239,12 @@ This method simply is a convenience wrapper for that use case. .desc = "Computes clearsky transmission of spectral radiances", .author = {"Richard Larsson"}, .methods = {"ray_pointBackground", - "atm_pathFromPath", - "freq_grid_pathFromPath", - "spectral_propmat_pathFromPath", + "spectral_propmat_and_atm_path_agendaExecute", "spectral_tramat_pathFromPath", "spectral_radCumulativeTransmission", "spectral_rad_jacFromBackground", "spectral_rad_jacAddPathPropagation"}, - .out = {"spectral_rad", "spectral_rad_jac"}, + .out = {"spectral_rad", "spectral_rad_jac", "ray_path"}, }); wsm_meta.push_back(WorkspaceMethodInternalMetaRecord{ diff --git a/src/workspace_method_class.cpp b/src/workspace_method_class.cpp index a825b8f2e6..a59983aa59 100644 --- a/src/workspace_method_class.cpp +++ b/src/workspace_method_class.cpp @@ -180,12 +180,14 @@ void Method::operator()(Workspace& ws) const try { Method outputs: {:B,} Method inputs: {:B,} -(NOTE: "_" and "@" prefixes are used for default and user inputs, respectively) +(NOTE: "{}" and "{}" prefixes are used for default and user inputs, respectively) {})", name, inargs, outargs, + internal_prefix, + named_input_prefix, std::string_view(e.what()))); } @@ -292,6 +294,11 @@ std::string Method::sphinx_list_item() const { return std::format("*{}*{}", name, setvals); } +void Method::change_default(Wsv value) { + assert(setval.has_value()); + setval = std::move(value); +} + void xml_io_stream::write(std::ostream& os, const Wsv& x, bofstream* pbofs, diff --git a/src/workspace_method_class.h b/src/workspace_method_class.h index 26d8839a00..db969eaec4 100644 --- a/src/workspace_method_class.h +++ b/src/workspace_method_class.h @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -46,6 +47,8 @@ class Method { [[nodiscard]] std::string sphinx_list_item() const; [[nodiscard]] bool is_callback() const; + + void change_default(Wsv value); }; template <> @@ -88,6 +91,10 @@ struct std::formatter { template FmtContext::iterator format(const Method& v, FmtContext& ctx) const { + if (tags.short_str) { + return tags.format(ctx, v.get_name()); + } + tags.add_if_bracket(ctx, '['); tags.add_if_bracket(ctx, '\n'); @@ -139,28 +146,21 @@ struct std::formatter { template FmtContext::iterator format(const Agenda& v, FmtContext& ctx) const { if (tags.short_str) { - tags.format(ctx, v.get_name(), '\n'); - - std::vector setvals; - for (auto& method : v.get_methods()) { - if (method.get_setval().has_value()) { - auto x = method.get_setval()->vformat("({})"sv); - x.replace(x.find('\n'), 1, " "sv); - if (x.size() > 50) x = x.substr(0, 47) + "..."; - setvals.push_back( - std::format("{} = {}", method.get_name().substr(1), x)); - } else { - tags.format(ctx, - " "sv, - method.get_name(), - '(', - std::format("{:,}", setvals), - ")\n"sv); - setvals.clear(); + tags.format(ctx, v.get_name(), ':', '\n'); + + for (const auto& method : v.get_methods()) { + const auto& method_name = method.get_name(); + if (method_name.starts_with(named_input_prefix) or + method_name.starts_with(internal_prefix)) { + continue; } + + tags.format(ctx, " - "sv, method.get_name(), '\n'); } - return tags.format(ctx, std::format("{:n}", setvals)); + if (v.get_methods().empty()) tags.format(ctx, " (no methods)", '\n'); + + return ctx.out(); } tags.add_if_bracket(ctx, '['); diff --git a/src/workspace_method_extra_doc.cpp b/src/workspace_method_extra_doc.cpp index 6c4db0fbfb..8dd8e6bb9e 100644 --- a/src/workspace_method_extra_doc.cpp +++ b/src/workspace_method_extra_doc.cpp @@ -158,7 +158,7 @@ Below follows a complete list of all single species tags that can be set using t x.FullName(), x.is_predefined() ? "Predefined Model"sv : x.is_normal() ? "Normal Isotopologue"sv - : "Joker"sv, + : "All Isotopologues"sv, x.is_normal() ? std::format("{}", x.mass) : ""s, x.is_normal() ? std::format("{}", x.gi) : ""s); } diff --git a/src/workspace_methods.cpp b/src/workspace_methods.cpp index 33ff30e29b..e2d37c886f 100644 --- a/src/workspace_methods.cpp +++ b/src/workspace_methods.cpp @@ -1145,18 +1145,6 @@ in your own control-flow. .in = {"ray_path", "atm_field"}, }; - wsm_data["freq_grid_pathFromPath"] = { - .desc = R"--(Gets the frequency grids along the path. - -The derivative transformation for wind parameters is also returned. - -See *spectral_propmat_jacWindFix* for use of the wind shift data. -)--", - .author = {"Richard Larsson"}, - .out = {"freq_grid_path", "freq_wind_shift_jac_path"}, - .in = {"freq_grid", "ray_path", "atm_path"}, - }; - wsm_data["spectral_propmat_pathFromPath"] = { .desc = R"--(Gets the propagation matrix and non-LTE source term along the path. @@ -1179,7 +1167,7 @@ Also outputs the *freq_grid_path* as a side effect (of wind). .pass_workspace = true, }; - wsm_data["spectral_propmat_pathAdaptiveHalfPath"] = { + wsm_data["spectral_propmat_pathAddAdaptiveHalfPath"] = { .desc = R"--(Same as *spectral_propmat_pathFromPath* but with adaptive path. @@ -1202,7 +1190,11 @@ Note that the input "freq_wind_shift_jac_path", "ray_path", "atm_path"}, - .in = {"freq_grid_path", + .in = {"spectral_propmat_path", + "spectral_nlte_srcvec_path", + "spectral_propmat_jac_path", + "spectral_nlte_srcvec_jac_path", + "freq_grid_path", "freq_wind_shift_jac_path", "ray_path", "atm_path", @@ -1221,6 +1213,34 @@ Note that the input .pass_workspace = true, }; + wsm_data["spectral_propmat_and_atm_path_agendaSetAdaptive"] = { + .desc = + R"--(Sets the *spectral_propmat_and_atm_path_agenda* to adaptive mode with the provided parameters. +)--", + .author = {"Richard Larsson"}, + .out = {"spectral_propmat_and_atm_path_agenda"}, + // Reuse parameters of the pathFromPath method + .gin = wsm_data.at("spectral_propmat_pathAddAdaptiveHalfPath").gin, + .gin_type = + wsm_data.at("spectral_propmat_pathAddAdaptiveHalfPath").gin_type, + .gin_value = + wsm_data.at("spectral_propmat_pathAddAdaptiveHalfPath").gin_value, + .gin_desc = + wsm_data.at("spectral_propmat_pathAddAdaptiveHalfPath").gin_desc, + }; + + wsm_data["freq_grid_pathFromPath"] = { + .desc = R"--(Gets the frequency grids along the path. + +The derivative transformation for wind parameters is also returned. + +See *spectral_propmat_jacWindFix* for use of the wind shift data. +)--", + .author = {"Richard Larsson"}, + .out = {"freq_grid_path", "freq_wind_shift_jac_path"}, + .in = {"freq_grid", "ray_path", "atm_path"}, + }; + wsm_data["spectral_propmat_path_species_splitFromPath"] = { .desc = R"--(As *spectral_propmat_pathFromPath* but the output is split between the species in the @@ -3631,6 +3651,19 @@ Points are added where the ray path crosses any of the three grids in pure geome {"The atmospheric field key for which the grid is expected if adding grid crossings is desired"}, }; + wsm_data["ray_pathAddGeometricAltitudeGridCrossings"] = { + .desc = + R"--(Fill the path with with points that crosses the *alt_grid* + +Also checks that only crossings of the *alt_grid* are in *ray_path* upon exit. + +The intent is to use this only with methods that work on altitude profiles. +)--", + .author = {"Richard Larsson"}, + .out = {"ray_path"}, + .in = {"ray_path", "alt_grid", "surf_field"}, + }; + wsm_data["ray_pathFillGeometricHalfStep"] = { .desc = R"--(Fill the path with geometric step points. @@ -4026,8 +4059,32 @@ the first 5 dimensions are computed in parallel. .pass_workspace = true, }; + wsm_data["atm_fieldWindIncludePlanetRotation"] = { + .desc = R"(Add planetary rotation to the wind field u-component. + +This is added as + +.. math:: + \delta u = \frac{2 \pi}{R} \left(\frac{a}{\sqrt{1 - e^2 \sin^2\theta}} + z\right) \cos\theta, + +where +:math:`R` is the rotation period of the planet, +:math:`\theta` is the latitude, +:math:`a` is the equatorial radius of the planet, +:math:`e` is the eccentricity of the planet, and +:math:`z` is the altitude. +)", + .author = {"Patrick Eriksson", "Richard Larsson"}, + .out = {"atm_field"}, + .in = {"atm_field", "surf_field"}, + .gin = {"planet_rotation_period"}, + .gin_type = {"Numeric"}, + .gin_value = {std::nullopt}, + .gin_desc = {"The rotation period of the planet in seconds"}, + }; + wsm_data["atm_fieldAppendBaseData"] = { - .desc = R"--(Append base data to the atmospheric field + .desc = R"--(Append base data to the atmospheric field This will look at the valid ``basename`` for files matching base data. The base data file names are of the form @@ -4053,28 +4110,35 @@ exists in the atmospheric field. The ``allow_missing_pressure`` and ``allow_missing_temperature`` are used to determine if the method should throw if the pressure or temperature is missing. )--", - .author = {"Richard Larsson"}, - .out = {"atm_field"}, - .in = {"atm_field"}, - .gin = {"basename", - "extrapolation", - "deal_with_field_component", - "replace_existing", - "allow_missing_pressure", - "allow_missing_temperature"}, - .gin_type = {"String", "String", "String", "Index", "Index", "Index"}, - .gin_value = {std::nullopt, - String{"Linear"}, - String{"Throw"}, - Index{1}, - Index{0}, - Index{0}}, - .gin_desc = {"The base name of the files", - "The extrapolation to use.", - "How to deal with the field component.", - "Whether or not to replace existing data", - "Whether or not to allow missing pressure data", - "Whether or not to allow missing temperature data"}, + .author = {"Richard Larsson"}, + .out = {"atm_field"}, + .in = {"atm_field"}, + .gin = {"basename", + "extrapolation", + "deal_with_field_component", + "replace_existing", + "allow_missing_pressure", + "allow_missing_temperature", + "pressure_log_interp"}, + .gin_type = + {"String", "String", "String", "Index", "Index", "Index", "Index"}, + .gin_value = + { + std::nullopt, + String{"Linear"}, + String{"Throw"}, + Index{1}, + Index{0}, + Index{0}, + Index{1}, + }, + .gin_desc = {"The base name of the files", + "The extrapolation to use.", + "How to deal with the field component.", + "Whether or not to replace existing data", + "Whether or not to allow missing pressure data", + "Whether or not to allow missing temperature data", + "Whether or not to use log interpolation for pressure data"}, }; wsm_data["atm_fieldAppendLineSpeciesData"] = { @@ -4944,7 +5008,7 @@ Hence, a temperature of 0 means 0s the edges of the *freq_grid*. "freq_grid", "atm_field", "surf_field", - "spectral_propmat_agenda", + "spectral_propmat_and_atm_path_agenda", "rte_option"}, .gin = {"depolarization_factor", "hse_derivative"}, .gin_type = {"Numeric", "Index"}, @@ -5975,6 +6039,54 @@ Additional work is requires if proper coverage of the limb is required .in = {"subsurf_field", "ray_path"}, }; + wsm_data["Profile2Path"] = { + .desc = + R"(Turns profile data into path data. + +Only valid when local spectral properties are identical in all directions. +Neither polarization nor wind calculations are possible with this method. +)", + .author = {"Richard Larsson"}, + .out = {"spectral_propmat_path", + "spectral_nlte_srcvec_path", + "spectral_propmat_jac_path", + "spectral_nlte_srcvec_jac_path", + "atm_path", + "freq_grid_path", + "freq_wind_shift_jac_path"}, + .in = {"freq_grid", + "alt_grid", + "ray_path", + "atm_profile", + "spectral_propmat_profile", + "spectral_propmat_jac_profile", + "spectral_nlte_srcvec_profile", + "spectral_nlte_srcvec_jac_profile"}, + }; + + wsm_data["ProfileFromAltitude"] = { + .desc = + R"(Turns profile data into path data. + +Only valid when local spectral properties are identical in all directions. +Neither polarization nor wind calculations are possible with this method. +)", + .author = {"Richard Larsson"}, + .out = {"atm_profile", + "spectral_propmat_profile", + "spectral_propmat_jac_profile", + "spectral_nlte_srcvec_profile", + "spectral_nlte_srcvec_jac_profile"}, + .in = {"spectral_propmat_agenda", + "jac_targets", + "alt_grid", + "atm_field", + "freq_grid", + "lat", + "lon"}, + .pass_workspace = true, + }; + /* MANUAL ENTRIES BELOW THIS POINT MESSES WITH THE LOGIC OF WSM_DATA */ @@ -5994,9 +6106,6 @@ Additional work is requires if proper coverage of the limb is required std::vector gin_desc{ "Whether or not to use the lookup table instead of pure line-by-line calculations"}; std::vector> gin_value{Index{0}}; - std::vector gout{}; - std::vector gout_type{}; - std::vector gout_desc{}; for (auto& method : {"spectral_propmatInit", "spectral_propmatAddCIA", @@ -6010,28 +6119,14 @@ Additional work is requires if proper coverage of the limb is required for (Size i = 0; i < x.gin.size(); i++) { if (auto ptr = stdr::find(gin, x.gin[i]); ptr != gin.end()) { - gin_desc[ptr - gin.begin()] += - std::format(", *{}*", std::string_view{method}); + gin_desc[ptr - gin.begin()] += std::format(", *{}*", method); } else { gin.emplace_back(x.gin[i]); gin_type.emplace_back(x.gin_type[i]); - gin_desc.emplace_back( - std::format("See *{}*", std::string_view{method})); + gin_desc.emplace_back(std::format("See *{}*", method)); gin_value.emplace_back(x.gin_value[i]); } } - - for (Size i = 0; i < x.gout.size(); i++) { - if (auto ptr = stdr::find(gout, x.gout[i]); ptr != gout.end()) { - gout_desc[ptr - gout.begin()] += - std::format(", *{}*", std::string_view{method}); - } else { - gout.emplace_back(x.gout[i]); - gout_type.emplace_back(x.gout_type[i]); - gout_desc.emplace_back( - std::format("See *{}*", std::string_view{method})); - } - } } catch (std::out_of_range&) { throw std::runtime_error(std::format( R"(Missing method: "{}", @@ -6044,7 +6139,7 @@ cannot generate automatic method signature for "{}" wsm_data[meta] = { .desc = - R"--(Sets the *spectral_propmat_agenda* automatically from absorption data and species tag meta information. + R"--(Sets the *spectral_propmat_agenda* automatically from absorption data on the workspace, and species tag meta information. The following methods are considered for addition to the agenda: @@ -6061,16 +6156,13 @@ Note that the signature of this method changes depending on the input methods. because several generic input parameters are used in the methods. Please see the individual methods for more information. )--", - .author = {"Richard Larsson"}, - .out = {"spectral_propmat_agenda"}, - .gout = gout, - .gout_type = gout_type, - .gout_desc = gout_desc, - .in = {"abs_species", "abs_bands"}, - .gin = gin, - .gin_type = gin_type, - .gin_value = gin_value, - .gin_desc = gin_desc, + .author = {"Richard Larsson"}, + .out = {"spectral_propmat_agenda"}, + .gin = gin, + .gin_type = gin_type, + .gin_value = gin_value, + .gin_desc = gin_desc, + .pass_workspace = true, }; } diff --git a/src/workspace_variables.cpp b/src/workspace_variables.cpp index 477f3291d3..5106d6da2f 100644 --- a/src/workspace_variables.cpp +++ b/src/workspace_variables.cpp @@ -372,12 +372,32 @@ The units are *spectral_rad_jac* per meter. .type = "ArrayOfStokvecMatrix", }; + wsv_data["spectral_nlte_srcvec_jac_profile"] = { + .desc = R"--(Additional non-LTE derivative in a propagation profile + +.. note:: + Polarization is considered but if the profile is polarized, + using it in anyways will yield invalid results. +)--", + .type = "ArrayOfStokvecMatrix", + }; + wsv_data["spectral_nlte_srcvec_path"] = { .desc = R"--(Additional non-LTE along the propagation path )--", .type = "ArrayOfStokvecVector", }; + wsv_data["spectral_nlte_srcvec_profile"] = { + .desc = R"--(Additional non-LTE in a propagation profile + +.. note:: + Polarization is considered but if the profile is polarized, + using it in anyways will yield invalid results. +)--", + .type = "ArrayOfStokvecVector", + }; + wsv_data["spectral_phamat_spectral"] = { .desc = R"--(The spectral phase matrix of totally random orientation particles at a single point along a path using spectral representation @@ -425,12 +445,32 @@ The units depend on what is set in *jac_targets* [1 / m / jacobian target's unit .type = "ArrayOfPropmatMatrix", }; + wsv_data["spectral_propmat_jac_profile"] = { + .desc = R"--(Propagation derivative matrices in a propagation profile + +.. note:: + Polarization is considered but if the profile is polarized, + using it in anyways will yield invalid results. +)--", + .type = "ArrayOfPropmatMatrix", + }; + wsv_data["spectral_propmat_path"] = { .desc = R"--(Propagation matrices along the propagation path )--", .type = "ArrayOfPropmatVector", }; + wsv_data["spectral_propmat_profile"] = { + .desc = R"--(Propagation matrices in a propagation profile + +.. note:: + Polarization is considered but if the profile is polarized, + using it in anyways will yield invalid results. +)--", + .type = "ArrayOfPropmatVector", + }; + wsv_data["spectral_propmat_scat"] = { .desc = R"--(This contains the propagation matrix for scattering for the current path point. diff --git a/tests/core/jac/full_arts_emission.py b/tests/core/jac/full_arts_emission.py index a5b1dbb2d8..dcd4b97650 100644 --- a/tests/core/jac/full_arts_emission.py +++ b/tests/core/jac/full_arts_emission.py @@ -76,4 +76,4 @@ dx2 = np.array(dx2) # Does this large discrepancy mean we need better code? -assert np.allclose(dx1 / dx2 - 1, 0, atol=0.02), "Should be within 2%" +assert np.allclose(dx1 / dx2 - 1, 0, atol=0.05), "Should be within 5%" diff --git a/tests/core/sensor/30cm-antenna.py b/tests/core/sensor/30cm-antenna.py index 72d626664a..ff5b4f0ee2 100644 --- a/tests/core/sensor/30cm-antenna.py +++ b/tests/core/sensor/30cm-antenna.py @@ -9,9 +9,9 @@ zen_grid = np.geomspace(1, 3, NZA) - 1 azi_grid = np.linspace(0, 360, NAZ + 1)[:-1] -ant1 = pa.arts.SensorGaussianAiryAntenna(zen_grid, azi_grid, 30e-2, "Iv") -ant2 = pa.arts.SensorGaussianAiryAntenna(zen_grid + 1e-3, azi_grid, 30e-2, "I") -ch = pa.arts.SensorBoxChannel(1e9, 1e12, 1001) +ant1 = pa.arts.sensor.GaussianAiryAntenna(zen_grid, 30e-2, NAZ, "Iv") +ant2 = pa.arts.sensor.GaussianAiryAntenna(zen_grid + 1e-3, 30e-2, NAZ, "I") +ch = pa.arts.sensor.BoxChannel(1e9, 1e12, 1001) obsel1 = ant1(ch, POS, LOS) obsel2 = ant2(ch, POS, LOS) @@ -30,7 +30,7 @@ assert len(poslos2) == NZA * NAZ, \ "Should have one poslos per non-zero DZA, but only one for zero DZA - has no zero DZA" -ant = pa.arts.SensorGaussianAiryAntenna(zen_grid, azi_grid, 30e-2, "I") +ant = pa.arts.sensor.GaussianAiryAntenna(zen_grid, 30e-2, NAZ, "I") pa.plots.plot(obsel1, point_spread=True) pa.plots.plot(obsel2, point_spread=True, frame='local', type='contour', levels=10) diff --git a/tests/core/sensor/600ghz-upper-sideband-sensor.py b/tests/core/sensor/600ghz-upper-sideband-sensor.py index 2e3f6fd2c8..12296f6bd8 100644 --- a/tests/core/sensor/600ghz-upper-sideband-sensor.py +++ b/tests/core/sensor/600ghz-upper-sideband-sensor.py @@ -18,27 +18,27 @@ if_offsets = pyarts.arts.AscendingGrid(np.linspace(IF_LOW, IF_HIGH, NCHANNELS)) print(f"if_offsets: {time.time() - start_time:.2f} seconds") -spectrometer = pyarts.arts.SensorSpectrometer( - pyarts.arts.SensorDiracChannel(), +spectrometer = pyarts.arts.sensor.Spectrometer( + pyarts.arts.sensor.DiracChannel(), if_offsets, ) print(f"spectrometer: {time.time() - start_time:.2f} seconds") -backend = pyarts.arts.SensorHeterodyneFrequencyRange( +backend = pyarts.arts.sensor.HeterodyneFrequencyRange( LO, np.array([LO + IF_LOW, LO + IF_HIGH]), ) print(f"backend: {time.time() - start_time:.2f} seconds") -antenna = pyarts.arts.SensorGaussianAiryAntenna( +antenna = pyarts.arts.sensor.GaussianAiryAntenna( ZEN, - AZI, APERTURE_DIAMETER, + 8, "I", ) print(f"antenna: {time.time() - start_time:.2f} seconds") -sensor = pyarts.arts.SensorBuilder(antenna, spectrometer, backend) +sensor = pyarts.arts.sensor.Builder(antenna, spectrometer, backend) print(f"sensor: {time.time() - start_time:.2f} seconds") ws = pyarts.Workspace() diff --git a/tests/core/sensor/heterodyne_frequency_response.py b/tests/core/sensor/heterodyne_frequency_response.py index a1c35b9e36..9a3400eea7 100644 --- a/tests/core/sensor/heterodyne_frequency_response.py +++ b/tests/core/sensor/heterodyne_frequency_response.py @@ -75,7 +75,7 @@ def make_asymmetric_sideband_filter(): def test_lowpass_then_lo_then_box_spectrometer(): - selector = pyarts.arts.SensorHeterodyneFrequencyRange() + selector = pyarts.arts.sensor.HeterodyneFrequencyRange() print("test_lowpass_then_lo_then_box_spectrometer") selector.lowpass(15.0) print("Apply lowpass keeping 0-15:", selector) @@ -90,7 +90,7 @@ def test_lowpass_then_lo_then_box_spectrometer(): expected_affine=[[16.0, -1.0]], ) - channel = pyarts.arts.SensorBoxChannel(2.6, 4.6, 5) + channel = pyarts.arts.sensor.BoxChannel(2.6, 4.6, 5) response = selector.channel_response(channel) assert_close( @@ -106,7 +106,7 @@ def test_lowpass_then_lo_then_box_spectrometer(): def test_default_positive_range_survives_mixes(): - selector = pyarts.arts.SensorHeterodyneFrequencyRange() + selector = pyarts.arts.sensor.HeterodyneFrequencyRange() selector.mix(10.0) assert_close( @@ -133,7 +133,7 @@ def test_default_positive_range_survives_mixes(): [[0.0, np.inf], [0.0, 1.0], [0.0, 9.0], [0.0, 1.0]], ) - response = selector.channel_response(pyarts.arts.SensorDiracChannel(0.25)) + response = selector.channel_response(pyarts.arts.sensor.DiracChannel(0.25)) assert_close( "default positive range dirac points after two LOs", @@ -148,7 +148,7 @@ def test_default_positive_range_survives_mixes(): def test_bandpass_lo_bandpass_lo_chain(): - selector = pyarts.arts.SensorHeterodyneFrequencyRange() + selector = pyarts.arts.sensor.HeterodyneFrequencyRange() selector.bandpass(np.array([5.0, 15.0])) selector.mix(3.0) selector.bandpass(np.array([3.0, 7.0])) @@ -162,7 +162,7 @@ def test_bandpass_lo_bandpass_lo_chain(): expected_affine=[[9.0, 1.0], [9.0, -1.0]], ) - channel = pyarts.arts.SensorBoxChannel(0.0, 1.0, 3) + channel = pyarts.arts.sensor.BoxChannel(0.0, 1.0, 3) response = selector.channel_response(channel) assert_close( @@ -178,7 +178,7 @@ def test_bandpass_lo_bandpass_lo_chain(): def test_ideal_split_about_lo(): - selector = pyarts.arts.SensorHeterodyneFrequencyRange() + selector = pyarts.arts.sensor.HeterodyneFrequencyRange() selector.bandpass(np.array([5.0, 15.0])) selector.mix(10.0) @@ -203,7 +203,7 @@ def test_ideal_split_about_lo(): def test_weighted_split_about_lo(): - selector = pyarts.arts.SensorHeterodyneFrequencyRange() + selector = pyarts.arts.sensor.HeterodyneFrequencyRange() selector.filter(make_triangular_filter()) selector.mix(10.0) @@ -226,7 +226,7 @@ def test_weighted_split_about_lo(): np.array([0.6]), ) - channel = pyarts.arts.SensorDiracChannel(2.0) + channel = pyarts.arts.sensor.DiracChannel(2.0) response = selector.channel_response(channel) assert_close( @@ -242,11 +242,11 @@ def test_weighted_split_about_lo(): def test_zero_frequency_is_not_double_counted(): - selector = pyarts.arts.SensorHeterodyneFrequencyRange() + selector = pyarts.arts.sensor.HeterodyneFrequencyRange() selector.bandpass(np.array([5.0, 15.0])) selector.mix(10.0) - response = selector.channel_response(pyarts.arts.SensorDiracChannel(0.0)) + response = selector.channel_response(pyarts.arts.sensor.DiracChannel(0.0)) assert_close( "zero-frequency channel point", @@ -261,7 +261,7 @@ def test_zero_frequency_is_not_double_counted(): def test_overlapping_sidebands_with_asymmetric_bandpass(): - selector = pyarts.arts.SensorHeterodyneFrequencyRange() + selector = pyarts.arts.sensor.HeterodyneFrequencyRange() selector.filter(make_asymmetric_sideband_filter()) selector.mix(10.0) @@ -289,9 +289,9 @@ def test_overlapping_sidebands_with_asymmetric_bandpass(): responses = selector.channel_responses( [ - pyarts.arts.SensorDiracChannel(1.0), - pyarts.arts.SensorDiracChannel(2.0), - pyarts.arts.SensorDiracChannel(3.0), + pyarts.arts.sensor.DiracChannel(1.0), + pyarts.arts.sensor.DiracChannel(2.0), + pyarts.arts.sensor.DiracChannel(3.0), ] ) @@ -327,4 +327,4 @@ def test_overlapping_sidebands_with_asymmetric_bandpass(): test_zero_frequency_is_not_double_counted() test_overlapping_sidebands_with_asymmetric_bandpass() -print("\nAll heterodyne frequency-response checks passed.") \ No newline at end of file +print("\nAll heterodyne frequency-response checks passed.") diff --git a/tests/core/sensor/sensor_builder.py b/tests/core/sensor/sensor_builder.py index d03804559c..bb16211953 100644 --- a/tests/core/sensor/sensor_builder.py +++ b/tests/core/sensor/sensor_builder.py @@ -11,12 +11,12 @@ def assert_close(name, got, expected, atol=1e-12): def test_sensor_builder_from_sequences(): - builder = pyarts.arts.SensorBuilder( + builder = pyarts.arts.sensor.Builder( [ - pyarts.arts.SensorBoxChannel(pyarts.arts.AscendingGrid([100.0, 101.0])), - pyarts.arts.SensorDiracChannel(200.0), + pyarts.arts.sensor.BoxChannel(pyarts.arts.AscendingGrid([100.0, 101.0])), + pyarts.arts.sensor.DiracChannel(200.0), ], - pyarts.arts.SensorPencilBeamAntenna(), + pyarts.arts.sensor.PencilBeamAntenna(), ) pos = [ @@ -47,9 +47,9 @@ def test_sensor_builder_from_sequences(): def test_sensor_builder_from_poslos_vector_and_length_guard(): - builder = pyarts.arts.SensorBuilder( - [pyarts.arts.SensorDiracChannel()], - pyarts.arts.SensorPencilBeamAntenna(), + builder = pyarts.arts.sensor.Builder( + [pyarts.arts.sensor.DiracChannel()], + pyarts.arts.sensor.PencilBeamAntenna(), ) poslos = pyarts.arts.SensorPosLosVector() @@ -77,8 +77,9 @@ def test_sensor_builder_from_poslos_vector_and_length_guard(): except ValueError: pass else: - raise AssertionError("SensorBuilder must reject mismatching position/LOS lengths") + raise AssertionError( + "Sensor builder must reject mismatching position/LOS lengths") test_sensor_builder_from_sequences() -test_sensor_builder_from_poslos_vector_and_length_guard() \ No newline at end of file +test_sensor_builder_from_poslos_vector_and_length_guard() diff --git a/tests/python/math/interp.py b/tests/python/math/interp.py index c217eb604d..0854e3e8c4 100644 --- a/tests/python/math/interp.py +++ b/tests/python/math/interp.py @@ -1,24 +1,30 @@ +import os + import matplotlib.pyplot as plt import numpy as np import pyarts3 as pyarts x = np.linspace(-160, 160, 1001) y = np.sin(np.deg2rad(x)) -xn = np.concatenate((np.linspace(-180, 180, 1000) - 360, np.linspace(-180, 180, 1000), np.linspace(-180, 180, 1000) + 360)) +xn = np.concatenate((np.linspace(-180, 180, 1000) - 360, + np.linspace(-180, 180, 1000), np.linspace(-180, 180, 1000) + 360)) # The test fails for N == 0 and N == 4 but the results look OK for N in range(1, 4): - lc = pyarts.arts.interp.ArrayOfLagrangeCyclic(x, xn, N, 0.0) - ll = pyarts.arts.interp.ArrayOfLagrange(x, xn, N, 0.0) - yc = pyarts.math.reinterp(y, lc) - yl = pyarts.math.reinterp(y, ll) + lc = pyarts.arts.interp.ArrayOfLagrangeCyclic(x, xn, N, 0.0) + ll = pyarts.arts.interp.ArrayOfLagrange(x, xn, N, 0.0) + yc = pyarts.math.reinterp(y, lc) + yl = pyarts.math.reinterp(y, ll) + + plt.plot(x, y, lw=2) + plt.plot(xn, yc, "--") + plt.plot(xn, yl, ":") + plt.ylim(-2, 2) - plt.plot(x, y, lw = 2) - plt.plot(xn, yc, "--") - plt.plot(xn, yl, ":") - plt.ylim(-2, 2) + assert np.allclose(yc[:1000], yc[-1000:]), \ + "cyclic grids should mostly overlap" + assert np.allclose(yl[:1000], -yl[-1000:][::-1]), \ + "linear grids should mostly negate" - assert np.allclose(yc[:1000], yc[-1000:]), \ - "cyclic grids should mostly overlap" - assert np.allclose(yl[:1000], -yl[-1000:][::-1]), \ - "linear grids should mostly negate" +if "ARTS_HEADLESS" not in os.environ: + plt.show()