ARTS  2.4.0(git:4fb77825)
m_cloudradar.cc
Go to the documentation of this file.
1 /* Copyright (C) 2012
2  Patrick Eriksson <Patrick.Eriksson@chalmers.se>
3  Stefan Buehler <sbuehler(at)ltu.se>
4 
5  This program is free software; you can redistribute it and/or modify it
6  under the terms of the GNU General Public License as published by the
7  Free Software Foundation; either version 2, or (at your option) any
8  later version.
9 
10  This program is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License
16  along with this program; if not, write to the Free Software
17  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
18  USA. */
19 
28 /*===========================================================================
29  === External declarations
30  ===========================================================================*/
31 
32 #include <cmath>
33 #include <stdexcept>
34 
35 #include "arts.h"
36 #include "auto_md.h"
37 #include "logic.h"
38 #include "messages.h"
39 #include "montecarlo.h"
40 #include "propagationmatrix.h"
41 #include "rte.h"
42 #include "sensor.h"
43 
44 extern const Numeric PI;
45 extern const Numeric SPEED_OF_LIGHT;
46 extern const String SCATSPECIES_MAINTAG;
47 
48 /* Workspace method: Doxygen documentation will be auto-generated */
50  Matrix& iy,
53  Vector& ppvar_p,
54  Vector& ppvar_t,
60  Matrix& ppvar_f,
62  const Index& stokes_dim,
63  const Vector& f_grid,
64  const Index& atmosphere_dim,
65  const Vector& p_grid,
66  const Tensor3& t_field,
68  const Tensor4& vmr_field,
70  const Tensor3& wind_u_field,
71  const Tensor3& wind_v_field,
72  const Tensor3& wind_w_field,
73  const Tensor3& mag_u_field,
74  const Tensor3& mag_v_field,
75  const Tensor3& mag_w_field,
76  const Index& cloudbox_on,
78  const Tensor4& pnd_field,
82  const Index& scat_data_checked,
84  const Index& jacobian_do,
86  const Ppath& ppath,
90  const Index& iy_agenda_call1,
91  const Tensor3& iy_transmission,
92  const Numeric& rte_alonglos_v,
93  const Index& trans_in_jacobian,
94  const Numeric& pext_scaling,
95  const Index& t_interp_order,
96  const Verbosity& verbosity _U_) {
97  // Some basic sizes
98  const Index nf = f_grid.nelem();
99  const Index ns = stokes_dim;
100  const Index np = ppath.np;
101  const Index ne = pnd_field.nbooks();
102 
103  // Radiative background index
104  const Index rbi = ppath_what_background(ppath);
105 
106  // Checks of input
107  // Throw error if unsupported features are requested
108  if (!iy_agenda_call1)
109  throw runtime_error(
110  "Recursive usage not possible (iy_agenda_call1 must be 1)");
111  if (!iy_transmission.empty())
112  throw runtime_error("*iy_transmission* must be empty");
113  if (rbi < 1 || rbi > 9)
114  throw runtime_error(
115  "ppath.background is invalid. Check your "
116  "calculation of *ppath*?");
117  if (rbi == 3 || rbi == 4)
118  throw runtime_error(
119  "The propagation path ends inside or at boundary of "
120  "the cloudbox.\nFor this method, *ppath* must be "
121  "calculated in this way:\n ppathCalc( cloudbox_on = 0 ).");
122  if (cloudbox_on) {
123  if (scat_data_checked != 1)
124  throw runtime_error(
125  "The scattering data must be flagged to have\n"
126  "passed a consistency check (scat_data_checked=1).");
127  if (ne != TotalNumberOfElements(scat_data))
128  throw runtime_error(
129  "*pnd_field* and *scat_data* inconsistent regarding total number of"
130  " scattering elements.");
131  }
132  if (jacobian_do) {
133  if (dpnd_field_dx.nelem() != jacobian_quantities.nelem())
134  throw runtime_error(
135  "*dpnd_field_dx* not properly initialized:\n"
136  "Number of elements in dpnd_field_dx must be equal number of jacobian"
137  " quantities.\n(Note: jacobians have to be defined BEFORE *pnd_field*"
138  " is calculated/set.");
139  }
140  // iy_aux_vars checked below
141  chk_if_in_range("pext_scaling", pext_scaling, 0, 2);
142 
143  // Transmitted signal
144  //
145  Matrix iy0;
146  //
148  iy0,
149  f_grid,
150  ppath.pos(np - 1, Range(0, atmosphere_dim)),
151  ppath.los(np - 1, joker),
153  if (iy0.ncols() != ns || iy0.nrows() != nf) {
154  ostringstream os;
155  os << "The size of *iy* returned from *iy_transmitter_agenda* is\n"
156  << "not correct:\n"
157  << " expected size = [" << nf << "," << stokes_dim << "]\n"
158  << " size of iy = [" << iy0.nrows() << "," << iy0.ncols() << "]\n";
159  throw runtime_error(os.str());
160  }
161  for (Index iv = 0; iv < nf; iv++) {
162  if (iy0(iv, 0) != 1)
163  throw runtime_error(
164  "The *iy* returned from *iy_transmitter_agenda* "
165  "must have the value 1 in the first column.");
166  }
167 
168  // Init Jacobian quantities
169  Index j_analytical_do = 0;
170  if (jacobian_do) {
171  FOR_ANALYTICAL_JACOBIANS_DO(j_analytical_do = 1;)
172  }
173  //
174  const Index nq = j_analytical_do ? jacobian_quantities.nelem() : 0;
175  ArrayOfTensor3 diy_dpath(nq);
176  ArrayOfIndex jac_species_i(nq), jac_scat_i(nq), jac_is_t(nq), jac_wind_i(nq);
177  ArrayOfIndex jac_mag_i(nq), jac_other(nq);
178 
179  if (j_analytical_do) {
180  rtmethods_jacobian_init(jac_species_i,
181  jac_scat_i,
182  jac_is_t,
183  jac_wind_i,
184  jac_mag_i,
185  jac_other,
186  diy_dx,
187  diy_dpath,
188  ns,
189  nf,
190  np,
191  nq,
192  abs_species,
193  cloudbox_on,
194  scat_species,
198  true);
199  }
200 
201  // Init iy_aux
202  const Index naux = iy_aux_vars.nelem();
203  iy_aux.resize(naux);
204  Index auxBackScat = -1;
205  Index auxOptDepth = -1;
206  Index auxPartAtte = -1;
207  //
208  for (Index i = 0; i < naux; i++) {
209  iy_aux[i].resize(nf * np, ns);
210  iy_aux[i] = 0;
211 
212  if (iy_aux_vars[i] == "Radiative background") {
213  iy_aux[i](joker, 0) = (Numeric)min((Index)2, rbi - 1);
214  } else if (iy_aux_vars[i] == "Backscattering") {
215  iy_aux[i] = 0;
216  auxBackScat = i;
217  } else if (iy_aux_vars[i] == "Optical depth") {
218  auxOptDepth = i;
219  } else if (iy_aux_vars[i] == "Particle extinction") {
220  iy_aux[i](Range(0, nf, np), joker) = 0;
221  auxPartAtte = i;
222  } else {
223  ostringstream os;
224  os << "The only allowed strings in *iy_aux_vars* are:\n"
225  << " \"Radiative background\"\n"
226  << " \"Backscattering\"\n"
227  << " \"Optical depth\"\n"
228  << " \"Particle extinction\"\n"
229  << "but you have selected: \"" << iy_aux_vars[i] << "\"";
230  throw runtime_error(os.str());
231  }
232  }
233 
234  // Get atmospheric and radiative variables along the propagation path
235  //
236  ppvar_trans_cumulat.resize(np, nf, ns, ns);
237  Tensor3 J(np, nf, ns);
238  Tensor4 trans_partial(np, nf, ns, ns);
239  Tensor5 dtrans_partial_dx_above(0, 0, 0, 0, 0);
240  Tensor5 dtrans_partial_dx_below(0, 0, 0, 0, 0);
241  Tensor5 Pe(ne, np, nf, ns, ns, 0);
242  ArrayOfMatrix ppvar_dpnd_dx(0);
243  ArrayOfIndex clear2cloudy;
244  Matrix scalar_ext(np, nf, 0); // Only used for iy_aux
245  Index nf_ssd = scat_data[0][0].pha_mat_data.nlibraries();
246  Index duplicate_freqs = ((nf == nf_ssd) ? 0 : 1);
247  Tensor6 pha_mat_1se(nf_ssd, 1, 1, 1, ns, ns);
248  Vector t_ok(1), t_array(1);
249  Matrix pdir(1, 2), idir(1, 2);
250  Index ptype;
251 
252  //
253  if (np == 1 && rbi == 1) // i.e. ppath is totally outside the atmosphere:
254  {
255  ppvar_p.resize(0);
256  ppvar_t.resize(0);
257  ppvar_vmr.resize(0, 0);
258  ppvar_wind.resize(0, 0);
259  ppvar_mag.resize(0, 0);
260  ppvar_pnd.resize(0, 0);
261  ppvar_f.resize(0, 0);
262  } else {
263  // Basic atmospheric variables
265  ppvar_t,
266  ppvar_nlte,
267  ppvar_vmr,
268  ppvar_wind,
269  ppvar_mag,
270  ppath,
272  p_grid,
273  t_field,
274  nlte_field,
275  vmr_field,
276  wind_u_field,
277  wind_v_field,
278  wind_w_field,
279  mag_u_field,
280  mag_v_field,
281  mag_w_field);
282 
283  get_ppath_f(
285 
286  // pnd_field
287  if (cloudbox_on) {
288  get_ppath_cloudvars(clear2cloudy,
289  ppvar_pnd,
290  ppvar_dpnd_dx,
291  ppath,
294  pnd_field,
295  dpnd_field_dx);
296  } else {
297  clear2cloudy.resize(np);
298  for (Index ip = 0; ip < np; ip++) {
299  clear2cloudy[ip] = -1;
300  }
301  }
302 
303  // Size radiative variables always used
304  PropagationMatrix K_this(nf, ns), K_past(nf, ns), Kp(nf, ns);
305  StokesVector a(nf, ns), S(nf, ns), Sp(nf, ns);
306  ArrayOfIndex lte(np);
307 
308  // Init variables only used if transmission part of jacobian
309  Vector dB_dT(0);
310  ArrayOfPropagationMatrix dK_this_dx(0), dK_past_dx(0), dKp_dx(0);
311  ArrayOfStokesVector da_dx(0), dS_dx(0), dSp_dx(0);
312  //
313  if (trans_in_jacobian && j_analytical_do) {
314  dtrans_partial_dx_above.resize(np, nq, nf, ns, ns);
315  dtrans_partial_dx_below.resize(np, nq, nf, ns, ns);
316  dtrans_partial_dx_above = 0; // Here all values are set to zero, to
317  dtrans_partial_dx_below = 0; // allow looping over non-defined values
318  dK_this_dx.resize(nq);
319  dK_past_dx.resize(nq);
320  dKp_dx.resize(nq);
321  da_dx.resize(nq);
322  dS_dx.resize(nq);
323  dSp_dx.resize(nq);
324  dB_dT.resize(nf);
325  FOR_ANALYTICAL_JACOBIANS_DO(dK_this_dx[iq] = PropagationMatrix(nf, ns);
326  dK_past_dx[iq] = PropagationMatrix(nf, ns);
327  dKp_dx[iq] = PropagationMatrix(nf, ns);
328  da_dx[iq] = StokesVector(nf, ns);
329  dS_dx[iq] = StokesVector(nf, ns);
330  dSp_dx[iq] = StokesVector(nf, ns);)
331  }
332 
333  // Loop ppath points and determine radiative properties
334  for (Index ip = 0; ip < np; ip++) {
336  K_this,
337  S,
338  lte[ip],
339  dK_this_dx,
340  dS_dx,
343  ppvar_f(joker, ip),
344  ppvar_mag(joker, ip),
345  ppath.los(ip, joker),
346  ppvar_nlte[ip],
347  ppvar_vmr(joker, ip),
348  ppvar_t[ip],
349  ppvar_p[ip],
350  jac_species_i,
351  trans_in_jacobian && j_analytical_do);
352 
353  if (trans_in_jacobian && j_analytical_do) {
355  dK_this_dx,
356  dS_dx,
358  ppvar_f(joker, ip),
359  ppath.los(ip, joker),
360  ppvar_vmr(joker, ip),
361  ppvar_t[ip],
362  ppvar_p[ip],
363  jac_species_i,
364  jac_wind_i,
365  lte[ip],
367  trans_in_jacobian && j_analytical_do);
368  }
369 
370  if (clear2cloudy[ip] + 1) {
372  Kp,
373  da_dx,
374  dKp_dx,
376  ppvar_pnd(joker, Range(ip, 1)),
377  ppvar_dpnd_dx,
378  ip,
379  scat_data,
380  ppath.los(ip, joker),
381  ppvar_t[Range(ip, 1)],
383  trans_in_jacobian && jacobian_do);
384 
385  if (abs(pext_scaling - 1) > 1e-6) {
386  Kp *= pext_scaling;
387  if (trans_in_jacobian && j_analytical_do) {
388  FOR_ANALYTICAL_JACOBIANS_DO(dKp_dx[iq] *= pext_scaling;)
389  }
390  }
391 
392  K_this += Kp;
393 
394  if (auxPartAtte >= 0) {
395  scalar_ext(ip, joker) = Kp.Kjj();
396  }
397 
398  if (trans_in_jacobian && j_analytical_do) {
399  FOR_ANALYTICAL_JACOBIANS_DO(dK_this_dx[iq] += dKp_dx[iq];)
400  }
401 
402  // Get back-scattering per particle, where relevant
403  {
404  // Direction of outgoing scattered radiation (which is reverse to LOS).
405  Vector los_sca;
406  mirror_los(los_sca, ppath.los(ip, joker), atmosphere_dim);
407  pdir(0, joker) = los_sca;
408 
409  // Obtain a length-2 vector for incoming direction
410  Vector los_inc;
411  if (atmosphere_dim == 3) {
412  los_inc = ppath.los(ip, joker);
413  } else // Mirror back to get a correct 3D LOS
414  {
415  mirror_los(los_inc, los_sca, 3);
416  }
417  idir(0, joker) = los_inc;
418 
419  Index i_se_flat = 0;
420  for (Index i_ss = 0; i_ss < scat_data.nelem(); i_ss++)
421  for (Index i_se = 0; i_se < scat_data[i_ss].nelem(); i_se++) {
422  // determine whether we have some valid pnd for this
423  // scatelem (in pnd or dpnd)
424  Index val_pnd = 0;
425  if (ppvar_pnd(i_se_flat, ip) != 0)
426  val_pnd = 1;
427  else if (j_analytical_do)
428  for (Index iq = 0; iq < nq && !val_pnd; iq++)
429  if (jac_scat_i[iq] >= 0)
430  if (ppvar_dpnd_dx[iq](i_se_flat, ip) != 0) {
431  val_pnd = 1;
432  break;
433  }
434  if (val_pnd) {
435  pha_mat_1ScatElem(pha_mat_1se,
436  ptype,
437  t_ok,
438  scat_data[i_ss][i_se],
439  ppvar_t[Range(ip, 1)],
440  pdir,
441  idir,
442  0,
443  t_interp_order);
444  if (t_ok[0] != 0)
445  if (duplicate_freqs)
446  for (Index iv = 0; iv < nf; iv++)
447  Pe(i_se_flat, ip, iv, joker, joker) =
448  pha_mat_1se(0, 0, 0, 0, joker, joker);
449  else
450  Pe(i_se_flat, ip, joker, joker, joker) =
451  pha_mat_1se(joker, 0, 0, 0, joker, joker);
452  else {
453  ostringstream os;
454  os << "Interpolation error for (flat-array) scattering"
455  << " element #" << i_se_flat << "\n"
456  << "at location/temperature point #" << ip << "\n";
457  throw runtime_error(os.str());
458  }
459  }
460  i_se_flat++;
461  }
462 
463  } // local scope
464  } // clear2cloudy
465 
466  if (trans_in_jacobian && j_analytical_do) {
469  trans_partial(ip, joker, joker, joker),
470  dtrans_partial_dx_above(ip, joker, joker, joker, joker),
471  dtrans_partial_dx_below(ip, joker, joker, joker, joker),
472  (ip > 0) ? ppvar_trans_cumulat(ip - 1, joker, joker, joker)
473  : Tensor3(0, 0, 0),
474  K_past,
475  K_this,
476  dK_past_dx,
477  dK_this_dx,
478  (ip > 0) ? ppath.lstep[ip - 1] : Numeric(1.0),
479  ip == 0);
480  } else {
483  trans_partial(ip, joker, joker, joker),
484  Tensor4(0, 0, 0, 0),
485  Tensor4(0, 0, 0, 0),
486  (ip > 0) ? ppvar_trans_cumulat(ip - 1, joker, joker, joker)
487  : Tensor3(0, 0, 0),
488  K_past,
489  K_this,
490  dK_past_dx,
491  dK_this_dx,
492  (ip > 0) ? ppath.lstep[ip - 1] : Numeric(1.0),
493  ip == 0);
494  }
495 
496  swap(K_past, K_this);
497  swap(dK_past_dx, dK_this_dx);
498  }
499  }
500 
501  // Transmission for reversed direction
502  //Tensor3 tr_rev( nf, ns, ns );
503  //for( Index iv=0; iv<nf; iv++ )
504  // { id_mat( tr_rev(iv,joker,joker) ); }
505 
506  // Size iy and set to zero
507  iy.resize(nf * np, ns);
508  iy = 0;
509 
510  // If neither cloudbox or aux variables, nothing more to do
511  if (!cloudbox_on && !naux) {
512  return;
513  }
514 
515  // Radiative transfer calculations
516  for (Index ip = 0; ip < np; ip++) {
517  // Radar return only possible if inside cloudbox
518  if (clear2cloudy[ip] >= 0) {
519  for (Index iv = 0; iv < nf; iv++) {
520  // Calculate bulk back-scattering
521  Matrix P(ns, ns, 0);
522  for (Index i = 0; i < ne; i++) {
523  if (ppvar_pnd(i, ip) != 0) {
524  for (Index is1 = 0; is1 < ns; is1++) {
525  for (Index is2 = 0; is2 < ns; is2++) {
526  P(is1, is2) += ppvar_pnd(i, ip) * Pe(i, ip, iv, is1, is2);
527  }
528  }
529  }
530  }
531 
532  // Combine iy0, double transmission and scattering matrix
533  Vector iy1(ns), iy2(ns);
534  const Index iout = iv * np + ip;
535  mult(iy1, ppvar_trans_cumulat(ip, iv, joker, joker), iy0(iv, joker));
536  mult(iy2, P, iy1);
537  mult(iy(iout, joker), ppvar_trans_cumulat(ip, iv, joker, joker), iy2);
538 
539  if (auxBackScat >= 0) {
540  mult(iy_aux[auxBackScat](iout, joker), P, iy0(iv, joker));
541  }
542  if (auxPartAtte >= 0 && ip > 0) {
543  iy_aux[auxPartAtte](iout, 0) =
544  iy_aux[auxPartAtte](iout - 1, 0) +
545  ppath.lstep[ip - 1] *
546  (scalar_ext(ip - 1, iv) + scalar_ext(ip - 1, iv));
547  }
548 
549  // Jacobians
550  if (j_analytical_do) {
551  // All this is to cover impact on back-scattering of
552  // scattering species
553  for (Index iq = 0; iq < nq; iq++) {
554  if (jacobian_quantities[iq].Analytical()) {
555  if (jac_scat_i[iq] >= 0) {
556  // Change of scattering scattering matrix
557  P = 0.0;
558  for (Index i = 0; i < ne; i++) {
559  if (ppvar_dpnd_dx[iq](i, ip) != 0) {
560  for (Index is1 = 0; is1 < ns; is1++) {
561  for (Index is2 = 0; is2 < ns; is2++) {
562  P(is1, is2) +=
563  ppvar_dpnd_dx[iq](i, ip) * Pe(i, ip, iv, is1, is2);
564  }
565  }
566  }
567  }
568 
569  // Apply transmissions as above
570  Vector iy_tmp(ns);
571  mult(iy_tmp, P, iy1);
572  mult(diy_dpath[iq](ip, iout, joker),
573  ppvar_trans_cumulat(ip, iv, joker, joker),
574  iy_tmp);
575  }
576  }
577  }
578  } // j_analytical_do
579  } // for iv
580  } // if cloudy
581  else {
582  if (auxPartAtte >= 0 && ip > 0) {
583  for (Index iv = 0; iv < nf; iv++) {
584  const Index iout = iv * np + ip;
585  iy_aux[auxPartAtte](iout, 0) = iy_aux[auxPartAtte](iout - nf, 0);
586  }
587  }
588  }
589 
590  if (auxOptDepth >= 0) {
591  for (Index iv = 0; iv < nf; iv++) {
592  const Index iout = iv * np + ip;
593  iy_aux[auxOptDepth](iout, 0) =
594  -2 * log(ppvar_trans_cumulat(ip, iv, 0, 0));
595  }
596  }
597 
598  // Update tr_rev
599  /*
600  if( ip<np-1 )
601  {
602  for( Index iv=0; iv<nf; iv++ )
603  {
604  Matrix tmp = tr_rev(iv,joker,joker);
605  mult( tr_rev(iv,joker,joker), trans_partial(ip+1,iv,joker,joker),
606  tmp );
607  }
608  }
609  */
610  } // for ip
611 
612  // Finalize analytical Jacobian
613  if (j_analytical_do) {
614  // Add impact on transmission, if selected to be included
615  if (trans_in_jacobian) {
616  Vector jterm(stokes_dim);
617  FOR_ANALYTICAL_JACOBIANS_DO(for (Index iv = 0; iv < nf; iv++) {
618  for (Index ip = 0; ip < np; ip++) // Index of x-point
619  {
620  for (Index j = ip; j < np; j++) // Index of back-scattering
621  {
622  // Nothing to do for j==0
623  if (j == 0) {
624  continue;
625  }
626 
627  const Index iout = iv * np + j;
628  // Impact on transmission towards the sensor
629  if (ip > 0) {
630  mult(jterm,
631  dtrans_partial_dx_above(ip, iq, iv, joker, joker),
632  iy(iout, joker));
633  Vector jnew(ns);
634  solve(jnew, trans_partial(ip, iv, joker, joker), jterm);
635  for (Index i = 0; i < jnew.nelem(); ++i) {
636  if (isnan(jnew[i])) jnew[i] = 0.0;
637  }
638  jnew *= 2;
639  diy_dpath[iq](ip, iout, joker) += jnew;
640  }
641  // Impact on transmission away from the sensor
642  if (j > ip && j < np - 1) {
643  mult(jterm,
644  dtrans_partial_dx_below(ip, iq, iv, joker, joker),
645  iy(iout, joker));
646  Vector jnew(ns);
647  solve(jnew, trans_partial(ip, iv, joker, joker), jterm);
648  for (Index i = 0; i < jnew.nelem(); ++i) {
649  if (isnan(jnew[i])) jnew[i] = 0.0;
650  }
651  jnew *= 2;
652  diy_dpath[iq](ip, iout, joker) += jnew;
653  }
654  }
655  }
656  })
657  }
658 
660  diy_dx,
661  diy_dpath,
662  ns,
663  nf,
664  np,
666  ppath,
667  ppvar_p,
668  ppvar_t,
669  ppvar_vmr,
674  jac_species_i,
675  jac_is_t);
676  }
677 }
678 
679 /* Workspace method: Doxygen documentation will be auto-generated */
681  Matrix& iy,
684  Vector& ppvar_p,
685  Vector& ppvar_t,
687  Matrix& ppvar_vmr,
689  Matrix& ppvar_mag,
690  Matrix& ppvar_pnd,
691  Matrix& ppvar_f,
693  const Index& stokes_dim,
694  const Vector& f_grid,
695  const Index& atmosphere_dim,
696  const Vector& p_grid,
697  const Tensor3& t_field,
698  const EnergyLevelMap& nlte_field,
699  const Tensor4& vmr_field,
701  const Tensor3& wind_u_field,
702  const Tensor3& wind_v_field,
703  const Tensor3& wind_w_field,
704  const Tensor3& mag_u_field,
705  const Tensor3& mag_v_field,
706  const Tensor3& mag_w_field,
707  const Index& cloudbox_on,
709  const Tensor4& pnd_field,
713  const Index& scat_data_checked,
714  const ArrayOfString& iy_aux_vars,
715  const Index& jacobian_do,
717  const Ppath& ppath,
719  const Agenda& water_p_eq_agenda,
721  const Index& iy_agenda_call1,
722  const Tensor3& iy_transmission,
723  const Numeric& rte_alonglos_v,
724  const Index& trans_in_jacobian,
725  const Numeric& pext_scaling,
726  const Index& t_interp_order,
727  const Verbosity& verbosity _U_) {
728  // Some basic sizes
729  const Index nf = f_grid.nelem();
730  const Index ns = stokes_dim;
731  const Index np = ppath.np;
732  const Index ne = pnd_field.nbooks();
733 
734  // Radiative background index
735  const Index rbi = ppath_what_background(ppath);
736 
737  // Checks of input
738  // Throw error if unsupported features are requested
739  if (!iy_agenda_call1)
740  throw runtime_error(
741  "Recursive usage not possible (iy_agenda_call1 must be 1)");
742  if (!iy_transmission.empty())
743  throw runtime_error("*iy_transmission* must be empty");
744  if (rbi < 1 || rbi > 9)
745  throw runtime_error(
746  "ppath.background is invalid. Check your "
747  "calculation of *ppath*?");
748  if (rbi == 3 || rbi == 4)
749  throw runtime_error(
750  "The propagation path ends inside or at boundary of "
751  "the cloudbox.\nFor this method, *ppath* must be "
752  "calculated in this way:\n ppathCalc( cloudbox_on = 0 ).");
753  if (cloudbox_on) {
754  if (scat_data_checked != 1)
755  throw runtime_error(
756  "The scattering data must be flagged to have\n"
757  "passed a consistency check (scat_data_checked=1).");
758  if (ne != TotalNumberOfElements(scat_data))
759  throw runtime_error(
760  "*pnd_field* and *scat_data* inconsistent regarding total number of"
761  " scattering elements.");
762  }
763  if (jacobian_do)
764  if (dpnd_field_dx.nelem() != jacobian_quantities.nelem())
765  throw runtime_error(
766  "*dpnd_field_dx* not properly initialized:\n"
767  "Number of elements in dpnd_field_dx must be equal number of jacobian"
768  " quantities.\n(Note: jacobians have to be defined BEFORE *pnd_field*"
769  " is calculated/set.");
770  // iy_aux_vars checked below
771  chk_if_in_range("pext_scaling", pext_scaling, 0, 2);
772 
773  // Transmitted signal
774  //
775  Matrix iy0;
776  //
778  iy0,
779  f_grid,
780  ppath.pos(np - 1, Range(0, atmosphere_dim)),
781  ppath.los(np - 1, joker),
783 
784  if (iy0.ncols() != ns || iy0.nrows() != nf) {
785  ostringstream os;
786  os << "The size of *iy* returned from *iy_transmitter_agenda* is\n"
787  << "not correct:\n"
788  << " expected size = [" << nf << "," << stokes_dim << "]\n"
789  << " size of iy = [" << iy0.nrows() << "," << iy0.ncols() << "]\n";
790  throw runtime_error(os.str());
791  }
792  for (Index iv = 0; iv < nf; iv++)
793  if (iy0(iv, 0) != 1)
794  throw runtime_error(
795  "The *iy* returned from *iy_transmitter_agenda* "
796  "must have the value 1 in the first column.");
797 
798  // Init Jacobian quantities
799  Index j_analytical_do = 0;
800  if (jacobian_do) FOR_ANALYTICAL_JACOBIANS_DO(j_analytical_do = 1;);
801  //
802  const Index nq = j_analytical_do ? jacobian_quantities.nelem() : 0;
803  ArrayOfTensor3 diy_dpath(nq);
804  ArrayOfIndex jac_species_i(nq), jac_scat_i(nq), jac_is_t(nq), jac_wind_i(nq);
805  ArrayOfIndex jac_mag_i(nq), jac_other(nq);
806 
807  if (j_analytical_do)
808  rtmethods_jacobian_init(jac_species_i,
809  jac_scat_i,
810  jac_is_t,
811  jac_wind_i,
812  jac_mag_i,
813  jac_other,
814  diy_dx,
815  diy_dpath,
816  ns,
817  nf,
818  np,
819  nq,
820  abs_species,
821  cloudbox_on,
822  scat_species,
826  true);
827 
828  // Init iy_aux
829  const Index naux = iy_aux_vars.nelem();
830  iy_aux.resize(naux);
831 
832  // Not implemented in this function yet
833  // Index auxBackScat = -1;
834  // Index auxOptDepth = -1;
835  // Index auxPartAtte = -1;
836  //
837  for (Index i = 0; i < naux; i++) {
838  iy_aux[i].resize(nf * np, ns);
839  iy_aux[i] = 0;
840 
841  if (iy_aux_vars[i] == "Radiative background")
842  iy_aux[i](joker, 0) = (Numeric)min((Index)2, rbi - 1);
843  // else if( iy_aux_vars[i] == "Backscattering" ) {
844  // iy_aux[i] = 0;
845  // auxBackScat = i;
846  // }
847  // else if( iy_aux_vars[i] == "Optical depth" )
848  // auxOptDepth = i;
849  // else if( iy_aux_vars[i] == "Particle extinction" ) {
850  // iy_aux[i](Range(0,nf,np),joker) = 0;
851  // auxPartAtte = i;
852  // }
853  else {
854  ostringstream os;
855  os << "The only allowed strings in *iy_aux_vars* are:\n"
856  << " \"Radiative background\"\n"
857  // << " \"Backscattering\"\n"
858  // << " \"Optical depth\"\n"
859  // << " \"Particle extinction\"\n"
860  << "but you have selected: \"" << iy_aux_vars[i] << "\"";
861  throw runtime_error(os.str());
862  }
863  }
864 
865  // Get atmospheric and radiative variables along the propagation path
866  ppvar_trans_cumulat.resize(np, nf, ns, ns);
867 
868  ArrayOfRadiationVector lvl_rad(np, RadiationVector(nf, ns));
870  np,
872  np, ArrayOfRadiationVector(nq, RadiationVector(nf, ns))));
873 
875  ArrayOfArrayOfTransmissionMatrix dlyr_tra_above(
877  ArrayOfArrayOfTransmissionMatrix dlyr_tra_below(
879 
880  // TEMPORARY VARIABLE, THIS COULD BE AN ArrayOfArrayOfPropagationMatrix most likely, to speedup calculations
881  Tensor5 Pe(ne, np, nf, ns, ns, 0);
882 
883  ArrayOfMatrix ppvar_dpnd_dx(0);
884  ArrayOfIndex clear2cloudy;
885  // Matrix scalar_ext(np,nf,0); // Only used for iy_aux
886  Index nf_ssd = scat_data[0][0].pha_mat_data.nlibraries();
887  Index duplicate_freqs = ((nf == nf_ssd) ? 0 : 1);
888  Tensor6 pha_mat_1se(nf_ssd, 1, 1, 1, ns, ns);
889  Vector t_ok(1), t_array(1);
890  Matrix pdir(1, 2), idir(1, 2);
891  Index ptype;
892 
893  if (np == 1 && rbi == 1) { // i.e. ppath is totally outside the atmosphere:
894  ppvar_p.resize(0);
895  ppvar_t.resize(0);
896  ppvar_vmr.resize(0, 0);
897  ppvar_wind.resize(0, 0);
898  ppvar_mag.resize(0, 0);
899  ppvar_pnd.resize(0, 0);
900  ppvar_f.resize(0, 0);
901  } else {
902  // Basic atmospheric variables
904  ppvar_t,
905  ppvar_nlte,
906  ppvar_vmr,
907  ppvar_wind,
908  ppvar_mag,
909  ppath,
911  p_grid,
912  t_field,
913  nlte_field,
914  vmr_field,
915  wind_u_field,
916  wind_v_field,
917  wind_w_field,
918  mag_u_field,
919  mag_v_field,
920  mag_w_field);
921 
922  get_ppath_f(
924 
925  // pnd_field
926  if (cloudbox_on)
927  get_ppath_cloudvars(clear2cloudy,
928  ppvar_pnd,
929  ppvar_dpnd_dx,
930  ppath,
933  pnd_field,
934  dpnd_field_dx);
935  else {
936  clear2cloudy.resize(np);
937  for (Index ip = 0; ip < np; ip++) clear2cloudy[ip] = -1;
938  }
939 
940  // Size radiative variables always used
941  PropagationMatrix K_this(nf, ns), K_past(nf, ns), Kp(nf, ns);
942  StokesVector a(nf, ns), S(nf, ns);
943  ArrayOfIndex lte(np);
944 
945  // Init variables only used if transmission part of jacobian
946  Vector dB_dT(0);
947  ArrayOfPropagationMatrix dK_this_dx(0), dK_past_dx(0), dKp_dx(0);
948  ArrayOfStokesVector da_dx(0), dS_dx(0);
949 
950  // HSE variables
951  Index temperature_derivative_position = -1;
952  bool do_hse = false;
953 
954  if (trans_in_jacobian && j_analytical_do) {
955  dK_this_dx.resize(nq);
956  dK_past_dx.resize(nq);
957  dKp_dx.resize(nq);
958  da_dx.resize(nq);
959  dS_dx.resize(nq);
960  dB_dT.resize(nf);
962  dK_this_dx[iq] = PropagationMatrix(nf, ns);
963  dK_past_dx[iq] = PropagationMatrix(nf, ns);
964  dKp_dx[iq] = PropagationMatrix(nf, ns);
965  da_dx[iq] = StokesVector(nf, ns);
966  dS_dx[iq] = StokesVector(nf, ns);
967  if (jacobian_quantities[iq] == JacPropMatType::Temperature) {
968  temperature_derivative_position = iq;
969  do_hse = jacobian_quantities[iq].Subtag() == "HSE on";
970  })
971  }
972 
973  // Loop ppath points and determine radiative properties
974  for (Index ip = 0; ip < np; ip++) {
976  K_this,
977  S,
978  lte[ip],
979  dK_this_dx,
980  dS_dx,
983  ppvar_f(joker, ip),
984  ppvar_mag(joker, ip),
985  ppath.los(ip, joker),
986  ppvar_nlte[ip],
987  ppvar_vmr(joker, ip),
988  ppvar_t[ip],
989  ppvar_p[ip],
990  jac_species_i,
991  trans_in_jacobian && j_analytical_do);
992 
993  if (trans_in_jacobian && j_analytical_do)
995  dK_this_dx,
996  dS_dx,
998  ppvar_f(joker, ip),
999  ppath.los(ip, joker),
1000  ppvar_vmr(joker, ip),
1001  ppvar_t[ip],
1002  ppvar_p[ip],
1003  jac_species_i,
1004  jac_wind_i,
1005  lte[ip],
1007  trans_in_jacobian && j_analytical_do);
1008 
1009  if (clear2cloudy[ip] + 1) {
1011  Kp,
1012  da_dx,
1013  dKp_dx,
1015  ppvar_pnd(joker, Range(ip, 1)),
1016  ppvar_dpnd_dx,
1017  ip,
1018  scat_data,
1019  ppath.los(ip, joker),
1020  ppvar_t[Range(ip, 1)],
1022  trans_in_jacobian && jacobian_do);
1023 
1024  if (abs(pext_scaling - 1) > 1e-6) {
1025  Kp *= pext_scaling;
1026  if (trans_in_jacobian && j_analytical_do) {
1027  FOR_ANALYTICAL_JACOBIANS_DO(dKp_dx[iq] *= pext_scaling;)
1028  }
1029  }
1030 
1031  K_this += Kp;
1032 
1033  // if( auxPartAtte >= 0 )
1034  // scalar_ext(ip,joker) = Kp.Kjj();
1035 
1036  if (trans_in_jacobian && j_analytical_do)
1037  FOR_ANALYTICAL_JACOBIANS_DO(dK_this_dx[iq] += dKp_dx[iq];);
1038 
1039  // Get back-scattering per particle, where relevant
1040  {
1041  // Direction of outgoing scattered radiation (which is reverse to LOS).
1042  Vector los_sca;
1043  mirror_los(los_sca, ppath.los(ip, joker), atmosphere_dim);
1044  pdir(0, joker) = los_sca;
1045 
1046  // Obtain a length-2 vector for incoming direction
1047  Vector los_inc;
1048  if (atmosphere_dim == 3) {
1049  los_inc = ppath.los(ip, joker);
1050  } else // Mirror back to get a correct 3D LOS
1051  {
1052  mirror_los(los_inc, los_sca, 3);
1053  }
1054  idir(0, joker) = los_inc;
1055 
1056  Index i_se_flat = 0;
1057  for (Index i_ss = 0; i_ss < scat_data.nelem(); i_ss++)
1058  for (Index i_se = 0; i_se < scat_data[i_ss].nelem(); i_se++) {
1059  // determine whether we have some valid pnd for this
1060  // scatelem (in pnd or dpnd)
1061  Index val_pnd = 0;
1062  if (ppvar_pnd(i_se_flat, ip) != 0)
1063  val_pnd = 1;
1064  else if (j_analytical_do)
1065  for (Index iq = 0; iq < nq && !val_pnd; iq++)
1066  if (jac_scat_i[iq] >= 0)
1067  if (ppvar_dpnd_dx[iq](i_se_flat, ip) != 0) {
1068  val_pnd = 1;
1069  break;
1070  }
1071  if (val_pnd) {
1072  pha_mat_1ScatElem(pha_mat_1se,
1073  ptype,
1074  t_ok,
1075  scat_data[i_ss][i_se],
1076  ppvar_t[Range(ip, 1)],
1077  pdir,
1078  idir,
1079  0,
1080  t_interp_order);
1081  if (t_ok[0] not_eq 0)
1082  if (duplicate_freqs)
1083  for (Index iv = 0; iv < nf; iv++)
1084  Pe(i_se_flat, ip, iv, joker, joker) =
1085  pha_mat_1se(0, 0, 0, 0, joker, joker);
1086  else
1087  Pe(i_se_flat, ip, joker, joker, joker) =
1088  pha_mat_1se(joker, 0, 0, 0, joker, joker);
1089  else {
1090  ostringstream os;
1091  os << "Interpolation error for (flat-array) scattering"
1092  << " element #" << i_se_flat << "\n"
1093  << "at location/temperature point #" << ip << "\n";
1094  throw runtime_error(os.str());
1095  }
1096  }
1097  i_se_flat++;
1098  }
1099 
1100  } // local scope
1101  } // clear2cloudy
1102 
1103  if (ip not_eq 0) {
1104  const Numeric dr_dT_past =
1105  do_hse ? ppath.lstep[ip - 1] / (2.0 * ppvar_t[ip - 1]) : 0;
1106  const Numeric dr_dT_this =
1107  do_hse ? ppath.lstep[ip - 1] / (2.0 * ppvar_t[ip]) : 0;
1108  stepwise_transmission(lyr_tra[ip],
1109  dlyr_tra_above[ip],
1110  dlyr_tra_below[ip],
1111  K_past,
1112  K_this,
1113  dK_past_dx,
1114  dK_this_dx,
1115  ppath.lstep[ip - 1],
1116  dr_dT_past,
1117  dr_dT_this,
1118  temperature_derivative_position);
1119  }
1120 
1121  swap(K_past, K_this);
1122  swap(dK_past_dx, dK_this_dx);
1123  }
1124  }
1125 
1126  const ArrayOfTransmissionMatrix tot_tra_forward =
1128  const ArrayOfTransmissionMatrix tot_tra_reflect =
1129  cumulative_transmission(lyr_tra, CumulativeTransmission::Forward);
1130  const ArrayOfTransmissionMatrix reflect_matrix =
1132  const ArrayOfArrayOfTransmissionMatrix dreflect_matrix =
1133  cumulative_backscatter_derivative(Pe, ppvar_dpnd_dx);
1134 
1135  lvl_rad[0] = iy0;
1136  RadiationVector rad_inc = RadiationVector(nf, ns);
1137  rad_inc = iy0;
1139  dlvl_rad,
1140  rad_inc,
1141  lyr_tra,
1142  tot_tra_forward,
1143  tot_tra_reflect,
1144  reflect_matrix,
1145  dlyr_tra_above,
1146  dlyr_tra_below,
1147  dreflect_matrix,
1148  BackscatterSolver::CommutativeTransmission);
1149 
1150 
1151  // Size iy and set to zero
1152  iy.resize(nf * np, ns); // iv*np + ip is the desired output order...
1153  for (Index ip = 0; ip < np; ip++) {
1154  for (Index iv = 0; iv < nf; iv++) {
1155  for (Index is = 0; is < stokes_dim; is++) {
1156  iy(iv * np + ip, is) = lvl_rad[ip](iv, is);
1157  if (j_analytical_do) {
1158  FOR_ANALYTICAL_JACOBIANS_DO(for (Index ip2 = 0; ip2 < np; ip2++)
1159  diy_dpath[iq](ip, iv * np + ip2, is) =
1160  dlvl_rad[ip][ip2][iq](iv, is););
1161  }
1162  }
1163  }
1164  }
1165  // FIXME: Add the aux-variables back
1166 
1167 
1168  // Finalize analytical Jacobian
1169  if (j_analytical_do)
1171  diy_dx,
1172  diy_dpath,
1173  ns,
1174  nf,
1175  np,
1177  ppath,
1178  ppvar_p,
1179  ppvar_t,
1180  ppvar_vmr,
1185  jac_species_i,
1186  jac_is_t);
1187 }
1188 
1189 /* Workspace method: Doxygen documentation will be auto-generated */
1191  Vector& y,
1192  Vector& y_f,
1194  Matrix& y_pos,
1195  Matrix& y_los,
1197  Matrix& y_geo,
1198  Matrix& jacobian,
1199  const Index& atmgeom_checked,
1200  const Index& atmfields_checked,
1201  const String& iy_unit,
1202  const ArrayOfString& iy_aux_vars,
1203  const Index& stokes_dim,
1204  const Vector& f_grid,
1205  const Index& atmosphere_dim,
1206  const EnergyLevelMap& nlte_field,
1207  const Index& cloudbox_on,
1208  const Index& cloudbox_checked,
1209  const Matrix& sensor_pos,
1210  const Matrix& sensor_los,
1211  const Index& sensor_checked,
1212  const Index& jacobian_do,
1214  const Agenda& iy_main_agenda,
1215  const Agenda& geo_pos_agenda,
1217  const Vector& range_bins,
1218  const Numeric& ze_tref,
1219  const Numeric& k2,
1220  const Numeric& dbze_min,
1221  const Verbosity&) {
1222  // Important sizes
1223  const Index npos = sensor_pos.nrows();
1224  const Index nbins = range_bins.nelem() - 1;
1225  const Index nf = f_grid.nelem();
1226  const Index naux = iy_aux_vars.nelem();
1227 
1228  //---------------------------------------------------------------------------
1229  // Input checks
1230  //---------------------------------------------------------------------------
1231 
1232  // Basics
1233  //
1234  chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
1235  //
1236  if (f_grid.empty()) {
1237  throw runtime_error("The frequency grid is empty.");
1238  }
1239  chk_if_increasing("f_grid", f_grid);
1240  if (f_grid[0] <= 0) {
1241 
1242  throw runtime_error("All frequencies in *f_grid* must be > 0.");
1243  }
1244  //
1245  chk_if_in_range("stokes_dim", stokes_dim, 1, 4);
1246  if (atmfields_checked != 1)
1247  throw runtime_error(
1248  "The atmospheric fields must be flagged to have "
1249  "passed a consistency check (atmfields_checked=1).");
1250  if (atmgeom_checked != 1)
1251  throw runtime_error(
1252  "The atmospheric geometry must be flagged to have "
1253  "passed a consistency check (atmgeom_checked=1).");
1254  if (cloudbox_checked != 1)
1255  throw runtime_error(
1256  "The cloudbox must be flagged to have "
1257  "passed a consistency check (cloudbox_checked=1).");
1258  if (sensor_checked != 1)
1259  throw runtime_error(
1260  "The sensor variables must be flagged to have "
1261  "passed a consistency check (sensor_checked=1).");
1262 
1263  // Method specific variables
1264  bool is_z = max(range_bins) > 1;
1265  if (!is_increasing(range_bins))
1266  throw runtime_error(
1267  "The vector *range_bins* must contain strictly "
1268  "increasing values.");
1269  if (!is_z && min(range_bins) < 0)
1270  throw runtime_error(
1271  "The vector *range_bins* is not allowed to contain "
1272  "negative times.");
1273  if (instrument_pol_array.nelem() != nf)
1274  throw runtime_error(
1275  "The main length of *instrument_pol_array* must match "
1276  "the number of frequencies.");
1277 
1278  // iy_unit and variables to handle conversion to Ze and dBZe
1279  Vector cfac(nf, 1.0);
1280  Numeric ze_min = 0;
1281  const Numeric jfac = 10 * log10(exp(1.0));
1282  if (iy_unit == "1") {
1283  } else if (iy_unit == "Ze") {
1284  ze_cfac(cfac, f_grid, ze_tref, k2);
1285  } else if (iy_unit == "dBZe") {
1286  ze_cfac(cfac, f_grid, ze_tref, k2);
1287  ze_min = pow(10.0, dbze_min / 10);
1288  } else {
1289  throw runtime_error(
1290  "For this method, *iy_unit* must be set to \"1\", "
1291  "\"Ze\" or \"dBZe\".");
1292  }
1293 
1294  //---------------------------------------------------------------------------
1295  // Various initializations
1296  //---------------------------------------------------------------------------
1297 
1298  // Variables to handle conversion from Stokes to instrument_pol
1299  ArrayOfIndex npolcum(nf + 1);
1300  npolcum[0] = 0;
1301  ArrayOfArrayOfVector W(nf);
1302  for (Index i = 0; i < nf; i++) {
1303  const Index ni = instrument_pol_array[i].nelem();
1304  npolcum[i + 1] = npolcum[i] + ni;
1305  W[i].resize(ni);
1306  for (Index j = 0; j < ni; j++) {
1307  W[i][j].resize(stokes_dim);
1308  stokes2pol(W[i][j], stokes_dim, instrument_pol_array[i][j], 0.5);
1309  }
1310  }
1311 
1312  // Resize and init *y* and *y_XXX*
1313  //
1314  const Index ntot = npos * npolcum[nf] * nbins;
1315  y.resize(ntot);
1316  y = NAN;
1317  y_f.resize(ntot);
1318  y_pol.resize(ntot);
1319  y_pos.resize(ntot, sensor_pos.ncols());
1320  y_los.resize(ntot, sensor_los.ncols());
1321  y_geo.resize(ntot, 5);
1322  y_geo = NAN; // Will be replaced if relavant data are provided
1323 
1324  // y_aux
1325  y_aux.resize(naux);
1326  for (Index i = 0; i < naux; i++) {
1327  y_aux[i].resize(ntot);
1328  y_aux[i] = NAN;
1329  }
1330 
1331  // Jacobian variables
1332  //
1333  Index j_analytical_do = 0;
1334  Index njq = 0;
1335  ArrayOfArrayOfIndex jacobian_indices;
1336  //
1337  if (jacobian_do) {
1338  bool any_affine;
1339  jac_ranges_indices(jacobian_indices, any_affine, jacobian_quantities, true);
1340 
1341  jacobian.resize(ntot,
1342  jacobian_indices[jacobian_indices.nelem() - 1][1] + 1);
1343  jacobian = 0;
1344  njq = jacobian_quantities.nelem();
1345  //
1346  FOR_ANALYTICAL_JACOBIANS_DO(j_analytical_do = 1;)
1347  } else {
1348  jacobian.resize(0, 0);
1349  }
1350 
1351  //---------------------------------------------------------------------------
1352  // The calculations
1353  //---------------------------------------------------------------------------
1354 
1355  // Loop positions
1356  for (Index p = 0; p < npos; p++) {
1357  // RT part
1358  Tensor3 iy_transmission(0, 0, 0);
1360  Vector rte_pos2(0);
1361  Matrix iy;
1362  Ppath ppath;
1364  const Index iy_id = (Index)1e6 * p;
1365  //
1367  iy,
1368  iy_aux,
1369  ppath,
1370  diy_dx,
1371  1,
1373  iy_aux_vars,
1374  iy_id,
1375  iy_unit,
1376  cloudbox_on,
1377  jacobian_do,
1378  f_grid,
1379  nlte_field,
1380  sensor_pos(p, joker),
1381  sensor_los(p, joker),
1382  rte_pos2,
1383  iy_main_agenda);
1384 
1385  // Check if path and size OK
1386  const Index np = ppath.np;
1387  if (np == 1)
1388  throw runtime_error(
1389  "A path consisting of a single point found. "
1390  "This is not allowed.");
1392  if (iy.nrows() != nf * np)
1393  throw runtime_error(
1394  "The size of *iy* returned from *iy_main_agenda* "
1395  "is not correct (for this method).");
1396 
1397  // Range of ppath, in altitude or time
1398  Vector range(np);
1399  if (is_z) {
1400  range = ppath.pos(joker, 0);
1401  } else { // Calculate round-trip time
1402  range[0] = 2 * ppath.end_lstep / SPEED_OF_LIGHT;
1403  for (Index i = 1; i < np; i++) {
1404  range[i] = range[i - 1] + ppath.lstep[i - 1] *
1405  (ppath.ngroup[i - 1] + ppath.ngroup[i]) /
1407  }
1408  }
1409  const Numeric range_end1 = min(range[0], range[np - 1]);
1410  const Numeric range_end2 = max(range[0], range[np - 1]);
1411 
1412  // Loop radar bins
1413  for (Index b = 0; b < nbins; b++) {
1414  if (!(range_bins[b] >= range_end2 || // Otherwise bin totally outside
1415  range_bins[b + 1] <= range_end1)) // range of ppath
1416  {
1417  // Bin limits
1418  Numeric blim1 = max(range_bins[b], range_end1);
1419  Numeric blim2 = min(range_bins[b + 1], range_end2);
1420 
1421  // Determine weight vector to obtain mean inside bin
1422  Vector hbin(np);
1423  integration_bin_by_vecmult(hbin, range, blim1, blim2);
1424  // The function above handles integration over the bin, while we
1425  // want the average, so divide weights with bin width
1426  hbin /= (blim2 - blim1);
1427 
1428  for (Index iv = 0; iv < nf; iv++) {
1429  // Pick out part of iy for frequency
1430  Matrix I = iy(Range(iv * np, np), joker);
1431  ArrayOfTensor3 dI(njq);
1432  if (j_analytical_do) {
1434  dI[iq] = diy_dx[iq](joker, Range(iv * np, np), joker);)
1435  }
1436  ArrayOfMatrix A(naux);
1437  for (Index a = 0; a < naux; a++) {
1438  A[a] = iy_aux[a](Range(iv * np, np), joker);
1439  }
1440 
1441  // Variables to hold data for one freq and one pol
1442  Vector refl(np);
1443  ArrayOfMatrix drefl(njq);
1444  if (j_analytical_do) {
1445  FOR_ANALYTICAL_JACOBIANS_DO(drefl[iq].resize(dI[iq].npages(), np);)
1446  }
1447  ArrayOfVector auxvar(naux);
1448  for (Index a = 0; a < naux; a++) {
1449  auxvar[a].resize(np);
1450  }
1451 
1452  for (Index ip = 0; ip < instrument_pol_array[iv].nelem(); ip++) {
1453  // Apply weights on each Stokes element
1454  mult(refl, I, W[iv][ip]);
1455  if (j_analytical_do) {
1457  k < drefl[iq].nrows();
1458  k++) {
1459  mult(drefl[iq](k, joker), dI[iq](k, joker, joker), W[iv][ip]);
1460  })
1461  }
1462  for (Index a = 0; a < naux; a++) {
1463  if (iy_aux_vars[a] == "Backscattering") {
1464  mult(auxvar[a], A[a], W[iv][ip]);
1465  } else {
1466  for (Index j = 0; j < np; j++) {
1467  auxvar[a][j] = A[a](j, 0);
1468  }
1469  }
1470  }
1471 
1472  // Apply bin weight vector to get final values.
1473  Index iout = nbins * (p * npolcum[nf] + npolcum[iv] + ip) + b;
1474  y[iout] = cfac[iv] * (hbin * refl);
1475  //
1476  if (j_analytical_do) {
1478  for (Index k = 0; k < drefl[iq].nrows(); k++) {
1479  jacobian(iout, jacobian_indices[iq][0] + k) =
1480  cfac[iv] * (hbin * drefl[iq](k, joker));
1481  if (iy_unit == "dBZe") {
1482  jacobian(iout, jacobian_indices[iq][0] + k) *=
1483  jfac / max(y[iout], ze_min);
1484  }
1485  })
1486  }
1487 
1488  if (iy_unit == "dBZe") {
1489  y[iout] = y[iout] <= ze_min ? dbze_min : 10 * log10(y[iout]);
1490  }
1491 
1492  // Same for aux variables
1493  for (Index a = 0; a < naux; a++) {
1494  if (iy_aux_vars[a] == "Backscattering") {
1495  y_aux[a][iout] = cfac[iv] * (hbin * auxvar[a]);
1496  if (iy_unit == "dBZe") {
1497  y_aux[a][iout] = y_aux[a][iout] <= ze_min
1498  ? dbze_min
1499  : 10 * log10(y_aux[a][iout]);
1500  }
1501  } else {
1502  y_aux[a][iout] = hbin * auxvar[a];
1503  }
1504  }
1505  }
1506  } // Frequency
1507  }
1508  }
1509 
1510  // Other aux variables
1511  //
1512  Vector geo_pos;
1514  if (geo_pos.nelem() && geo_pos.nelem() != atmosphere_dim) {
1515  throw runtime_error(
1516  "Wrong size of *geo_pos* obtained from "
1517  "*geo_pos_agenda*.\nThe length of *geo_pos* must "
1518  "be zero or equal to *atmosphere_dim*.");
1519  }
1520  //
1521  for (Index b = 0; b < nbins; b++) {
1522  for (Index iv = 0; iv < nf; iv++) {
1523  for (Index ip = 0; ip < instrument_pol_array[iv].nelem(); ip++) {
1524  const Index iout = nbins * (p * npolcum[nf] + npolcum[iv] + ip) + b;
1525  y_f[iout] = f_grid[iv];
1526  y_pol[iout] = instrument_pol_array[iv][ip];
1527  y_pos(iout, joker) = sensor_pos(p, joker);
1528  y_los(iout, joker) = sensor_los(p, joker);
1529  if (geo_pos.nelem()) {
1530  y_geo(iout, joker) = geo_pos;
1531  }
1532  }
1533  }
1534  }
1535  }
1536 }
ARTS::Var::ppvar_mag
Matrix ppvar_mag(Workspace &ws) noexcept
Definition: autoarts.h:5259
Matrix
The Matrix class.
Definition: matpackI.h:1193
ARTS::Var::y_pos
Matrix y_pos(Workspace &ws) noexcept
Definition: autoarts.h:7531
ARTS::Var::atmosphere_dim
Index atmosphere_dim(Workspace &ws) noexcept
Definition: autoarts.h:2510
ARTS::Var::Ppath::pos
std::size_t pos() const noexcept
Definition: autoarts.h:1300
ARTS::Var::wind_v_field
Tensor3 wind_v_field(Workspace &ws) noexcept
Definition: autoarts.h:7246
adapt_stepwise_partial_derivatives
void adapt_stepwise_partial_derivatives(ArrayOfPropagationMatrix &dK_dx, ArrayOfStokesVector &dS_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, ConstVectorView ppath_f_grid, ConstVectorView ppath_line_of_sight, ConstVectorView ppath_vmrs, const Numeric &ppath_temperature, const Numeric &ppath_pressure, const ArrayOfIndex &jacobian_species, const ArrayOfIndex &jacobian_wind, const Index &lte, const Index &atmosphere_dim, const bool &jacobian_do)
Adapts clearsky partial derivatives.
Definition: rte.cc:59
pha_mat_1ScatElem
void pha_mat_1ScatElem(Tensor6View pha_mat, Index &ptype, VectorView t_ok, const SingleScatteringData &ssd, const Vector &T_array, const Matrix &pdir_array, const Matrix &idir_array, const Index &f_start, const Index &t_interp_order)
Preparing phase matrix from one scattering element.
Definition: optproperties.cc:881
ARTS::Var::y_aux
ArrayOfVector y_aux(Workspace &ws) noexcept
Definition: autoarts.h:7423
ARTS::Var::geo_pos
Vector geo_pos(Workspace &ws) noexcept
Definition: autoarts.h:3544
cumulative_backscatter
ArrayOfTransmissionMatrix cumulative_backscatter(ConstTensor5View t, ConstMatrixView m)
Accumulated backscatter (???)
Definition: transmissionmatrix.cc:1714
ARTS::Group::PropagationMatrix
PropagationMatrix PropagationMatrix
Definition: autoarts.h:93
geo_pos_agendaExecute
void geo_pos_agendaExecute(Workspace &ws, Vector &geo_pos, const Ppath &ppath, const Agenda &input_agenda)
Definition: auto_md.cc:23915
StokesVector
Stokes vector is as Propagation matrix but only has 4 possible values.
Definition: propagationmatrix.h:1075
auto_md.h
ARTS::Var::atmgeom_checked
Index atmgeom_checked(Workspace &ws) noexcept
Definition: autoarts.h:2489
ARTS::Var::jacobian_quantities
ArrayOfRetrievalQuantity jacobian_quantities(Workspace &ws) noexcept
Definition: autoarts.h:3924
chk_if_increasing
void chk_if_increasing(const String &x_name, const ArrayOfIndex &x)
chk_if_increasing
Definition: check_input.cc:117
ARTS::Var::instrument_pol_array
ArrayOfArrayOfIndex instrument_pol_array(Workspace &ws) noexcept
Definition: autoarts.h:3628
Tensor3
The Tensor3 class.
Definition: matpackIII.h:339
iyActiveSingleScat2
void iyActiveSingleScat2(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, Tensor4 &ppvar_trans_cumulat, 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 Index &scat_data_checked, const ArrayOfString &iy_aux_vars, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Ppath &ppath, const Agenda &propmat_clearsky_agenda, const Agenda &water_p_eq_agenda, const Agenda &iy_transmitter_agenda, const Index &iy_agenda_call1, const Tensor3 &iy_transmission, const Numeric &rte_alonglos_v, const Index &trans_in_jacobian, const Numeric &pext_scaling, const Index &t_interp_order, const Verbosity &verbosity)
WORKSPACE METHOD: iyActiveSingleScat2.
Definition: m_cloudradar.cc:680
ARTS::Var::ppvar_vmr
Matrix ppvar_vmr(Workspace &ws) noexcept
Definition: autoarts.h:5372
ARTS::Var::iy_transmission
Tensor3 iy_transmission(Workspace &ws) noexcept
Definition: autoarts.h:3836
integration_bin_by_vecmult
void integration_bin_by_vecmult(VectorView h, ConstVectorView x_g_in, const Numeric &limit1, const Numeric &limit2)
integration_bin_by_vecmult
Definition: sensor.cc:1503
ARTS::Var::iy_transmitter_agenda
Agenda iy_transmitter_agenda(Workspace &ws) noexcept
Definition: autoarts.h:3843
joker
const Joker joker
ze_cfac
void ze_cfac(Vector &fac, const Vector &f_grid, const Numeric &ze_tref, const Numeric &k2)
Calculates factor to convert back-scattering to Ze.
Definition: rte.cc:2736
get_ppath_cloudvars
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.
Definition: rte.cc:1153
ARTS::Var::pnd_field
Tensor4 pnd_field(Workspace &ws) noexcept
Definition: autoarts.h:5081
ARTS::Var::iy_id
Index iy_id(Workspace &ws) noexcept
Definition: autoarts.h:3775
ARTS::Group::Tensor3
Tensor3 Tensor3
Definition: autoarts.h:104
abs
#define abs(x)
Definition: legacy_continua.cc:20626
ARTS::Var::verbosity
Verbosity verbosity(Workspace &ws) noexcept
Definition: autoarts.h:7112
ARTS::Group::Tensor4
Tensor4 Tensor4
Definition: autoarts.h:105
error_if_limb_ppath
void error_if_limb_ppath(const Ppath &ppath)
Throws an error if ppath altitudes not are strictly increasing or decreasing.
Definition: ppath.cc:555
mult
void mult(ComplexVectorView y, const ConstComplexMatrixView &M, const ConstComplexVectorView &x)
Matrix-Vector Multiplication.
Definition: complex.cc:1579
cumulative_backscatter_derivative
ArrayOfArrayOfTransmissionMatrix cumulative_backscatter_derivative(ConstTensor5View t, const ArrayOfMatrix &aom)
Accumulated backscatter derivative (???)
Definition: transmissionmatrix.cc:1755
swap
void swap(ComplexVector &v1, ComplexVector &v2)
Swaps two objects.
Definition: complex.cc:731
ARTS::Var::cloudbox_limits
ArrayOfIndex cloudbox_limits(Workspace &ws) noexcept
Definition: autoarts.h:2762
ARTS::Var::jacobian_do
Index jacobian_do(Workspace &ws) noexcept
Definition: autoarts.h:3912
ARTS::Var::wind_u_field
Tensor3 wind_u_field(Workspace &ws) noexcept
Definition: autoarts.h:7211
ARTS::Var::y
Vector y(Workspace &ws) noexcept
Definition: autoarts.h:7401
PropagationMatrix
Definition: propagationmatrix.h:87
Vector::resize
void resize(Index n)
Resize function.
Definition: matpackI.cc:404
ARTS::Var::diy_dx
ArrayOfTensor3 diy_dx(Workspace &ws) noexcept
Definition: autoarts.h:3007
ARTS::Var::y_los
Matrix y_los(Workspace &ws) noexcept
Definition: autoarts.h:7498
ARTS::Var::ppvar_f
Matrix ppvar_f(Workspace &ws) noexcept
Definition: autoarts.h:5233
ARTS::Var::mag_w_field
Tensor3 mag_w_field(Workspace &ws) noexcept
Definition: autoarts.h:4197
ARTS::Var::stokes_dim
Index stokes_dim(Workspace &ws) noexcept
Definition: autoarts.h:6650
ConstMatrixView::nrows
Index nrows() const
Returns the number of rows.
Definition: matpackI.cc:429
FOR_ANALYTICAL_JACOBIANS_DO
#define FOR_ANALYTICAL_JACOBIANS_DO(what_to_do)
Definition: jacobian.h:405
ARTS::Var::iy_unit
String iy_unit(Workspace &ws) noexcept
Definition: autoarts.h:3856
ARTS::Var::iy
Matrix iy(Workspace &ws) noexcept
Definition: autoarts.h:3690
CumulativeTransmission::Forward
@ Forward
montecarlo.h
TotalNumberOfElements
Index TotalNumberOfElements(const Array< Array< base > > &aa)
Determine total number of elements in an ArrayOfArray.
Definition: array.h:343
Tensor4
The Tensor4 class.
Definition: matpackIV.h:421
Ppath
The structure to describe a propagation path and releated quantities.
Definition: ppath.h:48
Agenda
The Agenda class.
Definition: agenda_class.h:44
ARTS::Var::range_bins
Vector range_bins(Workspace &ws) noexcept
Definition: autoarts.h:5493
SCATSPECIES_MAINTAG
const String SCATSPECIES_MAINTAG
rtmethods_jacobian_init
void rtmethods_jacobian_init(ArrayOfIndex &jac_species_i, ArrayOfIndex &jac_scat_i, ArrayOfIndex &jac_is_t, ArrayOfIndex &jac_wind_i, ArrayOfIndex &jac_mag_i, ArrayOfIndex &jac_other, ArrayOfTensor3 &diy_dx, ArrayOfTensor3 &diy_dpath, const Index &ns, const Index &nf, const Index &np, const Index &nq, const ArrayOfArrayOfSpeciesTag &abs_species, const Index &cloudbox_on, const ArrayOfString &scat_species, const ArrayOfTensor4 &dpnd_field_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, const Index &iy_agenda_call1, const bool is_active)
This function fixes the initial steps around Jacobian calculations, to be done inside radiative trans...
Definition: rte.cc:2348
ARTS::Var::ppvar_trans_cumulat
Tensor4 ppvar_trans_cumulat(Workspace &ws) noexcept
Definition: autoarts.h:5346
set_backscatter_radiation_vector
void set_backscatter_radiation_vector(ArrayOfRadiationVector &I, ArrayOfArrayOfArrayOfRadiationVector &dI, const RadiationVector &I_incoming, const ArrayOfTransmissionMatrix &T, const ArrayOfTransmissionMatrix &PiTf, const ArrayOfTransmissionMatrix &PiTr, const ArrayOfTransmissionMatrix &Z, const ArrayOfArrayOfTransmissionMatrix &dT1, const ArrayOfArrayOfTransmissionMatrix &dT2, const ArrayOfArrayOfTransmissionMatrix &dZ, const BackscatterSolver solver)
Set the backscatter radiation vector.
Definition: transmissionmatrix.cc:1548
get_stepwise_scattersky_propmat
void get_stepwise_scattersky_propmat(StokesVector &ap, PropagationMatrix &Kp, ArrayOfStokesVector &dap_dx, ArrayOfPropagationMatrix &dKp_dx, const ArrayOfRetrievalQuantity &jacobian_quantities, ConstMatrixView ppath_1p_pnd, const ArrayOfMatrix &ppath_dpnd_dx, const Index ppath_1p_id, const ArrayOfArrayOfSingleScatteringData &scat_data, ConstVectorView ppath_line_of_sight, ConstVectorView ppath_temperature, const Index &atmosphere_dim, const bool &jacobian_do)
Computes the contribution by scattering at propagation path point.
Definition: rte.cc:1591
ARTS::Var::rte_alonglos_v
Numeric rte_alonglos_v(Workspace &ws) noexcept
Definition: autoarts.h:5628
ARTS::Var::nlte_field
EnergyLevelMap nlte_field(Workspace &ws) noexcept
Definition: autoarts.h:4604
Array
This can be used to make arrays out of anything.
Definition: array.h:108
pow
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:628
yActive
void yActive(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 String &iy_unit, const ArrayOfString &iy_aux_vars, const Index &stokes_dim, const Vector &f_grid, const Index &atmosphere_dim, const EnergyLevelMap &nlte_field, const Index &cloudbox_on, const Index &cloudbox_checked, const Matrix &sensor_pos, const Matrix &sensor_los, const Index &sensor_checked, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Agenda &iy_main_agenda, const Agenda &geo_pos_agenda, const ArrayOfArrayOfIndex &instrument_pol_array, const Vector &range_bins, const Numeric &ze_tref, const Numeric &k2, const Numeric &dbze_min, const Verbosity &)
WORKSPACE METHOD: yActive.
Definition: m_cloudradar.cc:1190
sensor.h
Sensor modelling functions.
get_stepwise_clearsky_propmat
void get_stepwise_clearsky_propmat(Workspace &ws, PropagationMatrix &K, StokesVector &S, Index &lte, ArrayOfPropagationMatrix &dK_dx, ArrayOfStokesVector &dS_dx, const Agenda &propmat_clearsky_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, ConstVectorView ppath_f_grid, ConstVectorView ppath_magnetic_field, ConstVectorView ppath_line_of_sight, const EnergyLevelMap &ppath_nlte, ConstVectorView ppath_vmrs, const Numeric &ppath_temperature, const Numeric &ppath_pressure, const ArrayOfIndex &jacobian_species, const bool &jacobian_do)
Gets the clearsky propgation matrix and NLTE contributions.
Definition: rte.cc:1322
min
#define min(a, b)
Definition: legacy_continua.cc:20628
ArrayOfArrayOfRadiationVector
Array< ArrayOfRadiationVector > ArrayOfArrayOfRadiationVector
Definition: transmissionmatrix.h:900
Tensor5::resize
void resize(Index s, Index b, Index p, Index r, Index c)
Resize function.
Definition: matpackV.cc:1743
stepwise_transmission
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.
Definition: transmissionmatrix.cc:1353
messages.h
Declarations having to do with the four output streams.
ARTS::Var::y_f
Vector y_f(Workspace &ws) noexcept
Definition: autoarts.h:7466
ConstMatrixView::ncols
Index ncols() const
Returns the number of columns.
Definition: matpackI.cc:432
_U_
#define _U_
Definition: config.h:183
my_basic_string< char >
ARTS::Var::y_pol
ArrayOfIndex y_pol(Workspace &ws) noexcept
Definition: autoarts.h:7515
ARTS::Var::ppvar_pnd
Matrix ppvar_pnd(Workspace &ws) noexcept
Definition: autoarts.h:5317
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
SPEED_OF_LIGHT
const Numeric SPEED_OF_LIGHT
ArrayOfRadiationVector
Array< RadiationVector > ArrayOfRadiationVector
Definition: transmissionmatrix.h:899
ARTS::Var::vmr_field
Tensor4 vmr_field(Workspace &ws) noexcept
Definition: autoarts.h:7130
ARTS::Var::ppvar_nlte
EnergyLevelMap ppvar_nlte(Workspace &ws) noexcept
Definition: autoarts.h:5272
ARTS::Var::abs_species
ArrayOfArrayOfSpeciesTag abs_species(Workspace &ws) noexcept
Definition: autoarts.h:2157
ARTS::Var::ppvar_t
Vector ppvar_t(Workspace &ws) noexcept
Definition: autoarts.h:5330
ARTS::Var::atmfields_checked
Index atmfields_checked(Workspace &ws) noexcept
Definition: autoarts.h:2473
get_ppath_atmvars
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, ConstVectorView p_grid, ConstTensor3View t_field, const EnergyLevelMap &nlte_field, ConstTensor4View vmr_field, ConstTensor3View wind_u_field, ConstTensor3View wind_v_field, ConstTensor3View wind_w_field, ConstTensor3View mag_u_field, ConstTensor3View mag_v_field, ConstTensor3View mag_w_field)
Determines pressure, temperature, VMR, winds and magnetic field for each propgataion path point.
Definition: rte.cc:1034
ARTS::Var::p_grid
Vector p_grid(Workspace &ws) noexcept
Definition: autoarts.h:4763
ARTS::Var::ppath
Ppath ppath(Workspace &ws) noexcept
Definition: autoarts.h:5139
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
Verbosity
Definition: messages.h:49
ARTS::Var::iy_agenda_call1
Index iy_agenda_call1(Workspace &ws) noexcept
Definition: autoarts.h:3704
ARTS::Var::npages
Index npages(Workspace &ws) noexcept
Definition: autoarts.h:4675
ARTS::Var::iy_aux_vars
ArrayOfString iy_aux_vars(Workspace &ws) noexcept
Definition: autoarts.h:3735
ARTS::Var::f_grid
Vector f_grid(Workspace &ws) noexcept
Definition: autoarts.h:3449
ARTS::Var::sensor_pos
Matrix sensor_pos(Workspace &ws) noexcept
Definition: autoarts.h:6279
EnergyLevelMap
Definition: energylevelmap.h:60
max
#define max(a, b)
Definition: legacy_continua.cc:20629
iyActiveSingleScat
void iyActiveSingleScat(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, Tensor4 &ppvar_trans_cumulat, 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 Index &scat_data_checked, const ArrayOfString &iy_aux_vars, const Index &jacobian_do, const ArrayOfRetrievalQuantity &jacobian_quantities, const Ppath &ppath, const Agenda &propmat_clearsky_agenda, const Agenda &water_p_eq_agenda, const Agenda &iy_transmitter_agenda, const Index &iy_agenda_call1, const Tensor3 &iy_transmission, const Numeric &rte_alonglos_v, const Index &trans_in_jacobian, const Numeric &pext_scaling, const Index &t_interp_order, const Verbosity &verbosity)
WORKSPACE METHOD: iyActiveSingleScat.
Definition: m_cloudradar.cc:49
propagationmatrix.h
Stuff related to the propagation matrix.
ARTS::Var::ppvar_wind
Matrix ppvar_wind(Workspace &ws) noexcept
Definition: autoarts.h:5385
ARTS::Var::scat_data
ArrayOfArrayOfSingleScatteringData scat_data(Workspace &ws) noexcept
Definition: autoarts.h:5840
iy_transmitter_agendaExecute
void iy_transmitter_agendaExecute(Workspace &ws, Matrix &iy, const Vector &f_grid, const Vector &rtp_pos, const Vector &rtp_los, const Agenda &input_agenda)
Definition: auto_md.cc:24504
ARTS::Var::dpnd_field_dx
ArrayOfTensor4 dpnd_field_dx(Workspace &ws) noexcept
Definition: autoarts.h:3251
ArrayOfTransmissionMatrix
Array< TransmissionMatrix > ArrayOfTransmissionMatrix
Definition: transmissionmatrix.h:896
Tensor5
The Tensor5 class.
Definition: matpackV.h:506
cumulative_transmission
ArrayOfTransmissionMatrix cumulative_transmission(const ArrayOfTransmissionMatrix &T, const CumulativeTransmission type)
Accumulate the transmission matrix over all layers.
Definition: transmissionmatrix.cc:1527
ARTS::Var::scat_data_checked
Index scat_data_checked(Workspace &ws) noexcept
Definition: autoarts.h:5850
PI
const Numeric PI
Range
The range class.
Definition: matpackI.h:160
RadiationVector
Radiation Vector for Stokes dimension 1-4.
Definition: transmissionmatrix.h:395
ARTS::Var::iy_aux
ArrayOfMatrix iy_aux(Workspace &ws) noexcept
Definition: autoarts.h:3719
ARTS::Var::geo_pos_agenda
Agenda geo_pos_agenda(Workspace &ws) noexcept
Definition: autoarts.h:3551
is_increasing
bool is_increasing(ConstVectorView x)
Checks if a vector is sorted and strictly increasing.
Definition: logic.cc:215
logic.h
Header file for logic.cc.
ARTS::Var::iy_main_agenda
Agenda iy_main_agenda(Workspace &ws) noexcept
Definition: autoarts.h:3794
ARTS::Var::rte_pos2
Vector rte_pos2(Workspace &ws) noexcept
Definition: autoarts.h:5696
ARTS::Group::RadiationVector
RadiationVector RadiationVector
Definition: autoarts.h:95
stokes2pol
void stokes2pol(VectorView w, const Index &stokes_dim, const Index &ipol_1based, const Numeric nv)
stokes2pol
Definition: sensor.cc:1041
mirror_los
void mirror_los(Vector &los_mirrored, ConstVectorView los, const Index &atmosphere_dim)
Determines the backward direction for a given line-of-sight.
Definition: rte.cc:2290
rte.h
Declaration of functions in rte.cc.
ARTS::Var::mag_v_field
Tensor3 mag_v_field(Workspace &ws) noexcept
Definition: autoarts.h:4164
ARTS::Var::t_field
Tensor3 t_field(Workspace &ws) noexcept
Definition: autoarts.h:6947
Workspace
Workspace class.
Definition: workspace_ng.h:40
PropagationMatrix::Kjj
VectorView Kjj(const Index iz=0, const Index ia=0)
Vector view to diagonal elements.
Definition: propagationmatrix.h:799
ARTS::Var::wind_w_field
Tensor3 wind_w_field(Workspace &ws) noexcept
Definition: autoarts.h:7279
get_stepwise_transmission_matrix
void get_stepwise_transmission_matrix(Tensor3View cumulative_transmission, Tensor3View T, Tensor4View dT_close_dx, Tensor4View dT_far_dx, ConstTensor3View cumulative_transmission_close, const PropagationMatrix &K_close, const PropagationMatrix &K_far, const ArrayOfPropagationMatrix &dK_close_dx, const ArrayOfPropagationMatrix &dK_far_dx, const Numeric &ppath_distance, const bool &first_level, const Numeric &dppath_distance_dT_HSE_close, const Numeric &dppath_distance_dT_HSE_far, const Index &temperature_derivative_position_if_hse_is_active)
Computes layer transmission matrix and cumulative transmission.
Definition: rte.cc:1879
ARTS::Group::TransmissionMatrix
TransmissionMatrix TransmissionMatrix
Definition: autoarts.h:112
ARTS::Var::sensor_checked
Index sensor_checked(Workspace &ws) noexcept
Definition: autoarts.h:6144
ARTS::Var::propmat_clearsky_agenda
Agenda propmat_clearsky_agenda(Workspace &ws) noexcept
Definition: autoarts.h:5405
rtmethods_jacobian_finalisation
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_transmission, const Agenda &water_p_eq_agenda, const ArrayOfRetrievalQuantity &jacobian_quantities, const ArrayOfIndex jac_species_i, const ArrayOfIndex jac_is_t)
This function fixes the last steps to made on the Jacobian in some radiative transfer WSMs.
Definition: rte.cc:2416
iy_main_agendaExecute
void iy_main_agendaExecute(Workspace &ws, Matrix &iy, ArrayOfMatrix &iy_aux, Ppath &ppath, ArrayOfTensor3 &diy_dx, const Index iy_agenda_call1, const Tensor3 &iy_transmission, 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)
Definition: auto_md.cc:24221
ARTS::Var::jacobian
Matrix jacobian(Workspace &ws) noexcept
Definition: autoarts.h:3892
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
chk_if_in_range
void chk_if_in_range(const String &x_name, const Index &x, const Index &x_low, const Index &x_high)
chk_if_in_range
Definition: check_input.cc:89
ARTS::Var::scat_species
ArrayOfString scat_species(Workspace &ws) noexcept
Definition: autoarts.h:6086
Tensor6
The Tensor6 class.
Definition: matpackVI.h:1088
ARTS::Var::mag_u_field
Tensor3 mag_u_field(Workspace &ws) noexcept
Definition: autoarts.h:4130
ppath_what_background
Index ppath_what_background(const Ppath &ppath)
Returns the case number for the radiative background.
Definition: ppath.cc:1494
Vector
The Vector class.
Definition: matpackI.h:860
ARTS::Var::ppvar_p
Vector ppvar_p(Workspace &ws) noexcept
Definition: autoarts.h:5304
jac_ranges_indices
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.
Definition: jacobian.cc:58
ARTS::Var::water_p_eq_agenda
Agenda water_p_eq_agenda(Workspace &ws) noexcept
Definition: autoarts.h:7165
Array::nelem
Index nelem() const
Number of elements.
Definition: array.h:195
ARTS::Var::sensor_los
Matrix sensor_los(Workspace &ws) noexcept
Definition: autoarts.h:6201
ARTS::Var::y_geo
Matrix y_geo(Workspace &ws) noexcept
Definition: autoarts.h:7480
ARTS::Var::cloudbox_checked
Index cloudbox_checked(Workspace &ws) noexcept
Definition: autoarts.h:2650
ARTS::Group::StokesVector
StokesVector StokesVector
Definition: autoarts.h:101
get_ppath_f
void get_ppath_f(Matrix &ppath_f, const Ppath &ppath, ConstVectorView f_grid, const Index &atmosphere_dim, const Numeric &rte_alonglos_v, ConstMatrixView ppath_wind)
Determines the Doppler shifted frequencies along the propagation path.
Definition: rte.cc:1257
ns
#define ns
Definition: legacy_continua.cc:22183
arts.h
The global header file for ARTS.
ARTS::Var::cloudbox_on
Index cloudbox_on(Workspace &ws) noexcept
Definition: autoarts.h:2782
solve
void solve(VectorView w, const CovarianceMatrix &A, ConstVectorView v)
Definition: covariance_matrix.cc:594