ARTS  2.3.1285(git:92a29ea9-dirty)
lineshapes.cc
Go to the documentation of this file.
1 /* Copyright (C) 2000-2012
2  Axel von Engeln <engeln@uni-bremen.de>
3  Stefan Buehler <sbuehler@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 
34 #include <cmath>
35 #include <Faddeeva/Faddeeva.hh>
36 #include "absorption.h"
37 #include "array.h"
38 #include "arts.h"
39 #include "complex.h"
40 #include "matpackI.h"
41 
61 void lineshape_no_shape(Vector&,
62  Vector&,
63  Vector&,
64  Vector&,
65  Vector&,
66  Vector&,
67  Vector&,
68  const Numeric,
69  const Numeric,
70  const Numeric,
71  const Numeric,
72  const Numeric,
73  const Numeric,
74  const Numeric,
75  const Numeric,
76  ConstVectorView,
77  const bool,
78  const bool) {
79  // This function should never be called so throw an error here:
80  throw std::runtime_error(
81  "The no_shape lineshape is only a placeholder, but you tried\n"
82  "to use it like a real lineshape.");
83 }
84 
97 void lineshape_lorentz(Vector& ls_attenuation,
98  Vector& ls_phase,
99  Vector&, //ls_dattenuation_dfrequency_term
100  Vector&, //ls_dphase_dfrequency_term
101  Vector&, //ls_dattenuation_dpressure_term
102  Vector&, //ls_dphase_dpressure_term
103  Vector&, //X
104  const Numeric f0,
105  const Numeric gamma,
106  const Numeric,
107  const Numeric,
108  const Numeric,
109  const Numeric,
110  const Numeric, //sigma
111  const Numeric,
112  ConstVectorView f_grid,
113  const bool do_phase,
114  const bool) {
115  const Index nf = f_grid.nelem();
116 
117  // PI:
118  extern const Numeric PI;
119 
120  // assert( ls.nelem() == nf );
121 
122  Numeric gamma2 = gamma * gamma;
123  Numeric fac = gamma / PI;
124 
125  if (do_phase) {
126  for (Index i = 0; i < nf; ++i) {
127  ls_attenuation[i] = fac / ((f_grid[i] - f0) * (f_grid[i] - f0) + gamma2);
128  ls_phase[i] = (f_grid[i] - f0) / PI /
129  ((f_grid[i] - f0) * (f_grid[i] - f0) + gamma2);
130  }
131  } else {
132  for (Index i = 0; i < nf; ++i) {
133  ls_attenuation[i] = fac / ((f_grid[i] - f0) * (f_grid[i] - f0) + gamma2);
134  }
135  }
136 }
137 
148 void lineshape_mirrored_lorentz(Vector& ls_attenuation,
149  Vector& ls_phase,
150  Vector&, //ls_dattenuation_dfrequency_term
151  Vector&, //ls_dphase_dfrequency_term
152  Vector&, //ls_dattenuation_dpressure_term
153  Vector&, //ls_dphase_dpressure_term
154  Vector&,
155  const Numeric f0,
156  const Numeric gamma,
157  const Numeric,
158  const Numeric,
159  const Numeric,
160  const Numeric,
161  const Numeric,
162  const Numeric,
163  ConstVectorView f_grid,
164  const bool do_phase,
165  const bool) {
166  const Index nf = f_grid.nelem();
167 
168  // PI:
169  extern const Numeric PI;
170 
171  // assert( ls.nelem() == nf );
172 
173  Numeric gamma2 = gamma * gamma;
174  Numeric fac = gamma / PI;
175 
176  if (do_phase) {
177  for (Index i = 0; i < nf; ++i) {
178  ls_attenuation[i] = fac / ((f_grid[i] - f0) * (f_grid[i] - f0) + gamma2) +
179  fac / ((f_grid[i] + f0) * (f_grid[i] + f0) + gamma2);
180  ls_phase[i] = (f_grid[i] - f0) / PI /
181  ((f_grid[i] - f0) * (f_grid[i] - f0) + gamma2) +
182  (f_grid[i] + f0) / PI /
183  ((f_grid[i] + f0) * (f_grid[i] + f0) +
184  gamma2); //uncertain on this part
185  }
186  } else {
187  for (Index i = 0; i < nf; ++i) {
188  ls_attenuation[i] = fac / ((f_grid[i] - f0) * (f_grid[i] - f0) + gamma2) +
189  fac / ((f_grid[i] + f0) * (f_grid[i] + f0) + gamma2);
190  }
191  }
192 }
193 
206 void lineshape_doppler(Vector& ls_attenuation,
207  Vector&, //ls_phase
208  Vector&, //ls_dattenuation_dfrequency_term
209  Vector&, //ls_dphase_dfrequency_term
210  Vector&, //ls_dattenuation_dpressure_term
211  Vector&, //ls_dphase_dpressure_term
212  Vector&, //x
213  const Numeric f0,
214  const Numeric, //gamma
215  const Numeric,
216  const Numeric,
217  const Numeric,
218  const Numeric,
219  const Numeric sigma,
220  const Numeric,
221  ConstVectorView f_grid,
222  const bool,
223  const bool) {
224  const Index nf = f_grid.nelem();
225 
226  // SQRT(PI):
227  extern const Numeric PI;
228  const Numeric sqrtPI = sqrt(PI);
229 
230  // assert( ls_attenuation.nelem() == nf );
231 
232  Numeric sigma2 = sigma * sigma;
233  Numeric fac = 1.0 / (sqrtPI * sigma);
234 
235  for (Index i = 0; i < nf; ++i) {
236  ls_attenuation[i] = fac * exp(-pow(f_grid[i] - f0, 2) / sigma2);
237  }
238 }
239 
252 void deprecated(Vector&,
253  Vector&,
254  Vector&,
255  Vector&,
256  Vector&,
257  Vector&,
258  Vector&,
259  const Numeric,
260  const Numeric,
261  const Numeric,
262  const Numeric,
263  const Numeric,
264  const Numeric,
265  const Numeric,
266  const Numeric,
267  ConstVectorView,
268  const bool,
269  const bool) {
270  throw std::runtime_error(
271  "You are running a deprecated lineshape.\n"
272  "The deprecated line shapes are:\n"
273  "\t Voigt_Kuntz3 --- please replace with Voigt_Kuntz6\n"
274  "\t Voigt_Kuntz4 --- please replace with Voigt_Kuntz6\n");
275 }
276 //------------------------------------------------------------------------
277 
278 // help function for lineshape_voigt_kuntz1
279 long bfun6_(Numeric y, Numeric x) {
280  /* System generated locals */
281  long int ret_val;
282 
283  /* Local variables */
284  Numeric s = 0;
285 
286  /* -------------------------------------------------------------------- */
287  s = fabs(x) + y;
288  if (s >= 15.f) {
289  ret_val = 1;
290  } else if (s >= 5.5f) {
291  ret_val = 2;
292  } else if (y >= fabs(x) * .195f - .176f) {
293  ret_val = 3;
294  } else {
295  ret_val = 4;
296  }
297  return ret_val;
298 } /* bfun6_ */
299 
359 void lineshape_voigt_kuntz6(Vector& ls_attenuation,
360  Vector&, //ls_phase
361  Vector& ls_dattenuation_dfrequency_term,
362  Vector&, //ls_dphase_dfrequency_term
363  Vector& ls_dattenuation_dpressure_term,
364  Vector&, //ls_dphase_dpressure_term
365  Vector& x,
366  const Numeric f0,
367  const Numeric gamma,
368  const Numeric,
369  const Numeric,
370  const Numeric,
371  const Numeric,
372  const Numeric sigma,
373  const Numeric,
374  ConstVectorView f_grid,
375  const bool,
376  const bool do_partials)
377 
378 {
379  const Index nf = f_grid.nelem();
380 
381  // seems not necessary for Doppler correction
382  // extern const Numeric SQRT_NAT_LOG_2;
383 
384  // PI
385  extern const Numeric PI;
386 
387  // constant sqrt(1/pi)
388  const Numeric sqrt_invPI = sqrt(1 / PI);
389 
390  // constant normalization factor for voigt
391  Numeric fac = 1.0 / sigma * sqrt_invPI;
392 
393  /* Initialized data */
394 
395  Numeric yps1 = -1.0;
396  Numeric yps2 = -1.0;
397  Numeric yps3 = -1.0;
398  Numeric yps4 = -1.0;
399 
400  /* System generated locals */
401  long int i__1, i__2;
402  Numeric r__1;
403 
404  /* Local variables */
405  long int bmin = 0, lauf[16] = {0} /* was [4][4] */, bmax;
406  long int imin = 0, imax = 0, stack[80] = {0} /* was [20][4] */;
407  Numeric a1, a2, a3, a4, a5, a6, a8, b8, c8, d8, e8, f8, g8, h8, a7, b7, c7,
408  d7, e7, f7, o8, p8, q8, r8, s8, t8, g7, h7, o7, p7, q7, r7, s7, t7, b6,
409  c6, d6, e6, b5, c5, d5, e5, b4, c4, d4, b3, c3, d3, b1, y2;
410  a1 = a2 = a3 = a4 = a5 = a6 = a8 = b8 = c8 = d8 = e8 = f8 = g8 = h8 = a7 =
411  b7 = c7 = d7 = e7 = f7 = o8 = p8 = q8 = r8 = s8 = t8 = g7 = h7 = o7 = p7 =
412  q7 = r7 = s7 = t7 = b6 = c6 = d6 = e6 = b5 = c5 = d5 = e5 = b4 = c4 =
413  d4 = b3 = c3 = d3 = b1 = y2 = 0;
414  long int i2 = 0, i1 = 0;
415  Numeric x2 = 0, b2 = 0, c1 = 0;
416  long int stackp = 0, imitte = 0;
417  Numeric ym2 = 0;
418 
419  // variables needed in original c routine:
420 
421  // Ratio of the Lorentz halfwidth to the Doppler halfwidth
422  Numeric y = gamma / sigma;
423 
424  // Helper bools
425  bool do_x = false, do_y = false;
426 
427  // Forcing 0.0001 here, or 1/10000th of relevant broadening.
428  // This is from Wfuns agreeing on less than 1e-6 with this
429  // level of perturbation (faddeeva_algorithm_916 and Voigt_Kuntz6 were compared)
430  const Numeric dx = 0.0001, dy = 0.0001;
431 
432  // do_partials is put here because later on x will change...
433  if (do_partials) {
434  // pass variable for self-calling
435  Vector empty;
436 
437  // x = (f-f0)/sigma, so x+dx from f0+df0 means dx = - df0/sigma
438  const Numeric df0 = -dx * sigma;
439 
440  // Use this dx to see if reasonable to calculate dFa/dx.
441  // Too small dx and we will have constant
442  // Jacobian that are unreasonably large...
443  do_x = (f0 > (-df0 * 10));
444 
445  if (do_x) {
446  lineshape_voigt_kuntz6(ls_dattenuation_dfrequency_term,
447  empty,
448  empty,
449  empty,
450  empty,
451  empty,
452  x,
453  f0 + df0,
454  gamma,
455  0.0,
456  0.0,
457  0.0,
458  0.0,
459  sigma,
460  0.0,
461  f_grid,
462  false,
463  false);
464  }
465 
466  // y = gamma/sigma so y+dy from gamma+dgamma means dy=dgamma/sigma
467  const Numeric dgamma = dy * sigma;
468 
469  // Use this dy to see if reasonable to calculate dFa/dy.
470  // Too small dy and we will have constant
471  // Jacobian that are unreasonably large
472  do_y = (y > (dy * 10));
473 
474  if (do_y) {
475  lineshape_voigt_kuntz6(ls_dattenuation_dpressure_term,
476  empty,
477  empty,
478  empty,
479  empty,
480  empty,
481  x,
482  f0,
483  gamma + dgamma,
484  0.0,
485  0.0,
486  0.0,
487  0.0,
488  sigma,
489  0.0,
490  f_grid,
491  false,
492  false);
493  }
494  }
495 
496  // frequency in units of Doppler
497  for (i1 = 0; i1 < (int)nf; i1++) {
498  x[i1] = (f_grid[i1] - f0) / sigma;
499  }
500 
501  /* Parameter adjustments */
502  // this does not work for variables type Vector, corrected by
503  // adjusting the index to array ls and x
504  //--ls;
505  //--x;
506 
507  /* Function Body */
508  y2 = y * y;
509  if (y >= 15.0 || x[0] >= 15.0 || x[nf - 1] <= -15.0) {
510  lauf[0] = 1;
511  lauf[4] = nf;
512  lauf[8] = nf;
513  lauf[12] = 0;
514  goto L7;
515  }
516  for (i2 = 1; i2 <= 4; ++i2) {
517  for (i1 = 1; i1 <= 4; ++i1) {
518  lauf[i1 + (i2 << 2) - 5] = i2 % 2 * (nf + 1);
519  /* L1: */
520  }
521  }
522  stackp = 1;
523  stack[stackp - 1] = 1;
524  stack[stackp + 19] = nf;
525  stack[stackp + 39] = bfun6_(y, x[0]);
526  stack[stackp + 59] = bfun6_(y, x[nf - 1]);
527 L2:
528  imin = stack[stackp - 1];
529  imax = stack[stackp + 19];
530  bmin = stack[stackp + 39];
531  bmax = stack[stackp + 59];
532  if (bmin == bmax) {
533  if (x[imax - 1] < 0.f) {
534  /* Computing MIN */
535  i__1 = imin, i__2 = lauf[bmin - 1];
536  lauf[bmin - 1] = min(i__1, i__2);
537  /* Computing MAX */
538  i__1 = imax, i__2 = lauf[bmax + 3];
539  lauf[bmax + 3] = max(i__1, i__2);
540  --stackp;
541  goto L3;
542  } else if (x[imin - 1] >= 0.f) {
543  /* Computing MIN */
544  i__1 = imin, i__2 = lauf[bmin + 7];
545  lauf[bmin + 7] = min(i__1, i__2);
546  /* Computing MAX */
547  i__1 = imax, i__2 = lauf[bmax + 11];
548  lauf[bmax + 11] = max(i__1, i__2);
549  --stackp;
550  goto L3;
551  }
552  }
553  imitte = (imax + imin) / 2;
554  stack[stackp - 1] = imitte + 1;
555  stack[stackp + 19] = imax;
556  stack[stackp + 39] = bfun6_(y, x[imitte]);
557  stack[stackp + 59] = bmax;
558  ++stackp;
559  stack[stackp - 1] = imin;
560  stack[stackp + 19] = imitte;
561  stack[stackp + 39] = bmin;
562  stack[stackp + 59] = bfun6_(y, x[imitte - 1]);
563 L3:
564  if (stackp > 0) {
565  goto L2;
566  }
567  /* ---- Region 4 */
568  /* -------------------------------------------------------------------- */
569  if (lauf[7] >= lauf[3] || lauf[15] >= lauf[11]) {
570  if ((r__1 = y - yps4, fabs(r__1)) > 1e-8f) {
571  //yps4 = y;
572  a7 =
573  y *
574  (y2 *
575  (y2 *
576  (y2 *
577  (y2 *
578  (y2 *
579  (y2 *
580  (y2 *
581  (y2 *
582  (y2 *
583  (y2 *
584  (y2 *
585  (y2 *
586  (2.35944f -
587  y2 *
588  .56419f) -
589  72.9359f) +
590  571.687f) -
591  5860.68f) +
592  40649.2f) -
593  320772.f) +
594  1684100.f) -
595  9694630.f) +
596  40816800.f) -
597  1.53575e8f) +
598  4.56662e8f) -
599  9.86604e8f) +
600  1.16028e9f);
601  b7 =
602  y *
603  (y2 *
604  (y2 *
605  (y2 *
606  (y2 *
607  (y2 *
608  (y2 *
609  (y2 *
610  (y2 *
611  (y2 *
612  (y2 *
613  (y2 *
614  (23.0312f -
615  y2 *
616  7.33447f) -
617  234.143f) -
618  2269.19f) +
619  45251.3f) -
620  234417.f) +
621  3599150.f) -
622  7723590.f) +
623  86482900.f) -
624  2.91876e8f) +
625  8.06985e8f) -
626  9.85386e8f) -
627  5.60505e8f);
628  c7 =
629  y *
630  (y2 *
631  (y2 *
632  (y2 *
633  (y2 *
634  (y2 *
635  (y2 *
636  (y2 *
637  (y2 * (y2 * (y2 * (97.6203f -
638  y2 * 44.0068f) +
639  1097.77f) -
640  25338.3f) +
641  98079.1f) +
642  576054.f) -
643  2.3818e7f) +
644  22930200.f) -
645  2.04467e8f) +
646  2.94262e8f) +
647  2.47157e8f) -
648  6.51523e8f);
649  d7 =
650  y *
651  (y2 *
652  (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (228.563f -
653  y2 * 161.358f) +
654  8381.97f) -
655  66431.2f) -
656  303569.f) +
657  2240400.f) +
658  38311200.f) -
659  41501300.f) -
660  99622400.f) +
661  2.70167e8f) -
662  2.63894e8f);
663  e7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (296.38f -
664  y2 * 403.396f) +
665  23507.6f) -
666  66212.1f) -
667  1.003e6f) +
668  468142.f) +
669  24620100.f) +
670  5569650.f) +
671  1.40677e8f) -
672  63177100.f);
673  f7 =
674  y *
675  (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (125.591f - y2 * 726.113f) +
676  37544.8f) +
677  8820.94f) -
678  934717.f) -
679  1931140.f) -
680  33289600.f) +
681  4073820.f) -
682  16984600.f);
683  g7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (-260.198f - y2 * 968.15f) +
684  37371.9f) +
685  79902.5f) -
686  186682.f) -
687  900010.f) +
688  7528830.f) -
689  1231650.f);
690  h7 =
691  y *
692  (y2 * (y2 * (y2 * (y2 * (y2 * (-571.645f - y2 * 968.15f) + 23137.1f) +
693  72520.9f) +
694  153468.f) +
695  86407.6f) -
696  610622.f);
697  o7 = y * (y2 * (y2 * (y2 * (y2 * (-575.164f - y2 * 726.113f) + 8073.15f) +
698  26538.5f) +
699  49883.8f) -
700  23586.5f);
701  p7 = y * (y2 * (y2 * (y2 * (-352.467f - y2 * 403.396f) + 953.655f) +
702  2198.86f) -
703  8009.1f);
704  q7 = y * (y2 * (y2 * (-134.792f - y2 * 161.358f) - 271.202f) - 622.056f);
705  r7 = y * (y2 * (-29.7896f - y2 * 44.0068f) - 77.0535f);
706  s7 = y * (-2.92264f - y2 * 7.33447f);
707  t7 = y * -.56419f;
708  a8 =
709  y2 *
710  (y2 *
711  (y2 *
712  (y2 *
713  (y2 *
714  (y2 *
715  (y2 *
716  (y2 *
717  (y2 *
718  (y2 *
719  (y2 *
720  (y2 *
721  (y2 *
722  (y2 -
723  3.68288f) +
724  126.532f) -
725  955.194f) +
726  9504.65f) -
727  70946.1f) +
728  483737.f) -
729  2857210.f) +
730  14464700.f) -
731  61114800.f) +
732  2.11107e8f) -
733  5.79099e8f) +
734  1.17022e9f) -
735  1.5599e9f) +
736  1.02827e9f;
737  b8 =
738  y2 *
739  (y2 *
740  (y2 *
741  (y2 *
742  (y2 *
743  (y2 *
744  (y2 *
745  (y2 *
746  (y2 *
747  (y2 *
748  (y2 *
749  (y2 *
750  (y2 *
751  14.f -
752  40.5117f) +
753  533.254f) +
754  3058.26f) -
755  55600.f) +
756  498334.f) -
757  2849540.f) +
758  13946500.f) -
759  70135800.f) +
760  2.89676e8f) -
761  7.53828e8f) +
762  1.66421e9f) -
763  2.28855e9f) +
764  1.5599e9f;
765  c8 =
766  y2 *
767  (y2 *
768  (y2 *
769  (y2 *
770  (y2 *
771  (y2 *
772  (y2 *
773  (y2 * (y2 * (y2 * (y2 * (y2 * 91 -
774  198.876f) -
775  1500.17f) +
776  48153.3f) -
777  217801.f) -
778  1063520.f) +
779  1.4841e7f) -
780  46039600.f) +
781  63349600.f) -
782  6.60078e8f) +
783  1.06002e9f) -
784  1.66421e9f) +
785  1.17022e9f;
786  d8 =
787  y2 *
788  (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 364 -
789  567.164f) -
790  16493.7f) +
791  161461.f) +
792  280428.f) -
793  6890020.f) -
794  6876560.f) +
795  1.99846e8f) +
796  54036700.f) +
797  6.60078e8f) -
798  7.53828e8f) +
799  5.79099e8f;
800  e8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 1001 -
801  1012.79f) -
802  55582.f) +
803  240373.f) +
804  1954700.f) -
805  5257220.f) -
806  50101700.f) -
807  1.99846e8f) +
808  63349600.f) -
809  2.89676e8f) +
810  2.11107e8f;
811  f8 =
812  y2 *
813  (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 2002 - 1093.82f) -
814  106663.f) +
815  123052.f) +
816  3043160.f) +
817  5257220.f) -
818  6876560.f) +
819  46039600.f) -
820  70135800.f) +
821  61114800.f;
822  g8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 3003 - 486.14f) -
823  131337.f) -
824  123052.f) +
825  1954700.f) +
826  6890020.f) +
827  1.4841e7f) -
828  13946500.f) +
829  14464700.f;
830  h8 =
831  y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 3432 + 486.14f) - 106663.f) -
832  240373.f) +
833  280428.f) +
834  1063520.f) -
835  2849540.f) +
836  2857210.f;
837  o8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 3003 + 1093.82f) - 55582.f) -
838  161461.f) -
839  217801.f) -
840  498334.f) +
841  483737.f;
842  p8 =
843  y2 *
844  (y2 * (y2 * (y2 * (y2 * 2002 + 1012.79f) - 16493.7f) - 48153.3f) -
845  55600.f) +
846  70946.1f;
847  q8 = y2 * (y2 * (y2 * (y2 * 1001.f + 567.164f) - 1500.17f) - 3058.26f) +
848  9504.65f;
849  r8 = y2 * (y2 * (y2 * 364 + 198.876f) + 533.254f) + 955.194f;
850  s8 = y2 * (y2 * 91.f + 40.5117f) + 126.532f;
851  t8 = y2 * 14.f + 3.68288f;
852  }
853  ym2 = y * 2;
854  for (i2 = 1; i2 <= 3; i2 += 2) {
855  i__1 = lauf[((i2 + 1) << 2) - 1];
856  for (i1 = lauf[(i2 << 2) - 1]; i1 <= i__1; ++i1) {
857  x2 = x[i1 - 1] * x[i1 - 1];
858  ls_attenuation[i1 - 1] =
859  fac *
860  (exp(y2 - x2) * cos(x[i1 - 1] * ym2) -
861  (a7 +
862  x2 *
863  (b7 +
864  x2 *
865  (c7 +
866  x2 *
867  (d7 +
868  x2 *
869  (e7 +
870  x2 *
871  (f7 +
872  x2 *
873  (g7 +
874  x2 *
875  (h7 +
876  x2 *
877  (o7 +
878  x2 *
879  (p7 +
880  x2 *
881  (q7 +
882  x2 *
883  (r7 +
884  x2 *
885  (s7 +
886  x2 *
887  t7))))))))))))) /
888  (a8 +
889  x2 *
890  (b8 +
891  x2 *
892  (c8 +
893  x2 *
894  (d8 +
895  x2 *
896  (e8 +
897  x2 *
898  (f8 +
899  x2 *
900  (g8 +
901  x2 *
902  (h8 +
903  x2 *
904  (o8 +
905  x2 *
906  (p8 +
907  x2 *
908  (q8 +
909  x2 *
910  (r8 +
911  x2 *
912  (s8 +
913  x2 *
914  (t8 +
915  x2)))))))))))))));
916  /* L4: */
917  }
918  }
919  }
920  /* ---- Region 3 */
921  /* -------------------------------------------------------------------- */
922  if (lauf[6] >= lauf[2] || lauf[14] >= lauf[10]) {
923  if ((r__1 = y - yps3, fabs(r__1)) > 1e-8f) {
924  //yps3 = y;
925  a5 = y * (y * (y * (y * (y * (y * (y * (y * (y * .564224f + 7.55895f) +
926  49.5213f) +
927  204.501f) +
928  581.746f) +
929  1174.8f) +
930  1678.33f) +
931  1629.76f) +
932  973.778f) +
933  272.102f;
934  b5 = y * (y * (y * (y * (y * (y * (y * 2.25689f + 22.6778f) + 100.705f) +
935  247.198f) +
936  336.364f) +
937  220.843f) -
938  2.34403f) -
939  60.5644f;
940  c5 =
941  y * (y * (y * (y * (y * 3.38534f + 22.6798f) + 52.8454f) + 42.5683f) +
942  18.546f) +
943  4.58029f;
944  d5 = y * (y * (y * 2.25689f + 7.56186f) + 1.66203f) - .128922f;
945  e5 = y * .564224f + 9.71457e-4f;
946  a6 = y * (y * (y * (y * (y * (y * (y * (y * (y * (y + 13.3988f) +
947  88.2674f) +
948  369.199f) +
949  1074.41f) +
950  2256.98f) +
951  3447.63f) +
952  3764.97f) +
953  2802.87f) +
954  1280.83f) +
955  272.102f;
956  b6 = y * (y * (y * (y * (y * (y * (y * (y * 5.f + 53.5952f) + 266.299f) +
957  793.427f) +
958  1549.68f) +
959  2037.31f) +
960  1758.34f) +
961  902.306f) +
962  211.678f;
963  c6 = y * (y * (y * (y * (y * (y * 10.f + 80.3928f) + 269.292f) +
964  479.258f) +
965  497.302f) +
966  308.186f) +
967  78.866f;
968  d6 = y * (y * (y * (y * 10.f + 53.5952f) + 92.7568f) + 55.0293f) +
969  22.0353f;
970  e6 = y * (y * 5.f + 13.3988f) + 1.49645f;
971  }
972  for (i2 = 1; i2 <= 3; i2 += 2) {
973  i__1 = lauf[((i2 + 1) << 2) - 2];
974  for (i1 = lauf[(i2 << 2) - 2]; i1 <= i__1; ++i1) {
975  x2 = x[i1 - 1] * x[i1 - 1];
976  ls_attenuation[i1 - 1] =
977  fac * (a5 + x2 * (b5 + x2 * (c5 + x2 * (d5 + x2 * e5)))) /
978  (a6 + x2 * (b6 + x2 * (c6 + x2 * (d6 + x2 * (e6 + x2)))));
979  /* L5: */
980  }
981  }
982  }
983  /* ---- Region 2 */
984  /* -------------------------------------------------------------------- */
985  if (lauf[5] >= lauf[1] || lauf[13] >= lauf[9]) {
986  if ((r__1 = y - yps2, fabs(r__1)) > 1e-8f) {
987  //yps2 = y;
988  a3 = y * (y2 * (y2 * (y2 * .56419f + 3.10304f) + 4.65456f) + 1.05786f);
989  b3 = y * (y2 * (y2 * 1.69257f + .56419f) + 2.962f);
990  c3 = y * (y2 * 1.69257f - 2.53885f);
991  d3 = y * .56419f;
992  a4 = y2 * (y2 * (y2 * (y2 + 6.f) + 10.5f) + 4.5f) + .5625f;
993  b4 = y2 * (y2 * (y2 * 4.f + 6.f) + 9.f) - 4.5f;
994  c4 = y2 * (y2 * 6.f - 6.f) + 10.5f;
995  d4 = y2 * 4.f - 6.f;
996  }
997  for (i2 = 1; i2 <= 3; i2 += 2) {
998  i__1 = lauf[((i2 + 1) << 2) - 3];
999  for (i1 = lauf[(i2 << 2) - 3]; i1 <= i__1; ++i1) {
1000  x2 = x[i1 - 1] * x[i1 - 1];
1001  ls_attenuation[i1 - 1] = fac * (a3 + x2 * (b3 + x2 * (c3 + x2 * d3))) /
1002  (a4 + x2 * (b4 + x2 * (c4 + x2 * (d4 + x2))));
1003  /* L6: */
1004  }
1005  }
1006  }
1007  /* ---- Region 1 */
1008  /* -------------------------------------------------------------------- */
1009 L7:
1010  if (lauf[4] >= lauf[0] || lauf[12] >= lauf[8]) {
1011  if ((r__1 = y - yps1, fabs(r__1)) > 1e-8f) {
1012  //yps1 = y;
1013  a1 = y * .5641896f;
1014  b1 = y2 + .5f;
1015  a2 = y2 * 4;
1016  }
1017 
1018  c1 = fac * a1;
1019  for (i2 = 1; i2 <= 3; i2 += 2) {
1020  i__1 = lauf[((i2 + 1) << 2) - 4];
1021  for (i1 = lauf[(i2 << 2) - 4]; i1 <= i__1; ++i1) {
1022  x2 = x[i1 - 1] * x[i1 - 1];
1023  b2 = b1 - x2;
1024  ls_attenuation[i1 - 1] = c1 * (b1 + x2) / (b2 * b2 + a2 * x2);
1025  /* L8: */
1026  }
1027  }
1028  }
1029 
1030  if (do_partials) {
1031  // Set things right depending on if possible
1032  if (do_y && do_x) {
1033  for (Index iv = 0; iv < nf; iv++) {
1034  ls_dattenuation_dfrequency_term[iv] -= ls_attenuation[iv];
1035  ls_dattenuation_dfrequency_term[iv] /= dx;
1036  ls_dattenuation_dpressure_term[iv] -= ls_attenuation[iv];
1037  ls_dattenuation_dpressure_term[iv] /= dy;
1038  }
1039  } else if (do_x) {
1040  for (Index iv = 0; iv < nf; iv++) {
1041  ls_dattenuation_dfrequency_term[iv] -= ls_attenuation[iv];
1042  ls_dattenuation_dfrequency_term[iv] /= dx;
1043  ls_dattenuation_dpressure_term[iv] = 0.0;
1044  }
1045  } else if (do_y) {
1046  for (Index iv = 0; iv < nf; iv++) {
1047  ls_dattenuation_dfrequency_term[iv] = 0.0;
1048  ls_dattenuation_dpressure_term[iv] -= ls_attenuation[iv];
1049  ls_dattenuation_dpressure_term[iv] /= dy;
1050  }
1051  } else {
1052  for (Index iv = 0; iv < nf; iv++) {
1053  ls_dattenuation_dfrequency_term[iv] = 0.0;
1054  ls_dattenuation_dpressure_term[iv] = 0.0;
1055  }
1056  }
1057  }
1058 }
1059 
1060 //---------------------------------------------------------------
1061 // help function for lineshape_voigt_kuntz3
1062 
1063 long int bfun3_(Numeric y, Numeric x) {
1064  /* System generated locals */
1065  long int ret_val;
1066 
1067  /* Local variables */
1068  Numeric x2 = 0, y2 = 0;
1069 
1070  /* -------------------------------------------------------------------- */
1071  x2 = x * x;
1072  y2 = y * y;
1073  if (x2 * .4081676f + y2 > 21.159543f) {
1074  if (x2 * .7019639f + y2 > 1123.14221f) {
1075  ret_val = 0;
1076  } else {
1077  ret_val = 1;
1078  }
1079  } else {
1080  if (x2 * .20753051f + y2 > 4.20249292f) {
1081  ret_val = 2;
1082  } else if (y >= fabs(x) * .08f - .12f) {
1083  ret_val = 3;
1084  } else {
1085  ret_val = 4;
1086  }
1087  }
1088  return ret_val;
1089 } /* bfun3_ */
1090 
1133 /*** ROUTINE COMPUTES THE VOIGT FUNCTION Y/PI*INTEGRAL FROM ***/
1134 /*** - TO + INFINITY OF EXP(-T*T)/(Y*Y+(X-T)*(X-T)) DT ***/
1135 void lineshape_voigt_drayson(Vector& ls_attenuation,
1136  Vector&, //ls_phase
1137  Vector&, //ls_dattenuation_dfrequency_term
1138  Vector&, //ls_dphase_dfrequency_term
1139  Vector&, //ls_dattenuation_dpressure_term
1140  Vector&, //ls_dphase_dpressure_term
1141  Vector& x,
1142  const Numeric f0,
1143  const Numeric gamma,
1144  const Numeric,
1145  const Numeric,
1146  const Numeric,
1147  const Numeric,
1148  const Numeric sigma,
1149  const Numeric,
1150  ConstVectorView f_grid,
1151  const bool,
1152  const bool)
1153 
1154 {
1155  const Index nf = f_grid.nelem();
1156 
1157  // seems not necessary for Doppler correction
1158  // extern const Numeric SQRT_NAT_LOG_2;
1159 
1160  // PI
1161  extern const Numeric PI;
1162 
1163  // constant sqrt(1/pi)
1164  const Numeric sqrt_invPI = sqrt(1 / PI);
1165 
1166  // constant normalization factor for voigt
1167  Numeric fac = 1.0 / sigma * sqrt_invPI;
1168 
1169  Numeric B[22 + 1] = {0., 0., .7093602e-7};
1170  Numeric RI[15 + 1] = {0};
1171  const Numeric XN[15 + 1] = {
1172  0., 10., 9., 8., 8., 7., 6., 5., 4., 3., 3., 3., 3., 3., 3., 3.};
1173  const Numeric YN[15 + 1] = {
1174  0., .6, .6, .6, .5, .4, .4, .3, .3, .3, .3, 1., .9, .8, .7, .7};
1175  Numeric D0[25 + 1] = {0}, D1[25 + 1] = {0}, D2[25 + 1] = {0},
1176  D3[25 + 1] = {0}, D4[25 + 1] = {0};
1177  Numeric HN[25 + 1] = {0};
1178  Numeric H = .201;
1179  const Numeric XX[3 + 1] = {0., .5246476, 1.65068, .7071068};
1180  const Numeric HH[3 + 1] = {0., .2562121, .2588268e-1, .2820948};
1181  const Numeric NBY2[19 + 1] = {0., 9.5, 9., 8.5, 8., 7.5, 7., 6.5, 6., 5.5,
1182  5., 4.5, 4., 3.5, 3., 2.5, 2., 1.5, 1., .5};
1183  const Numeric C[21 + 1] = {
1184  0., .7093602e-7, -.2518434e-6, .856687e-6,
1185  -.2787638e-5, .866074e-5, -.2565551e-4, .7228775e-4,
1186  -.1933631e-3, .4899520e-3, -.1173267e-2, .2648762e-2,
1187  -.5623190e-2, .1119601e-1, -.2084976e-1, .3621573e-1,
1188  -.5851412e-1, .8770816e-1, -.121664, .15584,
1189  -.184, .2};
1190  Numeric CO = 0;
1191  Numeric U, V, UU, VV, Y2, DX;
1192  int I, J, K, MAX, MIN, N, i1;
1193  //int TRU = 0;
1194 
1195  // variables needed in original c routine:
1196 
1197  // Ratio of the Lorentz halfwidth to the Doppler halfwidth
1198  Numeric Y = gamma / sigma;
1199 
1200  // frequency in units of Doppler, Note: Drayson accepts only
1201  // positive values
1202  for (i1 = 0; i1 < (int)nf; i1++) {
1203  x[i1] = fabs((f_grid[i1] - f0)) / sigma;
1204  }
1205 
1206  //if (TRU) goto L104;
1207  /*** REGION I. COMPUTE DAWSON'S FUNCTION AT MESH POINTS ***/
1208  //TRU = 1;
1209  for (I = 1; I <= 15; I++) RI[I] = -I / 2.;
1210  for (I = 1; I <= 25; I++) {
1211  HN[I] = H * (I - .5);
1212  CO = 4. * HN[I] * HN[I] / 25. - 2.;
1213  for (J = 2; J <= 21; J++) B[J + 1] = CO * B[J] - B[J - 1] + C[J];
1214  D0[I] = HN[I] * (B[22] - B[21]) / 5.;
1215  D1[I] = 1. - 2. * HN[I] * D0[I];
1216  D2[I] = (HN[I] * D1[I] + D0[I]) / RI[2];
1217  D3[I] = (HN[I] * D2[I] + D1[I]) / RI[3];
1218  D4[I] = (HN[I] * D3[I] + D2[I]) / RI[4];
1219  }
1220  //L104:
1221  for (K = 0; K < (int)nf; K++) {
1222  if ((x[K] - 5.) < 0.)
1223  goto L105;
1224  else
1225  goto L112;
1226  L105:
1227  if ((Y - 1.) <= 0.)
1228  goto L110;
1229  else
1230  goto L106;
1231  L106:
1232  if (x[K] > 1.85 * (3.6 - Y)) goto L112;
1233  /*** REGION II: CONTINUED FRACTION. COMPUTE NUMBER OF TERMS NEEDED. ***/
1234  if (Y < 1.45) goto L107;
1235  I = (int)(Y + Y);
1236  goto L108;
1237  L107:
1238  I = (int)(11. * Y);
1239  L108:
1240  J = (int)(x[K] + x[K] + 1.85);
1241  MAX = (int)(XN[J] * YN[I] + .46);
1242  MIN = (int)((16 < 21 - 2 * MAX) ? 16 : 21 - 2 * MAX);
1243  /*** EVALUATE CONTINUED FRACTION ***/
1244  UU = Y;
1245  VV = x[K];
1246  for (J = MIN; J <= 19; J++) {
1247  U = NBY2[J] / (UU * UU + VV * VV);
1248  UU = Y + U * UU;
1249  VV = x[K] - U * VV;
1250  }
1251  ls_attenuation[K] = UU / (UU * UU + VV * VV) / 1.772454 * fac;
1252  continue;
1253  L110:
1254  Y2 = Y * Y;
1255  if (x[K] + Y >= 5.) goto L113;
1256  /*** REGION I. COMPUTE DAWSON'S FUNCTION AT X FROM TAYLOR SERIES. ***/
1257  N = (int)(x[K] / H);
1258  DX = x[K] - HN[N + 1];
1259  U = (((D4[N + 1] * DX + D3[N + 1]) * DX + D2[N + 1]) * DX + D1[N + 1]) *
1260  DX +
1261  D0[N + 1];
1262  V = 1. - 2. * x[K] * U;
1263  /*** TAYLOR SERIES EXPANSION ABOUT Y=0.0 ***/
1264  VV = exp(Y2 - x[K] * x[K]) * cos(2. * x[K] * Y) / 1.128379 - Y * V;
1265  UU = -Y;
1266  MAX = (int)(5. + (12.5 - x[K]) * .8 * Y);
1267  for (I = 2; I <= MAX; I += 2) {
1268  U = (x[K] * V + U) / RI[I];
1269  V = (x[K] * U + V) / RI[I + 1];
1270  UU = -UU * Y2;
1271  VV = VV + V * UU;
1272  }
1273  ls_attenuation[K] = 1.128379 * VV * fac;
1274  continue;
1275  L112:
1276  Y2 = Y * Y;
1277  if (Y < 11. - .6875 * x[K]) goto L113;
1278  /*** REGION IIIB: 2-POINT GAUSS-HERMITE QUADRATURE. ***/
1279  U = x[K] - XX[3];
1280  V = x[K] + XX[3];
1281  ls_attenuation[K] = Y * (HH[3] / (Y2 + U * U) + HH[3] / (Y2 + V * V)) * fac;
1282  continue;
1283  /*** REGION IIIA: 4-POINT GAUSS-HERMITE QUADRATURE. ***/
1284  L113:
1285  U = x[K] - XX[1];
1286  V = x[K] + XX[1];
1287  UU = x[K] - XX[2];
1288  VV = x[K] + XX[2];
1289  ls_attenuation[K] = Y *
1290  (HH[1] / (Y2 + U * U) + HH[1] / (Y2 + V * V) +
1291  HH[2] / (Y2 + UU * UU) + HH[2] / (Y2 + VV * VV)) *
1292  fac;
1293  continue;
1294  }
1295 }
1296 
1308 void chi_cousin(Numeric& chi, const Numeric& df) {
1309  // Conversion factor from Hz to cm-1
1310  extern const Numeric HZ2CM;
1311 
1312  const Numeric n2 = 0.79;
1313  const Numeric o2 = 0.21;
1314 
1315  const Numeric df_cm = df * HZ2CM;
1316  const Numeric df_cm_abs = fabs(df_cm);
1317 
1318  chi = 0;
1319 
1320  // N2 term
1321  if (df_cm_abs <= 5) {
1322  chi += n2;
1323  } else if (df_cm_abs <= 22) {
1324  chi += n2 * 1.968 * exp(-0.1354 * df_cm_abs);
1325  } else if (df_cm_abs <= 50) {
1326  chi += n2 * 0.160 * exp(-0.0214 * df_cm_abs);
1327  } else {
1328  chi += n2 * 0.162 * exp(-0.0216 * df_cm_abs);
1329  }
1330 
1331  // O2 term
1332  if (df_cm_abs <= 5) {
1333  chi += o2;
1334  } else if (df_cm_abs <= 22) {
1335  chi += o2 * 1.968 * exp(-0.1354 * df_cm_abs);
1336  } else if (df_cm_abs <= 50) {
1337  chi += o2 * 0.160 * exp(-0.0214 * df_cm_abs);
1338  } else {
1339  chi += o2 * 0.162 * exp(-0.0216 * df_cm_abs);
1340  }
1341 }
1342 
1352 void lineshape_CO2_lorentz(Vector& ls_attenuation,
1353  Vector&, //ls_phase
1354  Vector&, //ls_dattenuation_dfrequency_term
1355  Vector&, //ls_dattenuation_dfrequency_term
1356  Vector&, //ls_dphase_dpressure_term
1357  Vector&, //ls_dphase_dpressure_term
1358  Vector&, //X
1359  const Numeric f0,
1360  const Numeric gamma,
1361  const Numeric,
1362  const Numeric,
1363  const Numeric,
1364  const Numeric,
1365  const Numeric, //sigma
1366  const Numeric,
1367  ConstVectorView f_grid,
1368  const bool,
1369  const bool) {
1370  const Index nf = f_grid.nelem();
1371 
1372  // PI:
1373  extern const Numeric PI;
1374 
1375  // Some constant variables
1376  const Numeric gamma2 = gamma * gamma;
1377  const Numeric fac = gamma / PI;
1378 
1379  for (Index i = 0; i < nf; ++i) {
1380  // f-f0
1381  const Numeric df = f_grid[i] - f0;
1382 
1383  // The chi factor
1384  Numeric chi;
1385  chi_cousin(chi, df);
1386 
1387  // chi * Lorentz
1388  ls_attenuation[i] = chi * fac / (df * df + gamma2);
1389  }
1390 }
1391 
1403 void lineshape_CO2_drayson(Vector& ls_attenuation,
1404  Vector&, //ls_phase,
1405  Vector&, //ls_dattenuation_dfrequency_term
1406  Vector&, //ls_dphase_dfrequency_term
1407  Vector&, //ls_dattenuation_dpressure_term
1408  Vector&, //ls_dphase_dpressure_term
1409  Vector& X,
1410  const Numeric f0,
1411  const Numeric gamma,
1412  const Numeric,
1413  const Numeric,
1414  const Numeric,
1415  const Numeric,
1416  const Numeric sigma,
1417  const Numeric,
1418  ConstVectorView f_grid,
1419  const bool,
1420  const bool) {
1421  Vector tmp(0);
1422  lineshape_voigt_drayson(ls_attenuation,
1423  tmp,
1424  tmp,
1425  tmp,
1426  tmp,
1427  tmp,
1428  X,
1429  f0,
1430  gamma,
1431  0.,
1432  1.,
1433  0.,
1434  0.,
1435  sigma,
1436  0.,
1437  f_grid,
1438  0,
1439  0);
1440 
1441  const Index nf = f_grid.nelem();
1442  for (Index i = 0; i < nf; ++i) {
1443  // f-f0
1444  const Numeric df = f_grid[i] - f0;
1445 
1446  // The chi factor
1447  Numeric chi;
1448  chi_cousin(chi, df);
1449 
1450  // chi * Lorentz
1451  ls_attenuation[i] *= chi;
1452  }
1453 }
1454 
1476 void faddeeva_algorithm_916(Vector& ls_attenuation,
1477  Vector& ls_phase,
1478  Vector& ls_dattenuation_dfrequency_term,
1479  Vector& ls_dphase_dfrequency_term,
1480  Vector& ls_dattenuation_dpressure_term,
1481  Vector& ls_dphase_dpressure_term,
1482  Vector&,
1483  const Numeric f0,
1484  const Numeric gamma,
1485  const Numeric,
1486  const Numeric,
1487  const Numeric,
1488  const Numeric,
1489  const Numeric sigma,
1490  const Numeric,
1491  ConstVectorView f_grid,
1492  const bool do_phase,
1493  const bool do_partials)
1494 
1495 {
1496  const Index nf = f_grid.nelem();
1497 
1498  // PI
1499  extern const Numeric PI;
1500 
1501  // constant sqrt(1/pi)
1502  static const Numeric sqrt_invPI = sqrt(1 / PI);
1503 
1504  // constant normalization factor for voigt
1505  const Numeric fac = sqrt_invPI / sigma;
1506 
1507  // Ratio of the Lorentz halfwidth to the Doppler halfwidth
1508  const Numeric y = gamma / (sigma);
1509 
1510  // frequency in units of Doppler
1511  for (Index ii = 0; ii < nf; ii++) {
1512  const Numeric x = (f_grid[ii] - f0) / sigma;
1513  const Complex w = Faddeeva::w(Complex(x, y));
1514 
1515  ls_attenuation[ii] = fac * w.real();
1516  if (do_phase || do_partials) ls_phase[ii] = fac * w.imag();
1517  if (do_partials) {
1518  // Derivatives are from a paper that gives w(y-ix) = Fa + iFb but our formalism use
1519  // w(x+iy) = Fa + iFb. Thus signs on Fb-derivatives are difficult...
1520 
1521  // dFa/dx
1522  ls_dattenuation_dfrequency_term[ii] =
1523  2. * (y * ls_phase[ii] - x * ls_attenuation[ii]);
1524 
1525  // dFa/dy
1526  ls_dattenuation_dpressure_term[ii] =
1527  2. * (y * ls_attenuation[ii] + x * ls_phase[ii] - fac * sqrt_invPI);
1528 
1529  if (do_phase) {
1530  // dFb/dx
1531  ls_dphase_dfrequency_term[ii] = -ls_dattenuation_dpressure_term[ii];
1532 
1533  // dFb/dy
1534  ls_dphase_dpressure_term[ii] = ls_dattenuation_dfrequency_term[ii];
1535  }
1536  }
1537  }
1538 }
1539 
1540 // Helpers connected to internal partials.
1541 void w_x_plus_iy_dT(
1542  Vector& dx_dT,
1543  Numeric& dy_dT,
1544  Numeric& dFu_dT,
1545  ConstVectorView f,
1546  const Numeric&
1547  f0, // Note that this is NOT the line center, but the shifted center
1548  const Numeric& sigma,
1549  const Numeric& dPF_dT, // pressure shift temperature derivative
1550  const Numeric& dDF_dT, // line mixing shift temperature derivative
1551  const Numeric& dsigma_dT, // derivative of sigma with temperature
1552  const Numeric& gamma,
1553  const Numeric& dgamma_dT) // derivative of gamma with temperature
1554 {
1555  // Function is w(x+iy), and dw/dx and dw/dy are both known already.
1556  // This function serves to find dx/dT
1557 
1558  // x = (f-f0) / sigma;
1559 
1560  // f0 contains many different terms:
1561  // 1) Magnetic field splitting (not T-dependent)
1562  // 2) Pressure frequency shift (T-dependent)
1563  // 3) Pressure line mixing shift (T-dependent)
1564  // f contains the frequency grid influenced by wind
1565  // sigma is the halfwidth Doppler broadening
1566 
1567  // Since both denominator and nominator depends on temperature, we need to split the derivative into two terms
1568  // d(f-f0)/dT * 1/sigma and (f-f0) * d(1/sigma)/dT
1569 
1570  for (Index iv = 0; iv < f.nelem(); iv++)
1571  dx_dT[iv] = (-dPF_dT - dDF_dT - (f[iv] - f0) / sigma * dsigma_dT) / sigma;
1572 
1573  // y = gamma / sigma
1574 
1575  // Since both denominator and nominator depends on temperature, we need to split the derivative into two terms
1576  // d(gamma)/dT * 1/sigma and (gamma) * d(1/sigma)/dT
1577 
1578  dy_dT = (dgamma_dT - gamma / sigma * dsigma_dT) / sigma;
1579 
1580  // This line shape is normalization depends on temperature and is
1581  dFu_dT = -1.0 / sigma * dsigma_dT;
1582  // NOTE: The factor is 1/sqrt(pi)/sigma, its derivative is then - 1 / sqrt(pi) / sigma^2 * dsigma_dT, but
1583  // the factor is already in the lineshape that this will be multiplied to, so
1584  // derivative divided by the factor and only the above remains! Which should only be -1/2T for this case.
1585 }
1586 void w_x_plus_iy_dF(Numeric& dx_dF, const Numeric& sigma) {
1587  // Function is w(x+iy), and dw/dx and dw/dy are both known already.
1588  // This function serves to find dx/dT
1589 
1590  // x = (f-f0) / sigma;
1591 
1592  dx_dF = 1.0 / sigma;
1593 
1594  // No other terms depend on f
1595 }
1596 void w_x_plus_iy_dF0(Numeric& dx_dF0, const Numeric& sigma) {
1597  // Function is w(x+iy), and dw/dx and dw/dy are both known already.
1598  // This function serves to find dx/dT
1599 
1600  // x = (f-f0) / sigma;
1601 
1602  dx_dF0 = -1.0 / sigma;
1603 
1604  // No other terms depend on f
1605 }
1606 void w_x_plus_iy_dgamma(Numeric& dy_dgamma, const Numeric& sigma) {
1607  // Function is w(x+iy), and dw/dx and dw/dy are both known already.
1608  // This function serves to find dy/dgamma
1609 
1610  // y = gamma / sigma;
1611 
1612  dy_dgamma = 1.0 / sigma;
1613 
1614  // No other terms depend on gamma
1615 }
1616 void w_x_plus_iy_dH(Numeric& dx_dH,
1617  const Numeric& sigma,
1618  const Numeric& df0_dH) {
1619  // Function is w(x+iy), and dw/dx and dw/dy are both known already.
1620  // This function serves to find dx/dmag
1621 
1622  // x = (f-f0) / sigma;
1623 
1624  // f0 contains many different terms:
1625  // 1) Magnetic field splitting (mag-dependent)
1626  // 2) Pressure frequency shift (not mag-dependent)
1627  // 3) Pressure line mixing shift (not mag-dependent)
1628  // f contains the frequency grid influenced by wind
1629  // sigma is the halfwidth Doppler broadening
1630 
1631  dx_dH = -df0_dH / sigma;
1632 
1633  // No other terms depend on mag
1634 }
1635 void w_x_plus_iy_dDF(Numeric& dx_dDF, const Numeric& sigma) {
1636  // Function is w(x+iy), and dw/dx and dw/dy are both known already.
1637  // This function serves to find dx/dDF
1638 
1639  // x = (f-f0) / sigma;
1640 
1641  // f0 contains many different terms:
1642  // 1) Magnetic field splitting (not DF-dependent)
1643  // 2) Pressure frequency shift (not DF-dependent)
1644  // 3) Pressure line mixing shift (DF-dependent)
1645  // f contains the frequency grid influenced by wind
1646  // sigma is the halfwidth Doppler broadening
1647 
1648  dx_dDF = -1.0 / sigma;
1649 
1650  // No other terms depend on DF
1651 }
1652 
1682 void hui_etal_1978_lineshape(Vector& ls_attenuation,
1683  Vector& ls_phase,
1684  Vector&, //ls_dattenuation_dfrequency_term
1685  Vector&, //ls_dphase_dfrequency_term
1686  Vector&, //ls_dattenuation_dpressure_term
1687  Vector&, //ls_dphase_dpressure_term
1688  Vector& xvector,
1689  const Numeric f0,
1690  const Numeric gamma,
1691  const Numeric,
1692  const Numeric,
1693  const Numeric,
1694  const Numeric,
1695  const Numeric sigma,
1696  const Numeric,
1697  ConstVectorView f_grid,
1698  const bool do_phase,
1699  const bool)
1700 
1701 {
1702  const Index nf = f_grid.nelem();
1703 
1704  // PI
1705  extern const Numeric PI;
1706 
1707  // constant sqrt(1/pi)
1708  const Numeric sqrt_invPI = sqrt(1 / PI);
1709 
1710  // constant normalization factor for voigt
1711  const Numeric fac = sqrt_invPI / sigma;
1712 
1713  // Ratio of the Lorentz halfwidth to the Doppler halfwidth
1714  const Numeric y = gamma / (sigma);
1715 
1716  // frequency in units of Doppler
1717  for (Index ii = 0; ii < nf; ii++) {
1718  xvector[ii] = (f_grid[ii] - f0) / sigma;
1719 
1720  // Note that this works but I don't know why there is a difference
1721  // between the theory described above and this practical implementation.
1722  const Complex z(y, -xvector[ii]);
1723 
1724  const Complex A =
1725  (((((.5641896 * z + 5.912626) * z + 30.18014) * z + 93.15558) * z +
1726  181.9285) *
1727  z +
1728  214.3824) *
1729  z +
1730  122.6079;
1731  const Complex B =
1732  ((((((z + 10.47986) * z + 53.99291) * z + 170.3540) * z + 348.7039) *
1733  z +
1734  457.3345) *
1735  z +
1736  352.7306) *
1737  z +
1738  122.6079;
1739  const Complex C = A / B;
1740 
1741  ls_attenuation[ii] = fac * C.real();
1742  if (do_phase) ls_phase[ii] = fac * C.imag();
1743  }
1744 }
1745 
1771 void hartmann_tran_lineshape(Vector& ls_attenuation,
1772  Vector& ls_phase,
1773  Vector& ls_dattenuation_dfrequency_term,
1774  Vector& ls_dphase_dfrequency_term,
1775  Vector& ls_dattenuation_dpressure_term,
1776  Vector& ls_dphase_dpressure_term,
1777  Vector& xvector,
1778  const Numeric f0,
1779  const Numeric gamma_0,
1780  const Numeric gamma_2, //untested
1781  const Numeric eta, //untested
1782  const Numeric df_0, //untested
1783  const Numeric df_2, //untested
1784  const Numeric gamma_D,
1785  const Numeric f_VC, //untested
1786  ConstVectorView f_grid,
1787  const bool do_phase,
1788  const bool do_partials)
1789 
1790 {
1791  // NOTE: f_VC is the only non-scaled variable in the reference document.
1792  // If given as something other than cm-1 it might behave strange. Still not tested,
1793  // just a future "fixme". The others are probably OK, but I am very unsure on this
1794  // as well.
1795 
1796  // NOTE: need more outputs for the partial derivatives? This requires redesign so that
1797  // there is no longer such abstract derivatives. That would slow the other methods down.
1798  // So a better solution is to just allow complex vector outputs, and then save:
1799  // w(iZ-) and w(iZ+). Their partials can be found from knowing these two --- Y and X has to
1800  // be re-derived in code. The Forthomme etal 2015 method of finding the outputs with
1801  // regards to the other parameters of the equations can be done. See their supplementary
1802  // material for how to calculate the many partial derivatives (or redo the work yourself).
1803  // This does mean that, internally, this function will output derivatives that have
1804  // a very different meaning from all the others.
1805  // This might not work out, since Line mixing makes mixing of these parameters excruciating.
1806 
1807  const bool test1 =
1808  gamma_2 == df_2 && gamma_2 == f_VC && gamma_2 == eta && gamma_2 == 0;
1809 
1810  if (gamma_0 == df_0 && gamma_0 == gamma_2 && test1) {
1811  throw std::runtime_error(
1812  "Doppler line shape not supported in ARTS using HTP.\n");
1813  } else if (
1814  test1) // Standard Voigt function so use standard Voigt function call. NOTE: Partial derivations must know this.
1815  {
1816  faddeeva_algorithm_916(ls_attenuation,
1817  ls_phase,
1818  ls_dattenuation_dfrequency_term,
1819  ls_dphase_dfrequency_term,
1820  ls_dattenuation_dpressure_term,
1821  ls_dphase_dpressure_term,
1822  xvector,
1823  f0,
1824  gamma_0,
1825  gamma_2,
1826  eta,
1827  df_0,
1828  df_2,
1829  gamma_D,
1830  f_VC,
1831  f_grid,
1832  do_phase,
1833  do_partials);
1834  return;
1835  } else //This is mostly untested and there are likely errors. Like for c_2 == 0 and for eta == 1
1836  {
1837  if (do_partials)
1838  throw std::runtime_error(
1839  "It is not yet supported to do partial derivation with\n"
1840  "HTP line shape parameters.");
1841 
1842  // Constants
1843  extern const Numeric PI;
1844  extern const Numeric SPEED_OF_LIGHT;
1845  const Complex i(0., 1.);
1846 
1847  // Direct pressure term
1848  const Complex c_0(gamma_0, df_0);
1849 
1850  // Indirect pressure term
1851  const Complex c_2(gamma_2, df_2);
1852 
1853  // Constant c_0 helper
1854  const Complex c_0_tilde = (1. - eta) * (c_0 - 1.5 * c_2) + f_VC;
1855 
1856  // Constant c_2 helper
1857  const Complex c_2_tilde = (1. - eta) * c_2;
1858 
1859  // Average speed of molecule
1860  const Numeric v_a0 = SPEED_OF_LIGHT * gamma_D / f0;
1861 
1862  // Constant sqrt(pi)
1863  const Numeric sqrtPI = sqrt(PI);
1864 
1865  // Some calculation constants
1866  const Complex c1 = sqrtPI / gamma_D;
1867  const Complex v_a02 = v_a0 * v_a0;
1868  const Complex c4 = f_VC - eta * (c_0 - 1.5 * c_2);
1869  const Complex c5 = eta * c_2 / v_a02;
1870 
1871  // Special case when no speed dependency or no correlation
1872  if (c_2_tilde == Complex(0., 0.)) {
1873  for (Index nf = 0; nf < f_grid.nelem(); nf++) {
1874  const Complex Z1 = Complex(0., (f0 - f_grid[nf]) / gamma_D) + c_0_tilde;
1875 
1876  // NB. Z1 → infinity is treated specially in Tran etal. This I think already happens in
1877  // the Faddeeeva function as well. So I will do nothing special
1878 
1879  const Complex w1 = Faddeeva::w(i * Z1);
1880 
1881  const Complex A = c1 * w1;
1882 
1883  const Complex B = c1 * v_a02 * ((1. - Z1 * Z1) * w1 + Z1 / sqrtPI);
1884 
1885  // Hartmann-Tran line shape
1886  const Complex HTP = 1. / PI * A / (1. - c4 * A + c5 * B);
1887 
1888  // Output
1889  ls_attenuation[nf] = HTP.real();
1890  if (do_phase || do_partials) ls_phase[nf] = HTP.imag();
1891  }
1892  return;
1893  }
1894 
1895  // Y variable otherwise
1896  const Complex sqrtY = gamma_D / (2. * c_2_tilde);
1897  const Complex Y = sqrtY * sqrtY;
1898  const Numeric absY = abs(Y);
1899  const Complex sqrtY_t2 = 2. * sqrtY;
1900  const Complex c2 = v_a02 / c_2_tilde / c_2_tilde;
1901  const Complex c3 = sqrtPI / sqrtY_t2;
1902 
1903  for (Index nf = 0; nf < f_grid.nelem(); nf++) {
1904  const Numeric df = f0 - f_grid[nf];
1905 
1906  // X variable
1907  const Complex X = (Complex(0., df) + c_0_tilde) / c_2_tilde;
1908 
1909  // Tran etal says this must be between 3e-8 and 10e15 for numerical stability in their routine.
1910  // This is not the case for us since we have a different code, but is arbitrarily used anyways.
1911  const Numeric numerical_test = abs(X) / absY;
1912 
1913  // Same variables in all cases
1914  Complex A, B;
1915  if (numerical_test < 3e-8) {
1916  // Z_plus
1917  // Note that sqrt(X+Y)+sqrt(Y) → 2*sqrt(Y) in this case
1918  // This simplification is ignored below but is faster. Should it be used?
1919  const Complex Z_plus = sqrt(X + Y) + sqrtY;
1920 
1921  // Z_minus
1922  const Complex Z_minus = Complex(0., df / gamma_D) + c_0_tilde;
1923 
1924  // w(Z_minus)
1925  const Complex w_plus = Faddeeva::w(i * Z_plus);
1926 
1927  // w(Z_minus)
1928  const Complex w_minus = Faddeeva::w(i * Z_minus);
1929 
1930  // A
1931  A = c1 * (w_minus - w_plus);
1932 
1933  // B
1934  B = c2 * (-1. + c3 * ((1. - Z_minus * Z_minus) * w_minus -
1935  (1. - Z_plus * Z_plus) * w_plus));
1936  } else if (numerical_test > 1e15) {
1937  // Note that the Z:s are very similar but
1938  // Tran etal still makes a difference between them
1939  const Complex sqrtX = sqrt(X);
1940  const Complex Z2 = sqrt(X + Y);
1941 
1942  const Complex w1 = Faddeeva::w(i * sqrtX);
1943  const Complex w2 = Faddeeva::w(i * Z2);
1944 
1945  // Note that there is a warning for these equations as X → infinity that is ignored for now.
1946 
1947  Complex Aconst = (1. / sqrtPI - sqrtX * w1);
1948 
1949  // A
1950  A = 2. * c1 * Aconst;
1951 
1952  // B
1953  B = c2 * (-1. + 2. * sqrtPI * (1. - X - 2. * Y) * Aconst +
1954  2. * sqrtPI * Z2 * w2);
1955 
1956  } else {
1957  // Z_plus
1958  const Complex Z_plus = sqrt(X + Y) + sqrtY;
1959 
1960  // Z_minus
1961  const Complex Z_minus = Z_plus - sqrtY_t2;
1962 
1963  // w(Z_minus)
1964  const Complex w_plus = Faddeeva::w(i * Z_plus);
1965 
1966  // w(Z_minus)
1967  const Complex w_minus = Faddeeva::w(i * Z_minus);
1968 
1969  // A
1970  A = c1 * (w_minus - w_plus);
1971 
1972  // B
1973  B = c2 * (-1. + c3 * ((1. - Z_minus * Z_minus) * w_minus -
1974  (1. - Z_plus * Z_plus) * w_plus));
1975  }
1976 
1977  // Hartmann-Tran line shape
1978  const Complex HTP = 1. / PI * A / (1. - c4 * A + c5 * B);
1979 
1980  // Output
1981  ls_attenuation[nf] = HTP.real();
1982  if (do_phase || do_partials) ls_phase[nf] = HTP.imag();
1983  }
1984  }
1985 }
1986 
1999 void lineshape_o2nonresonant(Vector& ls_attenuation,
2000  Vector&, //ls_phase _U_,
2001  Vector&, //ls_dattenuation_dfrequency_term
2002  Vector&, //ls_dphase_dfrequency_term
2003  Vector&, //ls_dattenuation_dpressure_term
2004  Vector&, //ls_dphase_dpressure_term
2005  Vector&, //X,
2006  const Numeric f0,
2007  const Numeric gamma,
2008  const Numeric,
2009  const Numeric,
2010  const Numeric,
2011  const Numeric,
2012  const Numeric, //sigma
2013  const Numeric,
2014  ConstVectorView f_grid,
2015  const bool,
2016  const bool) {
2017  const Index nf = f_grid.nelem();
2018 
2019  // PI:
2020  extern const Numeric PI;
2021 
2022  // assert( ls.nelem() == nf );
2023 
2024  Numeric fac = gamma / PI;
2025 
2026  for (Index i = 0; i < nf; ++i) {
2027  ls_attenuation[i] = fac / ((f_grid[i] - f0) * (f_grid[i] - f0));
2028  }
2029 }
2030 
2031 //------------------------------------------------------------------------
2032 // Normalization Functions
2033 //------------------------------------------------------------------------
2034 
2043 void lineshape_norm_no_norm(Vector& fac,
2044  const Numeric, //f0
2045  ConstVectorView, //f_grid
2046  const Numeric) //T
2047 {
2048  fac = 1.0;
2049 }
2050 void lineshape_norm_no_norm_dT(Vector& fac,
2051  const Numeric, //f0
2052  ConstVectorView, //f_grid
2053  const Numeric) //T
2054 {
2055  fac = 0.0;
2056 }
2057 void lineshape_norm_no_norm_dF(Vector& fac,
2058  const Numeric, //f0
2059  ConstVectorView, //f_grid
2060  const Numeric) //T
2061 {
2062  fac = 0.0;
2063 }
2064 void lineshape_norm_no_norm_dF0(Vector& fac,
2065  const Numeric, //f0
2066  ConstVectorView, //f_grid
2067  const Numeric) //T
2068 {
2069  fac = 0.0;
2070 }
2071 
2080 void lineshape_norm_linear(Vector& fac,
2081  const Numeric f0,
2082  ConstVectorView f_grid,
2083  const Numeric) //T
2084 {
2085  const Index nf = f_grid.nelem();
2086 
2087  // Abs(f0) is constant in the loop:
2088  const Numeric abs_f0 = abs(f0);
2089 
2090  for (Index i = 0; i < nf; ++i) {
2091  fac[i] = f_grid[i] / abs_f0;
2092  }
2093 }
2094 void lineshape_norm_linear_dT(Vector& fac,
2095  const Numeric, // f0
2096  ConstVectorView, //f_grid
2097  const Numeric) //T
2098 {
2099  fac = 0.0;
2100 }
2101 void lineshape_norm_linear_dF(Vector& fac,
2102  const Numeric f0,
2103  ConstVectorView f_grid,
2104  const Numeric) //T
2105 {
2106  const Index nf = f_grid.nelem();
2107 
2108  // Abs(f0) is constant in the loop:
2109  const Numeric abs_f0 = abs(f0);
2110 
2111  for (Index i = 0; i < nf; ++i) {
2112  fac[i] = 1.0 / abs_f0;
2113  }
2114 }
2115 void lineshape_norm_linear_dF0(Vector& fac,
2116  const Numeric f0,
2117  ConstVectorView f_grid,
2118  const Numeric) //T
2119 {
2120  const Index nf = f_grid.nelem();
2121 
2122  // Abs(f0) is constant in the loop:
2123  const Numeric dabs_f0 =
2124  -f0 / abs(f0) / abs(f0) / abs(f0); //-sign(f0)/abs(f0)^2 Why is this abs?
2125 
2126  for (Index i = 0; i < nf; ++i) {
2127  fac[i] = f_grid[i] * dabs_f0;
2128  }
2129 }
2130 
2145 void lineshape_norm_quadratic_Rosenkranz(Vector& fac,
2146  const Numeric f0,
2147  ConstVectorView f_grid,
2148  const Numeric T) {
2149  extern const Numeric PLANCK_CONST;
2150  extern const Numeric BOLTZMAN_CONST;
2151  const Numeric mafac =
2152  (PLANCK_CONST * f0) / (2.000e0 * BOLTZMAN_CONST * T) /
2153  sinh((PLANCK_CONST * f0) / (2.000e0 * BOLTZMAN_CONST * T));
2154  const Index nf = f_grid.nelem();
2155  // don't do this for the whole loop
2156  const Numeric f0_2 = f0 * f0;
2157 
2158  for (Index i = 0; i < nf; ++i) {
2159  fac[i] = mafac * (f_grid[i] * f_grid[i]) / f0_2;
2160  }
2161 }
2162 void lineshape_norm_quadratic_Rosenkranz_dT(Vector& fac,
2163  const Numeric f0,
2164  ConstVectorView f_grid,
2165  const Numeric T) {
2166  extern const Numeric PLANCK_CONST;
2167  extern const Numeric BOLTZMAN_CONST;
2168  const Numeric hf0 = PLANCK_CONST * f0, kT = T * BOLTZMAN_CONST;
2169  const Numeric sinh_term = sinh((hf0) / (2 * kT));
2170  const Numeric cosh_term = cosh((hf0) / (2.0 * kT));
2171  const Numeric dmafac_dT = -(hf0 * (2.0 * kT * sinh_term - hf0 * cosh_term)) /
2172  (4.0 * T * kT * kT * sinh_term * sinh_term);
2173  const Index nf = f_grid.nelem();
2174  // don't do this for the whole loop
2175  const Numeric f0_2 = f0 * f0;
2176 
2177  for (Index i = 0; i < nf; ++i) {
2178  fac[i] = dmafac_dT * (f_grid[i] * f_grid[i]) / f0_2;
2179  }
2180 }
2181 void lineshape_norm_quadratic_Rosenkranz_dF(Vector& fac,
2182  const Numeric f0,
2183  ConstVectorView f_grid,
2184  const Numeric T) {
2185  extern const Numeric PLANCK_CONST;
2186  extern const Numeric BOLTZMAN_CONST;
2187  const Numeric mafac =
2188  (PLANCK_CONST * f0) / (2.000e0 * BOLTZMAN_CONST * T) /
2189  sinh((PLANCK_CONST * f0) / (2.000e0 * BOLTZMAN_CONST * T));
2190  const Index nf = f_grid.nelem();
2191  // don't do this for the whole loop
2192  const Numeric f0_2 = f0 * f0;
2193 
2194  for (Index i = 0; i < nf; ++i) {
2195  fac[i] = mafac * (2.0 * f_grid[i]) / f0_2;
2196  }
2197 }
2198 void lineshape_norm_quadratic_Rosenkranz_dF0(Vector& fac,
2199  const Numeric f0,
2200  ConstVectorView f_grid,
2201  const Numeric T) {
2202  extern const Numeric PLANCK_CONST;
2203  extern const Numeric BOLTZMAN_CONST;
2204  const Numeric kT = 2.0 * BOLTZMAN_CONST * T;
2205  const Index nf = f_grid.nelem();
2206  // don't do this for the whole loop
2207  const Numeric term1 = sinh((f0 * PLANCK_CONST) / kT);
2208 
2209  const Numeric ddenom =
2210  -(PLANCK_CONST *
2211  (kT * term1 + f0 * PLANCK_CONST * cosh((f0 * PLANCK_CONST) / kT))) /
2212  (f0 * f0 * kT * kT * term1 * term1);
2213 
2214  for (Index i = 0; i < nf; ++i) {
2215  fac[i] = (f_grid[i] * f_grid[i]) * ddenom;
2216  }
2217 }
2218 
2233 void lineshape_norm_VVH(Vector& fac,
2234  const Numeric f0,
2235  ConstVectorView f_grid,
2236  const Numeric T) {
2237  extern const Numeric PLANCK_CONST;
2238  extern const Numeric BOLTZMAN_CONST;
2239 
2240  const Index nf = f_grid.nelem();
2241 
2242  // 2kT is constant for the loop
2243  const Numeric kT = 2.0 * BOLTZMAN_CONST * T;
2244 
2245  // denominator is constant for the loop
2246  const Numeric denom = abs(f0) * tanh(PLANCK_CONST * abs(f0) / kT);
2247 
2248  for (Index i = 0; i < nf; ++i) {
2249  fac[i] = f_grid[i] * tanh(PLANCK_CONST * f_grid[i] / kT) / denom;
2250  }
2251 }
2252 void lineshape_norm_VVH_dT(Vector& fac,
2253  const Numeric f0,
2254  ConstVectorView f_grid,
2255  const Numeric T) {
2256  extern const Numeric PLANCK_CONST;
2257  extern const Numeric BOLTZMAN_CONST;
2258 
2259  const Index nf = f_grid.nelem();
2260 
2261  // Various constants for the loop
2262  const Numeric kT = BOLTZMAN_CONST * T;
2263  const Numeric hf0 = PLANCK_CONST * f0;
2264  const Numeric tanh_term2 = tanh((hf0) / (2.0 * kT));
2265  const Numeric tanh_term2_squared = tanh_term2 * tanh_term2;
2266 
2267  // denominator is constant for the loop
2268  const Numeric denom1 = 2.0 * T * kT * f0 * tanh_term2,
2269  denom2 = 2.0 * T * kT * tanh_term2_squared;
2270 
2271  for (Index i = 0; i < nf; ++i) {
2272  const Numeric hf = PLANCK_CONST * f_grid[i],
2273  tanh_term1 = tanh((hf) / (2.0 * kT)),
2274  tanh_term1_squared = tanh_term1 * tanh_term1;
2275 
2276  fac[i] = (f_grid[i] * hf * (tanh_term1_squared - 1.0)) / denom1 -
2277  (hf * tanh_term1 * (tanh_term2_squared - 1.0)) / denom2;
2278  }
2279 }
2280 void lineshape_norm_VVH_dF(Vector& fac,
2281  const Numeric f0,
2282  ConstVectorView f_grid,
2283  const Numeric T) {
2284  extern const Numeric PLANCK_CONST;
2285  extern const Numeric BOLTZMAN_CONST;
2286 
2287  const Index nf = f_grid.nelem();
2288 
2289  // Various constants for the loop
2290  const Numeric kT = BOLTZMAN_CONST * T, hf0 = PLANCK_CONST * f0,
2291  tanh_term2 = tanh((hf0) / (2.0 * kT)),
2292  denom = (2 * kT * f0 * tanh_term2);
2293 
2294  for (Index i = 0; i < nf; ++i) {
2295  const Numeric hf = PLANCK_CONST * f_grid[i],
2296  tanh_term1 = tanh((hf) / (2.0 * kT)),
2297  tanh_term1_squared = tanh_term1 * tanh_term1;
2298  fac[i] = (hf + 2 * kT * tanh_term1 - hf * tanh_term1_squared) / denom;
2299  }
2300 }
2301 void lineshape_norm_VVH_dF0(Vector& fac,
2302  const Numeric f0,
2303  ConstVectorView f_grid,
2304  const Numeric T) {
2305  extern const Numeric PLANCK_CONST;
2306  extern const Numeric BOLTZMAN_CONST;
2307 
2308  const Index nf = f_grid.nelem();
2309 
2310  // 2kT is constant for the loop
2311  const Numeric kT = 2.0 * BOLTZMAN_CONST * T,
2312  term1 = tanh((PLANCK_CONST * abs(f0)) / kT);
2313 
2314  // denominator is constant for the loop
2315  const Numeric ddenom = (PLANCK_CONST * f0 / abs(f0) * (term1 * term1 - 1)) /
2316  (kT * abs(f0) * term1 * term1) -
2317  f0 / abs(f0) / (abs(f0) * abs(f0) * term1);
2318 
2319  for (Index i = 0; i < nf; ++i) {
2320  fac[i] = f_grid[i] * tanh(PLANCK_CONST * f_grid[i] / kT) * ddenom;
2321  }
2322 }
2323 
2338 void lineshape_norm_VVW(Vector& fac,
2339  const Numeric f0,
2340  ConstVectorView f_grid,
2341  const Numeric) {
2342  const Index nf = f_grid.nelem();
2343 
2344  // denominator is constant for the loop
2345  const Numeric denom = f0 * f0;
2346 
2347  for (Index i = 0; i < nf; ++i) {
2348  fac[i] = f_grid[i] * f_grid[i] / denom;
2349  }
2350 }
2351 void lineshape_norm_VVW_dT(Vector& fac,
2352  const Numeric, // f0,
2353  ConstVectorView, // f_grid,
2354  const Numeric) {
2355  fac = 0.0;
2356 }
2357 void lineshape_norm_VVW_dF(Vector& fac,
2358  const Numeric f0,
2359  ConstVectorView f_grid,
2360  const Numeric) {
2361  const Index nf = f_grid.nelem();
2362 
2363  // denominator is constant for the loop
2364  const Numeric denom = f0 * f0;
2365 
2366  for (Index i = 0; i < nf; ++i) {
2367  fac[i] = 2.0 * f_grid[i] / denom;
2368  }
2369 }
2370 void lineshape_norm_VVW_dF0(Vector& fac,
2371  const Numeric f0,
2372  ConstVectorView f_grid,
2373  const Numeric) {
2374  const Index nf = f_grid.nelem();
2375 
2376  // denominator is constant for the loop
2377  const Numeric ddenom = -0.5 * f0 * f0 * f0;
2378 
2379  for (Index i = 0; i < nf; ++i) {
2380  fac[i] = f_grid[i] * f_grid[i] / ddenom;
2381  }
2382 }
2383 
2384 //------------------------------------------------------------------------
2385 // Available Lineshapes
2386 //------------------------------------------------------------------------
2387 
2389 namespace global_data {
2390 Array<LineshapeRecord> lineshape_data;
2391 }
2392 
2393 void define_lineshape_data() {
2394  using global_data::lineshape_data;
2395 
2396  const bool PHASE = true;
2397  const bool NO_PHASE = false;
2398  const bool PARTIALS = true;
2399  const bool NO_PARTIALS = false;
2400 
2401  // Initialize to empty, just in case.
2402  lineshape_data.resize(0);
2403 
2404  lineshape_data.push_back(LineshapeRecord(
2405  "no_shape",
2406  "This lineshape does nothing. It only exists, because formally\n"
2407  "you have to specify a lineshape also for continuum tags.",
2408  lineshape_no_shape,
2409  NO_PHASE,
2410  NO_PARTIALS));
2411 
2412  lineshape_data.push_back(LineshapeRecord("Lorentz",
2413  "The Lorentz line shape.",
2414  lineshape_lorentz,
2415  PHASE,
2416  NO_PARTIALS));
2417 
2418  lineshape_data.push_back(LineshapeRecord("Mirrored Lorentz",
2419  "The mirrored Lorentz line shape.",
2420  lineshape_mirrored_lorentz,
2421  PHASE,
2422  NO_PARTIALS));
2423 
2424  lineshape_data.push_back(LineshapeRecord("Doppler",
2425  "The Doppler line shape.",
2426  lineshape_doppler,
2427  NO_PHASE,
2428  NO_PARTIALS));
2429 
2430  lineshape_data.push_back(LineshapeRecord(
2431  "Voigt_Kuntz6",
2432  "The Voigt line shape. Approximation by Kuntz: Accuracy 2*10-6",
2433  lineshape_voigt_kuntz6,
2434  w_x_plus_iy_dT,
2435  w_x_plus_iy_dF,
2436  w_x_plus_iy_dF0,
2437  w_x_plus_iy_dgamma,
2438  w_x_plus_iy_dH,
2439  w_x_plus_iy_dDF,
2440  NO_PHASE,
2441  PARTIALS));
2442 
2443  lineshape_data.push_back(LineshapeRecord(
2444  "Voigt_Kuntz3",
2445  "The Voigt line shape. Approximation by Kuntz: Accuracy 2*10-3",
2446  deprecated,
2447  NO_PHASE,
2448  NO_PARTIALS));
2449 
2450  lineshape_data.push_back(LineshapeRecord(
2451  "Voigt_Kuntz4",
2452  "The Voigt line shape. Approximation by Kuntz: Accuracy 2*10-4",
2453  deprecated,
2454  NO_PHASE,
2455  NO_PARTIALS));
2456 
2457  lineshape_data.push_back(
2458  LineshapeRecord("Voigt_Drayson",
2459  "The Voigt line shape. Approximation by Drayson.",
2460  lineshape_voigt_drayson,
2461  NO_PHASE,
2462  NO_PARTIALS));
2463 
2464  lineshape_data.push_back(LineshapeRecord(
2465  "CO2_Lorentz",
2466  "Special line shape for CO2 in the infrared, neglecting Doppler\n"
2467  "broadening and details of line mixing. The line shape can be\n"
2468  "expressed as\n"
2469  " chi(f,f0) * Lorentz(f,f0) \n"
2470  "\n"
2471  "The chi-factor follows Cousin et al. 1985. Broadening by N2 and O2\n"
2472  "is considered, while self-broadening (CO2-CO2) is neglected."
2473  "\n"
2474  "NOTE: Temperature dependency is not yet included. The chi factor is\n"
2475  "valid for 238 K.",
2476  lineshape_CO2_lorentz,
2477  NO_PHASE,
2478  NO_PARTIALS));
2479 
2480  lineshape_data.push_back(LineshapeRecord(
2481  "CO2_Drayson",
2482  "Special line shape for CO2 in the infrared, neglecting details of\n"
2483  "line mixing. The line shape can be expressed as\n"
2484  " chi(f,f0) * Drayson(f,f0) \n"
2485  "\n"
2486  "The chi-factor follows Cousin et al. 1985. Broadening by N2 and O2\n"
2487  "is considered, while self-broadening (CO2-CO2) is neglected.\n"
2488  "\n"
2489  "NOTE: Temperature dependency is not yet included. The chi factor is\n"
2490  "valid for 238 K.",
2491  lineshape_CO2_drayson,
2492  NO_PHASE,
2493  NO_PARTIALS));
2494 
2495  lineshape_data.push_back(LineshapeRecord(
2496  "Faddeeva_Algorithm_916",
2497  "Voigt and Faraday-Voigt function as per Faddeeva function solution by JPL.\n"
2498  "Line shape is considered from w(z)=exp(-z^2)*erfc(-iz) where z=v'+ia, and \n"
2499  "v' is a Doppler weighted freqeuncy parameter and a is a Doppler weighted \n"
2500  "pressure parameter.",
2501  faddeeva_algorithm_916,
2502  w_x_plus_iy_dT,
2503  w_x_plus_iy_dF,
2504  w_x_plus_iy_dF0,
2505  w_x_plus_iy_dgamma,
2506  w_x_plus_iy_dH,
2507  w_x_plus_iy_dDF,
2508  PHASE,
2509  PARTIALS));
2510 
2511  lineshape_data.push_back(LineshapeRecord(
2512  "Hartmann-Tran",
2513  "Line shape is considered as described by the Hartmann-Tran profile.",
2514  hartmann_tran_lineshape,
2515  PHASE,
2516  NO_PARTIALS));
2517 
2518  lineshape_data.push_back(LineshapeRecord(
2519  "Hui_etal_1978",
2520  "Classic line shape. Solving the complex error function returns both parts\n"
2521  "of the refractive index.",
2522  hui_etal_1978_lineshape,
2523  PHASE,
2524  NO_PARTIALS));
2525 
2526  lineshape_data.push_back(
2527  LineshapeRecord("O2NonResonant",
2528  "Special line shape. Only use for non-resonant O2.",
2529  lineshape_o2nonresonant,
2530  NO_PHASE,
2531  NO_PARTIALS));
2532 }
2533 
2536 namespace global_data {
2537 Array<LineshapeNormRecord> lineshape_norm_data;
2538 }
2539 
2540 void define_lineshape_norm_data() {
2541  using global_data::lineshape_norm_data;
2542 
2543  // Initialize to empty, just in case.
2544  lineshape_norm_data.resize(0);
2545 
2546  lineshape_norm_data.push_back(
2547  LineshapeNormRecord("no_norm",
2548  "No normalization of the lineshape.",
2549  lineshape_norm_no_norm,
2550  lineshape_norm_no_norm_dT,
2551  lineshape_norm_no_norm_dF,
2552  lineshape_norm_no_norm_dF0));
2553 
2554  lineshape_norm_data.push_back(LineshapeNormRecord(
2555  "Rosenkranz_quadratic",
2556  "Quadratic normalization of the lineshape with (f/f0)^2*h*f0/(2*k*T)/sinh(h*f0/(2*k*T)).",
2557  lineshape_norm_quadratic_Rosenkranz,
2558  lineshape_norm_quadratic_Rosenkranz_dT,
2559  lineshape_norm_quadratic_Rosenkranz_dF,
2560  lineshape_norm_quadratic_Rosenkranz_dF0));
2561 
2562  lineshape_norm_data.push_back(LineshapeNormRecord(
2563  "VVH",
2564  "Van Vleck Huber normalization of the lineshape with\n"
2565  " (f*tanh(h*f/(2*k*T))) / (f0*tanh(h*f0/(2*k*T))).\n"
2566  " The denominator is a result of catalogue intensities.",
2567  lineshape_norm_VVH,
2568  lineshape_norm_VVH_dT,
2569  lineshape_norm_VVH_dF,
2570  lineshape_norm_VVH_dF0));
2571 
2572  lineshape_norm_data.push_back(LineshapeNormRecord(
2573  "VVW",
2574  "Van Vleck Weiskopf normalization of the lineshape with (f*f) / (f0*f0).\n",
2575  lineshape_norm_VVW,
2576  lineshape_norm_VVW_dT,
2577  lineshape_norm_VVW_dF,
2578  lineshape_norm_VVW_dF0));
2579 }
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].
N
#define N
Definition: rng.cc:164
lineshape_norm_VVW_dF
void lineshape_norm_VVW_dF(Vector &fac, const Numeric f0, ConstVectorView f_grid, const Numeric)
Definition: lineshapes.cc:2357
global_data::lineshape_norm_data
Array< LineshapeNormRecord > lineshape_norm_data
Definition: lineshapes.cc:2537
lineshape_no_shape
void lineshape_no_shape(Vector &, Vector &, Vector &, Vector &, Vector &, Vector &, Vector &, const Numeric, const Numeric, const Numeric, const Numeric, const Numeric, const Numeric, const Numeric, const Numeric, ConstVectorView, const bool, const bool)
Definition: lineshapes.cc:61
complex.h
A class implementing complex numbers for ARTS.
x2
#define x2
Definition: linefunctiondata.h:254
C
#define C(a, b)
Definition: Faddeeva.cc:255
lineshape_CO2_drayson
void lineshape_CO2_drayson(Vector &ls_attenuation, Vector &, Vector &, Vector &, Vector &, Vector &, Vector &X, const Numeric f0, const Numeric gamma, const Numeric, const Numeric, const Numeric, const Numeric, const Numeric sigma, const Numeric, ConstVectorView f_grid, const bool, const bool)
Definition: lineshapes.cc:1403
faddeeva_algorithm_916
void faddeeva_algorithm_916(Vector &ls_attenuation, Vector &ls_phase, Vector &ls_dattenuation_dfrequency_term, Vector &ls_dphase_dfrequency_term, Vector &ls_dattenuation_dpressure_term, Vector &ls_dphase_dpressure_term, Vector &, const Numeric f0, const Numeric gamma, const Numeric, const Numeric, const Numeric, const Numeric, const Numeric sigma, const Numeric, ConstVectorView f_grid, const bool do_phase, const bool do_partials)
Definition: lineshapes.cc:1476
lineshape_norm_linear_dF
void lineshape_norm_linear_dF(Vector &fac, const Numeric f0, ConstVectorView f_grid, const Numeric)
Definition: lineshapes.cc:2101
lineshape_o2nonresonant
void lineshape_o2nonresonant(Vector &ls_attenuation, Vector &, Vector &, Vector &, Vector &, Vector &, Vector &, const Numeric f0, const Numeric gamma, const Numeric, const Numeric, const Numeric, const Numeric, const Numeric, const Numeric, ConstVectorView f_grid, const bool, const bool)
Definition: lineshapes.cc:1999
BOLTZMAN_CONST
const Numeric BOLTZMAN_CONST
Global constant, the Boltzmann constant [J/K].
deprecated
void deprecated(Vector &, Vector &, Vector &, Vector &, Vector &, Vector &, Vector &, const Numeric, const Numeric, const Numeric, const Numeric, const Numeric, const Numeric, const Numeric, const Numeric, ConstVectorView, const bool, const bool)
Definition: lineshapes.cc:252
lineshape_norm_VVW_dT
void lineshape_norm_VVW_dT(Vector &fac, const Numeric, ConstVectorView, const Numeric)
Definition: lineshapes.cc:2351
Vector
The Vector class.
Definition: matpackI.h:860
lineshape_norm_quadratic_Rosenkranz
void lineshape_norm_quadratic_Rosenkranz(Vector &fac, const Numeric f0, ConstVectorView f_grid, const Numeric T)
Definition: lineshapes.cc:2145
abs
#define abs(x)
Definition: legacy_continua.cc:20626
fac
Numeric fac(const Index n)
fac
Definition: math_funcs.cc:63
lineshape_norm_VVH_dF0
void lineshape_norm_VVH_dF0(Vector &fac, const Numeric f0, ConstVectorView f_grid, const Numeric T)
Definition: lineshapes.cc:2301
test1
int test1()
Definition: test_matpack.cc:47
lineshape_norm_no_norm_dF
void lineshape_norm_no_norm_dF(Vector &fac, const Numeric, ConstVectorView, const Numeric)
Definition: lineshapes.cc:2057
w
cmplx FADDEEVA() w(cmplx z, double relerr)
Definition: Faddeeva.cc:680
bfun6_
long bfun6_(Numeric y, Numeric x)
Definition: lineshapes.cc:279
min
#define min(a, b)
Definition: legacy_continua.cc:20628
lineshape_norm_VVH
void lineshape_norm_VVH(Vector &fac, const Numeric f0, ConstVectorView f_grid, const Numeric T)
Definition: lineshapes.cc:2233
lineshape_norm_VVH_dF
void lineshape_norm_VVH_dF(Vector &fac, const Numeric f0, ConstVectorView f_grid, const Numeric T)
Definition: lineshapes.cc:2280
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:51
define_lineshape_data
void define_lineshape_data()
Definition: lineshapes.cc:2393
b2
#define b2
Definition: complex.h:59
array.h
This file contains the definition of Array.
w_x_plus_iy_dF
void w_x_plus_iy_dF(Numeric &dx_dF, const Numeric &sigma)
Definition: lineshapes.cc:1586
lineshape_norm_no_norm
void lineshape_norm_no_norm(Vector &fac, const Numeric, ConstVectorView, const Numeric)
Definition: lineshapes.cc:2043
lineshape_norm_quadratic_Rosenkranz_dF0
void lineshape_norm_quadratic_Rosenkranz_dF0(Vector &fac, const Numeric f0, ConstVectorView f_grid, const Numeric T)
Definition: lineshapes.cc:2198
arts.h
The global header file for ARTS.
hartmann_tran_lineshape
void hartmann_tran_lineshape(Vector &ls_attenuation, Vector &ls_phase, Vector &ls_dattenuation_dfrequency_term, Vector &ls_dphase_dfrequency_term, Vector &ls_dattenuation_dpressure_term, Vector &ls_dphase_dpressure_term, Vector &xvector, const Numeric f0, const Numeric gamma_0, const Numeric gamma_2, const Numeric eta, const Numeric df_0, const Numeric df_2, const Numeric gamma_D, const Numeric f_VC, ConstVectorView f_grid, const bool do_phase, const bool do_partials)
Definition: lineshapes.cc:1771
lineshape_voigt_drayson
void lineshape_voigt_drayson(Vector &ls_attenuation, Vector &, Vector &, Vector &, Vector &, Vector &, Vector &x, const Numeric f0, const Numeric gamma, const Numeric, const Numeric, const Numeric, const Numeric, const Numeric sigma, const Numeric, ConstVectorView f_grid, const bool, const bool)
Definition: lineshapes.cc:1135
w_x_plus_iy_dDF
void w_x_plus_iy_dDF(Numeric &dx_dDF, const Numeric &sigma)
Definition: lineshapes.cc:1635
Complex
std::complex< Numeric > Complex
Definition: complex.h:33
lineshape_norm_no_norm_dF0
void lineshape_norm_no_norm_dF0(Vector &fac, const Numeric, ConstVectorView, const Numeric)
Definition: lineshapes.cc:2064
lineshape_doppler
void lineshape_doppler(Vector &ls_attenuation, Vector &, Vector &, Vector &, Vector &, Vector &, Vector &, const Numeric f0, const Numeric, const Numeric, const Numeric, const Numeric, const Numeric, const Numeric sigma, const Numeric, ConstVectorView f_grid, const bool, const bool)
Definition: lineshapes.cc:206
HZ2CM
const Numeric HZ2CM
Global constant, converts Hz to cm-1.
lineshape_CO2_lorentz
void lineshape_CO2_lorentz(Vector &ls_attenuation, Vector &, Vector &, Vector &, Vector &, Vector &, Vector &, const Numeric f0, const Numeric gamma, const Numeric, const Numeric, const Numeric, const Numeric, const Numeric, const Numeric, ConstVectorView f_grid, const bool, const bool)
Definition: lineshapes.cc:1352
a1
#define a1
Definition: complex.h:56
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
w_x_plus_iy_dH
void w_x_plus_iy_dH(Numeric &dx_dH, const Numeric &sigma, const Numeric &df0_dH)
Definition: lineshapes.cc:1616
absorption.h
Declarations required for the calculation of absorption coefficients.
matpackI.h
Implementation of Matrix, Vector, and such stuff.
dx
#define dx
Definition: legacy_continua.cc:22180
pow
Numeric pow(const Rational base, Numeric exp)
Power of.
Definition: rational.h:628
w_x_plus_iy_dgamma
void w_x_plus_iy_dgamma(Numeric &dy_dgamma, const Numeric &sigma)
Definition: lineshapes.cc:1606
lineshape_norm_linear_dF0
void lineshape_norm_linear_dF0(Vector &fac, const Numeric f0, ConstVectorView f_grid, const Numeric)
Definition: lineshapes.cc:2115
Array
This can be used to make arrays out of anything.
Definition: array.h:40
lineshape_lorentz
void lineshape_lorentz(Vector &ls_attenuation, Vector &ls_phase, Vector &, Vector &, Vector &, Vector &, Vector &, const Numeric f0, const Numeric gamma, const Numeric, const Numeric, const Numeric, const Numeric, const Numeric, const Numeric, ConstVectorView f_grid, const bool do_phase, const bool)
Definition: lineshapes.cc:97
PI
const Numeric PI
Global constant, pi.
lineshape_norm_no_norm_dT
void lineshape_norm_no_norm_dT(Vector &fac, const Numeric, ConstVectorView, const Numeric)
Definition: lineshapes.cc:2050
global_data::lineshape_data
Array< LineshapeRecord > lineshape_data
Definition: lineshapes.cc:2390
lineshape_norm_linear_dT
void lineshape_norm_linear_dT(Vector &fac, const Numeric, ConstVectorView, const Numeric)
Definition: lineshapes.cc:2094
lineshape_mirrored_lorentz
void lineshape_mirrored_lorentz(Vector &ls_attenuation, Vector &ls_phase, Vector &, Vector &, Vector &, Vector &, Vector &, const Numeric f0, const Numeric gamma, const Numeric, const Numeric, const Numeric, const Numeric, const Numeric, const Numeric, ConstVectorView f_grid, const bool do_phase, const bool)
Definition: lineshapes.cc:148
max
#define max(a, b)
Definition: legacy_continua.cc:20629
global_data
Definition: agenda_record.cc:33
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:476
lineshape_norm_VVH_dT
void lineshape_norm_VVH_dT(Vector &fac, const Numeric f0, ConstVectorView f_grid, const Numeric T)
Definition: lineshapes.cc:2252
define_lineshape_norm_data
void define_lineshape_norm_data()
Definition: lineshapes.cc:2540
a2
#define a2
Definition: complex.h:58
w_x_plus_iy_dT
void w_x_plus_iy_dT(Vector &dx_dT, Numeric &dy_dT, Numeric &dFu_dT, ConstVectorView f, const Numeric &f0, const Numeric &sigma, const Numeric &dPF_dT, const Numeric &dDF_dT, const Numeric &dsigma_dT, const Numeric &gamma, const Numeric &dgamma_dT)
Definition: lineshapes.cc:1541
lineshape_norm_quadratic_Rosenkranz_dF
void lineshape_norm_quadratic_Rosenkranz_dF(Vector &fac, const Numeric f0, ConstVectorView f_grid, const Numeric T)
Definition: lineshapes.cc:2181
lineshape_norm_quadratic_Rosenkranz_dT
void lineshape_norm_quadratic_Rosenkranz_dT(Vector &fac, const Numeric f0, ConstVectorView f_grid, const Numeric T)
Definition: lineshapes.cc:2162
lineshape_norm_VVW
void lineshape_norm_VVW(Vector &fac, const Numeric f0, ConstVectorView f_grid, const Numeric)
Definition: lineshapes.cc:2338
b1
#define b1
Definition: complex.h:57
hui_etal_1978_lineshape
void hui_etal_1978_lineshape(Vector &ls_attenuation, Vector &ls_phase, Vector &, Vector &, Vector &, Vector &, Vector &xvector, const Numeric f0, const Numeric gamma, const Numeric, const Numeric, const Numeric, const Numeric, const Numeric sigma, const Numeric, ConstVectorView f_grid, const bool do_phase, const bool)
Definition: lineshapes.cc:1682
bfun3_
long int bfun3_(Numeric y, Numeric x)
Definition: lineshapes.cc:1063
lineshape_norm_VVW_dF0
void lineshape_norm_VVW_dF0(Vector &fac, const Numeric f0, ConstVectorView f_grid, const Numeric)
Definition: lineshapes.cc:2370
SPEED_OF_LIGHT
const Numeric SPEED_OF_LIGHT
lineshape_norm_linear
void lineshape_norm_linear(Vector &fac, const Numeric f0, ConstVectorView f_grid, const Numeric)
Definition: lineshapes.cc:2080
chi_cousin
void chi_cousin(Numeric &chi, const Numeric &df)
Definition: lineshapes.cc:1308
lineshape_voigt_kuntz6
void lineshape_voigt_kuntz6(Vector &ls_attenuation, Vector &, Vector &ls_dattenuation_dfrequency_term, Vector &, Vector &ls_dattenuation_dpressure_term, Vector &, Vector &x, const Numeric f0, const Numeric gamma, const Numeric, const Numeric, const Numeric, const Numeric, const Numeric sigma, const Numeric, ConstVectorView f_grid, const bool, const bool do_partials)
Definition: lineshapes.cc:359
w_x_plus_iy_dF0
void w_x_plus_iy_dF0(Numeric &dx_dF0, const Numeric &sigma)
Definition: lineshapes.cc:1596
sqrt
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:620