ARTS  2.3.1285(git:92a29ea9-dirty)
lineshapesdata.h
Go to the documentation of this file.
1 /* Copyright (C) 2017
2  * Richard Larsson <ric.larsson@gmail.com>
3  *
4  * This program is free software; you can redistribute it and/or modify it
5  * under the terms of the GNU General Public License as published by the
6  * Free Software Foundation; either version 2, or (at your option) any
7  * later version.
8  *
9  * This program is distributed in the hope that it will be useful,
10  * but WITHOUT ANY WARRANTY; without even the implied warranty of
11  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  * GNU General Public License for more details.
13  *
14  * You should have received a copy of the GNU General Public License
15  * along with this program; if not, write to the Free Software
16  * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
17  * USA. */
18 
30 #ifndef lineshapedata_h
31 #define lineshapedata_h
32 
33 #include "complex.h"
34 #include "linerecord.h"
35 #include "Faddeeva.hh"
36 #include "absorption.h"
37 #include "partial_derivatives.h"
38 #include <Eigen/Dense>
39 
40 /*
41  * Class to solve the problem
42  *
43  * cross-section of a line equals line strength times line shape
44  *
45  * or
46  *
47  * sigma = S x F,
48  *
49  * which means
50  *
51  * dsigma = dS x F + S x dF
52  *
53  * TODO: Add QuantumIdentifier input to test for line catalog-parameters derivatives in relevant places
54  *
55  * TODO: Add NLTE line strength calculator
56  *
57  * TODO: Better combination with Zeeman calculations
58  *
59  * TODO: Find work-around for incomplete line-shapes like "Voigt Kuntz"
60  */
61 class LineFunctions
62 {
63 public:
64 
65  enum LineShapeType {
66  LST_NONE, // Holder that should always fail if evoked X
67  LST_O2_NON_RESONANT, // Debye lineshape X
68  LST_DOPPLER, // Doppler lineshape X
69  LST_LORENTZ, // Lorentz lineshape X
70  LST_MIRRORED_LORENTZ, // Mirrored Lorentz lineshape X
71  LST_FADDEEVA_ALGORITHM916, // Faddeeva lineshape X
72  LST_HUI_ETAL_1978, // Faddeeva lineshape X
73  LST_HTP // Hartmann-Tran lineshape X
74  };
75 
76  enum LineShapeNorm {
77  LSN_NONE, // No normalization X
78  LSN_ROSENKRANZ_QUADRATIC, // Quadratic normalization (f/f0)^2*h*f0/(2*k*T)/sinh(h*f0/(2*k*T)) X
79  LSN_VVW, // Van Vleck Weiskopf normalization (f*f) / (f0*f0) X
80  LSN_VVH // Van Vleck Huber normalization f*tanh(h*f/(2*k*T))) / (f0*tanh(h*f0/(2*k*T))) X
81  };
82 
83  void set_lorentz(ComplexVector& F, // Sets the full complex line shape without line mixing
84  ArrayOfComplexVector& dF,
85  const Vector& f_grid,
86  const Numeric& zeeman_df,
87  const Numeric& magnetic_magnitude,
88  const Numeric& F0_noshift,
89  const Numeric& G0,
90  const Numeric& L0,
91  const PropmatPartialsData& derivatives_data=PropmatPartialsData(),
92  const Numeric& dG0_dT=0.0,
93  const Numeric& dL0_dT=0.0) const
94  {
95  const Index nf = f_grid.nelem(), nppd = derivatives_data.nelem();
96 
97  extern const Numeric PI;
98 
99  const Numeric invPI = 1.0 / PI;
100 
101  const Numeric F0 = F0_noshift + L0 + zeeman_df * magnetic_magnitude;
102 
103  // Signa change of F and F0?
104  const Complex denom0 = Complex(G0, F0);
105 
106  F = invPI;
107 
108  for(Index ii = 0; ii < nf; ii++)
109  {
110  F[ii] /= (denom0 - Complex(0.0, f_grid[ii]));
111  }
112 
113  if(nppd > 0)
114  {
115  dF[nppd-1] = F;
116  dF[nppd-1] *= F;
117  dF[nppd-1] *= PI;
118  }
119 
120  for(Index jj = 0; jj < nppd; jj++)
121  {
122 
123  if(derivatives_data(jj) == JQT_temperature)
124  {
125  dF[jj] = dF[nppd-1];
126  dF[jj] *= Complex(dG0_dT, dL0_dT);
127  }
128  else if(derivatives_data(jj) == JQT_frequency or
129  derivatives_data(jj) == JQT_wind_magnitude or
130  derivatives_data(jj) == JQT_wind_u or
131  derivatives_data(jj) == JQT_wind_v or
132  derivatives_data(jj) == JQT_wind_w)
133  {
134  dF[jj] = dF[nppd-1];
135  dF[jj] *= Complex(0.0, -1.0);
136  }
137  else if(derivatives_data(jj) == JQT_line_center)
138  {
139  dF[jj] = dF[nppd-1];
140  dF[jj] *= Complex(0.0, 1.0);
141  }
142  else if(derivatives_data(jj) == JQT_line_gamma_self or
143  derivatives_data(jj) == JQT_line_gamma_selfexponent or
144  derivatives_data(jj) == JQT_line_pressureshift_self or
145  derivatives_data(jj) == JQT_line_gamma_foreign or
146  derivatives_data(jj) == JQT_line_gamma_foreignexponent or
147  derivatives_data(jj) == JQT_line_pressureshift_foreign or
148  derivatives_data(jj) == JQT_line_gamma_water or
149  derivatives_data(jj) == JQT_line_gamma_waterexponent or
150  derivatives_data(jj) == JQT_line_pressureshift_water) // Only the zeroth order terms --- the derivative with respect to these have to happen later
151  {
152  dF[jj] = dF[nppd-1];
153  dF[jj] *= Complex(1.0, 1.0);
154  }
155  else if(derivatives_data(jj) == JQT_magnetic_magntitude)// No external inputs --- errors because of frequency shift when Zeeman is used?
156  {
157  dF[jj] = dF[nppd-1];
158  dF[jj] *= Complex(0.0, zeeman_df);
159  }
160  }
161  }
162 
163  void set_mirrored_lorentz(ComplexVector& F, // Sets the full complex line shape without line mixing
164  ArrayOfComplexVector& dF,
165  const Vector& f_grid,
166  const Numeric& zeeman_df,
167  const Numeric& magnetic_magnitude,
168  const Numeric& F0_noshift,
169  const Numeric& G0,
170  const Numeric& L0,
171  const PropmatPartialsData& derivatives_data=PropmatPartialsData(),
172  const Numeric& dG0_dT=0.0,
173  const Numeric& dL0_dT=0.0) const
174  {
175  Complex Fplus, Fminus;
176 
177  const Index nf = f_grid.nelem(), nppd = derivatives_data.nelem();
178 
179  extern const Numeric PI;
180 
181  const Numeric invPI = 1.0 / PI;
182 
183  const Numeric F0 = F0_noshift + L0 + zeeman_df * magnetic_magnitude;
184 
185  // Signa change of F and F0?
186  const Complex denom0 = Complex(G0, F0), denom1 = Complex(G0, -F0);
187 
188  for(Index ii = 0; ii < nf; ii++)
189  {
190  Fplus = invPI / (denom0 - Complex(0.0, f_grid[ii]));
191  Fminus = invPI / (denom1 - Complex(0.0, f_grid[ii]));
192  F[ii] = Fplus + Fminus;
193 
194  for(Index jj = 0; jj < nppd; jj++)
195  {
196  if(derivatives_data(jj) == JQT_temperature)
197  {
198  dF[jj][ii] = (Fplus * Fplus * Complex(dG0_dT, dL0_dT) + Fminus * Fminus * Complex(dG0_dT, -dL0_dT)) * PI;
199  }
200  else if(derivatives_data(jj) == JQT_frequency or
201  derivatives_data(jj) == JQT_wind_magnitude or
202  derivatives_data(jj) == JQT_wind_u or
203  derivatives_data(jj) == JQT_wind_v or
204  derivatives_data(jj) == JQT_wind_w)
205  {
206  dF[jj][ii] = (Fplus * Fplus + Fminus * Fminus) * PI;
207  dF[jj][ii] *= Complex(0.0, -1.0);
208  }
209  else if(derivatives_data(jj) == JQT_line_center)
210  {
211  dF[jj][ii] = (Fplus * Fplus * Complex(0.0, 1.0) + Fminus * Fminus * Complex(0.0, -1.0)) * PI;
212  }
213  else if(derivatives_data(jj) == JQT_line_gamma_self or
214  derivatives_data(jj) == JQT_line_gamma_selfexponent or
215  derivatives_data(jj) == JQT_line_pressureshift_self or
216  derivatives_data(jj) == JQT_line_gamma_foreign or
217  derivatives_data(jj) == JQT_line_gamma_foreignexponent or
218  derivatives_data(jj) == JQT_line_pressureshift_foreign or
219  derivatives_data(jj) == JQT_line_gamma_water or
220  derivatives_data(jj) == JQT_line_gamma_waterexponent or
221  derivatives_data(jj) == JQT_line_pressureshift_water) // Only the zeroth order terms --- the derivative with respect to these have to happen later
222  {
223  dF[jj][ii] = (Fplus * Fplus * Complex(1.0, 1.0) + Fminus * Fminus * Complex(1.0, -1.0)) * PI;
224  }
225  else if(derivatives_data(jj) == JQT_magnetic_magntitude)// No external inputs --- errors because of frequency shift when Zeeman is used?
226  {
227  dF[jj][ii] = (Fplus * Fplus * Complex(0.0, zeeman_df) + Fminus * Fminus * Complex(0.0, -zeeman_df)) * PI;
228  }
229  }
230  }
231  }
232 
233  void set_doppler(ComplexVector& F, // Sets the full complex line shape without line mixing
234  ArrayOfComplexVector& dF,
235  const Vector& f_grid,
236  const Numeric& zeeman_df,
237  const Numeric& magnetic_magnitude,
238  const Numeric& F0_noshift,
239  const Numeric& GD_div_F0,
240  const PropmatPartialsData& derivatives_data=PropmatPartialsData(),
241  const Numeric& dGD_div_F0_dT=0.0) const
242  {
243  set_faddeeva_algorithm916(F,
244  dF,
245  f_grid,
246  zeeman_df,
247  magnetic_magnitude,
248  F0_noshift,
249  GD_div_F0, 0.0, 0.0,
250  derivatives_data,
251  dGD_div_F0_dT);
252  }
253 
254  void set_htp(ComplexVector& F, // Sets the full complex line shape without line mixing
255  ArrayOfComplexVector& dF,
256  const Vector& f_grid,
257  const Numeric& zeeman_df,
258  const Numeric& magnetic_magnitude,
259  const Numeric& F0_noshift,
260  const Numeric& GD_div_F0,
261  const Numeric& G0,
262  const Numeric& L0,
263  const Numeric& L2,
264  const Numeric& G2,
265  const Numeric& eta,
266  const Numeric& FVC,
267  const PropmatPartialsData& derivatives_data=PropmatPartialsData(),
268  const Numeric& dGD_div_F0_dT=0.0,
269  const Numeric& dG0_dT=0.0,
270  const Numeric& dL0_dT=0.0,
271  const Numeric& dG2_dT=0.0,
272  const Numeric& dL2_dT=0.0,
273  const Numeric& deta_dT=0.0,
274  const Numeric& dFVC_dT=0.0) const
275  {
276  /*
277  *
278  * Function is meant to compute the Hartmann-Tran lineshape
279  *
280  * The assumptions on input are described in accompanying pdf-documents
281  *
282  * Note that Hartmann-Tran lineshape does not support line mixing in this version
283  *
284  */
285  extern const Numeric PI;
286  static const Complex i(0.0, 1.0), one_plus_one_i(1.0, 1.0);
287  static const Numeric ln2 = log(2.0), sqrtLN2 = sqrt(ln2);
288  static const Numeric sqrtPI = sqrt(PI), invPI = 1.0 / PI, invSqrtPI = sqrt(invPI);
289 
290  // Main lineshape parameters
291  Complex A, B, Zp, Zm, Zp2, Zm2, wiZp, wiZm, X, sqrtXY, invG;
292 
293  // Derivatives parameters
294  Complex dA, dB, dC0, dC0t, dC2, dC2t, dZp, dZm, dwiZp, dwiZm, dX, dY, dg;
295 
296  // Number of frequencies
297  const Index nf = f_grid.nelem();
298 
299  // Number of derivatives
300  const Index nppd = derivatives_data.nelem();
301 
302  // Zeeman effect and other line center shifts means that the true line center is shifted
303  const Numeric F0 = F0_noshift + zeeman_df * magnetic_magnitude;
304  const Numeric F0_ratio = F0_noshift / F0;
305  const Numeric GD = GD_div_F0 * F0;
306  const Numeric one_minus_eta = 1.0 - eta;
307  const Numeric dGD_dT = dGD_div_F0_dT * F0;
308 
309  // Index of derivative for first occurence of similarly natured partial derivatives
310  const Index first_frequency = derivatives_data.get_first_frequency(),
311  first_pressure_broadening = derivatives_data.get_first_pressure_term();
312 
313  // Pressure broadening terms
314  const Complex C0 = G0 + i*L0;
315  const Complex C2 = G2 + i*L2;
316 
317  // Pressure terms adjusted by collisions and correlations
318  const Complex C0_m1p5_C2 = (C0 - 1.5 * C2);
319  const Complex C0t = one_minus_eta * C0_m1p5_C2 + FVC;
320  const Complex invC2t = 1.0 / (one_minus_eta * C2);
321 
322  // Relative pressure broadening terms to be used in the shape function
323  const Complex sqrtY = 0.5 * invC2t * GD / sqrtLN2;
324  const Complex Y = sqrtY * sqrtY;
325 
326  // Scale factor
327  const Numeric const1 = ((sqrtLN2 * sqrtPI) / GD);
328 
329  // Loop over frequencies that cannot be parallelized
330  for(Index ii = 0; ii < nf; ii++)
331  {
332  // Relative frequency
333  X = (C0t - i*(F0 - f_grid[ii])) * invC2t;
334 
335  // Setting up the Z terms
336  sqrtXY = sqrt(X+Y);
337  Zm = Zp = sqrtXY;
338  Zm -= sqrtY;
339  Zp += sqrtY;
340 
341  // The shape functions
342  wiZm = Faddeeva::w(i*Zm);
343  wiZp = Faddeeva::w(i*Zp);
344 
345  // Main line shape is computed here --- WARNING when wiZM and wiZp are too close, this could lead to numerical errors --- should be caught before frequency loop
346  A = const1 * (wiZm - wiZp);
347 
348  // If there are correlations or if there is a temperature dependency
349  if(eta not_eq 0.0 or deta_dT not_eq 0.0)
350  {
351  Zm2 = Zm * Zm;
352  Zp2 = Zp * Zp;
353  B = (1.0 / one_minus_eta) * (-1.0 + 0.5 * sqrtPI / sqrtY * ((1.0 - Zm2)*wiZm - (1.0 - Zp2)*wiZp));
354  invG = 1.0 / (1.0 - (FVC - eta * C0_m1p5_C2)*A + eta * B);
355  }
356  else
357  {
358  invG = 1.0 / (1.0 - FVC*A); // WARNING if FVC x A == 1 then this fail --- no way around it, though A should be very small and FVC should not be too large
359  }
360  F[ii] = A * invG * invPI;
361 
362  for(Index jj = 0; jj < nppd; jj++)
363  {
364  if(derivatives_data(jj) == JQT_frequency or
365  derivatives_data(jj) == JQT_wind_magnitude or
366  derivatives_data(jj) == JQT_wind_u or
367  derivatives_data(jj) == JQT_wind_v or
368  derivatives_data(jj) == JQT_wind_w) // No external inputs
369  {
370  // If this is the first time it is calculated this frequency bin, do the full calculation
371  if(first_frequency == jj)
372  {
373  dX = invC2t * i;
374  dZp = dZm = 0.5 * dX / sqrtXY;
375 
376  dwiZm = dwiZp = i * invSqrtPI;
377  dwiZm -= Zm * wiZm;
378  dwiZp -= Zp * wiZp;
379  dwiZm *= 2.0 * dZm;
380  dwiZp *= 2.0 * dZp;
381 
382  dA = const1 * (dwiZm - dwiZp);
383 
384  if(eta not_eq 0.0)
385  {
386  dB = 0.5 * sqrtPI / (sqrtY * one_minus_eta) * ((1.0 - Zm2) * dwiZm - 2.0 * Zm * dZm * wiZm + (1.0 - Zp2) * dwiZp + 2.0 * Zp * dZp * wiZp);
387  dg = eta * (C0_m1p5_C2 * dA - dB) - FVC * dA;
388  }
389  else
390  {
391  dg = - FVC * dA;
392  }
393 
394  dF[jj][ii] = invG * (invPI * dA - F[ii] * dg);
395  }
396  else // copy for repeated occurences
397  {
398  dF[jj][ii] = dF[first_frequency][ii];
399  }
400  }
401  else if(derivatives_data(jj) == JQT_temperature)
402  {
403  dC0 = dG0_dT + i*dL0_dT;
404  dC2 = dG2_dT + i*dL2_dT;
405 
406  dC0t = one_minus_eta * (dC0 - 1.5 * dC2) - deta_dT * C0_m1p5_C2 + dFVC_dT;
407  dC2t = one_minus_eta * dC2 - deta_dT * C2;
408 
409  dY = GD / (2.0 * ln2) * invC2t*invC2t * (dGD_dT - GD * invC2t * dC2t);
410  dX = invC2t * (dC0t - X*dC2t);
411 
412  dZm = dZp = 0.5 * (dX + dY) / sqrtXY;
413  dZm -= 0.5 / sqrtY * dY;
414  dZp += 0.5 / sqrtY * dY;
415 
416  dwiZm = dwiZp = i * invSqrtPI;
417  dwiZm -= Zm * wiZm;
418  dwiZp -= Zp * wiZp;
419  dwiZm *= 2.0 * dZm;
420  dwiZp *= 2.0 * dZp;
421 
422  dA = const1 * (dwiZm - dwiZp - A * GD_div_F0);
423 
424  if(eta not_eq 0)
425  {
426  dB = 1.0 / one_minus_eta * (1.0 / one_minus_eta * deta_dT + 0.5 * sqrtPI / sqrtY * ( - 0.5 * ((1.0 - Zm2) * wiZm - (1.0 - Zp2) * wiZp) / Y * dY + (1.0 - Zm2) * dwiZm - 2.0 * Zm * dZm * wiZm + (1.0 - Zp2) * dwiZp + 2.0 * Zp * dZp * wiZp));
427  dg = (dFVC_dT - deta_dT * C0_m1p5_C2 - eta * (dC0 - 1.5 * dC2)) * A - (FVC - eta * C0_m1p5_C2) * dA + deta_dT * B + eta * dB;
428  }
429  else
430  {
431  dg = (dFVC_dT - deta_dT * C0_m1p5_C2) * A - FVC * dA + deta_dT * B;
432  }
433 
434  dF[jj][ii] = invG * (invPI * dA - F[ii] * dg);
435  }
436  else if(derivatives_data(jj) == JQT_line_center) // No external inputs --- errors because of frequency shift when Zeeman is used?
437  {
438  dY = GD / (2.0 * ln2) * invC2t*invC2t * GD_div_F0 * F0_ratio;
439  dX = - i * invC2t * F0_ratio;
440 
441  dZm = dZp = 0.5 * (dX + dY) / sqrtXY;
442  dZm -= 0.5 / sqrtY * dY;
443  dZp += 0.5 / sqrtY * dY;
444 
445  dwiZm = dwiZp = i * invSqrtPI;
446  dwiZm -= Zm * wiZm;
447  dwiZp -= Zp * wiZp;
448  dwiZm *= 2.0 * dZm;
449  dwiZp *= 2.0 * dZp;
450 
451  dA = const1 * (dwiZm - dwiZp - A * GD_div_F0);
452 
453  if(eta not_eq 0.0)
454  {
455  dB = 0.5 * sqrtPI / (sqrtY * one_minus_eta) * (- 0.5 * ((1.0 - Zm2) * wiZm - (1.0 - Zp2) * wiZp) / Y * dY + (1.0 - Zm2) * dwiZm - 2.0 * Zm * dZm * wiZm + (1.0 - Zp2) * dwiZp + 2.0 * Zp * dZp * wiZp);
456  dg = eta * (C0_m1p5_C2 * dA - dB) - FVC * dA;
457  }
458  else
459  {
460  dg = - FVC * dA;
461  }
462 
463  dF[jj][ii] = invG * (invPI * dA - F[ii] * dg);
464  }
465  else if(derivatives_data(jj) == JQT_line_gamma_self or
466  derivatives_data(jj) == JQT_line_gamma_selfexponent or
467  derivatives_data(jj) == JQT_line_pressureshift_self or
468  derivatives_data(jj) == JQT_line_gamma_foreign or
469  derivatives_data(jj) == JQT_line_gamma_foreignexponent or
470  derivatives_data(jj) == JQT_line_pressureshift_foreign or
471  derivatives_data(jj) == JQT_line_gamma_water or
472  derivatives_data(jj) == JQT_line_gamma_waterexponent or
473  derivatives_data(jj) == JQT_line_pressureshift_water) // Only the zeroth order terms --- the derivative with respect to these have to happen later
474  {
475  // Note that if the species vmr partial derivative is necessary here is where it goes
476  if(first_pressure_broadening == jj)
477  {
478  dC0t = one_minus_eta * one_plus_one_i;
479  dX = invC2t * dC0t;
480  dZm = dZp = 0.5 * dX / sqrtXY;
481 
482  dwiZm = dwiZp = i * invSqrtPI;
483  dwiZm -= Zm * wiZm;
484  dwiZp -= Zp * wiZp;
485  dwiZm *= 2.0 * dZm;
486  dwiZp *= 2.0 * dZp;
487 
488  dA = const1 * (dwiZm - dwiZp);
489 
490  if(eta not_eq 0.0)
491  {
492  dB = 0.5 * sqrtPI / (sqrtY * one_minus_eta) * ((1.0 - Zm2) * dwiZm - 2.0 * Zm * dZm * wiZm + (1.0 - Zp2) * dwiZp + 2.0 * Zp * dZp * wiZp);
493  dg = eta * (C0_m1p5_C2 * dA - dB - dC0 * A) - FVC * dA;
494  }
495  else
496  {
497  dg = - FVC * dA;
498  }
499 
500  dF[jj][ii] = invG * (invPI * dA - F[ii] * dg);
501  }
502  else // copy for repeated occurences
503  {
504  dF[jj][ii] = dF[first_frequency][ii];
505  }
506  }
507  else if(derivatives_data(jj) == JQT_magnetic_magntitude)// No external inputs --- errors because of frequency shift when Zeeman is used?
508  {
509  dY = GD / (2.0 * ln2) * invC2t*invC2t * GD_div_F0 * (1.0 - F0_ratio) / magnetic_magnitude;
510  dX = - i * invC2t * (1.0 - F0_ratio) / magnetic_magnitude;
511 
512  dZm = dZp = 0.5 * (dX + dY) / sqrtXY;
513  dZm -= 0.5 / sqrtY * dY;
514  dZp += 0.5 / sqrtY * dY;
515 
516  dwiZm = dwiZp = i * invSqrtPI;
517  dwiZm -= Zm * wiZm;
518  dwiZp -= Zp * wiZp;
519  dwiZm *= 2.0 * dZm;
520  dwiZp *= 2.0 * dZp;
521 
522  dA = const1 * (dwiZm - dwiZp - A * GD_div_F0);
523 
524  if(eta not_eq 0.0)
525  {
526  dB = 0.5 * sqrtPI / (sqrtY * one_minus_eta) *
527  (- 0.5 * ((1.0 - Zm2) * wiZm - (1.0 - Zp2) * wiZp) / Y * dY +
528  (1.0 - Zm2) * dwiZm - 2.0 * Zm * dZm * wiZm +
529  (1.0 - Zp2) * dwiZp + 2.0 * Zp * dZp * wiZp);
530  dg = eta * (C0_m1p5_C2 * dA - dB) - FVC * dA;
531  }
532  else
533  {
534  dg = - FVC * dA;
535  }
536 
537  dF[jj][ii] = invG * (invPI * dA - F[ii] * dg);
538  }
539  }
540  }
541  }
542 
543  void set_faddeeva_algorithm916(ComplexVector& F, // Sets the full complex line shape without line mixing
544  ArrayOfComplexVector& dF,
545  const Vector& f_grid,
546  const Numeric& zeeman_df,
547  const Numeric& magnetic_magnitude,
548  const Numeric& F0_noshift,
549  const Numeric& GD_div_F0,
550  const Numeric& G0,
551  const Numeric& L0,
552  const PropmatPartialsData& derivatives_data=PropmatPartialsData(),
553  const Numeric& dGD_div_F0_dT=0.0,
554  const Numeric& dG0_dT=0.0,
555  const Numeric& dL0_dT=0.0) const
556  {
557  // For calculations
558  Numeric dx;
559  Complex w, z, dw_over_dz, dz;
560 
561  const Index nf = f_grid.nelem(), nppd = derivatives_data.nelem();
562 
563  // PI
564  extern const Numeric PI;
565 
566  // constant sqrt(1/pi)
567  static const Numeric sqrt_invPI = sqrt(1.0/PI);
568 
569  // Doppler broadening and line center
570  const Numeric F0 = F0_noshift + zeeman_df * magnetic_magnitude + L0;
571  const Numeric GD = GD_div_F0 * F0;
572  const Numeric invGD = 1.0 / GD;
573  const Numeric dGD_dT = dGD_div_F0_dT * F0;
574 
575  // constant normalization factor for voigt
576  const Numeric fac = sqrt_invPI * invGD;
577  F = fac;
578 
579  // Ratio of the Lorentz halfwidth to the Doppler halfwidth
580  const Complex z0 = Complex(-F0, G0) * invGD;
581 
582  const Index first_pressure_broadening = derivatives_data.get_first_pressure_term(),
583  first_frequency = derivatives_data.get_first_frequency();
584 
585  // frequency in units of Doppler
586  for (Index ii=0; ii<nf; ii++)
587  {
588  dx = f_grid[ii] * invGD;
589  z = z0 + dx;
590  w = Faddeeva::w(z);
591 
592  F[ii] *= w;
593 
594  for(Index jj = 0; jj < nppd; jj++)
595  {
596  if(jj==0)
597  dw_over_dz = 2.0 * (z * w - sqrt_invPI);
598 
599  if(derivatives_data(jj) == JQT_frequency or
600  derivatives_data(jj) == JQT_wind_magnitude or
601  derivatives_data(jj) == JQT_wind_u or
602  derivatives_data(jj) == JQT_wind_v or
603  derivatives_data(jj) == JQT_wind_w) // No external inputs
604  {
605  // If this is the first time it is calculated this frequency bin, do the full calculation
606  if(first_frequency == jj)
607  {
608  //dz = Complex(invGD, 0.0);
609 
610  dF[jj][ii] = fac * dw_over_dz * invGD; //dz;
611  }
612  else // copy for repeated occurences
613  {
614  dF[jj][ii] = dF[first_frequency][ii];
615  }
616  }
617  else if(derivatives_data(jj) == JQT_temperature)
618  {
619  dz = Complex(-dL0_dT, dG0_dT) - z * dGD_dT;
620 
621  dF[jj][ii] = -F[ii] * dGD_dT;
622  dF[jj][ii] += fac * dw_over_dz * dz;
623  dF[jj][ii] *= invGD;
624  }
625  else if(derivatives_data(jj) == JQT_line_center) // No external inputs --- errors because of frequency shift when Zeeman is used?
626  {
627  dz = -z * GD_div_F0 - 1.0;
628 
629  dF[jj][ii] = -F[ii] * GD_div_F0;
630  dF[jj][ii] += dw_over_dz * dz;
631  dF[jj][ii] *= fac * invGD;
632  }
633  else if(derivatives_data(jj) == JQT_line_gamma_self or
634  derivatives_data(jj) == JQT_line_gamma_selfexponent or
635  derivatives_data(jj) == JQT_line_pressureshift_self or
636  derivatives_data(jj) == JQT_line_gamma_foreign or
637  derivatives_data(jj) == JQT_line_gamma_foreignexponent or
638  derivatives_data(jj) == JQT_line_pressureshift_foreign or
639  derivatives_data(jj) == JQT_line_gamma_water or
640  derivatives_data(jj) == JQT_line_gamma_waterexponent or
641  derivatives_data(jj) == JQT_line_pressureshift_water) // Only the zeroth order terms --- the derivative with respect to these have to happen later
642  {
643  // Note that if the species vmr partial derivative is necessary here is where it goes
644  if(first_pressure_broadening == jj)
645  {
646  dz = Complex(-1.0, 1.0) * invGD;
647  dF[jj][ii] = fac * dw_over_dz * dz;
648  }
649  else
650  {
651  dF[jj][ii] = dF[first_frequency][ii];
652  }
653  }
654  else if(derivatives_data(jj) == JQT_magnetic_magntitude)// No external inputs --- errors because of frequency shift when Zeeman is used?
655  {
656  // dz = Complex(- zeeman_df * invGD, 0.0);
657 
658  dF[jj][ii] = fac * dw_over_dz * (- zeeman_df * invGD); //* dz;
659  }
660  else if(derivatives_data(jj) == JQT_line_mixing_DF or
661  derivatives_data(jj) == JQT_line_mixing_DF0 or
662  derivatives_data(jj) == JQT_line_mixing_DF1 or
663  derivatives_data(jj) == JQT_line_mixing_DFexp)
664  {
665  // dz = Complex(-invGD, 0.0);
666 
667  dF[jj][ii] = fac * dw_over_dz * (-invGD); //* dz;
668  }
669  }
670 
671  }
672  }
673 
674  void set_faddeeva_from_full_linemixing(ComplexVector& F, // Sets the full complex line shape without line mixing
675  ArrayOfComplexVector& dF,
676  const Vector& f_grid,
677  const Complex& eigenvalue_no_shift,
678  const Numeric& GD_div_F0,
679  const Numeric& L0,
680  const PropmatPartialsData& derivatives_data=PropmatPartialsData(),
681  const Numeric& dGD_div_F0_dT=0.0,
682  const Complex& deigenvalue_dT=0.0,
683  const Numeric& dL0_dT=0.0) const
684  {
685  // For calculations
686  Numeric dx;
687  Complex w, z, dw_over_dz, dz;
688 
689  const Index nf = f_grid.nelem(), nppd = derivatives_data.nelem();
690 
691  // PI
692  extern const Numeric PI;
693 
694  // constant sqrt(1/pi)
695  static const Numeric sqrt_invPI = sqrt(1.0/PI);
696 
697  // Doppler broadening and line center
698  const Complex eigenvalue = eigenvalue_no_shift + L0;
699  const Numeric F0 = eigenvalue.real() + L0;
700  const Numeric GD = GD_div_F0 * F0;
701  const Numeric invGD = 1.0 / GD;
702  const Numeric dGD_dT = dGD_div_F0_dT * F0;
703 
704  // constant normalization factor for voigt
705  const Numeric fac = sqrt_invPI * invGD;
706  F = fac;
707 
708  // Ratio of the Lorentz halfwidth to the Doppler halfwidth
709  const Complex z0 = -eigenvalue * invGD;
710 
711  const Index first_pressure_broadening = derivatives_data.get_first_pressure_term(),
712  first_frequency = derivatives_data.get_first_frequency();
713 
714  // frequency in units of Doppler
715  for (Index ii=0; ii<nf; ii++)
716  {
717  dx = f_grid[ii] * invGD;
718  z = z0 + dx;
719  w = Faddeeva::w(z);
720 
721  F[ii] *= w;
722 
723  for(Index jj = 0; jj < nppd; jj++)
724  {
725  if(jj==0)
726  dw_over_dz = 2.0 * (z * w - sqrt_invPI);
727 
728  if(derivatives_data(jj) == JQT_frequency or
729  derivatives_data(jj) == JQT_wind_magnitude or
730  derivatives_data(jj) == JQT_wind_u or
731  derivatives_data(jj) == JQT_wind_v or
732  derivatives_data(jj) == JQT_wind_w) // No external inputs
733  {
734  // If this is the first time it is calculated this frequency bin, do the full calculation
735  if(first_frequency == jj)
736  {
737  //dz = Complex(invGD, 0.0);
738 
739  dF[jj][ii] = fac * dw_over_dz * invGD; //dz;
740  }
741  else // copy for repeated occurences
742  {
743  dF[jj][ii] = dF[first_frequency][ii];
744  }
745  }
746  else if(derivatives_data(jj) == JQT_temperature)
747  {
748  dz = (deigenvalue_dT - dL0_dT) - z * dGD_dT;
749 
750  dF[jj][ii] = -F[ii] * dGD_dT;
751  dF[jj][ii] += fac * dw_over_dz * dz;
752  dF[jj][ii] *= invGD;
753  }
754  else if(derivatives_data(jj) == JQT_line_center) // No external inputs --- errors because of frequency shift when Zeeman is used?
755  {
756  dz = -z * GD_div_F0 - 1.0;
757 
758  dF[jj][ii] = -F[ii] * GD_div_F0;
759  dF[jj][ii] += dw_over_dz * dz;
760  dF[jj][ii] *= fac * invGD;
761  }
762  else if(derivatives_data(jj) == JQT_line_gamma_self or
763  derivatives_data(jj) == JQT_line_gamma_selfexponent or
764  derivatives_data(jj) == JQT_line_pressureshift_self or
765  derivatives_data(jj) == JQT_line_gamma_foreign or
766  derivatives_data(jj) == JQT_line_gamma_foreignexponent or
767  derivatives_data(jj) == JQT_line_pressureshift_foreign or
768  derivatives_data(jj) == JQT_line_gamma_water or
769  derivatives_data(jj) == JQT_line_gamma_waterexponent or
770  derivatives_data(jj) == JQT_line_pressureshift_water) // Only the zeroth order terms --- the derivative with respect to these have to happen later
771  {
772  // Note that if the species vmr partial derivative is necessary here is where it goes
773  if(first_pressure_broadening == jj)
774  {
775  dz = Complex(-1.0, 1.0) * invGD;
776  dF[jj][ii] = fac * dw_over_dz * dz;
777  }
778  else
779  {
780  dF[jj][ii] = dF[first_frequency][ii];
781  }
782  }
783  else if(derivatives_data(jj) == JQT_line_mixing_DF or
784  derivatives_data(jj) == JQT_line_mixing_DF0 or
785  derivatives_data(jj) == JQT_line_mixing_DF1 or
786  derivatives_data(jj) == JQT_line_mixing_DFexp)
787  {
788  // dz = Complex(-invGD, 0.0);
789 
790  dF[jj][ii] = fac * dw_over_dz * (-invGD); //* dz;
791  }
792  }
793 
794  }
795  }
796 
797  void set_hui_etal_1978(ComplexVector& F, // Sets the full complex line shape without line mixing
798  ArrayOfComplexVector& dF,
799  const Vector& f_grid,
800  const Numeric& zeeman_df,
801  const Numeric& magnetic_magnitude,
802  const Numeric& F0_noshift,
803  const Numeric& GD_div_F0,
804  const Numeric& G0,
805  const Numeric& L0,
806  const PropmatPartialsData& derivatives_data=PropmatPartialsData(),
807  const Numeric& dGD_div_F0_dT=0.0,
808  const Numeric& dG0_dT=0.0,
809  const Numeric& dL0_dT=0.0) const
810  {
811  /*
812  * Exact copy of get_faddeeva_algorithm916 but with the change of A and B calculations from empirical numbers
813  *
814  * For future extensions and corrections, make sure get_faddeeva_algorithm916 first works and then update this
815  */
816 
817  // For calculations
818  Numeric dx;
819  Complex w, z, dw_over_dz, dz, A, B;
820 
821  const Index nf = f_grid.nelem(), nppd = derivatives_data.nelem();
822 
823  // PI
824  extern const Numeric PI;
825 
826  // constant sqrt(1/pi)
827  static const Numeric sqrt_invPI = sqrt(1/PI);
828 
829  // Doppler broadening and line center
830  const Numeric F0 = F0_noshift + zeeman_df * magnetic_magnitude + L0;
831  const Numeric GD = GD_div_F0 * F0;
832  const Numeric invGD = 1.0 / GD;
833  const Numeric dGD_dT = dGD_div_F0_dT * F0;
834 
835  // constant normalization factor for voigt
836  const Numeric fac = sqrt_invPI * invGD;
837 
838  // Ratio of the Lorentz halfwidth to the Doppler halfwidth
839  const Complex z0 = Complex(-F0, G0) * invGD;
840 
841  const Index first_pressure_broadening = derivatives_data.get_first_pressure_term(),
842  first_frequency = derivatives_data.get_first_frequency();
843 
844  // frequency in units of Doppler
845  for (Index ii=0; ii<nf; ii++)
846  {
847  dx = f_grid[ii] * invGD;
848  z = z0 + dx;
849  /* Since this is ported from old FORTRAN code, intention must be that A and B coeffs are floats?
850  A = (((((.5641896*z+5.912626)*z+30.18014)*z+
851  93.15558)*z+181.9285)*z+214.3824)*z+122.6079;
852  B = ((((((z+10.47986)*z+53.99291)*z+170.3540)*z+
853  348.7039)*z+457.3345)*z+352.7306)*z+122.6079;
854  */
855  A = (((((.5641896f*z+5.912626f)*z+30.18014f)*z+
856  93.15558f)*z+181.9285f)*z+214.3824f)*z+122.6079f;
857  B = ((((((z+10.47986f)*z+53.99291f)*z+170.3540f)*z+
858  348.7039f)*z+457.3345f)*z+352.7306f)*z+122.6079f;
859  w = A / B;
860 
861  F[ii] = fac * w;
862 
863  for(Index jj = 0; jj < nppd; jj++)
864  {
865  if(jj==0)
866  dw_over_dz = 2.0 * (z * w - sqrt_invPI);
867 
868  if(derivatives_data(jj) == JQT_frequency or
869  derivatives_data(jj) == JQT_wind_magnitude or
870  derivatives_data(jj) == JQT_wind_u or
871  derivatives_data(jj) == JQT_wind_v or
872  derivatives_data(jj) == JQT_wind_w) // No external inputs
873  {
874  // If this is the first time it is calculated this frequency bin, do the full calculation
875  if(first_frequency == jj)
876  {
877  //dz = Complex(invGD, 0.0);
878 
879  dF[jj][ii] = fac * dw_over_dz * invGD; //dz;
880  }
881  else // copy for repeated occurences
882  {
883  dF[jj][ii] = dF[first_frequency][ii];
884  }
885  }
886  else if(derivatives_data(jj) == JQT_temperature)
887  {
888  dz = Complex(-dL0_dT, dG0_dT) - z * dGD_dT;
889 
890  dF[jj][ii] = -F[ii] * dGD_dT;
891  dF[jj][ii] += fac * dw_over_dz * dz;
892  dF[jj][ii] *= invGD;
893  }
894  else if(derivatives_data(jj) == JQT_line_center) // No external inputs --- errors because of frequency shift when Zeeman is used?
895  {
896  dz = -z * GD_div_F0 - 1.0;
897 
898  dF[jj][ii] = -F[ii] * GD_div_F0;
899  dF[jj][ii] += dw_over_dz * dz;
900  dF[jj][ii] *= fac * invGD;
901  }
902  else if(derivatives_data(jj) == JQT_line_gamma_self or
903  derivatives_data(jj) == JQT_line_gamma_selfexponent or
904  derivatives_data(jj) == JQT_line_pressureshift_self or
905  derivatives_data(jj) == JQT_line_gamma_foreign or
906  derivatives_data(jj) == JQT_line_gamma_foreignexponent or
907  derivatives_data(jj) == JQT_line_pressureshift_foreign or
908  derivatives_data(jj) == JQT_line_gamma_water or
909  derivatives_data(jj) == JQT_line_gamma_waterexponent or
910  derivatives_data(jj) == JQT_line_pressureshift_water) // Only the zeroth order terms --- the derivative with respect to these have to happen later
911  {
912  // Note that if the species vmr partial derivative is necessary here is where it goes
913  if(first_pressure_broadening == jj)
914  {
915  dz = Complex(-1.0, 1.0) * invGD;
916  dF[jj][ii] = fac * dw_over_dz * dz;
917  }
918  else
919  {
920  dF[jj][ii] = dF[first_frequency][ii];
921  }
922  }
923  else if(derivatives_data(jj) == JQT_magnetic_magntitude)// No external inputs --- errors because of frequency shift when Zeeman is used?
924  {
925  // dz = Complex(- zeeman_df * invGD, 0.0);
926 
927  dF[jj][ii] = fac * dw_over_dz * (- zeeman_df * invGD); //* dz;
928  }
929  else if(derivatives_data(jj) == JQT_line_mixing_DF or
930  derivatives_data(jj) == JQT_line_mixing_DF0 or
931  derivatives_data(jj) == JQT_line_mixing_DF1 or
932  derivatives_data(jj) == JQT_line_mixing_DFexp)
933  {
934  // dz = Complex(-invGD, 0.0);
935 
936  dF[jj][ii] = fac * dw_over_dz * (-invGD); //* dz;
937  }
938  }
939  }
940  }
941 
942  void set_o2_non_resonant(ComplexVector& F, // Sets the full complex line shape without line mixing
943  ArrayOfComplexVector& dF,
944  const Vector& f_grid,
945  const Numeric& F0,
946  const Numeric& G0,
947  const PropmatPartialsData& derivatives_data=PropmatPartialsData(),
948  const Numeric& dG0_dT=0.0)
949  {
950  const Index nf = f_grid.nelem(), nppd = derivatives_data.nelem();
951 
952  extern const Numeric PI;
953 
954  const Numeric fac = G0/PI, invG0 = 1.0/G0;
955 
956  for(Index ii = 0; ii < nf; ii++)
957  {
958  F[ii] = fac / ((f_grid[ii]-F0) * (f_grid[ii]-F0));
959  }
960 
961  for(Index jj = 0; jj < nppd; jj++)
962  {
963  if(derivatives_data(jj) == JQT_temperature)
964  {
965  dF[jj] = F;
966  dF[jj] *= dG0_dT * invG0;
967  }
968  else if(derivatives_data(jj) == JQT_frequency or
969  derivatives_data(jj) == JQT_wind_magnitude or
970  derivatives_data(jj) == JQT_wind_u or
971  derivatives_data(jj) == JQT_wind_v or
972  derivatives_data(jj) == JQT_wind_w)
973  {
974  dF[jj] = F;
975  for(Index ii = 0; ii < nf; ii++)
976  {
977  dF[jj][ii] /= -2.0 / (f_grid[ii] - F0);
978  }
979  }
980  else if(derivatives_data(jj) == JQT_line_center)
981  {
982  dF[jj] = F;
983  for(Index ii = 0; ii < nf; ii++)
984  {
985  dF[jj][ii] *= 2.0 / (f_grid[ii] - F0);
986  }
987  }
988  else if(derivatives_data(jj) == JQT_line_gamma_self or
989  derivatives_data(jj) == JQT_line_gamma_selfexponent or
990  derivatives_data(jj) == JQT_line_gamma_foreign or
991  derivatives_data(jj) == JQT_line_gamma_foreignexponent or
992  derivatives_data(jj) == JQT_line_gamma_water or
993  derivatives_data(jj) == JQT_line_gamma_waterexponent)
994  {
995  dF[jj] = F;
996  dF[jj] *= invG0;
997  }
998  }
999  }
1000 
1001  void apply_linemixing(ComplexVector& F, // Returns the full complex or normalized line shape with line mixing
1002  ArrayOfComplexVector& dF,
1003  const Numeric& Y,
1004  const Numeric& G,
1005  const PropmatPartialsData& derivatives_data=PropmatPartialsData(),
1006  const Numeric& dY_dT=0.0,
1007  const Numeric& dG_dT=0.0) const
1008  {
1009  const Index nf = F.nelem(), nppd = derivatives_data.nelem();
1010 
1011  const Complex LM = Complex(1.0 + G, Y);
1012  const Complex dLM_dT = Complex(dG_dT, dY_dT);
1013 
1014  for(Index iq = 0; iq < nppd; iq++)
1015  {
1016  if(derivatives_data(iq) == JQT_temperature)
1017  {
1018  dF[iq] *= LM;
1019  for(Index ii = 0; ii < nf; ii++)
1020  {
1021  dF[iq][ii] += F[ii] * dLM_dT;
1022  }
1023  }
1024  else if(derivatives_data(iq) == JQT_line_mixing_Y or
1025  derivatives_data(iq) == JQT_line_mixing_Y0 or
1026  derivatives_data(iq) == JQT_line_mixing_Y1 or
1027  derivatives_data(iq) == JQT_line_mixing_Yexp)
1028  {
1029  dF[iq] = F;
1030  dF[iq] *= Complex(0.0, 1.0);
1031  }
1032  else if(derivatives_data(iq) == JQT_line_mixing_G or
1033  derivatives_data(iq) == JQT_line_mixing_G0 or
1034  derivatives_data(iq) == JQT_line_mixing_G1 or
1035  derivatives_data(iq) == JQT_line_mixing_Gexp)
1036  {
1037  dF[iq] = F;
1038  }
1039  else
1040  {
1041  dF[iq] *= LM;
1042  }
1043  }
1044 
1045  F *= LM;
1046  }
1047 
1048  void apply_rosenkranz_quadratic(ComplexVector& F, // Returns as normalized complex line shape with or without line mixing
1049  ArrayOfComplexVector& dF,
1050  const Vector& f_grid,
1051  const Numeric& F0, // Only line center without any shifts
1052  const Numeric& T,
1053  const PropmatPartialsData& derivatives_data=PropmatPartialsData()) const
1054  {
1055  const Index nf = f_grid.nelem(), nppd = derivatives_data.nelem();
1056 
1057  extern const Numeric PLANCK_CONST;
1058  extern const Numeric BOLTZMAN_CONST;
1059 
1060  const Numeric invF0 = 1.0/F0;
1061  const Numeric mafac = (PLANCK_CONST) / (2.000e0 * BOLTZMAN_CONST * T) /
1062  sinh((PLANCK_CONST * F0) / (2.000e0 * BOLTZMAN_CONST * T)) * invF0;
1063 
1064  Numeric dmafac_dF0_div_fun, dmafac_dT_div_fun;
1065  if(derivatives_data.do_line_center())
1066  {
1067  dmafac_dF0_div_fun = -invF0 -
1068  PLANCK_CONST/(2.0*BOLTZMAN_CONST*T*tanh(F0*PLANCK_CONST/(2.0*BOLTZMAN_CONST*T)));
1069  }
1070  if(derivatives_data.do_temperature())
1071  {
1072  dmafac_dT_div_fun = -(BOLTZMAN_CONST*T - F0*PLANCK_CONST/
1073  (2.0*tanh(F0*PLANCK_CONST/(2.0*BOLTZMAN_CONST*T))))/(BOLTZMAN_CONST*T*T);
1074  }
1075 
1076  Numeric fun;
1077 
1078  for (Index ii=0; ii < nf; ii++)
1079  {
1080  fun = mafac * (f_grid[ii] * f_grid[ii]);
1081  F[ii] *= fun;
1082 
1083  for(Index jj = 0; jj < nppd; jj++)
1084  {
1085  dF[jj][ii] *= fun;
1086  if(derivatives_data(jj) == JQT_temperature)
1087  {
1088  dF[jj][ii] += dmafac_dT_div_fun * F[ii];
1089  }
1090  else if(derivatives_data(jj) == JQT_line_center)
1091  {
1092  dF[jj][ii] += dmafac_dF0_div_fun * F[ii];
1093  }
1094  else if(derivatives_data(jj) == JQT_frequency or
1095  derivatives_data(jj) == JQT_wind_magnitude or
1096  derivatives_data(jj) == JQT_wind_u or
1097  derivatives_data(jj) == JQT_wind_v or
1098  derivatives_data(jj) == JQT_wind_w)
1099  {
1100  dF[jj][ii] += (2.0 / f_grid[ii]) * F[ii];
1101  }
1102  }
1103  }
1104  }
1105 
1106  void apply_VVH(ComplexVector& F, // Returns as normalized complex line shape with or without line mixing
1107  ArrayOfComplexVector& dF,
1108  const Vector& f_grid,
1109  const Numeric& F0, // Only line center without any shifts
1110  const Numeric& T,
1111  const PropmatPartialsData& derivatives_data=PropmatPartialsData()) const
1112  {
1113  extern const Numeric PLANCK_CONST;
1114  extern const Numeric BOLTZMAN_CONST;
1115 
1116  const Index nf = f_grid.nelem(), nppd = derivatives_data.nelem();
1117 
1118  // 2kT is constant for the loop
1119  const Numeric kT = 2.0 * BOLTZMAN_CONST * T;
1120 
1121  // denominator is constant for the loop
1122  const Numeric absF0 = abs(F0);
1123  const Numeric tanh_f0part = tanh(PLANCK_CONST * absF0 / kT);
1124  const Numeric denom = absF0 * tanh_f0part;
1125 
1126  Numeric fac_df0;
1127  if(derivatives_data.do_line_center())
1128  fac_df0 = (-1.0/absF0 + PLANCK_CONST*tanh_f0part/(kT) - PLANCK_CONST/(kT*tanh_f0part)) * F0/absF0;
1129 
1130  for(Index ii=0; ii < nf; ii++)
1131  {
1132  const Numeric tanh_fpart = tanh( PLANCK_CONST * f_grid[ii] / kT );
1133  const Numeric fun = f_grid[ii] * tanh_fpart / denom;
1134  F[ii] *= fun;
1135 
1136  for(Index jj = 0; jj < nppd; jj++)
1137  {
1138  dF[jj][ii] *= fun;
1139  if(derivatives_data(jj) == JQT_temperature)
1140  {
1141  dF[jj][ii] += (-PLANCK_CONST*(denom - absF0/tanh_f0part -
1142  f_grid[ii]*tanh_fpart + f_grid[ii]/tanh_fpart)/(kT*T)) * F[ii];
1143  }
1144  else if(derivatives_data(jj) == JQT_line_center)
1145  {
1146  dF[jj][ii] += fac_df0 * F[ii];
1147  }
1148  else if(derivatives_data(jj) == JQT_frequency or
1149  derivatives_data(jj) == JQT_wind_magnitude or
1150  derivatives_data(jj) == JQT_wind_u or
1151  derivatives_data(jj) == JQT_wind_v or
1152  derivatives_data(jj) == JQT_wind_w)
1153  {
1154  dF[jj][ii] += (1.0/f_grid[ii] -PLANCK_CONST*tanh_fpart/kT + PLANCK_CONST/(kT*tanh_fpart)) * F[ii];
1155  }
1156  }
1157  }
1158  }
1159 
1160  void apply_VVW(ComplexVector& F, // Returns as normalized line shape with or without line mixing
1161  ArrayOfComplexVector& dF,
1162  const Vector& f_grid,
1163  const Numeric& F0, // Only line center without any shifts
1164  const PropmatPartialsData& derivatives_data=PropmatPartialsData()) const
1165  {
1166  const Index nf = f_grid.nelem(), nppd = derivatives_data.nelem();
1167 
1168  // denominator is constant for the loop
1169  const Numeric invF0 = 1.0 / F0;
1170  const Numeric invF02 = invF0 * invF0;
1171 
1172  for(Index ii = 0; ii < nf; ii++)
1173  {
1174  const Numeric fun = f_grid[ii] * f_grid[ii] * invF02;
1175  F[ii] *= fun;
1176 
1177  for(Index jj = 0; jj < nppd; jj++)
1178  {
1179  dF[jj][ii] *= fun;
1180  if(derivatives_data(jj) == JQT_line_center)
1181  {
1182  dF[jj][ii] -= 2.0 * invF0 * F[ii];
1183  }
1184  else if(derivatives_data(jj) == JQT_frequency or
1185  derivatives_data(jj) == JQT_wind_magnitude or
1186  derivatives_data(jj) == JQT_wind_u or
1187  derivatives_data(jj) == JQT_wind_v or
1188  derivatives_data(jj) == JQT_wind_w)
1189  {
1190  dF[jj][ii] += 2.0 * f_grid[ii] * F[ii];
1191  }
1192  }
1193  }
1194  }
1195 
1196  void apply_linestrength(ComplexVector& F, // Returns the full complex line shape with or without line mixing
1197  ArrayOfComplexVector& dF,
1198  const Numeric& S0,
1199  const Numeric& isotopic_ratio,
1200  const Numeric& QT,
1201  const Numeric& QT0,
1202  const Numeric& K1,
1203  const Numeric& K2,
1204  //const Numeric& K3, Add again when NLTE is considered
1205  //const Numeric& K4, Add again when NLTE is considered
1206  const PropmatPartialsData& derivatives_data=PropmatPartialsData(),
1207  const Numeric& dQT_dT=0.0,
1208  const Numeric& dK1_dT=0.0,
1209  const Numeric& dK2_dT=0.0,
1210  const Numeric& dK2_dF0=0.0) const
1211  {
1212  const Index nppd = derivatives_data.nelem();
1213 
1214  const Numeric invQT = 1.0/QT;
1215  const Numeric QT_ratio = QT0 * invQT;
1216  const Numeric S = S0 * isotopic_ratio * QT_ratio * K1 * K2;
1217 
1218  for(Index jj = 0; jj < nppd; jj++)
1219  {
1220  if(derivatives_data(jj) == JQT_temperature)
1221  {
1222  Eigen::VectorXcd eig_dF = MapToEigen(dF[jj]);
1223 
1224  eig_dF *= S;
1225  eig_dF += MapToEigen(F) * (S0 * isotopic_ratio * QT_ratio *
1226  (K1 * dK2_dT +
1227  dK1_dT * K2 -
1228  invQT * dQT_dT * K1 * K2));
1229  }
1230  else if(derivatives_data(jj) == JQT_line_strength)
1231  {
1232  dF[jj] = F;
1233  dF[jj] *= isotopic_ratio;
1234  }
1235  else if(derivatives_data(jj) == JQT_line_center)
1236  {
1237  Eigen::VectorXcd eig_dF = MapToEigen(dF[jj]);
1238 
1239  eig_dF *= S;
1240  eig_dF += MapToEigen(F) * (S0 * isotopic_ratio * QT_ratio * K1 * dK2_dF0);
1241  }
1242  else
1243  {
1244  dF[jj] *= S;
1245  }
1246  }
1247  F *= S;
1248  }
1249 
1250  void apply_linestrength_from_full_linemixing(ComplexVector& F, // Returns the full complex line shape with line mixing
1251  ArrayOfComplexVector& dF,
1252  const Numeric& F0,
1253  const Numeric& T,
1254  const Complex& S_LM,
1255  const Numeric& isotopic_ratio,
1256  const PropmatPartialsData& derivatives_data=PropmatPartialsData(),
1257  const Complex& dS_LM_dT=0.0) const
1258  {
1259  extern const Numeric PLANCK_CONST;
1260  extern const Numeric BOLTZMAN_CONST;
1261  static const Numeric C1 = - PLANCK_CONST / BOLTZMAN_CONST;
1262 
1263  const Index nppd = derivatives_data.nelem();
1264 
1265  const Numeric invT = 1.0 / T,
1266  F0_invT = F0 * invT,
1267  exp_factor = exp(C1 * F0_invT),
1268  f0_factor = F0 * (1.0 - exp_factor);
1269 
1270  const Complex S = S_LM * f0_factor * isotopic_ratio;
1271 
1272  for(Index jj = 0; jj < nppd; jj++)
1273  {
1274  if(derivatives_data(jj) == JQT_temperature)
1275  {
1276  Eigen::VectorXcd eig_dF = MapToEigen(dF[jj]);
1277 
1278  eig_dF *= S_LM;
1279  eig_dF += MapToEigen(F) * (dS_LM_dT * f0_factor +
1280  S_LM * C1 * F0_invT * F0_invT * exp_factor) * isotopic_ratio;
1281  }
1282  else if(derivatives_data(jj) == JQT_line_strength)
1283  {
1284  throw std::runtime_error("Not working yet");
1285  }
1286  else if(derivatives_data(jj) == JQT_line_center)
1287  {
1288  throw std::runtime_error("Not working yet");
1289  }
1290  else
1291  {
1292  dF[jj] *= S;
1293  }
1294  }
1295 
1296  F *= S;
1297  }
1298 
1299  void apply_dipole(ComplexVector& F, // Returns the full complex line shape without line mixing
1300  ArrayOfComplexVector& dF, // Returns the derivatives of the full line shape for list_of_derivatives
1301  const Numeric& F0,
1302  const Numeric& T,
1303  const Numeric& d0,
1304  const Numeric& rho,
1305  const Numeric& isotopic_ratio,
1306  const PropmatPartialsData& derivatives_data=PropmatPartialsData(),
1307  const Numeric& drho_dT=0.0) const
1308  {
1309  // Output is d0^2 * rho * F * isotopic_ratio * F0 * (1-e^(hF0/kT))
1310 
1311  extern const Numeric PLANCK_CONST;
1312  extern const Numeric BOLTZMAN_CONST;
1313  static const Numeric C1 = - PLANCK_CONST / BOLTZMAN_CONST;
1314 
1315  const Index nppd = derivatives_data.nelem();
1316 
1317  const Numeric S = d0 * d0 * rho * isotopic_ratio,
1318  invT = 1.0 / T,
1319  F0_invT = F0 * invT,
1320  exp_factor = exp(C1 * F0_invT),
1321  f0_factor = F0 * (1.0 - exp_factor);
1322 
1323  for(Index jj = 0; jj < nppd; jj++)
1324  {
1325  if(derivatives_data(jj) == JQT_temperature)
1326  {
1327  ComplexMatrixViewMap eig_dF = MapToEigen(dF[jj]);
1328 
1329  eig_dF *= S * f0_factor;
1330  eig_dF += MapToEigen(F) * (d0 * d0 * (drho_dT * f0_factor +
1331  rho * C1 * F0_invT * F0_invT * exp_factor) * isotopic_ratio);
1332  }
1333  else if(derivatives_data(jj) == JQT_line_center)
1334  {
1335  throw std::runtime_error("Not working yet");
1336  }
1337  else
1338  {
1339  dF[jj] *= S * f0_factor;
1340  }
1341  }
1342 
1343  F *= S * f0_factor;
1344  }
1345 
1346  void apply_pressurebroadening_jacobian(ArrayOfComplexVector& dF,
1347  const PropmatPartialsData& derivatives_data,
1348  const Numeric& dgamma_dline_gamma_self=0.0,
1349  const Numeric& dgamma_dline_gamma_foreign=0.0,
1350  const Numeric& dgamma_dline_gamma_water=0.0,
1351  const Numeric& dgamma_dline_pressureshift_self=0.0,
1352  const Numeric& dgamma_dline_pressureshift_foreign=0.0,
1353  const Numeric& dgamma_dline_pressureshift_water=0.0,
1354  const Complex& dgamma_dline_gamma_selfexponent=0.0,
1355  const Complex& dgamma_dline_gamma_foreignexponent=0.0,
1356  const Complex& dgamma_dline_gamma_waterexponent=0.0)
1357  {
1358  const Index nppd = derivatives_data.nelem();
1359 
1360  for(Index jj = 0; jj < nppd; jj++)
1361  {
1362  if(derivatives_data(jj) == JQT_line_gamma_self)
1363  {
1364  dF[jj] *= dgamma_dline_gamma_self;
1365  }
1366  else if(derivatives_data(jj) == JQT_line_gamma_foreign)
1367  {
1368  dF[jj] *= dgamma_dline_gamma_foreign;
1369  }
1370  else if(derivatives_data(jj) == JQT_line_gamma_water)
1371  {
1372  dF[jj] *= dgamma_dline_gamma_water;
1373  }
1374  else if(derivatives_data(jj) == JQT_line_pressureshift_self)
1375  {
1376  dF[jj] *= Complex(0.0, dgamma_dline_pressureshift_self);
1377  }
1378  else if(derivatives_data(jj) == JQT_line_pressureshift_foreign)
1379  {
1380  dF[jj] *= Complex(0.0, dgamma_dline_pressureshift_foreign);
1381  }
1382  else if(derivatives_data(jj) == JQT_line_pressureshift_water)
1383  {
1384  dF[jj] *= Complex(0.0, dgamma_dline_pressureshift_water);
1385  }
1386  else if(derivatives_data(jj) == JQT_line_gamma_selfexponent)
1387  {
1388  dF[jj] *= dgamma_dline_gamma_selfexponent;
1389  }
1390  else if(derivatives_data(jj) == JQT_line_gamma_foreignexponent)
1391  {
1392  dF[jj] *= dgamma_dline_gamma_foreignexponent;
1393  }
1394  else if(derivatives_data(jj) == JQT_line_gamma_waterexponent)
1395  {
1396  dF[jj] *= dgamma_dline_gamma_waterexponent;
1397  }
1398  }
1399  }
1400 
1401 private:
1402  LineShapeType mtype;
1403  LineShapeNorm mnorm;
1404  Numeric cutoff;
1405 };
1406 
1407 typedef Array<LineFunctions> ArrayOfLineFunctions;
1408 
1409 #endif //lineshapedata_h
1410 
1411 
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
PLANCK_CONST
const Numeric PLANCK_CONST
Global constant, the Planck constant [Js].
complex.h
A class implementing complex numbers for ARTS.
LineFunctions::set_hui_etal_1978
void set_hui_etal_1978(ComplexVector &F, ArrayOfComplexVector &dF, const Vector &f_grid, const Numeric &zeeman_df, const Numeric &magnetic_magnitude, const Numeric &F0_noshift, const Numeric &GD_div_F0, const Numeric &G0, const Numeric &L0, const PropmatPartialsData &derivatives_data=PropmatPartialsData(), const Numeric &dGD_div_F0_dT=0.0, const Numeric &dG0_dT=0.0, const Numeric &dL0_dT=0.0) const
Definition: lineshapesdata.h:797
LineFunctions::apply_linestrength
void apply_linestrength(ComplexVector &F, ArrayOfComplexVector &dF, const Numeric &S0, const Numeric &isotopic_ratio, const Numeric &QT, const Numeric &QT0, const Numeric &K1, const Numeric &K2, const PropmatPartialsData &derivatives_data=PropmatPartialsData(), const Numeric &dQT_dT=0.0, const Numeric &dK1_dT=0.0, const Numeric &dK2_dT=0.0, const Numeric &dK2_dF0=0.0) const
Definition: lineshapesdata.h:1196
BOLTZMAN_CONST
const Numeric BOLTZMAN_CONST
Global constant, the Boltzmann constant [J/K].
LineFunctions::apply_VVH
void apply_VVH(ComplexVector &F, ArrayOfComplexVector &dF, const Vector &f_grid, const Numeric &F0, const Numeric &T, const PropmatPartialsData &derivatives_data=PropmatPartialsData()) const
Definition: lineshapesdata.h:1106
Vector
The Vector class.
Definition: matpackI.h:860
abs
#define abs(x)
Definition: legacy_continua.cc:20626
fac
Numeric fac(const Index n)
fac
Definition: math_funcs.cc:63
w
cmplx FADDEEVA() w(cmplx z, double relerr)
Definition: Faddeeva.cc:680
ArrayOfLineFunctions
Array< LineFunctions > ArrayOfLineFunctions
Definition: lineshapesdata.h:1407
LineFunctions::apply_pressurebroadening_jacobian
void apply_pressurebroadening_jacobian(ArrayOfComplexVector &dF, const PropmatPartialsData &derivatives_data, const Numeric &dgamma_dline_gamma_self=0.0, const Numeric &dgamma_dline_gamma_foreign=0.0, const Numeric &dgamma_dline_gamma_water=0.0, const Numeric &dgamma_dline_pressureshift_self=0.0, const Numeric &dgamma_dline_pressureshift_foreign=0.0, const Numeric &dgamma_dline_pressureshift_water=0.0, const Complex &dgamma_dline_gamma_selfexponent=0.0, const Complex &dgamma_dline_gamma_foreignexponent=0.0, const Complex &dgamma_dline_gamma_waterexponent=0.0)
Definition: lineshapesdata.h:1346
LineFunctions::apply_linestrength_from_full_linemixing
void apply_linestrength_from_full_linemixing(ComplexVector &F, ArrayOfComplexVector &dF, const Numeric &F0, const Numeric &T, const Complex &S_LM, const Numeric &isotopic_ratio, const PropmatPartialsData &derivatives_data=PropmatPartialsData(), const Complex &dS_LM_dT=0.0) const
Definition: lineshapesdata.h:1250
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
LineFunctions::apply_linemixing
void apply_linemixing(ComplexVector &F, ArrayOfComplexVector &dF, const Numeric &Y, const Numeric &G, const PropmatPartialsData &derivatives_data=PropmatPartialsData(), const Numeric &dY_dT=0.0, const Numeric &dG_dT=0.0) const
Definition: lineshapesdata.h:1001
d0
#define d0
Definition: legacy_continua.cc:21715
Complex
std::complex< Numeric > Complex
Definition: complex.h:33
ConstComplexVectorView::nelem
Index nelem() const
Definition: complex.h:280
oem::LM
invlib::LevenbergMarquardt< Numeric, CovarianceMatrix, Std > LM
Levenberg-Marquardt (LM) optimization using normed ARTS QR solver.
Definition: oem.h:173
LineFunctions::mnorm
LineShapeNorm mnorm
Definition: lineshapesdata.h:1403
LineFunctions::set_doppler
void set_doppler(ComplexVector &F, ArrayOfComplexVector &dF, const Vector &f_grid, const Numeric &zeeman_df, const Numeric &magnetic_magnitude, const Numeric &F0_noshift, const Numeric &GD_div_F0, const PropmatPartialsData &derivatives_data=PropmatPartialsData(), const Numeric &dGD_div_F0_dT=0.0) const
Definition: lineshapesdata.h:233
ComplexMatrixViewMap
Eigen::Map< ComplexMatrixType, 0, ComplexStrideType > ComplexMatrixViewMap
Definition: complex.h:170
LineFunctions::apply_dipole
void apply_dipole(ComplexVector &F, ArrayOfComplexVector &dF, const Numeric &F0, const Numeric &T, const Numeric &d0, const Numeric &rho, const Numeric &isotopic_ratio, const PropmatPartialsData &derivatives_data=PropmatPartialsData(), const Numeric &drho_dT=0.0) const
Definition: lineshapesdata.h:1299
LineFunctions::mtype
LineShapeType mtype
Definition: lineshapesdata.h:1402
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
absorption.h
Declarations required for the calculation of absorption coefficients.
dx
#define dx
Definition: legacy_continua.cc:22180
LineFunctions::set_lorentz
void set_lorentz(ComplexVector &F, ArrayOfComplexVector &dF, const Vector &f_grid, const Numeric &zeeman_df, const Numeric &magnetic_magnitude, const Numeric &F0_noshift, const Numeric &G0, const Numeric &L0, const PropmatPartialsData &derivatives_data=PropmatPartialsData(), const Numeric &dG0_dT=0.0, const Numeric &dL0_dT=0.0) const
Definition: lineshapesdata.h:83
ComplexVector
The ComplexVector class.
Definition: complex.h:573
partial_derivatives.h
Computes partial derivatives (old method)
LineFunctions::set_faddeeva_from_full_linemixing
void set_faddeeva_from_full_linemixing(ComplexVector &F, ArrayOfComplexVector &dF, const Vector &f_grid, const Complex &eigenvalue_no_shift, const Numeric &GD_div_F0, const Numeric &L0, const PropmatPartialsData &derivatives_data=PropmatPartialsData(), const Numeric &dGD_div_F0_dT=0.0, const Complex &deigenvalue_dT=0.0, const Numeric &dL0_dT=0.0) const
Definition: lineshapesdata.h:674
Array
This can be used to make arrays out of anything.
Definition: array.h:40
LineFunctions::set_o2_non_resonant
void set_o2_non_resonant(ComplexVector &F, ArrayOfComplexVector &dF, const Vector &f_grid, const Numeric &F0, const Numeric &G0, const PropmatPartialsData &derivatives_data=PropmatPartialsData(), const Numeric &dG0_dT=0.0)
Definition: lineshapesdata.h:942
PI
const Numeric PI
Global constant, pi.
LineFunctions::apply_rosenkranz_quadratic
void apply_rosenkranz_quadratic(ComplexVector &F, ArrayOfComplexVector &dF, const Vector &f_grid, const Numeric &F0, const Numeric &T, const PropmatPartialsData &derivatives_data=PropmatPartialsData()) const
Definition: lineshapesdata.h:1048
F0
G0 G2 FVC Y DV F0
Definition: arts_api_classes.cc:218
linerecord.h
LineRecord class for managing line catalog data.
LineFunctions::set_faddeeva_algorithm916
void set_faddeeva_algorithm916(ComplexVector &F, ArrayOfComplexVector &dF, const Vector &f_grid, const Numeric &zeeman_df, const Numeric &magnetic_magnitude, const Numeric &F0_noshift, const Numeric &GD_div_F0, const Numeric &G0, const Numeric &L0, const PropmatPartialsData &derivatives_data=PropmatPartialsData(), const Numeric &dGD_div_F0_dT=0.0, const Numeric &dG0_dT=0.0, const Numeric &dL0_dT=0.0) const
Definition: lineshapesdata.h:543
LineFunctions::set_htp
void set_htp(ComplexVector &F, ArrayOfComplexVector &dF, const Vector &f_grid, const Numeric &zeeman_df, const Numeric &magnetic_magnitude, const Numeric &F0_noshift, const Numeric &GD_div_F0, const Numeric &G0, const Numeric &L0, const Numeric &L2, const Numeric &G2, const Numeric &eta, const Numeric &FVC, const PropmatPartialsData &derivatives_data=PropmatPartialsData(), const Numeric &dGD_div_F0_dT=0.0, const Numeric &dG0_dT=0.0, const Numeric &dL0_dT=0.0, const Numeric &dG2_dT=0.0, const Numeric &dL2_dT=0.0, const Numeric &deta_dT=0.0, const Numeric &dFVC_dT=0.0) const
Definition: lineshapesdata.h:254
MapToEigen
ComplexConstMatrixViewMap MapToEigen(const ConstComplexMatrixView &A)
Definition: complex.cc:1655
LineFunctions::apply_VVW
void apply_VVW(ComplexVector &F, ArrayOfComplexVector &dF, const Vector &f_grid, const Numeric &F0, const PropmatPartialsData &derivatives_data=PropmatPartialsData()) const
Definition: lineshapesdata.h:1160
LineFunctions::set_mirrored_lorentz
void set_mirrored_lorentz(ComplexVector &F, ArrayOfComplexVector &dF, const Vector &f_grid, const Numeric &zeeman_df, const Numeric &magnetic_magnitude, const Numeric &F0_noshift, const Numeric &G0, const Numeric &L0, const PropmatPartialsData &derivatives_data=PropmatPartialsData(), const Numeric &dG0_dT=0.0, const Numeric &dL0_dT=0.0) const
Definition: lineshapesdata.h:163
sqrt
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:620