64 const Index& stokes_dim,
70 "No need to use this method with *iy_unit* = \"1\".");
73 "The spectrum matrix *iy* is required to have original radiance\n"
74 "unit, but this seems not to be the case. This as a value above\n"
75 "1e-3 is found in *iy*.")
79 for (
Index is = 0; is < stokes_dim; is++) {
85 for (
Index i = 0; i < iy_aux_vars.
nelem(); i++) {
86 if (iy_aux_vars[i] ==
"iy" || iy_aux_vars[i] ==
"Error" ||
87 iy_aux_vars[i] ==
"Error (uncorrelated)") {
98 const Index& atmfields_checked,
99 const Index& atmgeom_checked,
102 const Index& cloudbox_on,
103 const Index& cloudbox_checked,
104 const Index& scat_data_checked,
111 const Agenda& iy_main_agenda,
116 "The atmospheric fields must be flagged to have\n"
117 "passed a consistency check (atmfields_checked=1).");
119 "The atmospheric geometry must be flagged to have\n"
120 "passed a consistency check (atmgeom_checked=1).");
122 "The cloudbox must be flagged to have\n"
123 "passed a consistency check (cloudbox_checked=1).");
126 "The scattering data must be flagged to have\n"
127 "passed a consistency check (scat_data_checked=1).");
130 Tensor3 iy_transmittance(0, 0, 0);
155 "One or several NaNs found in *iy*.");
176 const Index& stokes_dim,
178 const Index& atmosphere_dim,
190 const Index& cloudbox_on,
198 const Index& jacobian_do,
200 const Agenda& propmat_clearsky_agenda,
201 const Agenda& water_p_eq_agenda,
202 const String& rt_integration_option,
203 const Agenda& iy_main_agenda,
204 const Agenda& iy_space_agenda,
205 const Agenda& iy_surface_agenda,
206 const Agenda& iy_cloudbox_agenda,
207 const Index& iy_agenda_call1,
208 const Tensor3& iy_transmittance,
212 const Tensor3& surface_props_data,
216 const Index& t_interp_order,
257 propmat_clearsky_agenda,
259 rt_integration_option,
272 const Index j_analytical_do = jacobian_do ? do_analytical_jacobian<1>(jacobian_quantities) : 0;
278 const Index nq = j_analytical_do ? jacobian_quantities.
nelem() : 0;
284 if (atmosphere_dim != 1)
286 "With cloudbox on, this method handles only 1D calculations.");
287 if (Naa < 3)
throw runtime_error(
"Naa must be > 2.");
289 if (dpnd_field_dx.
nelem() != jacobian_quantities.
nelem())
291 "*dpnd_field_dx* not properly initialized:\n"
292 "Number of elements in dpnd_field_dx must be equal number of jacobian"
293 " quantities.\n(Note: jacobians have to be defined BEFORE *pnd_field*"
294 " is calculated/set.");
295 if (rbi < 1 || rbi > 9)
297 "ppath.background is invalid. Check your "
298 "calculation of *ppath*?");
299 if (rbi == 3 || rbi == 4)
301 "The propagation path ends inside or at boundary of "
302 "the cloudbox.\nFor this method, *ppath* must be "
303 "calculated in this way:\n ppathCalc( cloudbox_on = 0 ).");
306 if (cloudbox_field.
ncols() != stokes_dim)
308 "Obtained *cloudbox_field* number of Stokes elements inconsistent with "
310 if (cloudbox_field.
nrows() != 1)
312 "Obtained *cloudbox_field* has wrong number of azimuth angles.");
315 "Obtained *cloudbox_field* number of zenith angles inconsistent with "
317 if (cloudbox_field.
nbooks() != 1)
319 "Obtained *cloudbox_field* has wrong number of longitude points.");
322 "Obtained *cloudbox_field* has wrong number of latitude points.");
323 if (cloudbox_field.
nvitrines() != cloudbox_limits[1] - cloudbox_limits[0] + 1)
325 "Obtained *cloudbox_field* number of pressure points inconsistent with "
326 "*cloudbox_limits*.");
329 "Obtained *cloudbox_field* number of frequency points inconsistent with "
348 for (
Index i = 0; i < naux; i++) {
349 iy_aux[i].resize(nf,
ns);
351 if (iy_aux_vars[i] ==
"Optical depth") {
353 else if (iy_aux_vars[i] ==
"Radiative background")
357 os <<
"The only allowed strings in *iy_aux_vars* are:\n"
358 <<
" \"Radiative background\"\n"
359 <<
" \"Optical depth\"\n"
360 <<
"but you have selected: \"" << iy_aux_vars[i] <<
"\"";
361 throw runtime_error(os.str());
385 if (np == 1 && rbi == 1) {
392 ppvar_trans_cumulat = 0;
393 ppvar_trans_partial = 0;
394 for (
Index iv = 0; iv < nf; iv++) {
395 for (
Index is = 0; is <
ns; is++) {
396 ppvar_trans_cumulat(0,iv,is,is) = 1;
397 ppvar_trans_partial(0,iv,is,is) = 1;
422 ppvar_f, ppath, f_grid, atmosphere_dim, rte_alonglos_v, ppvar_wind);
449 Index temperature_derivative_position = -1;
452 if (j_analytical_do) {
462 temperature_derivative_position = iq;
463 do_hse = jacobian_quantities[iq].Subtag() ==
"HSE on";
466 const bool temperature_jacobian =
470 for (
Index ip = 0; ip < np; ip++) {
472 B, dB_dT, ppvar_f(
joker, ip), ppvar_t[ip], temperature_jacobian);
480 propmat_clearsky_agenda,
483 ppvar_mag(
joker, ip),
486 ppvar_vmr(
joker, ip),
501 if (clear2cloudy[ip] + 1) {
512 ppvar_t[
Range(ip, 1)],
520 dK_this_dx[iq] += dKp_dx[iq];)
528 ppvar_pnd(
joker, ip),
537 ppvar_t[
Range(ip, 1)],
553 do_hse ? ppath.
lstep[ip - 1] / (2.0 * ppvar_t[ip - 1]) : 0;
555 do_hse ? ppath.
lstep[ip - 1] / (2.0 * ppvar_t[ip]) : 0;
566 temperature_derivative_position);
582 swap(K_past, K_this);
583 swap(dK_past_dx, dK_this_dx);
593 iy_trans_new = tot_tra[np - 1];
598 for (
Index i = 0; i < naux; i++)
599 if (iy_aux_vars[i] ==
"Optical depth")
600 for (
Index iv = 0; iv < nf; iv++)
601 iy_aux[i](iv,
joker) = -log(ppvar_trans_cumulat(np - 1, iv, 0, 0));
627 lvl_rad[np - 1] = iy;
630 for (
Index ip = np - 2; ip >= 0; ip--) {
631 lvl_rad[ip] = lvl_rad[ip + 1];
641 dlyr_tra_above[ip + 1],
642 dlyr_tra_below[ip + 1],
657 for (
Index ip = 0; ip < lvl_rad.
nelem(); ip++) {
716 const Index& stokes_dim,
718 const Index& atmosphere_dim,
730 const Index& cloudbox_on,
733 const Index& jacobian_do,
737 const Agenda& propmat_clearsky_agenda,
738 const Agenda& water_p_eq_agenda,
739 const String& rt_integration_option,
740 const Agenda& iy_main_agenda,
741 const Agenda& iy_space_agenda,
742 const Agenda& iy_surface_agenda,
743 const Agenda& iy_cloudbox_agenda,
744 const Index& iy_agenda_call1,
745 const Tensor3& iy_transmittance,
747 const Tensor3& surface_props_data,
750 const Index j_analytical_do = jacobian_do ? do_analytical_jacobian<2>(jacobian_quantities) : 0;
756 const Index nq = j_analytical_do ? jacobian_quantities.
nelem() : 0;
763 "ppath.background is invalid. Check your "
764 "calculation of *ppath*?");
766 "A secondary propagation path starting at the "
767 "surface and is going directly into the surface "
768 "is found. This is not allowed.");
783 Index auxOptDepth = -1;
785 for (
Index i = 0; i < naux; i++) {
786 iy_aux[i].resize(nf,
ns);
789 if (iy_aux_vars[i] ==
"Radiative background")
791 else if (iy_aux_vars[i] ==
"Optical depth")
795 "The only allowed strings in *iy_aux_vars* are:\n"
796 " \"Radiative background\"\n"
797 " \"Optical depth\"\n"
798 "but you have selected: \"", iy_aux_vars[i],
"\"")
827 if (np == 1 && rbi == 1) {
834 ppvar_trans_cumulat = 0;
835 ppvar_trans_partial = 0;
836 for (
Index iv = 0; iv < nf; iv++) {
837 for (
Index is = 0; is <
ns; is++) {
838 ppvar_trans_cumulat(0,iv,is,is) = 1;
839 ppvar_trans_partial(0,iv,is,is) = 1;
865 ppvar_f, ppath, f_grid, atmosphere_dim, rte_alonglos_v, ppvar_wind);
867 const bool temperature_jacobian =
875 Vector dB_dT(temperature_jacobian ? nf : 0);
879 Index temperature_derivative_position = -1;
882 if (j_analytical_do) {
883 for (
Index ip = 0; ip < np; ip++) {
884 dK_dx[ip].resize(nq);
890 temperature_derivative_position = iq;
891 do_hse = jacobian_quantities[iq].Subtag() ==
"HSE on";
895 Agenda l_propmat_clearsky_agenda(propmat_clearsky_agenda);
898 bool do_abort =
false;
901#pragma omp parallel for if (!arts_omp_in_parallel()) \
902 firstprivate(l_ws, l_propmat_clearsky_agenda, a, B, dB_dT, S, da_dx, dS_dx)
903 for (
Index ip = 0; ip < np; ip++) {
904 if (do_abort)
continue;
907 B, dB_dT, ppvar_f(
joker, ip), ppvar_t[ip], temperature_jacobian);
916 l_propmat_clearsky_agenda,
919 ppvar_mag(
joker, ip),
922 ppvar_vmr(
joker, ip),
954 }
catch (
const std::runtime_error& e) {
956 os <<
"Runtime-error in source calculation at index " << ip
959#pragma omp critical(iyEmissionStandard_source)
962 fail_msg.push_back(os.str());
967#pragma omp parallel for if (!arts_omp_in_parallel())
968 for (
Index ip = 1; ip < np; ip++) {
969 if (do_abort)
continue;
972 do_hse ? ppath.
lstep[ip - 1] / (2.0 * ppvar_t[ip - 1]) : 0;
974 do_hse ? ppath.
lstep[ip - 1] / (2.0 * ppvar_t[ip]) : 0;
985 temperature_derivative_position);
987 r[ip - 1] = ppath.
lstep[ip - 1];
988 if (temperature_derivative_position >= 0){
989 dr_below[ip][temperature_derivative_position] = dr_dT_past;
990 dr_above[ip][temperature_derivative_position] = dr_dT_this;
992 }
catch (
const std::runtime_error& e) {
994 os <<
"Runtime-error in transmission calculation at index " << ip
997#pragma omp critical(iyEmissionStandard_transmission)
1000 fail_msg.push_back(os.str());
1006 "Error messages from failed cases:\n", fail_msg)
1014 if (iy_agenda_call1)
1015 iy_trans_new = tot_tra[np - 1];
1020 if (auxOptDepth >= 0)
1021 for (
Index iv = 0; iv < nf; iv++)
1022 iy_aux[auxOptDepth](iv, 0) = -std::log(tot_tra[np - 1](iv, 0, 0));
1031 jacobian_quantities,
1048 lvl_rad[np - 1] = iy;
1051 if (rt_integration_option ==
"first order" || rt_integration_option ==
"default") {
1052 for (
Index ip = np - 2; ip >= 0; ip--) {
1053 lvl_rad[ip] = lvl_rad[ip + 1];
1063 dlyr_tra_above[ip + 1],
1064 dlyr_tra_below[ip + 1],
1076 }
else if (rt_integration_option ==
"second order") {
1077 for (
Index ip = np - 2; ip >= 0; ip--) {
1078 lvl_rad[ip] = lvl_rad[ip + 1];
1088 dlyr_tra_above[ip + 1],
1089 dlyr_tra_below[ip + 1],
1108 for (
Index ip = 0; ip < lvl_rad.
nelem(); ip++) {
1112 if (j_analytical_do)
1118 if (j_analytical_do) {
1133 jacobian_quantities,
1138 if (iy_agenda_call1) {
1146 jacobian_quantities,
1161 const Index& atmosphere_dim,
1177 const Index& cloudbox_on,
1180 const Matrix& particle_masses,
1181 const Agenda& ppath_agenda,
1183 const Numeric& ppath_lraytrace,
1184 const Index& iy_agenda_call1,
1186 const Tensor3& iy_transmittance,
1190 const Index& jacobian_do,
1192 const Agenda& iy_independent_beam_approx_agenda,
1193 const Index& return_atm1d,
1194 const Index& skip_vmr,
1195 const Index& skip_pnd,
1196 const Index& return_masses,
1200 "Jacobians not provided by the method, *jacobian_do* "
1203 "This method does not yet support non-empty *nlte_field*.");
1205 "This method does not yet support non-empty *wind_u_field*.");
1207 "This method does not yet support non-empty *wind_v_field*.");
1209 "This method does not yet support non-empty *wind_w_field*.");
1211 "This method does not yet support non-empty *mag_u_field*.");
1213 "This method does not yet support non-empty *mag_v_field*.");
1215 "This method does not yet support non-empty *mag_w_field*.");
1217 if (return_masses) {
1219 "Sizes of *pnd_field* and *particle_masses* "
1220 "are inconsistent.");
1253 if (cloudbox_on && ppath.
end_lstep == 0) {
1254 Vector los_tmp = rte_los;
1255 if (
abs(rte_los[0]) < 90) {
1277 const Index np = ppath.
np + ppath2.
np - 1;
1283 for (
Index i = 0; i < ppath.
np; i++) {
1284 const Index ip = ppath.
np - i - 1;
1285 gp_p[i] = ppath.
gp_p[ip];
1286 if (atmosphere_dim > 1) {
1287 gp_lat[i] = ppath.
gp_lat[ip];
1288 if (atmosphere_dim == 3) {
1289 gp_lon[i] = ppath.
gp_lon[ip];
1294 for (
Index i = ppath.
np; i < np; i++) {
1295 const Index ip = i - ppath.
np + 1;
1296 gp_p[i] = ppath2.
gp_p[ip];
1297 if (atmosphere_dim > 1) {
1298 gp_lat[i] = ppath2.
gp_lat[ip];
1299 if (atmosphere_dim == 3) {
1300 gp_lon[i] = ppath2.
gp_lon[ip];
1306 for (
Index i = 0; i < ppath2.
np - 1; i++) {
1307 const Index ip = ppath2.
np - i - 1;
1308 gp_p[i] = ppath2.
gp_p[ip];
1309 if (atmosphere_dim > 1) {
1310 gp_lat[i] = ppath2.
gp_lat[ip];
1311 if (atmosphere_dim == 3) {
1312 gp_lon[i] = ppath2.
gp_lon[ip];
1317 for (
Index i = ppath2.
np - 1; i < np; i++) {
1318 const Index ip = i - ppath2.
np + 1;
1319 gp_p[i] = ppath.
gp_p[ip];
1320 if (atmosphere_dim > 1) {
1321 gp_lat[i] = ppath.
gp_lat[ip];
1322 if (atmosphere_dim == 3) {
1323 gp_lon[i] = ppath.
gp_lon[ip];
1334 itw2p(p1, p_grid, gp_p, itw);
1338 Vector lat_true1(1), lon_true1(1);
1339 if (atmosphere_dim == 3) {
1342 interp(lat_true1, itw, lat_grid, gp1);
1345 interp(lon_true1, itw, lon_grid, gp1);
1346 }
else if (atmosphere_dim == 2) {
1349 interp(lat_true1, itw, lat_true, gp1);
1350 interp(lon_true1, itw, lon_true, gp1);
1352 lat_true1[0] = lat_true[0];
1353 lon_true1[0] = lon_true[0];
1362 t1(
joker, 0, 0), atmosphere_dim, t_field, gp_p, gp_lat, gp_lon, itw);
1367 z1(
joker, 0, 0), atmosphere_dim, z_field, gp_p, gp_lat, gp_lon, itw);
1371 for (
Index is = 0; is < vmr_field.
nbooks(); is++) {
1383 zsurf1(0, 0) = z1(0, 0, 0);
1387 pos1[0] = rte_pos[0];
1389 los1[0] =
abs(rte_los[0]);
1391 if (rte_pos2.
nelem()) {
1397 Index cbox_on1 = cloudbox_on;
1405 for (
Index i = 0; i < np; i++) {
1429 const Index extra_bot = ifirst == 0 ? 0 : 1;
1430 const Index extra_top = ilast == np - 1 ? 0 : 1;
1432 cbox_lims1.resize(2);
1433 cbox_lims1[0] = ifirst - extra_bot;
1434 cbox_lims1[1] = ilast + extra_top;
1438 pnd1.
resize(pnd_field.
nbooks(), cbox_lims1[1] - cbox_lims1[0] + 1, 1, 1);
1443 for (
Index i = extra_bot; i < pnd1.
npages() - extra_top; i++) {
1444 const Index i0 = cbox_lims1[0] + i;
1471 const Index adim1 = 1;
1505 iy_independent_beam_approx_agenda);
1513 const Index nmass = return_masses ? particle_masses.
ncols() : 0;
1514 const Index ntot = 2 + nvmr + npnd + nmass;
1516 atm_fields_compact.
resize(ntot, np, 1, 1);
1519 field_names[0] =
"Geometric altitudes";
1523 field_names[1] =
"Temperature";
1528 for (
Index i = 0; i < nvmr; i++) {
1529 const Index iout = 2 + i;
1531 sstr <<
"VMR species " << i;
1532 field_names[iout] = sstr.str();
1533 atm_fields_compact.
data(iout,
joker, 0, 0) = vmr1(i,
joker, 0, 0);
1539 for (
Index i = 0; i < npnd; i++) {
1540 const Index iout = 2 + nvmr + i;
1542 sstr <<
"Scattering element " << i;
1543 field_names[iout] = sstr.str();
1544 atm_fields_compact.
data(iout,
joker, 0, 0) = 0;
1545 atm_fields_compact.
data(
1546 iout,
Range(cbox_lims1[0], pnd1.
npages()), 0, 0) =
1547 pnd1(i,
joker, 0, 0);
1553 for (
Index i = 0; i < nmass; i++) {
1554 const Index iout = 2 + nvmr + npnd + i;
1556 sstr <<
"Mass category " << i;
1557 field_names[iout] = sstr.str();
1558 atm_fields_compact.
data(iout,
joker, 0, 0) = 0;
1559 for (
Index ip = cbox_lims1[0]; ip < pnd1.
npages(); ip++) {
1561 atm_fields_compact.
data(iout, ip, 0, 0) +=
1562 particle_masses(is, i) * pnd1(is, ip, 0, 0);
1571 "Data created by *iyIndependentBeamApproximation*");
1574 "Atmospheric quantity");
1592 const Index& iy_agenda_call1,
1593 const Tensor3& iy_transmittance,
1597 const Index& stokes_dim,
1599 const Agenda& iy_loop_freqs_agenda,
1603 "Recursive usage not possible (iy_agenda_call1 must be 1).");
1605 "*iy_transmittance* must be empty.");
1609 for (
Index i = 0; i < nf; i++) {
1628 iy_loop_freqs_agenda);
1632 iy.
resize(nf, stokes_dim);
1634 iy_aux.resize(iy_aux1.
nelem());
1636 iy_aux[
q].resize(nf, stokes_dim);
1639 diy_dx.resize(diy_dx1.
nelem());
1641 diy_dx[
q].resize(diy_dx1[
q].npages(), nf, stokes_dim);
1661 const Index& iy_agenda_call1,
1662 const Tensor3& iy_transmittance,
1666 const Index& jacobian_do,
1667 const Index& atmosphere_dim,
1674 const Vector& refellipsoid,
1676 const Index& cloudbox_on,
1678 const Index& stokes_dim,
1681 const Agenda& iy_space_agenda,
1682 const Agenda& surface_rtprop_agenda,
1683 const Agenda& propmat_clearsky_agenda,
1684 const Agenda& ppath_step_agenda,
1686 const Numeric& ppath_lraytrace,
1690 const Index& mc_max_time,
1691 const Index& mc_max_iter,
1692 const Index& mc_min_iter,
1693 const Numeric& mc_taustep_limit,
1694 const Index& t_interp_order,
1698 "Only 3D atmospheres are allowed (atmosphere_dim must be 3)");
1700 "The cloudbox must be activated (cloudbox_on must be 1)");
1702 "This method does not provide any jacobians (jacobian_do must be 0)");
1704 "Recursive usage not possible (iy_agenda_call1 must be 1)");
1706 "*iy_transmittance* must be empty");
1712 iy.
resize(nf, stokes_dim);
1716 Index auxError = -1;
1719 iy_aux.resize(naux);
1721 for (
Index i = 0; i < naux; i++) {
1722 if (iy_aux_vars[i] ==
"Error (uncorrelated)") {
1724 iy_aux[i].resize(nf, stokes_dim);
1727 "In *iy_aux_vars* you have included: \"", iy_aux_vars[i],
1728 "\"\nThis choice is not recognised.")
1740 Matrix pos(1, 3), los(1, 2);
1742 pos(0,
joker) = rte_pos;
1743 los(0,
joker) = rte_los;
1746 Agenda l_ppath_step_agenda(ppath_step_agenda);
1747 Agenda l_iy_space_agenda(iy_space_agenda);
1748 Agenda l_propmat_clearsky_agenda(propmat_clearsky_agenda);
1749 Agenda l_surface_rtprop_agenda(surface_rtprop_agenda);
1752 bool failed =
false;
1755#pragma omp parallel for if (!arts_omp_in_parallel() && nf > 1) \
1756 firstprivate(l_ws, \
1757 l_ppath_step_agenda, \
1758 l_iy_space_agenda, \
1759 l_propmat_clearsky_agenda, \
1760 l_surface_rtprop_agenda)
1761 for (
Index f_index = 0; f_index < nf; f_index++) {
1762 if (failed)
continue;
1771 Index mc_iteration_count;
1789 l_ppath_step_agenda,
1793 l_surface_rtprop_agenda,
1794 l_propmat_clearsky_agenda,
1824 iy(f_index,
joker) = y;
1826 if (auxError >= 0) {
1827 iy_aux[auxError](f_index,
joker) = mc_error;
1829 }
catch (
const std::exception& e) {
1831 os <<
"Error for f_index = " << f_index <<
" (" << f_grid[f_index]
1834#pragma omp critical(iyMC_fail)
1837 fail_msg = os.str();
1850 const Index& jacobian_do,
1854 "*iy_aux* and *iy_aux_vars* must have the same "
1855 "number of elements.");
1858 "This method can not provide any jacobians and "
1859 "*jacobian_do* must be 0.");
1863 for (
Index i = 0; i < iy_aux.
nelem() && !ready; i++) {
1864 if (iy_aux_vars[i] == aux_var) {
1871 "The selected auxiliary variable to insert in *iy* "
1872 "is either not defined at all or is not set.");
1877 Matrix& ppvar_optical_depth,
1878 const Tensor4& ppvar_trans_cumulat,
1880 ppvar_optical_depth = ppvar_trans_cumulat(
joker,
joker, 0, 0);
1881 transform(ppvar_optical_depth, log, ppvar_optical_depth);
1882 ppvar_optical_depth *= -1;
1895 const Index& atmgeom_checked,
1896 const Index& atmfields_checked,
1897 const Index& atmosphere_dim,
1899 const Index& cloudbox_on,
1900 const Index& cloudbox_checked,
1901 const Index& scat_data_checked,
1902 const Index& sensor_checked,
1903 const Index& stokes_dim,
1905 const Matrix& sensor_pos,
1906 const Matrix& sensor_los,
1907 const Matrix& transmitter_pos,
1908 const Matrix& mblock_dlos_grid,
1909 const Sparse& sensor_response,
1910 const Vector& sensor_response_f,
1912 const Matrix& sensor_response_dlos,
1914 const Agenda& iy_main_agenda,
1915 const Agenda& geo_pos_agenda,
1916 const Agenda& jacobian_agenda,
1917 const Index& jacobian_do,
1928 "The frequency grid is empty.");
1931 "All frequencies in *f_grid* must be > 0.");
1934 "The atmospheric fields must be flagged to have\n"
1935 "passed a consistency check (atmfields_checked=1).");
1937 "The atmospheric geometry must be flagged to have\n"
1938 "passed a consistency check (atmgeom_checked=1).");
1940 "The cloudbox must be flagged to have\n"
1941 "passed a consistency check (cloudbox_checked=1).");
1944 "The scattering data must be flagged to have\n"
1945 "passed a consistency check (scat_data_checked=1).");
1947 "The sensor variables must be flagged to have\n"
1948 "passed a consistency check (sensor_checked=1).");
1955 const Index niyb = nf * nlos * stokes_dim;
1964 y_f.
resize(nmblock * n1y);
1965 y_pol.resize(nmblock * n1y);
1968 y_geo.
resize(nmblock * n1y, 5);
1977 Index j_analytical_do = 0;
1983 jacobian.
resize(nmblock * n1y,
1984 jacobian_indices[jacobian_indices.
nelem() - 1][1] + 1);
1997 bool failed =
false;
2000 (nf <= nmblock && nmblock >= nlos)) {
2001 out3 <<
" Parallelizing mblock loop (" << nmblock <<
" iterations)\n";
2006 Agenda l_jacobian_agenda(jacobian_agenda);
2007 Agenda l_iy_main_agenda(iy_main_agenda);
2008 Agenda l_geo_pos_agenda(geo_pos_agenda);
2010#pragma omp parallel for firstprivate( \
2011 l_ws, l_jacobian_agenda, l_iy_main_agenda, l_geo_pos_agenda)
2012 for (
Index mblock_index = 0; mblock_index < nmblock; mblock_index++) {
2014 if (failed)
continue;
2038 sensor_response_pol,
2039 sensor_response_dlos,
2045 jacobian_quantities,
2054 out3 <<
" Not parallelizing mblock loop (" << nmblock <<
" iterations)\n";
2056 for (
Index mblock_index = 0; mblock_index < nmblock; mblock_index++) {
2058 if (failed)
continue;
2082 sensor_response_pol,
2083 sensor_response_dlos,
2089 jacobian_quantities,
2108 y_aux[
q].resize(nmblock * n1y);
2110 for (
Index mblock_index = 0; mblock_index < nmblock; mblock_index++) {
2111 const Range rowind =
2117 if (iy_aux_vars[
q] ==
"Error (uncorrelated)") {
2118 for (
Index i = 0; i < n1y; i++) {
2119 const Index row = row0 + i;
2121 for (
Index j = 0; j < niyb; j++) {
2123 pow(sensor_response(i, j) * iyb_aux_array[mblock_index][
q][j],
2126 y_aux[
q][row] =
sqrt(y_aux[
q][row]);
2129 mult(y_aux[
q][rowind], sensor_response, iyb_aux_array[mblock_index][
q]);
2146 const Index& atmfields_checked,
2147 const Index& atmgeom_checked,
2148 const Index& atmosphere_dim,
2150 const Index& cloudbox_on,
2151 const Index& cloudbox_checked,
2152 const Index& scat_data_checked,
2153 const Index& sensor_checked,
2154 const Index& stokes_dim,
2156 const Matrix& sensor_pos,
2157 const Matrix& sensor_los,
2158 const Matrix& transmitter_pos,
2159 const Matrix& mblock_dlos_grid,
2160 const Sparse& sensor_response,
2161 const Vector& sensor_response_f,
2163 const Matrix& sensor_response_dlos,
2165 const Agenda& iy_main_agenda,
2166 const Agenda& geo_pos_agenda,
2167 const Agenda& jacobian_agenda,
2168 const Index& jacobian_do,
2171 const Index& append_instrument_wfs,
2178 jacobian_indices_copy, any_affine, jacobian_quantities_copy,
true);
2186 "Input *y* is empty. Use *yCalc*");
2188 "Lengths of input *y* and *y_f* are inconsistent.");
2190 "Lengths of input *y* and *y_pol* are inconsistent.");
2192 "Sizes of input *y* and *y_pos* are inconsistent.");
2194 "Sizes of input *y* and *y_los* are inconsistent.");
2196 "Sizes of input *y* and *y_geo* are inconsistent.");
2198 nrq1 = jacobian_quantities_copy.
nelem();
2200 "Sizes of *y* and *jacobian* are inconsistent.");
2202 "Size of input *jacobian* and size implied "
2203 "*jacobian_quantities_copy* are inconsistent.");
2209 Matrix y_pos2, y_los2, y_geo2, jacobian2;
2238 sensor_response_pol,
2239 sensor_response_dlos,
2245 jacobian_quantities,
2251 "Different number of columns in *y_pos* between the measurements.");
2253 "Different number of columns in *y_los* between the measurements.");
2261 const Vector y1 = y, y_f1 = y_f;
2262 const Matrix y_pos1 = y_pos, y_los1 = y_los, y_geo1 = y_geo;
2267 y[
Range(0, n1)] = y1;
2268 y[
Range(n1, n2)] = y2;
2271 y_f[
Range(0, n1)] = y_f1;
2272 y_f[
Range(n1, n2)] = y_f2;
2278 y_los.
resize(n1 + n2, y_los1.ncols());
2282 y_geo.
resize(n1 + n2, y_geo1.ncols());
2286 y_pol.resize(n1 + n2);
2287 for (
Index i = 0; i < n1; i++) {
2288 y_pol[i] = y_pol1[i];
2290 for (
Index i = 0; i < n2; i++) {
2291 y_pol[n1 + i] = y_pol2[i];
2302 y_aux[
a].resize(n1 + n2);
2304 y_aux[
a][
Range(0, n1)] = y_aux1[
a];
2306 y_aux[
a][
Range(0, n1)] = 0;
2309 y_aux[
a][
Range(n1, n2)] = y_aux2[
a];
2311 y_aux[
a][
Range(n1, n2)] = 0;
2323 jacobian_quantities = jacobian_quantities_copy;
2324 jacobian_indices = jacobian_indices_copy;
2329 const Index nrq2 = jacobian_quantities2.
nelem();
2332 for (
Index q2 = 0; q2 < nrq2; q2++) {
2337 if (jacobian_quantities2[q2].Target().isSpeciesVMR() ||
2340 jacobian_quantities2[q2].Target().isWind() ||
2341 jacobian_quantities2[q2] == Jacobian::Special::SurfaceString ||
2342 append_instrument_wfs) {
2344 if (jacobian_quantities2[q2].Target().sameTargetType(jacobian_quantities_copy[
q1].Target())) {
2345 if (jacobian_quantities2[q2].Target().isSpeciesVMR()) {
2346 if (jacobian_quantities2[q2].Subtag() ==
2347 jacobian_quantities_copy[
q1].Subtag()) {
2348 if (jacobian_quantities2[q2].Mode() ==
2349 jacobian_quantities_copy[
q1].Mode()) {
2353 "Jacobians for ", jacobian_quantities2[q2],
2354 " shall be appended.\nThis requires "
2355 "that the same retrieval unit is used "
2356 "but it seems that this requirement is "
2361 if (jacobian_quantities2[q2].Subtag() ==
2362 jacobian_quantities_copy[
q1].Subtag()) {
2366 "Jacobians for ", jacobian_quantities2[q2],
2367 " shall be appended.\nThis requires "
2368 "that HSE is either ON or OFF for both "
2369 "parts but it seems that this requirement "
2373 if ((jacobian_quantities2[q2].Subtag() ==
2374 jacobian_quantities_copy[
q1].Subtag()) &&
2375 (jacobian_quantities2[q2].SubSubtag() ==
2376 jacobian_quantities_copy[
q1].SubSubtag())) {
2379 }
else if (jacobian_quantities2[q2].Subtag() == jacobian_quantities_copy[
q1].Subtag()) {
2388 map_table[q2] = jacobian_quantities.
nelem();
2389 jacobian_quantities.push_back(jacobian_quantities2[q2]);
2391 indices[0] = jacobian_indices[jacobian_indices.
nelem() - 1][1] + 1;
2393 indices[0] + jacobian_indices2[q2][1] - jacobian_indices2[q2][0];
2394 jacobian_indices.push_back(indices);
2398 map_table[q2] = pos;
2400 ArrayOfVector grids1 = jacobian_quantities_copy[pos].Grids();
2402 bool any_wrong =
false;
2407 if (grids1[g].
nelem() != grids2[g].
nelem()) {
2410 for (
Index e = 0; e < grids1[g].
nelem(); e++) {
2411 const Numeric v1 = grids1[g][e];
2412 const Numeric v2 = grids2[g][e];
2413 if ((v1 == 0 &&
abs(v2) > 1e-9) ||
abs(v1 - v2) / v1 > 1e-6) {
2422 "Jacobians for ", jacobian_quantities2[q2],
2423 " shall be appended.\nThis requires that the "
2424 "same grids are used for both measurements,\nbut "
2425 "it seems that this requirement is not met.")
2432 const Index nrq = jacobian_quantities.
nelem();
2433 const Matrix jacobian1 = jacobian;
2435 jacobian.
resize(n1 + n2, jacobian_indices[nrq - 1][1] + 1);
2439 jacobian(
Range(0, n1),
Range(0, jacobian_indices_copy[nrq1 - 1][1] + 1)) =
2442 for (
Index q2 = 0; q2 < nrq2; q2++) {
2443 jacobian(
Range(n1, n2),
2444 Range(jacobian_indices[map_table[q2]][0],
2445 jacobian_indices[map_table[q2]][1] -
2446 jacobian_indices[map_table[q2]][0] + 1)) =
2449 Range(jacobian_indices2[q2][0],
2450 jacobian_indices2[q2][1] - jacobian_indices2[q2][0] + 1));
2463 "No need to use this method with *iy_unit* = \"1\".");
2466 "The spectrum vector *y* is required to have original radiance\n"
2467 "unit, but this seems not to be the case. This as a value above\n"
2468 "1e-3 is found in *y*.")
2474 const bool do_j = jacobian.
nrows() == ny;
2478 "The method can not be used with jacobian quantities that are not\n"
2479 "obtained through radiative transfer calculations. One example on\n"
2480 "quantity that can not be handled is *jacobianAddPolyfit*.\n"
2481 "The maximum value of *jacobian* indicates that one or several\n"
2482 "such jacobian quantities are included.")
2486 if (iy_unit ==
"PlanckBT") {
2499 while (i0 + n < ny && y_f[i0] == y_f[i0 + n]) {
2505 bool any_quv =
false;
2507 for (
Index i = 0; i < n; i++) {
2508 const Index ix = i0 + i;
2510 i_pol[i] = y_pol[ix];
2511 if (i_pol[i] > 1 && i_pol[i] < 5) {
2521 "The conversion to PlanckBT, of the Jacobian and "
2522 "errors for Q, U and V, requires that I (first Stokes "
2523 "element) is at hand and that the data are sorted in "
2524 "such way that I comes first for each frequency.")
2535 y[ii] = yv(0,
joker);
2549 for (
Index i = 0; i < ny; i++) {
2551 i_pol[0] = y_pol[i];
Array< Index > ArrayOfIndex
An array of Index.
Array< Tensor3 > ArrayOfTensor3
An array of Tensor3.
The global header file for ARTS.
int arts_omp_get_max_threads()
Wrapper for omp_get_max_threads.
Header file for helper functions for OpenMP.
void iy_loop_freqs_agendaExecute(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, Ppath &ppath, ArrayOfTensor3 &diy_dx, const Index iy_agenda_call1, const Tensor3 &iy_transmittance, const ArrayOfString &iy_aux_vars, const Index iy_id, const Vector &f_grid, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Agenda &input_agenda)
void iy_main_agendaExecute(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, Ppath &ppath, ArrayOfTensor3 &diy_dx, const Index iy_agenda_call1, const Tensor3 &iy_transmittance, const ArrayOfString &iy_aux_vars, const Index iy_id, const String &iy_unit, const Index cloudbox_on, const Index jacobian_do, const Vector &f_grid, const EnergyLevelMap &nlte_field, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Agenda &input_agenda)
void ppath_agendaExecute(Workspace &ws, Ppath &ppath, const Numeric ppath_lmax, const Numeric ppath_lraytrace, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Index cloudbox_on, const Index ppath_inside_cloudbox_do, const Vector &f_grid, const Agenda &input_agenda)
void iy_independent_beam_approx_agendaExecute(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, Ppath &ppath, ArrayOfTensor3 &diy_dx, const Index iy_agenda_call1, const String &iy_unit, const Tensor3 &iy_transmittance, const ArrayOfString &iy_aux_vars, const Index iy_id, const Index atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &lat_true, const Vector &lon_true, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const Matrix &z_surface, const Numeric ppath_lmax, const Numeric ppath_lraytrace, const Index cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const Index jacobian_do, const Vector &f_grid, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Agenda &input_agenda)
This can be used to make arrays out of anything.
Index nelem() const ARTS_NOEXCEPT
Number of elements.
Index nrows() const noexcept
Index ncols() const noexcept
bool empty() const noexcept
Index ncols() const
Returns the number of columns.
bool empty() const noexcept
Index nbooks() const noexcept
Index npages() const noexcept
Index ncols() const noexcept
Index npages() const noexcept
Index nrows() const noexcept
Index nlibraries() const noexcept
Index nvitrines() const noexcept
Index nshelves() const noexcept
Index nbooks() const noexcept
Index nelem() const noexcept
Returns the number of elements.
bool empty() const noexcept
Returns true if variable size is zero.
void resize(const GriddedField4 &gf)
Make this GriddedField4 the same size as the given one.
void set_name(const String &s)
Set name of this gridded field.
void set_grid_name(Index i, const String &s)
Set grid name.
void set_grid(Index i, const Vector &g)
Set a numeric grid.
void resize(Index r, Index c)
Resize function.
constexpr Index get_start() const noexcept
Returns the start index of the range.
Stokes vector is as Propagation matrix but only has 4 possible values.
void resize(Index p, Index r, Index c)
Resize function.
void resize(Index b, Index p, Index r, Index c)
Resize function.
void resize(Index n)
Resize function.
bool is_gp_inside_cloudbox(const GridPos &gp_p, const GridPos &gp_lat, const GridPos &gp_lon, const ArrayOfIndex &cloudbox_limits, const bool &include_boundaries, const Index &atmosphere_dim)
void mult(MatrixView C, ConstMatrixView A, const Block &B)
#define ARTS_ASSERT(condition,...)
#define ARTS_USER_ERROR(...)
#define ARTS_USER_ERROR_IF(condition,...)
Numeric interp(ConstVectorView itw, ConstVectorView a, const GridPos &tc)
Red 1D Interpolate.
ArrayOfIndex get_pointers_for_analytical_species(const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfSpeciesTag &abs_species)
Help function for analytical jacobian calculations.
ArrayOfTensor3 get_standard_starting_diy_dx(const ArrayOfRetrievalQuantity &jacobian_quantities, Index np, Index nf, Index ns, bool active)
Help function for analytical jacobian calculations.
bool do_temperature_jacobian(const ArrayOfRetrievalQuantity &js) noexcept
Returns if the array wants the temperature derivative.
ArrayOfTensor3 get_standard_diy_dpath(const ArrayOfRetrievalQuantity &jacobian_quantities, Index np, Index nf, Index ns, bool active)
Help function for analytical jacobian calculations.
ArrayOfIndex get_pointers_for_scat_species(const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfString &scat_species, const bool cloudbox_on)
Help function for analytical jacobian calculations.
void jac_ranges_indices(ArrayOfArrayOfIndex &jis, bool &any_affine, const ArrayOfRetrievalQuantity &jqs, const bool &before_affine)
Determines the index range inside x and the Jacobian for each retrieval quantity.
Routines for setting up the jacobian.
#define FOR_ANALYTICAL_JACOBIANS_DO2(what_to_do)
#define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do)
Header file for logic.cc.
void MCSetSeedFromTime(Index &mc_seed, const Verbosity &)
WORKSPACE METHOD: MCSetSeedFromTime.
void MCGeneral(Workspace &ws, Vector &y, Index &mc_iteration_count, Vector &mc_error, Tensor3 &mc_points, ArrayOfIndex &mc_source_domain, ArrayOfIndex &mc_scat_order, const MCAntenna &mc_antenna, const Vector &f_grid, const Index &f_index, const Matrix &sensor_pos, const Matrix &sensor_los, const Index &stokes_dim, const Index &atmosphere_dim, const Agenda &ppath_step_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Agenda &iy_space_agenda, const Agenda &surface_rtprop_agenda, const Agenda &propmat_clearsky_agenda, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Vector &refellipsoid, const Matrix &z_surface, const Tensor3 &t_field, const Tensor4 &vmr_field, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfArrayOfSingleScatteringData &scat_data, const Index &atmfields_checked, const Index &atmgeom_checked, const Index &scat_data_checked, const Index &cloudbox_checked, const String &iy_unit, const Index &mc_seed, const Numeric &std_err, const Index &max_time, const Index &max_iter, const Index &min_iter, const Numeric &taustep_limit, const Index &l_mc_scat_order, const Index &t_interp_order, const Verbosity &verbosity)
WORKSPACE METHOD: MCGeneral.
const Index GFIELD4_P_GRID
void iyLoopFrequencies(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, Ppath &ppath, ArrayOfTensor3 &diy_dx, const ArrayOfString &iy_aux_vars, const Index &iy_agenda_call1, const Tensor3 &iy_transmittance, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Index &stokes_dim, const Vector &f_grid, const Agenda &iy_loop_freqs_agenda, const Verbosity &)
WORKSPACE METHOD: iyLoopFrequencies.
void iyApplyUnit(Matrix &iy, ArrayOfMatrix &iy_aux, const Index &stokes_dim, const Vector &f_grid, const ArrayOfString &iy_aux_vars, const String &iy_unit, const Verbosity &)
WORKSPACE METHOD: iyApplyUnit.
void iyCalc(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, Ppath &ppath, const Index &atmfields_checked, const Index &atmgeom_checked, const ArrayOfString &iy_aux_vars, const Index &iy_id, const Index &cloudbox_on, const Index &cloudbox_checked, const Index &scat_data_checked, const Vector &f_grid, const EnergyLevelMap &nlte_field, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const String &iy_unit, const Agenda &iy_main_agenda, const Verbosity &)
WORKSPACE METHOD: iyCalc.
void iyEmissionStandard(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, ArrayOfTensor3 &diy_dx, Vector &ppvar_p, Vector &ppvar_t, EnergyLevelMap &ppvar_nlte, Matrix &ppvar_vmr, Matrix &ppvar_wind, Matrix &ppvar_mag, Matrix &ppvar_f, Tensor3 &ppvar_iy, Tensor4 &ppvar_trans_cumulat, Tensor4 &ppvar_trans_partial, const Index &iy_id, const Index &stokes_dim, const Vector &f_grid, const Index &atmosphere_dim, const Vector &p_grid, const Tensor3 &t_field, const EnergyLevelMap &nlte_field, const Tensor4 &vmr_field, const ArrayOfArrayOfSpeciesTag &abs_species, const Tensor3 &wind_u_field, const Tensor3 &wind_v_field, const Tensor3 &wind_w_field, const Tensor3 &mag_u_field, const Tensor3 &mag_v_field, const Tensor3 &mag_w_field, const Index &cloudbox_on, const String &iy_unit, const ArrayOfString &iy_aux_vars, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Ppath &ppath, const Vector &rte_pos2, const Agenda &propmat_clearsky_agenda, const Agenda &water_p_eq_agenda, const String &rt_integration_option, const Agenda &iy_main_agenda, const Agenda &iy_space_agenda, const Agenda &iy_surface_agenda, const Agenda &iy_cloudbox_agenda, const Index &iy_agenda_call1, const Tensor3 &iy_transmittance, const Numeric &rte_alonglos_v, const Tensor3 &surface_props_data, const Verbosity &verbosity)
WORKSPACE METHOD: iyEmissionStandard.
void yCalcAppend(Workspace &ws, Vector &y, Vector &y_f, ArrayOfIndex &y_pol, Matrix &y_pos, Matrix &y_los, ArrayOfVector &y_aux, Matrix &y_geo, Matrix &jacobian, ArrayOfRetrievalQuantity &jacobian_quantities, const Index &atmfields_checked, const Index &atmgeom_checked, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &cloudbox_checked, const Index &scat_data_checked, const Index &sensor_checked, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_pos, const Matrix &sensor_los, const Matrix &transmitter_pos, const Matrix &mblock_dlos_grid, const Sparse &sensor_response, const Vector &sensor_response_f, const ArrayOfIndex &sensor_response_pol, const Matrix &sensor_response_dlos, const String &iy_unit, const Agenda &iy_main_agenda, const Agenda &geo_pos_agenda, const Agenda &jacobian_agenda, const Index &jacobian_do, const ArrayOfString &iy_aux_vars, const ArrayOfRetrievalQuantity &jacobian_quantities_copy, const Index &append_instrument_wfs, const Verbosity &verbosity)
WORKSPACE METHOD: yCalcAppend.
void iyIndependentBeamApproximation(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, Ppath &ppath, ArrayOfTensor3 &diy_dx, GriddedField4 &atm_fields_compact, const Index &iy_id, const Vector &f_grid, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Vector &lat_true, const Vector &lon_true, const Tensor3 &t_field, const Tensor3 &z_field, const Tensor4 &vmr_field, const EnergyLevelMap &nlte_field, const Tensor3 &wind_u_field, const Tensor3 &wind_v_field, const Tensor3 &wind_w_field, const Tensor3 &mag_u_field, const Tensor3 &mag_v_field, const Tensor3 &mag_w_field, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const Matrix &particle_masses, const Agenda &ppath_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Index &iy_agenda_call1, const String &iy_unit, const Tensor3 &iy_transmittance, const Vector &rte_pos, const Vector &rte_los, const Vector &rte_pos2, const Index &jacobian_do, const ArrayOfString &iy_aux_vars, const Agenda &iy_independent_beam_approx_agenda, const Index &return_atm1d, const Index &skip_vmr, const Index &skip_pnd, const Index &return_masses, const Verbosity &)
WORKSPACE METHOD: iyIndependentBeamApproximation.
void iyMC(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, ArrayOfTensor3 &diy_dx, const Index &iy_agenda_call1, const Tensor3 &iy_transmittance, const Vector &rte_pos, const Vector &rte_los, const ArrayOfString &iy_aux_vars, const Index &jacobian_do, const Index &atmosphere_dim, const Vector &p_grid, const Vector &lat_grid, const Vector &lon_grid, const Tensor3 &z_field, const Tensor3 &t_field, const Tensor4 &vmr_field, const Vector &refellipsoid, const Matrix &z_surface, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Index &stokes_dim, const Vector &f_grid, const ArrayOfArrayOfSingleScatteringData &scat_data, const Agenda &iy_space_agenda, const Agenda &surface_rtprop_agenda, const Agenda &propmat_clearsky_agenda, const Agenda &ppath_step_agenda, const Numeric &ppath_lmax, const Numeric &ppath_lraytrace, const Tensor4 &pnd_field, const String &iy_unit, const Numeric &mc_std_err, const Index &mc_max_time, const Index &mc_max_iter, const Index &mc_min_iter, const Numeric &mc_taustep_limit, const Index &t_interp_order, const Verbosity &verbosity)
WORKSPACE METHOD: iyMC.
const Index GFIELD4_FIELD_NAMES
void yApplyUnit(Vector &y, Matrix &jacobian, const Vector &y_f, const ArrayOfIndex &y_pol, const String &iy_unit, const Verbosity &)
WORKSPACE METHOD: yApplyUnit.
void iyReplaceFromAux(Matrix &iy, const ArrayOfMatrix &iy_aux, const ArrayOfString &iy_aux_vars, const Index &jacobian_do, const String &aux_var, const Verbosity &)
WORKSPACE METHOD: iyReplaceFromAux.
const Numeric SPEED_OF_LIGHT
void yCalc(Workspace &ws, Vector &y, Vector &y_f, ArrayOfIndex &y_pol, Matrix &y_pos, Matrix &y_los, ArrayOfVector &y_aux, Matrix &y_geo, Matrix &jacobian, const Index &atmgeom_checked, const Index &atmfields_checked, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &cloudbox_checked, const Index &scat_data_checked, const Index &sensor_checked, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_pos, const Matrix &sensor_los, const Matrix &transmitter_pos, const Matrix &mblock_dlos_grid, const Sparse &sensor_response, const Vector &sensor_response_f, const ArrayOfIndex &sensor_response_pol, const Matrix &sensor_response_dlos, const String &iy_unit, const Agenda &iy_main_agenda, const Agenda &geo_pos_agenda, const Agenda &jacobian_agenda, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfString &iy_aux_vars, const Verbosity &verbosity)
WORKSPACE METHOD: yCalc.
const Index GFIELD4_LAT_GRID
void ppvar_optical_depthFromPpvar_trans_cumulat(Matrix &ppvar_optical_depth, const Tensor4 &ppvar_trans_cumulat, const Verbosity &)
WORKSPACE METHOD: ppvar_optical_depthFromPpvar_trans_cumulat.
void iyEmissionHybrid(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, ArrayOfTensor3 &diy_dx, Vector &ppvar_p, Vector &ppvar_t, EnergyLevelMap &ppvar_nlte, Matrix &ppvar_vmr, Matrix &ppvar_wind, Matrix &ppvar_mag, Matrix &ppvar_pnd, Matrix &ppvar_f, Tensor3 &ppvar_iy, Tensor4 &ppvar_trans_cumulat, Tensor4 &ppvar_trans_partial, const Index &iy_id, const Index &stokes_dim, const Vector &f_grid, const Index &atmosphere_dim, const Vector &p_grid, const Tensor3 &t_field, const EnergyLevelMap &nlte_field, const Tensor4 &vmr_field, const ArrayOfArrayOfSpeciesTag &abs_species, const Tensor3 &wind_u_field, const Tensor3 &wind_v_field, const Tensor3 &wind_w_field, const Tensor3 &mag_u_field, const Tensor3 &mag_v_field, const Tensor3 &mag_w_field, const Index &cloudbox_on, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfTensor4 &dpnd_field_dx, const ArrayOfString &scat_species, const ArrayOfArrayOfSingleScatteringData &scat_data, const String &iy_unit, const ArrayOfString &iy_aux_vars, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Agenda &propmat_clearsky_agenda, const Agenda &water_p_eq_agenda, const String &rt_integration_option, const Agenda &iy_main_agenda, const Agenda &iy_space_agenda, const Agenda &iy_surface_agenda, const Agenda &iy_cloudbox_agenda, const Index &iy_agenda_call1, const Tensor3 &iy_transmittance, const Ppath &ppath, const Vector &rte_pos2, const Numeric &rte_alonglos_v, const Tensor3 &surface_props_data, const Tensor7 &cloudbox_field, const Vector &za_grid, const Index &Naa, const Index &t_interp_order, const Verbosity &verbosity)
WORKSPACE METHOD: iyEmissionHybrid.
const Index GFIELD4_LON_GRID
void nlinspace(Vector &x, const Numeric start, const Numeric stop, const Index n)
nlinspace
void transform(VectorView y, double(&my_func)(double), ConstVectorView x)
A generic transform function for vectors, which can be used to implement mathematical functions opera...
NUMERIC Numeric
The type to use for all floating point numbers.
INDEX Index
The type to use for all integer numbers and indices.
ConstComplexMatrixView transpose(ConstComplexMatrixView m)
Const version of transpose.
void swap(ComplexVector &v1, ComplexVector &v2)
Swaps two objects.
Declarations having to do with the four output streams.
Index nelem(const Lines &l)
Number of lines.
Particulates ENUMCLASS(Line, char, VMR, Strength, Center, ShapeG0X0, ShapeG0X1, ShapeG0X2, ShapeG0X3, ShapeD0X0, ShapeD0X1, ShapeD0X2, ShapeD0X3, ShapeG2X0, ShapeG2X1, ShapeG2X2, ShapeG2X3, ShapeD2X0, ShapeD2X1, ShapeD2X2, ShapeD2X3, ShapeFVCX0, ShapeFVCX1, ShapeFVCX2, ShapeFVCX3, ShapeETAX0, ShapeETAX1, ShapeETAX2, ShapeETAX3, ShapeYX0, ShapeYX1, ShapeYX2, ShapeYX3, ShapeGX0, ShapeGX1, ShapeGX2, ShapeGX3, ShapeDVX0, ShapeDVX1, ShapeDVX2, ShapeDVX3, ECS_SCALINGX0, ECS_SCALINGX1, ECS_SCALINGX2, ECS_SCALINGX3, ECS_BETAX0, ECS_BETAX1, ECS_BETAX2, ECS_BETAX3, ECS_LAMBDAX0, ECS_LAMBDAX1, ECS_LAMBDAX2, ECS_LAMBDAX3, ECS_DCX0, ECS_DCX1, ECS_DCX2, ECS_DCX3, NLTE) static_assert(Index(Line ScatteringString
constexpr bool isnan(double d) noexcept
invlib::Vector< ArtsVector > Vector
invlib wrapper type for ARTS vectors.
This file contains declerations of functions of physical character.
void error_if_limb_ppath(const Ppath &ppath)
Throws an error if ppath altitudes not are strictly increasing or decreasing.
Index ppath_what_background(const Ppath &ppath)
Returns the case number for the radiative background.
Propagation path structure and functions.
Array< PropagationMatrix > ArrayOfPropagationMatrix
Numeric pow(const Rational base, Numeric exp)
Power of.
Numeric sqrt(const Rational r)
Square root.
void get_stepwise_scattersky_propmat(StokesVector &ap, PropagationMatrix &Kp, ArrayOfStokesVector &dap_dx, ArrayOfPropagationMatrix &dKp_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const ConstMatrixView &ppath_1p_pnd, const ArrayOfMatrix &ppath_dpnd_dx, const Index ppath_1p_id, const ArrayOfArrayOfSingleScatteringData &scat_data, const ConstVectorView &ppath_line_of_sight, const ConstVectorView &ppath_temperature, const Index &atmosphere_dim, const bool &jacobian_do)
Computes the contribution by scattering at propagation path point.
void get_stepwise_clearsky_propmat(Workspace &ws, PropagationMatrix &K, StokesVector &S, Index <e, ArrayOfPropagationMatrix &dK_dx, ArrayOfStokesVector &dS_dx, const Agenda &propmat_clearsky_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const ConstVectorView &ppath_f_grid, const ConstVectorView &ppath_magnetic_field, const ConstVectorView &ppath_line_of_sight, const EnergyLevelMap &ppath_nlte, const ConstVectorView &ppath_vmrs, const Numeric &ppath_temperature, const Numeric &ppath_pressure, const bool &jacobian_do)
Gets the clearsky propgation matrix and NLTE contributions.
void get_iy_of_background(Workspace &ws, Matrix &iy, ArrayOfTensor3 &diy_dx, const ConstTensor3View &iy_transmittance, const Index &iy_id, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Ppath &ppath, const ConstVectorView &rte_pos2, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const ConstVectorView &f_grid, const String &iy_unit, const ConstTensor3View &surface_props_data, const Agenda &iy_main_agenda, const Agenda &iy_space_agenda, const Agenda &iy_surface_agenda, const Agenda &iy_cloudbox_agenda, const Index &iy_agenda_call1, const Verbosity &verbosity)
Determines iy of the "background" of a propgation path.
void get_stepwise_scattersky_source(StokesVector &Sp, ArrayOfStokesVector &dSp_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const ConstVectorView &ppath_1p_pnd, const ArrayOfMatrix &ppath_dpnd_dx, const Index ppath_1p_id, const ArrayOfArrayOfSingleScatteringData &scat_data, const ConstTensor7View &cloudbox_field, const ConstVectorView &za_grid, const ConstVectorView &aa_grid, const ConstMatrixView &ppath_line_of_sight, const GridPos &ppath_pressure, const Vector &temperature, const Index &atmosphere_dim, const bool &jacobian_do, const Index &t_interp_order)
Calculates the stepwise scattering source terms.
void get_ppath_atmvars(Vector &ppath_p, Vector &ppath_t, EnergyLevelMap &ppath_nlte, Matrix &ppath_vmr, Matrix &ppath_wind, Matrix &ppath_mag, const Ppath &ppath, const Index &atmosphere_dim, const ConstVectorView &p_grid, const ConstTensor3View &t_field, const EnergyLevelMap &nlte_field, const ConstTensor4View &vmr_field, const ConstTensor3View &wind_u_field, const ConstTensor3View &wind_v_field, const ConstTensor3View &wind_w_field, const ConstTensor3View &mag_u_field, const ConstTensor3View &mag_v_field, const ConstTensor3View &mag_w_field)
Determines pressure, temperature, VMR, winds and magnetic field for each propgataion path point.
void apply_iy_unit2(Tensor3View J, const ConstMatrixView &iy, const String &iy_unit, const ConstVectorView &f_grid, const Numeric &n, const ArrayOfIndex &i_pol)
Largely as apply_iy_unit but operates on jacobian data.
void yCalc_mblock_loop_body(bool &failed, String &fail_msg, ArrayOfArrayOfVector &iyb_aux_array, Workspace &ws, Vector &y, Vector &y_f, ArrayOfIndex &y_pol, Matrix &y_pos, Matrix &y_los, Matrix &y_geo, Matrix &jacobian, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &stokes_dim, const Vector &f_grid, const Matrix &sensor_pos, const Matrix &sensor_los, const Matrix &transmitter_pos, const Matrix &mblock_dlos_grid, const Sparse &sensor_response, const Vector &sensor_response_f, const ArrayOfIndex &sensor_response_pol, const Matrix &sensor_response_dlos, const String &iy_unit, const Agenda &iy_main_agenda, const Agenda &geo_pos_agenda, const Agenda &jacobian_agenda, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfArrayOfIndex &jacobian_indices, const ArrayOfString &iy_aux_vars, const Verbosity &verbosity, const Index &mblock_index, const Index &n1y, const Index &j_analytical_do)
Performs calculations for one measurement block, on y-level.
void iy_transmittance_mult(Tensor3 &iy_trans_total, const ConstTensor3View &iy_trans_old, const ConstTensor3View &iy_trans_new)
Multiplicates iy_transmittance with transmissions.
void apply_iy_unit(MatrixView iy, const String &iy_unit, const ConstVectorView &f_grid, const Numeric &n, const ArrayOfIndex &i_pol)
Performs conversion from radiance to other units, as well as applies refractive index to fulfill the ...
void get_ppath_f(Matrix &ppath_f, const Ppath &ppath, const ConstVectorView &f_grid, const Index &atmosphere_dim, const Numeric &rte_alonglos_v, const ConstMatrixView &ppath_wind)
Determines the Doppler shifted frequencies along the propagation path.
void rtmethods_unit_conversion(Matrix &iy, ArrayOfTensor3 &diy_dx, Tensor3 &ppvar_iy, const Index &ns, const Index &np, const Vector &f_grid, const Ppath &ppath, const ArrayOfRetrievalQuantity &jacobian_quantities, const Index &j_analytical_do, const String &iy_unit)
This function handles the unit conversion to be done at the end of some radiative transfer WSMs.
void rtmethods_jacobian_finalisation(Workspace &ws, ArrayOfTensor3 &diy_dx, ArrayOfTensor3 &diy_dpath, const Index &ns, const Index &nf, const Index &np, const Index &atmosphere_dim, const Ppath &ppath, const Vector &ppvar_p, const Vector &ppvar_t, const Matrix &ppvar_vmr, const Index &iy_agenda_call1, const Tensor3 &iy_transmittance, const Agenda &water_p_eq_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfIndex &jac_species_i)
This function fixes the last steps to made on the Jacobian in some radiative transfer WSMs.
Range get_rowindex_for_mblock(const Sparse &sensor_response, const Index &mblock_index)
Returns the "range" of y corresponding to a measurement block.
void adapt_stepwise_partial_derivatives(ArrayOfPropagationMatrix &dK_dx, ArrayOfStokesVector &dS_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const ConstVectorView &ppath_f_grid, const ConstVectorView &ppath_line_of_sight, const Index <e, const Index &atmosphere_dim, const bool &jacobian_do)
Adapts clearsky partial derivatives.
void get_stepwise_blackbody_radiation(VectorView B, VectorView dB_dT, const ConstVectorView &ppath_f_grid, const Numeric &ppath_temperature, const bool &do_temperature_derivative)
Get the blackbody radiation at propagation path point.
void get_ppath_cloudvars(ArrayOfIndex &clear2cloudy, Matrix &ppath_pnd, ArrayOfMatrix &ppath_dpnd_dx, const Ppath &ppath, const Index &atmosphere_dim, const ArrayOfIndex &cloudbox_limits, const Tensor4 &pnd_field, const ArrayOfTensor4 &dpnd_field_dx)
Determines the particle fields along a propagation path.
Declaration of functions in rte.cc.
void interp_atmfield_gp2itw(Matrix &itw, const Index &atmosphere_dim, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon)
Converts atmospheric grid positions to weights for interpolation of an atmospheric field.
void itw2p(VectorView p_values, ConstVectorView p_grid, const ArrayOfGridPos &gp, ConstMatrixView itw)
Converts interpolation weights to pressures.
void interp_cloudfield_gp2itw(VectorView itw, GridPos &gp_p_out, GridPos &gp_lat_out, GridPos &gp_lon_out, const GridPos &gp_p_in, const GridPos &gp_lat_in, const GridPos &gp_lon_in, const Index &atmosphere_dim, const ArrayOfIndex &cloudbox_limits)
Converts atmospheric a grid position to weights for interpolation of a field defined ONLY inside the ...
void interp_atmfield_by_itw(VectorView x, const Index &atmosphere_dim, ConstTensor3View x_field, const ArrayOfGridPos &gp_p, const ArrayOfGridPos &gp_lat, const ArrayOfGridPos &gp_lon, ConstMatrixView itw)
Interpolates an atmospheric field with pre-calculated weights by interp_atmfield_gp2itw.
Header file for special_interp.cc.
An Antenna object used by MCGeneral.
void set_pencil_beam()
set_pencil_beam
The structure to describe a propagation path and releated quantities.
Matrix los
Line-of-sight at each ppath point.
ArrayOfGridPos gp_lon
Index position with respect to the longitude grid.
Index np
Number of points describing the ppath.
Matrix pos
The distance between start pos and the last position in pos.
ArrayOfGridPos gp_lat
Index position with respect to the latitude grid.
Numeric end_lstep
The distance between end pos and the first position in pos.
Vector lstep
The length between ppath points.
ArrayOfGridPos gp_p
Index position with respect to the pressure grid.
Radiation Vector for Stokes dimension 1-4.
Index nrows() const
Returns the number of rows.
Class to keep track of Transmission Matrices for Stokes Dim 1-4.
void stepwise_source(RadiationVector &J, ArrayOfRadiationVector &dJ, const PropagationMatrix &K, const StokesVector &a, const StokesVector &S, const ArrayOfPropagationMatrix &dK, const ArrayOfStokesVector &da, const ArrayOfStokesVector &dS, const ConstVectorView &B, const ConstVectorView &dB_dT, const ArrayOfRetrievalQuantity &jacobian_quantities, const bool &jacobian_do)
Set the stepwise source.
ArrayOfTransmissionMatrix cumulative_transmission(const ArrayOfTransmissionMatrix &T, const CumulativeTransmission type)
Accumulate the transmission matrix over all layers.
void update_radiation_vector(RadiationVector &I, ArrayOfRadiationVector &dI1, ArrayOfRadiationVector &dI2, const RadiationVector &J1, const RadiationVector &J2, const ArrayOfRadiationVector &dJ1, const ArrayOfRadiationVector &dJ2, const TransmissionMatrix &T, const TransmissionMatrix &PiT, const ArrayOfTransmissionMatrix &dT1, const ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric r, const Vector &dr1, const Vector &dr2, const Index ia, const Index iz, const RadiativeTransferSolver solver)
Update the Radiation Vector.
void stepwise_transmission(TransmissionMatrix &T, ArrayOfTransmissionMatrix &dT1, ArrayOfTransmissionMatrix &dT2, const PropagationMatrix &K1, const PropagationMatrix &K2, const ArrayOfPropagationMatrix &dK1, const ArrayOfPropagationMatrix &dK2, const Numeric &r, const Numeric &dr_dtemp1, const Numeric &dr_dtemp2, const Index temp_deriv_pos)
Set the stepwise transmission matrix.
Stuff related to the transmission matrix.
Array< RadiationVector > ArrayOfRadiationVector
Array< TransmissionMatrix > ArrayOfTransmissionMatrix