ARTS  2.2.66
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 "arts.h"
36 #include "matpackI.h"
37 #include "array.h"
38 #include "absorption.h"
39 #include "Faddeeva.hh"
40 
60 void lineshape_no_shape( Vector& ls_attenuation _U_,
61  Vector& ls_phase _U_,
62  Vector& X _U_,
63  const Numeric f0 _U_,
64  const Numeric gamma _U_,
65  const Numeric sigma _U_,
66  ConstVectorView f_grid _U_)
67 {
68  // This function should never be called so throw an error here:
69  throw runtime_error("The no_shape lineshape is only a placeholder, but you tried\n"
70  "to use it like a real lineshape.");
71 }
72 
73 
86 void lineshape_lorentz(Vector& ls_attenuation,
87  Vector& ls_phase,
88  Vector& X _U_,
89  const Numeric f0,
90  const Numeric gamma,
91  const Numeric sigma _U_,
92  ConstVectorView f_grid)
93 {
94  const Index nf = f_grid.nelem();
95 
96  // PI:
97  extern const Numeric PI;
98 
99  // assert( ls.nelem() == nf );
100 
101  Numeric gamma2 = gamma * gamma;
102  Numeric fac = gamma/PI;
103 
104  for ( Index i=0; i<nf; ++i )
105  {
106  ls_attenuation[i] = fac / ( (f_grid[i]-f0) * (f_grid[i]-f0) + gamma2 );
107  ls_phase[i] = ( (f_grid[i]-f0) * (f_grid[i]-f0) + gamma2 ) * (f_grid[i]-f0) / PI ;
108  }
109 }
110 
123 void lineshape_doppler(Vector& ls_attenuation,
124  Vector& ls_phase _U_,
125  Vector& x _U_,
126  const Numeric f0,
127  const Numeric gamma _U_,
128  const Numeric sigma,
129  ConstVectorView f_grid)
130 {
131  const Index nf = f_grid.nelem();
132 
133  // SQRT(PI):
134  extern const Numeric PI;
135  const Numeric sqrtPI = sqrt(PI);
136 
137  // assert( ls_attenuation.nelem() == nf );
138 
139  Numeric sigma2 = sigma * sigma;
140  Numeric fac = 1.0 / (sqrtPI * sigma);
141 
142  for ( Index i=0; i<nf ; ++i )
143  {
144  ls_attenuation[i] = fac * exp( - pow( f_grid[i]-f0, 2) / sigma2 );
145  }
146 }
147 
148 
149 
150 //------------------------------------------------------------------------
151 
152 // help function for lineshape_voigt_kuntz1
153 long bfun6_(Numeric y, Numeric x)
154 {
155  /* System generated locals */
156  long int ret_val;
157 
158  /* Local variables */
159  Numeric s = 0;
160 
161  /* -------------------------------------------------------------------- */
162  s = fabs(x) + y;
163  if (s >= 15.f) {
164  ret_val = 1;
165  } else if (s >= 5.5f) {
166  ret_val = 2;
167  } else if (y >= fabs(x) * .195f - .176f) {
168  ret_val = 3;
169  } else {
170  ret_val = 4;
171  }
172  return ret_val;
173 } /* bfun6_ */
174 
175 
176 
236 void lineshape_voigt_kuntz6(Vector& ls_attenuation,
237  Vector& ls_phase _U_,
238  Vector& x,
239  const Numeric f0,
240  const Numeric gamma,
241  const Numeric sigma,
242  ConstVectorView f_grid)
243 
244 {
245  const Index nf = f_grid.nelem();
246 
247  // seems not necessary for Doppler correction
248  // extern const Numeric SQRT_NAT_LOG_2;
249 
250  // PI
251  extern const Numeric PI;
252 
253  // constant sqrt(1/pi)
254  const Numeric sqrt_invPI = sqrt(1/PI);
255 
256  // constant normalization factor for voigt
257  Numeric fac = 1.0 / sigma * sqrt_invPI;
258 
259 
260  /* Initialized data */
261 
262  Numeric yps1 = -1.0;
263  Numeric yps2 = -1.0;
264  Numeric yps3 = -1.0;
265  Numeric yps4 = -1.0;
266 
267  /* System generated locals */
268  long int i__1, i__2;
269  Numeric r__1;
270 
271  /* Local variables */
272  long int bmin = 0, lauf[16] = {0} /* was [4][4] */, bmax;
273  long int imin = 0, imax = 0, stack[80] = {0} /* was [20][4] */;
274  Numeric a1, a2, a3, a4, a5, a6, a8, b8, c8, d8, e8, f8, g8, h8, a7,
275  b7, c7, d7, e7, f7, o8, p8, q8, r8, s8, t8, g7, h7, o7, p7, q7,
276  r7, s7, t7, b6, c6, d6, e6, b5, c5, d5, e5, b4, c4, d4, b3, c3,
277  d3, b1, y2;
278  a1 = a2 = a3 = a4 = a5 = a6 = a8 = b8 = c8 = d8 = e8 = f8 = g8 = h8 = a7
279  = b7 = c7 = d7 = e7 = f7 = o8 = p8 = q8 = r8 = s8 = t8 = g7 = h7 = o7 = p7
280  = q7 = r7 = s7 = t7 = b6 = c6 = d6 = e6 = b5 = c5 = d5 = e5 = b4 = c4 = d4
281  = b3 = c3 = d3 = b1 = y2 = 0;
282  long int i2 = 0, i1 = 0;
283  Numeric x2 = 0, b2 = 0, c1 = 0;
284  long int stackp = 0, imitte = 0;
285  Numeric ym2 = 0;
286 
287 
288  // variables needed in original c routine:
289 
290  // Ratio of the Lorentz halfwidth to the Doppler halfwidth
291  Numeric y = gamma / sigma;
292 
293  // frequency in units of Doppler
294  for (i1=0; i1< (int) nf; i1++)
295  {
296  x[i1] = (f_grid[i1] - f0) / sigma;
297  }
298 
299  /* Parameter adjustments */
300  // this does not work for variables type Vector, corrected by
301  // adjusting the index to array ls and x
302  //--ls;
303  //--x;
304 
305  /* Function Body */
306  y2 = y * y;
307  if (y >= 15.0 || x[0] >= 15.0 || x[nf-1] <= -15.0) {
308  lauf[0] = 1;
309  lauf[4] = nf;
310  lauf[8] = nf;
311  lauf[12] = 0;
312  goto L7;
313  }
314  for (i2 = 1; i2 <= 4; ++i2) {
315  for (i1 = 1; i1 <= 4; ++i1) {
316  lauf[i1 + (i2 << 2) - 5] = i2 % 2 * (nf + 1);
317  /* L1: */
318  }
319  }
320  stackp = 1;
321  stack[stackp - 1] = 1;
322  stack[stackp + 19] = nf;
323  stack[stackp + 39] = bfun6_(y, x[0]);
324  stack[stackp + 59] = bfun6_(y, x[nf-1]);
325  L2:
326  imin = stack[stackp - 1];
327  imax = stack[stackp + 19];
328  bmin = stack[stackp + 39];
329  bmax = stack[stackp + 59];
330  if (bmin == bmax) {
331  if (x[imax-1] < 0.f) {
332  /* Computing MIN */
333  i__1 = imin, i__2 = lauf[bmin - 1];
334  lauf[bmin - 1] = min(i__1,i__2);
335  /* Computing MAX */
336  i__1 = imax, i__2 = lauf[bmax + 3];
337  lauf[bmax + 3] = max(i__1,i__2);
338  --stackp;
339  goto L3;
340  } else if (x[imin-1] >= 0.f) {
341  /* Computing MIN */
342  i__1 = imin, i__2 = lauf[bmin + 7];
343  lauf[bmin + 7] = min(i__1,i__2);
344  /* Computing MAX */
345  i__1 = imax, i__2 = lauf[bmax + 11];
346  lauf[bmax + 11] = max(i__1,i__2);
347  --stackp;
348  goto L3;
349  }
350  }
351  imitte = (imax + imin) / 2;
352  stack[stackp - 1] = imitte + 1;
353  stack[stackp + 19] = imax;
354  stack[stackp + 39] = bfun6_(y, x[imitte]);
355  stack[stackp + 59] = bmax;
356  ++stackp;
357  stack[stackp - 1] = imin;
358  stack[stackp + 19] = imitte;
359  stack[stackp + 39] = bmin;
360  stack[stackp + 59] = bfun6_(y, x[imitte-1]);
361  L3:
362  if (stackp > 0) {
363  goto L2;
364  }
365  /* ---- Region 4 */
366  /* -------------------------------------------------------------------- */
367  if (lauf[7] >= lauf[3] || lauf[15] >= lauf[11]) {
368  if ((r__1 = y - yps4, fabs(r__1)) > 1e-8f) {
369  //yps4 = y;
370  a7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
371  y2 * (y2 * (y2 * (2.35944f - y2 * .56419f) - 72.9359f) +
372  571.687f) - 5860.68f) + 40649.2f) - 320772.f) + 1684100.f)
373  - 9694630.f) + 40816800.f) - 1.53575e8f) + 4.56662e8f) -
374  9.86604e8f) + 1.16028e9f);
375  b7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
376  y2 * (y2 * (23.0312f - y2 * 7.33447f) - 234.143f) -
377  2269.19f) + 45251.3f) - 234417.f) + 3599150.f) -
378  7723590.f) + 86482900.f) - 2.91876e8f) + 8.06985e8f) -
379  9.85386e8f) - 5.60505e8f);
380  c7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
381  y2 * (97.6203f - y2 * 44.0068f) + 1097.77f) - 25338.3f) +
382  98079.1f) + 576054.f) - 2.3818e7f) + 22930200.f) -
383  2.04467e8f) + 2.94262e8f) + 2.47157e8f) - 6.51523e8f);
384  d7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
385  228.563f - y2 * 161.358f) + 8381.97f) - 66431.2f) -
386  303569.f) + 2240400.f) + 38311200.f) - 41501300.f) -
387  99622400.f) + 2.70167e8f) - 2.63894e8f);
388  e7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
389  296.38f - y2 * 403.396f) + 23507.6f) - 66212.1f) -
390  1.003e6f) + 468142.f) + 24620100.f) + 5569650.f) +
391  1.40677e8f) - 63177100.f);
392  f7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (125.591f -
393  y2 * 726.113f) + 37544.8f) + 8820.94f) - 934717.f) -
394  1931140.f) - 33289600.f) + 4073820.f) - 16984600.f);
395  g7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (-260.198f - y2 *
396  968.15f) + 37371.9f) + 79902.5f) - 186682.f) - 900010.f)
397  + 7528830.f) - 1231650.f);
398  h7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (-571.645f - y2 * 968.15f)
399  + 23137.1f) + 72520.9f) + 153468.f) + 86407.6f) -
400  610622.f);
401  o7 = y * (y2 * (y2 * (y2 * (y2 * (-575.164f - y2 * 726.113f) +
402  8073.15f) + 26538.5f) + 49883.8f) - 23586.5f);
403  p7 = y * (y2 * (y2 * (y2 * (-352.467f - y2 * 403.396f) +
404  953.655f) + 2198.86f) - 8009.1f);
405  q7 = y * (y2 * (y2 * (-134.792f - y2 * 161.358f) - 271.202f) -
406  622.056f);
407  r7 = y * (y2 * (-29.7896f - y2 * 44.0068f) - 77.0535f);
408  s7 = y * (-2.92264f - y2 * 7.33447f);
409  t7 = y * -.56419f;
410  a8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
411  y2 * (y2 * (y2 * (y2 - 3.68288f) + 126.532f) - 955.194f)
412  + 9504.65f) - 70946.1f) + 483737.f) - 2857210.f) +
413  14464700.f) - 61114800.f) + 2.11107e8f) - 5.79099e8f) +
414  1.17022e9f) - 1.5599e9f) + 1.02827e9f;
415  b8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
416  y2 * (y2 * (y2 * 14.f - 40.5117f) + 533.254f) + 3058.26f)
417  - 55600.f) + 498334.f) - 2849540.f) + 13946500.f) -
418  70135800.f) + 2.89676e8f) - 7.53828e8f) + 1.66421e9f) -
419  2.28855e9f) + 1.5599e9f;
420  c8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
421  y2 * (y2 * 91 - 198.876f) - 1500.17f) + 48153.3f) -
422  217801.f) - 1063520.f) + 1.4841e7f) - 46039600.f) +
423  63349600.f) - 6.60078e8f) + 1.06002e9f) - 1.66421e9f) +
424  1.17022e9f;
425  d8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
426  y2 * 364 - 567.164f) - 16493.7f) + 161461.f) + 280428.f)
427  - 6890020.f) - 6876560.f) + 1.99846e8f) + 54036700.f) +
428  6.60078e8f) - 7.53828e8f) + 5.79099e8f;
429  e8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 *
430  1001 - 1012.79f) - 55582.f) + 240373.f) + 1954700.f) -
431  5257220.f) - 50101700.f) - 1.99846e8f) + 63349600.f) -
432  2.89676e8f) + 2.11107e8f;
433  f8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 2002 -
434  1093.82f) - 106663.f) + 123052.f) + 3043160.f) +
435  5257220.f) - 6876560.f) + 46039600.f) - 70135800.f) +
436  61114800.f;
437  g8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 3003 -
438  486.14f) - 131337.f) - 123052.f) + 1954700.f) + 6890020.f)
439  + 1.4841e7f) - 13946500.f) + 14464700.f;
440  h8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 3432 + 486.14f) -
441  106663.f) - 240373.f) + 280428.f) + 1063520.f) -
442  2849540.f) + 2857210.f;
443  o8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 3003 + 1093.82f) -
444  55582.f) - 161461.f) - 217801.f) - 498334.f) + 483737.f;
445  p8 = y2 * (y2 * (y2 * (y2 * (y2 * 2002 + 1012.79f) - 16493.7f) -
446  48153.3f) - 55600.f) + 70946.1f;
447  q8 = y2 * (y2 * (y2 * (y2 * 1001.f + 567.164f) - 1500.17f) -
448  3058.26f) + 9504.65f;
449  r8 = y2 * (y2 * (y2 * 364 + 198.876f) + 533.254f) + 955.194f;
450  s8 = y2 * (y2 * 91.f + 40.5117f) + 126.532f;
451  t8 = y2 * 14.f + 3.68288f;
452  }
453  ym2 = y * 2;
454  for (i2 = 1; i2 <= 3; i2 += 2) {
455  i__1 = lauf[((i2 + 1) << 2) - 1];
456  for (i1 = lauf[(i2 << 2) - 1]; i1 <= i__1; ++i1) {
457  x2 = x[i1-1] * x[i1-1];
458  ls_attenuation[i1-1] = fac * (exp(y2 - x2) * cos(x[i1-1] * ym2) - (a7 + x2 *
459  (b7 + x2 * (c7 + x2 * (d7 + x2 * (e7 + x2 * (f7 + x2
460  * (g7 + x2 * (h7 + x2 * (o7 + x2 * (p7 + x2 * (q7 +
461  x2 * (r7 + x2 * (s7 + x2 * t7))))))))))))) / (a8 + x2
462  * (b8 + x2 * (c8 + x2 * (d8 + x2 * (e8 + x2 * (f8 +
463  x2 * (g8 + x2 * (h8 + x2 * (o8 + x2 * (p8 + x2 * (q8
464  + x2 * (r8 + x2 * (s8 + x2 * (t8 + x2)))))))))))))));
465  /* L4: */
466  }
467  }
468  }
469  /* ---- Region 3 */
470  /* -------------------------------------------------------------------- */
471  if (lauf[6] >= lauf[2] || lauf[14] >= lauf[10]) {
472  if ((r__1 = y - yps3, fabs(r__1)) > 1e-8f) {
473  //yps3 = y;
474  a5 = y * (y * (y * (y * (y * (y * (y * (y * (y *
475  .564224f + 7.55895f) + 49.5213f) + 204.501f) + 581.746f)
476  + 1174.8f) + 1678.33f) + 1629.76f) + 973.778f) + 272.102f;
477  b5 = y * (y * (y * (y * (y * (y * (y * 2.25689f + 22.6778f)
478  + 100.705f) + 247.198f) + 336.364f) + 220.843f) -
479  2.34403f) - 60.5644f;
480  c5 = y * (y * (y * (y * (y * 3.38534f + 22.6798f) + 52.8454f)
481  + 42.5683f) + 18.546f) + 4.58029f;
482  d5 = y * (y * (y * 2.25689f + 7.56186f) + 1.66203f) - .128922f;
483  e5 = y * .564224f + 9.71457e-4f;
484  a6 = y * (y * (y * (y * (y * (y * (y * (y * (y * (y +
485  13.3988f) + 88.2674f) + 369.199f) + 1074.41f) + 2256.98f)
486  + 3447.63f) + 3764.97f) + 2802.87f) + 1280.83f) +
487  272.102f;
488  b6 = y * (y * (y * (y * (y * (y * (y * (y * 5.f +
489  53.5952f) + 266.299f) + 793.427f) + 1549.68f) + 2037.31f)
490  + 1758.34f) + 902.306f) + 211.678f;
491  c6 = y * (y * (y * (y * (y * (y * 10.f + 80.3928f) +
492  269.292f) + 479.258f) + 497.302f) + 308.186f) + 78.866f;
493  d6 = y * (y * (y * (y * 10.f + 53.5952f) + 92.7568f) +
494  55.0293f) + 22.0353f;
495  e6 = y * (y * 5.f + 13.3988f) + 1.49645f;
496  }
497  for (i2 = 1; i2 <= 3; i2 += 2) {
498  i__1 = lauf[((i2 + 1) << 2) - 2];
499  for (i1 = lauf[(i2 << 2) - 2]; i1 <= i__1; ++i1) {
500  x2 = x[i1-1] * x[i1-1];
501  ls_attenuation[i1-1] = fac * (a5 + x2 * (b5 + x2 * (c5 + x2 * (d5 + x2 *
502  e5)))) / (a6 + x2 * (b6 + x2 * (c6 + x2 * (d6 + x2 * (
503  e6 + x2)))));
504  /* L5: */
505  }
506  }
507  }
508  /* ---- Region 2 */
509  /* -------------------------------------------------------------------- */
510  if (lauf[5] >= lauf[1] || lauf[13] >= lauf[9]) {
511  if ((r__1 = y - yps2, fabs(r__1)) > 1e-8f) {
512  //yps2 = y;
513  a3 = y * (y2 * (y2 * (y2 * .56419f + 3.10304f) + 4.65456f) +
514  1.05786f);
515  b3 = y * (y2 * (y2 * 1.69257f + .56419f) + 2.962f);
516  c3 = y * (y2 * 1.69257f - 2.53885f);
517  d3 = y * .56419f;
518  a4 = y2 * (y2 * (y2 * (y2 + 6.f) + 10.5f) + 4.5f) + .5625f;
519  b4 = y2 * (y2 * (y2 * 4.f + 6.f) + 9.f) - 4.5f;
520  c4 = y2 * (y2 * 6.f - 6.f) + 10.5f;
521  d4 = y2 * 4.f - 6.f;
522  }
523  for (i2 = 1; i2 <= 3; i2 += 2) {
524  i__1 = lauf[((i2 + 1) << 2) - 3];
525  for (i1 = lauf[(i2 << 2) - 3]; i1 <= i__1; ++i1) {
526  x2 = x[i1-1] * x[i1-1];
527  ls_attenuation[i1-1] = fac * (a3 + x2 * (b3 + x2 * (c3 + x2 * d3))) / (a4
528  + x2 * (b4 + x2 * (c4 + x2 * (d4 + x2))));
529  /* L6: */
530  }
531  }
532  }
533  /* ---- Region 1 */
534  /* -------------------------------------------------------------------- */
535  L7:
536  if (lauf[4] >= lauf[0] || lauf[12] >= lauf[8]) {
537  if ((r__1 = y - yps1, fabs(r__1)) > 1e-8f) {
538  //yps1 = y;
539  a1 = y * .5641896f;
540  b1 = y2 + .5f;
541  a2 = y2 * 4;
542  }
543 
544  c1 = fac * a1;
545  for (i2 = 1; i2 <= 3; i2 += 2) {
546  i__1 = lauf[((i2 + 1) << 2) - 4];
547  for (i1 = lauf[(i2 << 2) - 4]; i1 <= i__1; ++i1) {
548  x2 = x[i1-1] * x[i1-1];
549  b2 = b1 - x2;
550  ls_attenuation[i1-1] = c1 * (b1 + x2) / (b2 * b2 + a2 * x2);
551  /* L8: */
552  }
553  }
554  }
555 }
556 
557 
558 //---------------------------------------------------------------
559 // help function for lineshape_voigt_kuntz3
560 
561 long int bfun3_(Numeric y, Numeric x)
562 {
563  /* System generated locals */
564  long int ret_val;
565 
566  /* Local variables */
567  Numeric x2 = 0, y2 = 0;
568 
569 /* -------------------------------------------------------------------- */
570  x2 = x * x;
571  y2 = y * y;
572  if (x2 * .4081676f + y2 > 21.159543f) {
573  if (x2 * .7019639f + y2 > 1123.14221f) {
574  ret_val = 0;
575  } else {
576  ret_val = 1;
577  }
578  } else {
579  if (x2 * .20753051f + y2 > 4.20249292f) {
580  ret_val = 2;
581  } else if (y >= fabs(x) * .08f - .12f) {
582  ret_val = 3;
583  } else {
584  ret_val = 4;
585  }
586  }
587  return ret_val;
588 } /* bfun3_ */
589 
590 
591 
650 void lineshape_voigt_kuntz3(Vector& ls_attenuation,
651  Vector& ls_phase _U_,
652  Vector& x,
653  const Numeric f0,
654  const Numeric gamma,
655  const Numeric sigma,
656  ConstVectorView f_grid)
657 
658 {
659  const Index nf = f_grid.nelem();
660 
661  // seems not necessary for Doppler correction
662  // extern const Numeric SQRT_NAT_LOG_2;
663 
664  // PI
665  extern const Numeric PI;
666 
667  // constant sqrt(1/pi)
668  const Numeric sqrt_invPI = sqrt(1/PI);
669 
670  // constant normalization factor for voigt
671  Numeric fac = 1.0 / sigma * sqrt_invPI;
672 
673  /* Initialized data */
674 
675  Numeric yps0 = -1.0;
676  Numeric yps1 = -1.0;
677  Numeric yps2 = -1.0;
678  Numeric yps3 = -1.0;
679  Numeric yps4 = -1.0;
680 
681  /* System generated locals */
682  long i__1, i__2;
683  Numeric r__1;
684 
685  /* Local variables */
686  long bmin = 0, lauf[20] = {0} /* was [5][4] */, bmax, imin, imax;
687  long stack[80] = {0} /* was [20][4] */;
688  Numeric a0, a1, a2, a3, a4, a5, a6, a7, a8, b8, c8, d8, e8, f8, g8,
689  h8, b7, c7, d7, e7, f7, g7, o8, p8, q8, r8, s8, t8, h7, o7, p7,
690  q7, r7, s7, t7, b6, c6, d6, e6, b5, c5, d5, e5, b4, c4, d4, b3,
691  c3, d3, b1, y2;
692  a0 = a1 = a2 = a3 = a4 = a5 = a6 = a7 = a8 = b8 = c8 = d8 = e8 = f8 = g8
693  = h8 = b7 = c7 = d7 = e7 = f7 = g7 = o8 = p8 = q8 = r8 = s8 = t8 = h7 = o7
694  = p7 = q7 = r7 = s7 = t7 = b6 = c6 = d6 = e6 = b5 = c5 = d5 = e5 = b4 = c4
695  = d4 = b3 = c3 = d3 = b1 = y2 = 0;
696  long i2 = 0, i1 = 0;
697  Numeric x2 = 0, b2 = 0, b0 = 0, c1 = 0;
698  long stackp = 0, imitte = 0;
699  Numeric ym2 = 0;
700 
701  // variables needed in original c routine:
702 
703  // Ratio of the Lorentz halfwidth to the Doppler halfwidth
704  Numeric y = gamma / sigma;
705 
706  // frequency in units of Doppler
707  for (i1=0; i1< (int) nf; i1++)
708  {
709  x[i1] = (f_grid[i1] - f0) / sigma;
710  }
711 
712 
713  /* Parameter adjustments */
714  // this does not work for variables type Vector, corrected by
715  // adjusting the index to array ls and x
716  //--ls;
717  //--x;
718 
719  /* Function Body */
720  y2 = y * y;
721  if (y >= 23.0 || x[0] >= 39.0 || x[nf-1] <= -39.0) {
722  lauf[0] = 1;
723  lauf[5] = nf;
724  lauf[10] = nf;
725  lauf[15] = 0;
726  goto L8;
727  }
728  for (i2 = 1; i2 <= 4; ++i2) {
729  for (i1 = 0; i1 <= 4; ++i1) {
730  lauf[i1 + i2 * 5 - 5] = i2 % 2 * (nf + 1);
731  /* L1: */
732  }
733  }
734  stackp = 1;
735  stack[stackp - 1] = 1;
736  stack[stackp + 19] = nf;
737  stack[stackp + 39] = bfun3_(y, x[0]);
738  stack[stackp + 59] = bfun3_(y, x[nf-1]);
739  L2:
740  imin = stack[stackp - 1];
741  imax = stack[stackp + 19];
742  bmin = stack[stackp + 39];
743  bmax = stack[stackp + 59];
744  if (bmin == bmax) {
745  if (x[imax-1] < 0.f) {
746  /* Computing MIN */
747  i__1 = imin, i__2 = lauf[bmin];
748  lauf[bmin] = min(i__1,i__2);
749  /* Computing MAX */
750  i__1 = imax, i__2 = lauf[bmax + 5];
751  lauf[bmax + 5] = max(i__1,i__2);
752  --stackp;
753  goto L3;
754  } else if (x[imin-1] >= 0.f) {
755  /* Computing MIN */
756  i__1 = imin, i__2 = lauf[bmin + 10];
757  lauf[bmin + 10] = min(i__1,i__2);
758  /* Computing MAX */
759  i__1 = imax, i__2 = lauf[bmax + 15];
760  lauf[bmax + 15] = max(i__1,i__2);
761  --stackp;
762  goto L3;
763  }
764  }
765  imitte = (imax + imin) / 2;
766  stack[stackp - 1] = imitte + 1;
767  stack[stackp + 19] = imax;
768  stack[stackp + 39] = bfun3_(y, x[imitte]);
769  stack[stackp + 59] = bmax;
770  ++stackp;
771  stack[stackp - 1] = imin;
772  stack[stackp + 19] = imitte;
773  stack[stackp + 39] = bmin;
774  stack[stackp + 59] = bfun3_(y, x[imitte-1]);
775  L3:
776  if (stackp > 0) {
777  goto L2;
778  }
779  /* ---- Region 4 */
780  /* -------------------------------------------------------------------- */
781  if (lauf[9] >= lauf[4] || lauf[19] >= lauf[14]) {
782  if ((r__1 = y - yps4, fabs(r__1)) > 1e-8f) {
783  //yps4 = y;
784  a7 = y * (y2 * (y2 * 4.56662e8f - 9.86604e8f) + 1.16028e9f);
785  b7 = y * (y2 * (y2 * 8.06985e8f - 9.85386e8f) - 5.60505e8f);
786  c7 = y * (y2 * (y2 * 2.94262e8f + 2.47157e8f) - 6.51523e8f);
787  d7 = y * (y2 * (2.70167e8f - y2 * 99622400.f) - 2.63894e8f);
788  e7 = y * (y2 * (y2 * 5569650.f + 1.40677e8f) - 63177100.f);
789  f7 = y * (y2 * (4073820.f - y2 * 33289600.f) - 16984600.f);
790  g7 = y * (y2 * (7528830.f - y2 * 900010) - 1231650.f);
791  h7 = y * (y2 * (y2 * 153468 + 86407.6f) - 610622.f);
792  o7 = y * (y2 * (y2 * 26538.5f + 49883.8f) - 23586.5f);
793  p7 = y * (y2 * (y2 * 953.655f + 2198.86f) - 8009.1f);
794  q7 = y * (y2 * (-271.202f - y2 * 134.792f) - 622.056f);
795  r7 = y * (y2 * (-29.7896f - y2 * 44.0068f) - 77.0535f);
796  s7 = y * (-2.92264f - y2 * 7.33447f);
797  t7 = y * -.56419f;
798  a8 = y2 * (y2 * 1.17022e9f - 1.5599e9f) + 1.02827e9f;
799  b8 = y2 * (y2 * 1.66421e9f - 2.28855e9f) + 1.5599e9f;
800  c8 = y2 * (y2 * 1.06002e9f - 1.66421e9f) + 1.17022e9f;
801  d8 = y2 * (y2 * 6.60078e8f - 7.53828e8f) + 5.79099e8f;
802  e8 = y2 * (y2 * 63349600.f - 2.89676e8f) + 2.11107e8f;
803  f8 = y2 * (y2 * 46039600.f - 70135800.f) + 61114800.f;
804  g8 = y2 * (y2 * 1.4841e7f - 13946500.f) + 14464700.f;
805  h8 = y2 * (y2 * 1063520.f - 2849540.f) + 2857210.f;
806  o8 = y2 * (-498334.f - y2 * 217801.f) + 483737.f;
807  p8 = y2 * (-55600.f - y2 * 48153.3f) + 70946.1f;
808  q8 = y2 * (-3058.26f - y2 * 1500.17f) + 9504.65f;
809  r8 = y2 * (y2 * 198.876f + 533.254f) + 955.194f;
810  s8 = y2 * (y2 * 91.f + 40.5117f) + 126.532f;
811  t8 = y2 * 14.f + 3.68288f;
812  }
813  ym2 = y * 2;
814  for (i2 = 1; i2 <= 3; i2 += 2) {
815  i__1 = lauf[(i2 + 1) * 5 - 1];
816  for (i1 = lauf[i2 * 5 - 1]; i1 <= i__1; ++i1) {
817  x2 = x[i1-1] * x[i1-1];
818  ls_attenuation[i1-1] = fac * (exp(y2 - x2) * cos(x[i1-1] * ym2) - (a7 + x2 *
819  (b7 + x2 * (c7 + x2 * (d7 + x2 * (e7 + x2 * (f7 + x2
820  * (g7 + x2 * (h7 + x2 * (o7 + x2 * (p7 + x2 * (q7 +
821  x2 * (r7 + x2 * (s7 + x2 * t7))))))))))))) / (a8 + x2
822  * (b8 + x2 * (c8 + x2 * (d8 + x2 * (e8 + x2 * (f8 +
823  x2 * (g8 + x2 * (h8 + x2 * (o8 + x2 * (p8 + x2 * (q8
824  + x2 * (r8 + x2 * (s8 + x2 * (t8 + x2)))))))))))))));
825  /* L4: */
826  }
827  }
828  }
829  /* ---- Region 3 */
830  /* -------------------------------------------------------------------- */
831  if (lauf[8] >= lauf[3] || lauf[18] >= lauf[13]) {
832  if ((r__1 = y - yps3, fabs(r__1)) > 1e-8f) {
833  //yps3 = y;
834  a5 = y * (y * (y * (y * (y * (y * (y * (y * (y *
835  .564224f + 7.55895f) + 49.5213f) + 204.501f) + 581.746f)
836  + 1174.8f) + 1678.33f) + 1629.76f) + 973.778f) + 272.102f;
837  b5 = y * (y * (y * (y * (y * (y * (y * 2.25689f + 22.6778f)
838  + 100.705f) + 247.198f) + 336.364f) + 220.843f) -
839  2.34403f) - 60.5644f;
840  c5 = y * (y * (y * (y * (y * 3.38534f + 22.6798f) + 52.8454f)
841  + 42.5683f) + 18.546f) + 4.58029f;
842  d5 = y * (y * (y * 2.25689f + 7.56186f) + 1.66203f) - .128922f;
843  e5 = y * .564224f + 9.71457e-4f;
844  a6 = y * (y * (y * (y * (y * (y * (y * (y * (y * (y +
845  13.3988f) + 88.2674f) + 369.199f) + 1074.41f) + 2256.98f)
846  + 3447.63f) + 3764.97f) + 2802.87f) + 1280.83f) +
847  272.102f;
848  b6 = y * (y * (y * (y * (y * (y * (y * (y * 5.f +
849  53.5952f) + 266.299f) + 793.427f) + 1549.68f) + 2037.31f)
850  + 1758.34f) + 902.306f) + 211.678f;
851  c6 = y * (y * (y * (y * (y * (y * 10.f + 80.3928f) +
852  269.292f) + 479.258f) + 497.302f) + 308.186f) + 78.866f;
853  d6 = y * (y * (y * (y * 10.f + 53.5952f) + 92.7568f) +
854  55.0293f) + 22.0353f;
855  e6 = y * (y * 5.f + 13.3988f) + 1.49645f;
856  }
857  for (i2 = 1; i2 <= 3; i2 += 2) {
858  i__1 = lauf[(i2 + 1) * 5 - 2];
859  for (i1 = lauf[i2 * 5 - 2]; i1 <= i__1; ++i1) {
860  x2 = x[i1-1] * x[i1-1];
861  ls_attenuation[i1-1] = fac * (a5 + x2 * (b5 + x2 * (c5 + x2 * (d5 + x2 *
862  e5)))) / (a6 + x2 * (b6 + x2 * (c6 + x2 * (d6 + x2 * (
863  e6 + x2)))));
864  /* L5: */
865  }
866  }
867  }
868  /* ---- Region 2 */
869  /* -------------------------------------------------------------------- */
870  if (lauf[7] >= lauf[2] || lauf[17] >= lauf[12]) {
871  if ((r__1 = y - yps2, fabs(r__1)) > 1e-8f) {
872  //yps2 = y;
873  a3 = y * (y2 * (y2 * (y2 * .56419f + 3.10304f) + 4.65456f) +
874  1.05786f);
875  b3 = y * (y2 * (y2 * 1.69257f + .56419f) + 2.962f);
876  c3 = y * (y2 * 1.69257f - 2.53885f);
877  d3 = y * .56419f;
878  a4 = y2 * (y2 * (y2 * (y2 + 6.f) + 10.5f) + 4.5f) + .5625f;
879  b4 = y2 * (y2 * (y2 * 4.f + 6.f) + 9.f) - 4.5f;
880  c4 = y2 * (y2 * 6.f - 6.f) + 10.5f;
881  d4 = y2 * 4.f - 6.f;
882  }
883  for (i2 = 1; i2 <= 3; i2 += 2) {
884  i__1 = lauf[(i2 + 1) * 5 - 3];
885  for (i1 = lauf[i2 * 5 - 3]; i1 <= i__1; ++i1) {
886  x2 = x[i1-1] * x[i1-1];
887  ls_attenuation[i1-1] = fac * (a3 + x2 * (b3 + x2 * (c3 + x2 * d3))) / (a4
888  + x2 * (b4 + x2 * (c4 + x2 * (d4 + x2))));
889  /* L6: */
890  }
891  }
892  }
893  /* ---- Region 1 */
894  /* -------------------------------------------------------------------- */
895  if (lauf[6] >= lauf[1] || lauf[16] >= lauf[11]) {
896  if ((r__1 = y - yps1, fabs(r__1)) > 1e-8f) {
897  //yps1 = y;
898  a1 = y * .5641896f;
899  b1 = y2 + .5f;
900  a2 = y2 * 4;
901  }
902 
903  c1 = a1 * fac;
904  for (i2 = 1; i2 <= 3; i2 += 2) {
905  i__1 = lauf[(i2 + 1) * 5 - 4];
906  for (i1 = lauf[i2 * 5 - 4]; i1 <= i__1; ++i1) {
907  x2 = x[i1-1] * x[i1-1];
908  b2 = b1 - x2;
909  ls_attenuation[i1-1] = c1 * (b1 + x2) / (b2 * b2 + a2 * x2);
910  /* L7: */
911  }
912  }
913  }
914  /* ---- Region 0 (Lorentz) */
915  /* -------------------------------------------------------------------- */
916 L8:
917  if (lauf[5] >= lauf[0] || lauf[15] >= lauf[10]) {
918  if ((r__1 = y - yps0, fabs(r__1)) > 1e-8f) {
919  //yps0 = y;
920  a0 = y * .5641896f;
921  }
922 
923  b0 = a0 * fac;
924  for (i2 = 1; i2 <= 3; i2 += 2) {
925  i__1 = lauf[(i2 + 1) * 5 - 5];
926  for (i1 = lauf[i2 * 5 - 5]; i1 <= i__1; ++i1) {
927  ls_attenuation[i1-1] = b0 / (x[i1-1] * x[i1-1] + y2);
928  /* L9: */
929  }
930  }
931  }
932 }
933 
934 
935 //------------------------------------------------------------------
936 // help function for lineshape_voigt_kuntz4
937 
938 long bfun4_(Numeric y, Numeric x)
939 {
940  /* System generated locals */
941  long ret_val;
942 
943  /* Local variables */
944  Numeric x2 = 0, y2 = 0;
945 
946  x2 = x * x;
947  y2 = y * y;
948  if (x2 * .0062f + y2 * .01417f > 1.f) {
949  if (x2 * 6.2e-5f + y2 * 1.98373e-4f > 1.f) {
950  ret_val = 0;
951  } else {
952  ret_val = 1;
953  }
954  } else {
955  if (x2 * .041649f + y2 * .111111111f > 1.f) {
956  ret_val = 2;
957  } else if (y >= fabs(x) * .19487f - .1753846f) {
958  ret_val = 3;
959  } else {
960  ret_val = 4;
961  }
962  }
963  return ret_val;
964 }
965 
1024 void lineshape_voigt_kuntz4(Vector& ls_attenuation,
1025  Vector& ls_phase _U_,
1026  Vector& x,
1027  const Numeric f0,
1028  const Numeric gamma,
1029  const Numeric sigma,
1030  ConstVectorView f_grid)
1031 {
1032  const Index nf = f_grid.nelem();
1033 
1034  // seems not necessary for Doppler correction
1035  // extern const Numeric SQRT_NAT_LOG_2;
1036 
1037  // PI
1038  extern const Numeric PI;
1039 
1040  // constant sqrt(1/pi)
1041  const Numeric sqrt_invPI = sqrt(1/PI);
1042 
1043  // constant normalization factor for voigt
1044  Numeric fac = 1.0 / sigma * sqrt_invPI;
1045 
1046 
1047  /* Initialized data */
1048 
1049  float yps0 = -1.f;
1050  float yps1 = -1.f;
1051  float yps2 = -1.f;
1052  float yps3 = -1.f;
1053  float yps4 = -1.f;
1054 
1055  /* System generated locals */
1056  long i__1, i__2;
1057  float r__1;
1058 
1059  /* Local variables */
1060  long bmin = 0, lauf[20] = {0} /* was [5][4] */, bmax, imin, imax;
1061  long stack[80] = {0} /* was [20][4] */;
1062  Numeric a0, a1, a2, a3, a4, a5, a6, a7, a8, b8, c8, d8, e8, f8, g8,
1063  h8, b7, c7, d7, e7, f7, g7, o8, p8, q8, r8, s8, t8, h7, o7, p7,
1064  q7, r7, s7, t7, b6, c6, d6, e6, b5, c5, d5, e5, b4, c4, d4, b3,
1065  c3, d3, b1, y2;
1066  a0 = a1 = a2 = a3 = a4 = a5 = a6 = a7 = a8 = b8 = c8 = d8 = e8 = f8 = g8
1067  = h8 = b7 = c7 = d7 = e7 = f7 = g7 = o8 = p8 = q8 = r8 = s8 = t8 = h7
1068  = o7 = p7 = q7 = r7 = s7 = t7 = b6 = c6 = d6 = e6 = b5 = c5 = d5 = e5
1069  = b4 = c4 = d4 = b3 = c3 = d3 = b1 = y2 = 0;
1070  long i2 = 0, i1 = 0;
1071  Numeric x2 = 0, b2 = 0, b0 = 0, c1 = 0;
1072  long stackp = 0, imitte = 0;
1073  Numeric ym2 = 0;
1074 
1075  // variables needed in original c routine:
1076 
1077  // Ratio of the Lorentz halfwidth to the Doppler halfwidth
1078  Numeric y = gamma / sigma;
1079 
1080  // frequency in units of Doppler
1081  for (i1=0; i1< (int) nf; i1++)
1082  {
1083  x[i1] = (f_grid[i1] - f0) / sigma;
1084  }
1085 
1086 
1087  /* Parameter adjustments */
1088  // this does not work for variables type Vector, corrected by
1089  // adjusting the index to array ls and x
1090  //--ls;
1091  //--x;
1092 
1093  /* Function Body */
1094  y2 = y * y;
1095  if (y >= 71.f || x[0] >= 123.f || x[nf-1] <= -123.f) {
1096  lauf[0] = 1;
1097  lauf[5] = nf;
1098  lauf[10] = nf;
1099  lauf[15] = 0;
1100  goto L8;
1101  }
1102  for (i2 = 1; i2 <= 4; ++i2) {
1103  for (i1 = 0; i1 <= 4; ++i1) {
1104  lauf[i1 + i2 * 5 - 5] = i2 % 2 * (nf + 1);
1105  /* L1: */
1106  }
1107  }
1108  stackp = 1;
1109  stack[stackp - 1] = 1;
1110  stack[stackp + 19] = nf;
1111  stack[stackp + 39] = bfun4_(y, x[0]);
1112  stack[stackp + 59] = bfun4_(y, x[nf-1]);
1113  L2:
1114  imin = stack[stackp - 1];
1115  imax = stack[stackp + 19];
1116  bmin = stack[stackp + 39];
1117  bmax = stack[stackp + 59];
1118  if (bmin == bmax) {
1119  if (x[imax-1] < 0.f) {
1120  /* Computing MIN */
1121  i__1 = imin, i__2 = lauf[bmin];
1122  lauf[bmin] = min(i__1,i__2);
1123  /* Computing MAX */
1124  i__1 = imax, i__2 = lauf[bmax + 5];
1125  lauf[bmax + 5] = max(i__1,i__2);
1126  --stackp;
1127  goto L3;
1128  } else if (x[imin-1] >= 0.f) {
1129  /* Computing MIN */
1130  i__1 = imin, i__2 = lauf[bmin + 10];
1131  lauf[bmin + 10] = min(i__1,i__2);
1132  /* Computing MAX */
1133  i__1 = imax, i__2 = lauf[bmax + 15];
1134  lauf[bmax + 15] = max(i__1,i__2);
1135  --stackp;
1136  goto L3;
1137  }
1138  }
1139  imitte = (imax + imin) / 2;
1140  stack[stackp - 1] = imitte + 1;
1141  stack[stackp + 19] = imax;
1142  stack[stackp + 39] = bfun4_(y, x[imitte]);
1143  stack[stackp + 59] = bmax;
1144  ++stackp;
1145  stack[stackp - 1] = imin;
1146  stack[stackp + 19] = imitte;
1147  stack[stackp + 39] = bmin;
1148  stack[stackp + 59] = bfun4_(y, x[imitte-1]);
1149  L3:
1150  if (stackp > 0) {
1151  goto L2;
1152  }
1153  /* ---- Region 4 */
1154  /* -------------------------------------------------------------------- */
1155  if (lauf[9] >= lauf[4] || lauf[19] >= lauf[14]) {
1156  if ((r__1 = (float)(y - yps4), fabs(r__1)) > 1e-8f) {
1157  //yps4 = (float)y;
1158  a7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
1159  y2 * (y2 * (y2 * (2.35944f - y2 * .56419f) - 72.9359f) +
1160  571.687f) - 5860.68f) + 40649.2f) - 320772.f) + 1684100.f)
1161  - 9694630.f) + 40816800.f) - 1.53575e8f) + 4.56662e8f) -
1162  9.86604e8f) + 1.16028e9f);
1163  b7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
1164  y2 * (y2 * (23.0312f - y2 * 7.33447f) - 234.143f) -
1165  2269.19f) + 45251.3f) - 234417.f) + 3599150.f) -
1166  7723590.f) + 86482900.f) - 2.91876e8f) + 8.06985e8f) -
1167  9.85386e8f) - 5.60505e8f);
1168  c7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
1169  y2 * (97.6203f - y2 * 44.0068f) + 1097.77f) - 25338.3f) +
1170  98079.1f) + 576054.f) - 2.3818e7f) + 22930200.f) -
1171  2.04467e8f) + 2.94262e8f) + 2.47157e8f) - 6.51523e8f);
1172  d7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
1173  228.563f - y2 * 161.358f) + 8381.97f) - 66431.2f) -
1174  303569.f) + 2240400.f) + 38311200.f) - 41501300.f) -
1175  99622400.f) + 2.70167e8f) - 2.63894e8f);
1176  e7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
1177  296.38f - y2 * 403.396f) + 23507.6f) - 66212.1f) -
1178  1.003e6f) + 468142.f) + 24620100.f) + 5569650.f) +
1179  1.40677e8f) - 63177100.f);
1180  f7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (125.591f -
1181  y2 * 726.113f) + 37544.8f) + 8820.94f) - 934717.f) -
1182  1931140.f) - 33289600.f) + 4073820.f) - 16984600.f);
1183  g7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (-260.198f - y2 *
1184  968.15f) + 37371.9f) + 79902.5f) - 186682.f) - 900010.f)
1185  + 7528830.f) - 1231650.f);
1186  h7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (-571.645f - y2 * 968.15f)
1187  + 23137.1f) + 72520.9f) + 153468.f) + 86407.6f) -
1188  610622.f);
1189  o7 = y * (y2 * (y2 * (y2 * (y2 * (-575.164f - y2 * 726.113f) +
1190  8073.15f) + 26538.5f) + 49883.8f) - 23586.5f);
1191  p7 = y * (y2 * (y2 * (y2 * (-352.467f - y2 * 403.396f) +
1192  953.655f) + 2198.86f) - 8009.1f);
1193  q7 = y * (y2 * (y2 * (-134.792f - y2 * 161.358f) - 271.202f) -
1194  622.056f);
1195  r7 = y * (y2 * (-29.7896f - y2 * 44.0068f) - 77.0535f);
1196  s7 = y * (-2.92264f - y2 * 7.33447f);
1197  t7 = y * -.56419f;
1198  a8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
1199  y2 * (y2 * (y2 * (y2 - 3.68288f) + 126.532f) - 955.194f)
1200  + 9504.65f) - 70946.1f) + 483737.f) - 2857210.f) +
1201  14464700.f) - 61114800.f) + 2.11107e8f) - 5.79099e8f) +
1202  1.17022e9f) - 1.5599e9f) + 1.02827e9f;
1203  b8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
1204  y2 * (y2 * (y2 * 14.f - 40.5117f) + 533.254f) + 3058.26f)
1205  - 55600.f) + 498334.f) - 2849540.f) + 13946500.f) -
1206  70135800.f) + 2.89676e8f) - 7.53828e8f) + 1.66421e9f) -
1207  2.28855e9f) + 1.5599e9f;
1208  c8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
1209  y2 * (y2 * 91 - 198.876f) - 1500.17f) + 48153.3f) -
1210  217801.f) - 1063520.f) + 1.4841e7f) - 46039600.f) +
1211  63349600.f) - 6.60078e8f) + 1.06002e9f) - 1.66421e9f) +
1212  1.17022e9f;
1213  d8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
1214  y2 * 364 - 567.164f) - 16493.7f) + 161461.f) + 280428.f)
1215  - 6890020.f) - 6876560.f) + 1.99846e8f) + 54036700.f) +
1216  6.60078e8f) - 7.53828e8f) + 5.79099e8f;
1217  e8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 *
1218  1001 - 1012.79f) - 55582.f) + 240373.f) + 1954700.f) -
1219  5257220.f) - 50101700.f) - 1.99846e8f) + 63349600.f) -
1220  2.89676e8f) + 2.11107e8f;
1221  f8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 2002 -
1222  1093.82f) - 106663.f) + 123052.f) + 3043160.f) +
1223  5257220.f) - 6876560.f) + 46039600.f) - 70135800.f) +
1224  61114800.f;
1225  g8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 3003 -
1226  486.14f) - 131337.f) - 123052.f) + 1954700.f) + 6890020.f)
1227  + 1.4841e7f) - 13946500.f) + 14464700.f;
1228  h8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 3432 + 486.14f) -
1229  106663.f) - 240373.f) + 280428.f) + 1063520.f) -
1230  2849540.f) + 2857210.f;
1231  o8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 3003 + 1093.82f) -
1232  55582.f) - 161461.f) - 217801.f) - 498334.f) + 483737.f;
1233  p8 = y2 * (y2 * (y2 * (y2 * (y2 * 2002 + 1012.79f) - 16493.7f) -
1234  48153.3f) - 55600.f) + 70946.1f;
1235  q8 = y2 * (y2 * (y2 * (y2 * 1001.f + 567.164f) - 1500.17f) -
1236  3058.26f) + 9504.65f;
1237  r8 = y2 * (y2 * (y2 * 364 + 198.876f) + 533.254f) + 955.194f;
1238  s8 = y2 * (y2 * 91.f + 40.5117f) + 126.532f;
1239  t8 = y2 * 14.f + 3.68288f;
1240  }
1241  ym2 = y * 2;
1242  for (i2 = 1; i2 <= 3; i2 += 2) {
1243  i__1 = lauf[(i2 + 1) * 5 - 1];
1244  for (i1 = lauf[i2 * 5 - 1]; i1 <= i__1; ++i1) {
1245  x2 = x[i1-1] * x[i1-1];
1246  ls_attenuation[i1-1] = fac * (exp(y2 - x2) * cos(x[i1-1] * ym2) - (a7 + x2 *
1247  (b7 + x2 * (c7 + x2 * (d7 + x2 * (e7 + x2 * (f7 + x2
1248  * (g7 + x2 * (h7 + x2 * (o7 + x2 * (p7 + x2 * (q7 +
1249  x2 * (r7 + x2 * (s7 + x2 * t7))))))))))))) / (a8 + x2
1250  * (b8 + x2 * (c8 + x2 * (d8 + x2 * (e8 + x2 * (f8 +
1251  x2 * (g8 + x2 * (h8 + x2 * (o8 + x2 * (p8 + x2 * (q8
1252  + x2 * (r8 + x2 * (s8 + x2 * (t8 + x2)))))))))))))));
1253  /* L4: */
1254  }
1255  }
1256  }
1257  /* ---- Region 3 */
1258  /* -------------------------------------------------------------------- */
1259  if (lauf[8] >= lauf[3] || lauf[18] >= lauf[13]) {
1260  if ((r__1 = (float)(y - yps3), fabs(r__1)) > 1e-8f) {
1261  //yps3 = (float)y;
1262  a5 = y * (y * (y * (y * (y * (y * (y * (y * (y *
1263  .564224f + 7.55895f) + 49.5213f) + 204.501f) + 581.746f)
1264  + 1174.8f) + 1678.33f) + 1629.76f) + 973.778f) + 272.102f;
1265  b5 = y * (y * (y * (y * (y * (y * (y * 2.25689f + 22.6778f)
1266  + 100.705f) + 247.198f) + 336.364f) + 220.843f) -
1267  2.34403f) - 60.5644f;
1268  c5 = y * (y * (y * (y * (y * 3.38534f + 22.6798f) + 52.8454f)
1269  + 42.5683f) + 18.546f) + 4.58029f;
1270  d5 = y * (y * (y * 2.25689f + 7.56186f) + 1.66203f) - .128922f;
1271  e5 = y * .564224f + 9.71457e-4f;
1272  a6 = y * (y * (y * (y * (y * (y * (y * (y * (y * (y +
1273  13.3988f) + 88.2674f) + 369.199f) + 1074.41f) + 2256.98f)
1274  + 3447.63f) + 3764.97f) + 2802.87f) + 1280.83f) +
1275  272.102f;
1276  b6 = y * (y * (y * (y * (y * (y * (y * (y * 5.f +
1277  53.5952f) + 266.299f) + 793.427f) + 1549.68f) + 2037.31f)
1278  + 1758.34f) + 902.306f) + 211.678f;
1279  c6 = y * (y * (y * (y * (y * (y * 10.f + 80.3928f) +
1280  269.292f) + 479.258f) + 497.302f) + 308.186f) + 78.866f;
1281  d6 = y * (y * (y * (y * 10.f + 53.5952f) + 92.7568f) +
1282  55.0293f) + 22.0353f;
1283  e6 = y * (y * 5.f + 13.3988f) + 1.49645f;
1284  }
1285  for (i2 = 1; i2 <= 3; i2 += 2) {
1286  i__1 = lauf[(i2 + 1) * 5 - 2];
1287  for (i1 = lauf[i2 * 5 - 2]; i1 <= i__1; ++i1) {
1288  x2 = x[i1-1] * x[i1-1];
1289  ls_attenuation[i1-1] = fac * (a5 + x2 * (b5 + x2 * (c5 + x2 * (d5 + x2 *
1290  e5)))) / (a6 + x2 * (b6 + x2 * (c6 + x2 * (d6 + x2 * (
1291  e6 + x2)))));
1292  /* L5: */
1293  }
1294  }
1295  }
1296  /* ---- Region 2 */
1297  /* -------------------------------------------------------------------- */
1298  if (lauf[7] >= lauf[2] || lauf[17] >= lauf[12]) {
1299  if ((r__1 = (float)(y - yps2), fabs(r__1)) > 1e-8f) {
1300  //yps2 = (float)y;
1301  a3 = y * (y2 * (y2 * (y2 * .56419f + 3.10304f) + 4.65456f) +
1302  1.05786f);
1303  b3 = y * (y2 * (y2 * 1.69257f + .56419f) + 2.962f);
1304  c3 = y * (y2 * 1.69257f - 2.53885f);
1305  d3 = y * .56419f;
1306  a4 = y2 * (y2 * (y2 * (y2 + 6.f) + 10.5f) + 4.5f) + .5625f;
1307  b4 = y2 * (y2 * (y2 * 4.f + 6.f) + 9.f) - 4.5f;
1308  c4 = y2 * (y2 * 6.f - 6.f) + 10.5f;
1309  d4 = y2 * 4.f - 6.f;
1310  }
1311  for (i2 = 1; i2 <= 3; i2 += 2) {
1312  i__1 = lauf[(i2 + 1) * 5 - 3];
1313  for (i1 = lauf[i2 * 5 - 3]; i1 <= i__1; ++i1) {
1314  x2 = x[i1-1] * x[i1-1];
1315  ls_attenuation[i1-1] = fac * (a3 + x2 * (b3 + x2 * (c3 + x2 * d3))) / (a4
1316  + x2 * (b4 + x2 * (c4 + x2 * (d4 + x2))));
1317  /* L6: */
1318  }
1319  }
1320  }
1321  /* ---- Region 1 */
1322  /* -------------------------------------------------------------------- */
1323  if (lauf[6] >= lauf[1] || lauf[16] >= lauf[11]) {
1324  if ((r__1 = (float)(y - yps1), fabs(r__1)) > 1e-8f) {
1325  //yps1 = (float)y;
1326  a1 = y * .5641896f;
1327  b1 = y2 + .5f;
1328  a2 = y2 * 4;
1329  }
1330 
1331  c1 = a1 * fac;
1332  for (i2 = 1; i2 <= 3; i2 += 2) {
1333  i__1 = lauf[(i2 + 1) * 5 - 4];
1334  for (i1 = lauf[i2 * 5 - 4]; i1 <= i__1; ++i1) {
1335  x2 = x[i1-1] * x[i1-1];
1336  b2 = b1 - x2;
1337  ls_attenuation[i1-1] = c1 * (b1 + x2) / (b2 * b2 + a2 * x2);
1338  /* L7: */
1339  }
1340  }
1341  }
1342  /* ---- Region 0 (Lorentz) */
1343  /* -------------------------------------------------------------------- */
1344  L8:
1345  if (lauf[5] >= lauf[0] || lauf[15] >= lauf[10]) {
1346  if ((r__1 = (float)(y - yps0), fabs(r__1)) > 1e-8f) {
1347  //yps0 = (float)y;
1348  a0 = y * .5641896f;
1349  }
1350 
1351  b0 = a0 * fac;
1352  for (i2 = 1; i2 <= 3; i2 += 2) {
1353  i__1 = lauf[(i2 + 1) * 5 - 5];
1354  for (i1 = lauf[i2 * 5 - 5]; i1 <= i__1; ++i1) {
1355  ls_attenuation[i1-1] = b0 / (x[i1-1] * x[i1-1] + y2);
1356  /* L9: */
1357  }
1358  }
1359  }
1360 }
1361 
1362 
1363 
1406 /*** ROUTINE COMPUTES THE VOIGT FUNCTION Y/PI*INTEGRAL FROM ***/
1407 /*** - TO + INFINITY OF EXP(-T*T)/(Y*Y+(X-T)*(X-T)) DT ***/
1408 void lineshape_voigt_drayson(Vector& ls_attenuation,
1409  Vector& ls_phase _U_,
1410  Vector& x,
1411  const Numeric f0,
1412  const Numeric gamma,
1413  const Numeric sigma,
1414  ConstVectorView f_grid)
1415 
1416 {
1417  const Index nf = f_grid.nelem();
1418 
1419  // seems not necessary for Doppler correction
1420  // extern const Numeric SQRT_NAT_LOG_2;
1421 
1422  // PI
1423  extern const Numeric PI;
1424 
1425  // constant sqrt(1/pi)
1426  const Numeric sqrt_invPI = sqrt(1/PI);
1427 
1428  // constant normalization factor for voigt
1429  Numeric fac = 1.0 / sigma * sqrt_invPI;
1430 
1431  Numeric B[22+1] = {0.,0.,.7093602e-7};
1432  Numeric RI[15+1] = {0};
1433  const Numeric XN[15+1] = {0.,10.,9.,8.,8.,7.,6.,5.,4.,3.,3.,3.,3.,3.,3.,3.};
1434  const Numeric YN[15+1] = {0.,.6,.6,.6,.5,.4,.4,.3,.3,.3,.3,1.,.9,.8,.7,.7};
1435  Numeric D0[25+1] = {0}, D1[25+1] = {0}, D2[25+1] = {0}, D3[25+1] = {0}, D4[25+1] = {0};
1436  Numeric HN[25+1] = {0};
1437  Numeric H = .201;
1438  const Numeric XX[3+1] = {0.,.5246476,1.65068,.7071068};
1439  const Numeric HH[3+1] = {0.,.2562121,.2588268e-1,.2820948};
1440  const Numeric NBY2[19+1] = {0.,9.5,9.,8.5,8.,7.5,7.,6.5
1441  ,6.,5.5,5.,4.5,4.,3.5,3.,2.5,2.,1.5,1.,.5};
1442  const Numeric C[21+1] = {0.,.7093602e-7,-.2518434e-6,.856687e-6,
1443  -.2787638e-5,.866074e-5,-.2565551e-4,.7228775e-4,
1444  -.1933631e-3,.4899520e-3,-.1173267e-2,.2648762e-2,
1445  -.5623190e-2,.1119601e-1,-.2084976e-1,.3621573e-1,
1446  -.5851412e-1,.8770816e-1,-.121664,.15584,-.184,.2};
1447  Numeric CO = 0;
1448  Numeric U, V, UU, VV, Y2, DX;
1449  int I, J, K, MAX, MIN, N, i1;
1450  //int TRU = 0;
1451 
1452 
1453  // variables needed in original c routine:
1454 
1455  // Ratio of the Lorentz halfwidth to the Doppler halfwidth
1456  Numeric Y = gamma / sigma;
1457 
1458  // frequency in units of Doppler, Note: Drayson accepts only
1459  // positive values
1460  for (i1=0; i1< (int) nf; i1++)
1461  {
1462  x[i1] = fabs( (f_grid[i1] - f0) )/ sigma;
1463  }
1464 
1465 
1466  //if (TRU) goto L104;
1467 /*** REGION I. COMPUTE DAWSON'S FUNCTION AT MESH POINTS ***/
1468  //TRU = 1;
1469  for (I=1; I<=15; I++)
1470  RI[I] = -I/2.;
1471  for (I=1; I<=25; I++)
1472  {
1473  HN[I] = H*(I-.5);
1474  CO = 4.*HN[I]*HN[I]/25.-2.;
1475  for (J=2; J<=21; J++)
1476  B[J+1] = CO*B[J]-B[J-1]+C[J];
1477  D0[I] = HN[I]*(B[22]-B[21])/5.;
1478  D1[I] = 1.-2.*HN[I]*D0[I];
1479  D2[I] = (HN[I]*D1[I]+D0[I])/RI[2];
1480  D3[I] = (HN[I]*D2[I]+D1[I])/RI[3];
1481  D4[I] = (HN[I]*D3[I]+D2[I])/RI[4];
1482  }
1483 //L104:
1484  for (K=0; K<(int) nf; K++)
1485  {
1486  if ((x[K]-5.) < 0.) goto L105; else goto L112;
1487  L105: if ((Y-1.) <= 0.) goto L110; else goto L106;
1488  L106: if (x[K] > 1.85*(3.6-Y)) goto L112;
1489  /*** REGION II: CONTINUED FRACTION. COMPUTE NUMBER OF TERMS NEEDED. ***/
1490  if (Y < 1.45) goto L107;
1491  I = (int) (Y+Y);
1492  goto L108;
1493  L107: I = (int) (11.*Y);
1494  L108: J = (int) (x[K]+x[K]+1.85);
1495  MAX = (int) (XN[J]*YN[I]+.46);
1496  MIN = (int) ( (16 < 21-2*MAX) ? 16 : 21-2*MAX );
1497  /*** EVALUATE CONTINUED FRACTION ***/
1498  UU = Y;
1499  VV = x[K];
1500  for (J=MIN; J <= 19; J++)
1501  {
1502  U = NBY2[J]/(UU*UU+VV*VV);
1503  UU = Y+U*UU;
1504  VV = x[K]-U*VV;
1505  }
1506  ls_attenuation[K] = UU/(UU*UU+VV*VV)/1.772454*fac;
1507  continue;
1508  L110: Y2 = Y*Y;
1509  if (x[K]+Y >= 5.) goto L113;
1510  /*** REGION I. COMPUTE DAWSON'S FUNCTION AT X FROM TAYLOR SERIES. ***/
1511  N = (int) (x[K]/H);
1512  DX = x[K]-HN[N+1];
1513  U = (((D4[N+1]*DX+D3[N+1])*DX+D2[N+1])*DX+D1[N+1])*DX+D0[N+1];
1514  V = 1.-2.*x[K]*U;
1515  /*** TAYLOR SERIES EXPANSION ABOUT Y=0.0 ***/
1516  VV = exp(Y2-x[K]*x[K])*cos(2.*x[K]*Y)/1.128379-Y*V;
1517  UU = -Y;
1518  MAX = (int) (5.+(12.5-x[K])*.8*Y);
1519  for (I=2; I<=MAX; I+=2)
1520  {
1521  U = (x[K]*V+U)/RI[I];
1522  V = (x[K]*U+V)/RI[I+1];
1523  UU = -UU*Y2;
1524  VV = VV+V*UU;
1525  }
1526  ls_attenuation[K] = 1.128379*VV*fac;
1527  continue;
1528  L112: Y2 = Y*Y;
1529  if (Y < 11.-.6875*x[K]) goto L113;
1530  /*** REGION IIIB: 2-POINT GAUSS-HERMITE QUADRATURE. ***/
1531  U = x[K]-XX[3];
1532  V = x[K]+XX[3];
1533  ls_attenuation[K] = Y*(HH[3]/(Y2+U*U)+HH[3]/(Y2+V*V))*fac;
1534  continue;
1535  /*** REGION IIIA: 4-POINT GAUSS-HERMITE QUADRATURE. ***/
1536  L113: U = x[K]-XX[1];
1537  V = x[K]+XX[1];
1538  UU = x[K]-XX[2];
1539  VV = x[K]+XX[2];
1540  ls_attenuation[K] = Y*(HH[1]/(Y2+U*U)+HH[1]/(Y2+V*V)+HH[2]/(Y2+UU*UU)+HH[2]/(Y2+
1541  VV*VV))*fac;
1542  continue;
1543  }
1544 }
1545 
1546 
1558 void chi_cousin(
1559  Numeric& chi,
1560  const Numeric& df )
1561 {
1562  // Conversion factor from Hz to cm-1
1563  extern const Numeric HZ2CM;
1564 
1565  const Numeric n2 = 0.79;
1566  const Numeric o2 = 0.21;
1567 
1568  const Numeric df_cm = df * HZ2CM;
1569  const Numeric df_cm_abs = fabs( df_cm );
1570 
1571  chi = 0;
1572 
1573  // N2 term
1574  if( df_cm_abs <= 5 )
1575  { chi += n2; }
1576  else if( df_cm_abs <= 22 )
1577  { chi += n2 * 1.968 * exp( -0.1354 * df_cm_abs ); }
1578  else if( df_cm_abs <= 50 )
1579  { chi += n2 * 0.160 * exp( -0.0214 * df_cm_abs ); }
1580  else
1581  { chi += n2 * 0.162 * exp( -0.0216 * df_cm_abs ); }
1582 
1583  // O2 term
1584  if( df_cm_abs <= 5 )
1585  { chi += o2; }
1586  else if( df_cm_abs <= 22 )
1587  { chi += o2 * 1.968 * exp( -0.1354 * df_cm_abs ); }
1588  else if( df_cm_abs <= 50 )
1589  { chi += o2 * 0.160 * exp( -0.0214 * df_cm_abs ); }
1590  else
1591  { chi += o2 * 0.162 * exp( -0.0216 * df_cm_abs ); }
1592 }
1593 
1594 
1595 
1607 void lineshape_CO2_lorentz(Vector& ls_attenuation,
1608  Vector& ls_phase _U_,
1609  Vector& X _U_,
1610  const Numeric f0,
1611  const Numeric gamma,
1612  const Numeric sigma _U_,
1613  ConstVectorView f_grid)
1614 {
1615  const Index nf = f_grid.nelem();
1616 
1617  // PI:
1618  extern const Numeric PI;
1619 
1620  // Some constant variables
1621  const Numeric gamma2 = gamma * gamma;
1622  const Numeric fac = gamma/PI;
1623 
1624  for ( Index i=0; i<nf; ++i )
1625  {
1626  // f-f0
1627  const Numeric df = f_grid[i] - f0;
1628 
1629  // The chi factor
1630  Numeric chi;
1631  chi_cousin( chi, df );
1632 
1633  // chi * Lorentz
1634  ls_attenuation[i] = chi * fac / ( df * df + gamma2 );
1635  }
1636 }
1637 
1638 
1639 
1651 void lineshape_CO2_drayson(Vector& ls_attenuation,
1652  Vector& ls_phase _U_,
1653  Vector& X,
1654  const Numeric f0,
1655  const Numeric gamma,
1656  const Numeric sigma,
1657  ConstVectorView f_grid)
1658 {
1659  lineshape_voigt_drayson( ls_attenuation, ls_phase, X, f0, gamma, sigma, f_grid );
1660 
1661  const Index nf = f_grid.nelem();
1662  for ( Index i=0; i<nf; ++i )
1663  {
1664  // f-f0
1665  const Numeric df = f_grid[i] - f0;
1666 
1667  // The chi factor
1668  Numeric chi;
1669  chi_cousin( chi, df );
1670 
1671  // chi * Lorentz
1672  ls_attenuation[i] *= chi;
1673  }
1674 }
1675 
1676 
1699 void faddeeva_algorithm_916( Vector& ls_attenuation,
1700  Vector& ls_phase,
1701  Vector& xvector,
1702  const Numeric f0,
1703  const Numeric gamma,
1704  const Numeric sigma,
1705  ConstVectorView f_grid)
1706 
1707 {
1708  const Index nf = f_grid.nelem();
1709 
1710  // PI
1711  extern const Numeric PI;
1712 
1713  // constant sqrt(1/pi)
1714  const Numeric sqrt_invPI = sqrt(1/PI);
1715 
1716  // constant normalization factor for voigt
1717  const Numeric fac = sqrt_invPI / sigma;
1718 
1719  // Ratio of the Lorentz halfwidth to the Doppler halfwidth
1720  const Numeric y = gamma / (sigma);
1721 
1722  // frequency in units of Doppler
1723  for (Index ii=0; ii<nf; ii++)
1724  {
1725  xvector[ii] = (f_grid[ii] - f0) / sigma;
1726  std::complex<Numeric> z(xvector[ii], y);
1727 
1728  z = Faddeeva::w(z);
1729 
1730  ls_attenuation[ii] = fac * z.real();
1731  ls_phase[ii] = fac * z.imag();
1732  }
1733 }
1734 
1735 
1765 void hui_etal_1978_lineshape( Vector& ls_attenuation,
1766  Vector& ls_phase,
1767  Vector& xvector,
1768  const Numeric f0,
1769  const Numeric gamma,
1770  const Numeric sigma,
1771  ConstVectorView f_grid)
1772 
1773 {
1774  const Index nf = f_grid.nelem();
1775 
1776  // PI
1777  extern const Numeric PI;
1778 
1779  // constant sqrt(1/pi)
1780  const Numeric sqrt_invPI = sqrt(1/PI);
1781 
1782  // constant normalization factor for voigt
1783  const Numeric fac = sqrt_invPI / sigma;
1784 
1785  // Ratio of the Lorentz halfwidth to the Doppler halfwidth
1786  const Numeric y = gamma / (sigma);
1787 
1788  // frequency in units of Doppler
1789  for (Index ii=0; ii<nf; ii++)
1790  {
1791  xvector[ii] = (f_grid[ii] - f0) / sigma;
1792 
1793  // Note that this works but I don't know why there is a difference
1794  // between the theory described above and this practical implementation.
1795  const std::complex<Numeric> z(y , - xvector[ii]);
1796 
1797  const std::complex<Numeric> A = (((((.5641896*z+5.912626)*z+30.18014)*z+
1798  93.15558)*z+181.9285)*z+214.3824)*z+122.6079;
1799  const std::complex<Numeric> B = ((((((z+10.47986)*z+53.99291)*z+170.3540)*z+
1800  348.7039)*z+457.3345)*z+352.7306)*z+122.6079;
1801  const std::complex<Numeric> C = A / B;
1802 
1803  ls_attenuation[ii] = fac * C.real();
1804  ls_phase[ii] = fac * C.imag();
1805  }
1806 }
1807 
1808 //------------------------------------------------------------------------
1809 // Normalization Functions
1810 //------------------------------------------------------------------------
1811 
1812 
1821 void lineshape_norm_no_norm(Vector& fac,
1822  const Numeric f0 _U_,
1823  ConstVectorView f_grid _U_,
1824  const Numeric T _U_)
1825 {
1826  fac = 1.0;
1827 }
1828 
1829 
1830 
1839 void lineshape_norm_linear(Vector& fac,
1840  const Numeric f0,
1841  ConstVectorView f_grid,
1842  const Numeric T _U_)
1843 {
1844  const Index nf = f_grid.nelem();
1845 
1846  // Abs(f0) is constant in the loop:
1847  const Numeric abs_f0 = abs(f0);
1848 
1849  for ( Index i=0; i<nf; ++i )
1850  {
1851  fac[i] = f_grid[i] / abs_f0;
1852  }
1853 }
1854 
1863 void lineshape_norm_quadratic(Vector& fac,
1864  const Numeric f0,
1865  ConstVectorView f_grid,
1866  const Numeric T _U_)
1867 {
1868  const Index nf = f_grid.nelem();
1869 
1870  // don't do this for the whole loop
1871  const Numeric f0_2 = f0 * f0;
1872 
1873  for ( Index i=0; i<nf; ++i )
1874  {
1875  fac[i] = (f_grid[i] * f_grid[i]) / f0_2;
1876  }
1877 }
1878 
1893 void lineshape_norm_VVH(Vector& fac,
1894  const Numeric f0,
1895  ConstVectorView f_grid,
1896  const Numeric T)
1897 {
1898  extern const Numeric PLANCK_CONST;
1899  extern const Numeric BOLTZMAN_CONST;
1900 
1901  const Index nf = f_grid.nelem();
1902 
1903  // 2kT is constant for the loop
1904  const Numeric kT = 2.0 * BOLTZMAN_CONST * T;
1905 
1906  // denominator is constant for the loop
1907  const Numeric denom = abs(f0) * tanh( PLANCK_CONST * abs(f0) / kT );
1908 
1909  for ( Index i=0; i<nf; ++i )
1910  {
1911  fac[i] = f_grid[i] * tanh( PLANCK_CONST * f_grid[i] / kT ) /
1912  denom;
1913  }
1914 }
1915 
1916 
1917 
1918 
1919 //------------------------------------------------------------------------
1920 // Available Lineshapes
1921 //------------------------------------------------------------------------
1922 
1924 namespace global_data {
1925 Array<LineshapeRecord> lineshape_data;
1926 }
1927 
1928 void define_lineshape_data()
1929 {
1930  using global_data::lineshape_data;
1931 
1932  const bool PHASE = true;
1933  const bool NO_PHASE = false;
1934 
1935  // Initialize to empty, just in case.
1936  lineshape_data.resize(0);
1937 
1938  lineshape_data.push_back
1939  (LineshapeRecord
1940  ("no_shape",
1941  "This lineshape does nothing. It only exists, because formally\n"
1942  "you have to specify a lineshape also for continuum tags.",
1943  lineshape_no_shape, NO_PHASE));
1944 
1945  lineshape_data.push_back
1946  (LineshapeRecord
1947  ("Lorentz",
1948  "The Lorentz line shape.",
1949  lineshape_lorentz, PHASE));
1950 
1951  lineshape_data.push_back
1952  (LineshapeRecord
1953  ("Doppler",
1954  "The Doppler line shape.",
1955  lineshape_doppler, NO_PHASE));
1956 
1957  lineshape_data.push_back
1958  (LineshapeRecord
1959  ("Voigt_Kuntz6",
1960  "The Voigt line shape. Approximation by Kuntz: Accuracy 2*10-6",
1961  lineshape_voigt_kuntz6, NO_PHASE));
1962 
1963  lineshape_data.push_back
1964  (LineshapeRecord
1965  ("Voigt_Kuntz3",
1966  "The Voigt line shape. Approximation by Kuntz: Accuracy 2*10-3",
1967  lineshape_voigt_kuntz3, NO_PHASE));
1968 
1969  lineshape_data.push_back
1970  (LineshapeRecord
1971  ("Voigt_Kuntz4",
1972  "The Voigt line shape. Approximation by Kuntz: Accuracy 2*10-4",
1973  lineshape_voigt_kuntz4, NO_PHASE));
1974 
1975  lineshape_data.push_back
1976  (LineshapeRecord
1977  ("Voigt_Drayson",
1978  "The Voigt line shape. Approximation by Drayson.",
1979  lineshape_voigt_drayson, NO_PHASE));
1980 
1981  lineshape_data.push_back
1982  (LineshapeRecord
1983  ("CO2_Lorentz",
1984  "Special line shape for CO2 in the infrared, neglecting Doppler\n"
1985  "broadening and details of line mixing. The line shape can be\n"
1986  "expressed as\n"
1987  " chi(f,f0) * Lorentz(f,f0) \n"
1988  "\n"
1989  "The chi-factor follows Cousin et al. 1985. Broadening by N2 and O2\n"
1990  "is considered, while self-broadening (CO2-CO2) is neglected."
1991  "\n"
1992  "NOTE: Temperature dependency is not yet included. The chi factor is\n"
1993  "valid for 238 K.",
1994  lineshape_CO2_lorentz, NO_PHASE));
1995 
1996  lineshape_data.push_back
1997  (LineshapeRecord
1998  ("CO2_Drayson",
1999  "Special line shape for CO2 in the infrared, neglecting details of\n"
2000  "line mixing. The line shape can be expressed as\n"
2001  " chi(f,f0) * Drayson(f,f0) \n"
2002  "\n"
2003  "The chi-factor follows Cousin et al. 1985. Broadening by N2 and O2\n"
2004  "is considered, while self-broadening (CO2-CO2) is neglected.\n"
2005  "\n"
2006  "NOTE: Temperature dependency is not yet included. The chi factor is\n"
2007  "valid for 238 K.",
2008  lineshape_CO2_drayson, NO_PHASE));
2009 
2010  lineshape_data.push_back
2011  (LineshapeRecord
2012  ("Faddeeva_Algorithm_916",
2013  "Voigt and Faraday-Voigt function as per Faddeeva function solution by JPL.\n"
2014  "Line shape is considered from w(z)=exp(-z^2)*erfc(-iz) where z=v'+ia, and \n"
2015  "v' is a Doppler weighted freqeuncy parameter and a is a Doppler weighted \n"
2016  "pressure parameter.",
2017  faddeeva_algorithm_916, PHASE));
2018 
2019  lineshape_data.push_back
2020  (LineshapeRecord
2021  ("Hui_etal_1978",
2022  "Classic line shape. Solving the complex error function returns both parts\n"
2023  "of the refractive index.",
2024  hui_etal_1978_lineshape, PHASE));
2025 }
2026 
2029 namespace global_data {
2030 Array<LineshapeNormRecord> lineshape_norm_data;
2031 }
2032 
2033 void define_lineshape_norm_data()
2034 {
2035  using global_data::lineshape_norm_data;
2036 
2037  // Initialize to empty, just in case.
2038  lineshape_norm_data.resize(0);
2039 
2040  lineshape_norm_data.push_back
2041  (LineshapeNormRecord
2042  ("no_norm",
2043  "No normalization of the lineshape.",
2044  lineshape_norm_no_norm));
2045 
2046  lineshape_norm_data.push_back
2047  (LineshapeNormRecord
2048  ("quadratic",
2049  "Quadratic normalization of the lineshape with (f/f0)^2.",
2050  lineshape_norm_quadratic));
2051 
2052  lineshape_norm_data.push_back
2053  (LineshapeNormRecord
2054  ("VVH",
2055  "Van Vleck Huber normalization of the lineshape with\n"
2056  " (f*tanh(h*f/(2*k*T))) / (f0*tanh(h*f0/(2*k*T))).\n"
2057  " The denominator is a result of catalogue intensities.",
2058  lineshape_norm_VVH));
2059 }
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:35
PLANCK_CONST
const Numeric PLANCK_CONST
Global constant, the Planck constant [Js].
hui_etal_1978_lineshape
void hui_etal_1978_lineshape(Vector &ls_attenuation, Vector &ls_phase, Vector &xvector, const Numeric f0, const Numeric gamma, const Numeric sigma, ConstVectorView f_grid)
Definition: lineshapes.cc:1765
N
#define N
Definition: rng.cc:195
bfun4_
long bfun4_(Numeric y, Numeric x)
Definition: lineshapes.cc:938
C
#define C(a, b)
Definition: Faddeeva.cc:254
BOLTZMAN_CONST
const Numeric BOLTZMAN_CONST
Global constant, the Boltzmann constant [J/K].
lineshape_voigt_kuntz3
void lineshape_voigt_kuntz3(Vector &ls_attenuation, Vector &ls_phase, Vector &x, const Numeric f0, const Numeric gamma, const Numeric sigma, ConstVectorView f_grid)
Definition: lineshapes.cc:650
lineshape_norm_no_norm
void lineshape_norm_no_norm(Vector &fac, const Numeric f0, ConstVectorView f_grid, const Numeric T)
Definition: lineshapes.cc:1821
Vector
The Vector class.
Definition: matpackI.h:556
fac
Numeric fac(const Index n)
fac
Definition: math_funcs.cc:68
lineshape_norm_quadratic
void lineshape_norm_quadratic(Vector &fac, const Numeric f0, ConstVectorView f_grid, const Numeric T)
Definition: lineshapes.cc:1863
lineshape_doppler
void lineshape_doppler(Vector &ls_attenuation, Vector &ls_phase, Vector &x, const Numeric f0, const Numeric gamma, const Numeric sigma, ConstVectorView f_grid)
Definition: lineshapes.cc:123
w
cmplx FADDEEVA() w(cmplx z, double relerr)
Definition: Faddeeva.cc:679
bfun6_
long bfun6_(Numeric y, Numeric x)
Definition: lineshapes.cc:153
lineshape_norm_VVH
void lineshape_norm_VVH(Vector &fac, const Numeric f0, ConstVectorView f_grid, const Numeric T)
Definition: lineshapes.cc:1893
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.cc:180
define_lineshape_data
void define_lineshape_data()
Definition: lineshapes.cc:1928
global_data::lineshape_data
const Array< LineshapeRecord > lineshape_data
The lookup data for the different lineshapes.
Definition: lineshapes.cc:1925
array.h
This file contains the definition of Array.
lineshape_norm_linear
void lineshape_norm_linear(Vector &fac, const Numeric f0, ConstVectorView f_grid, const Numeric T)
Definition: lineshapes.cc:1839
arts.h
The global header file for ARTS.
LineshapeNormRecord
Lineshape related normalization function information.
Definition: absorption.h:107
PI
const Numeric PI
max
#define max(a, b)
Definition: continua.cc:20461
b0
#define b0
Definition: continua.cc:21471
lineshape_lorentz
void lineshape_lorentz(Vector &ls_attenuation, Vector &ls_phase, Vector &X, const Numeric f0, const Numeric gamma, const Numeric sigma, ConstVectorView f_grid)
Definition: lineshapes.cc:86
HZ2CM
const Numeric HZ2CM
Global constant, converts Hz to cm-1.
lineshape_voigt_kuntz4
void lineshape_voigt_kuntz4(Vector &ls_attenuation, Vector &ls_phase, Vector &x, const Numeric f0, const Numeric gamma, const Numeric sigma, ConstVectorView f_grid)
Definition: lineshapes.cc:1024
abs
#define abs(x)
Definition: continua.cc:20458
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:29
absorption.h
Declarations required for the calculation of absorption coefficients.
global_data::lineshape_norm_data
const Array< LineshapeNormRecord > lineshape_norm_data
The lookup data for the different normalization factors to the lineshapes.
Definition: lineshapes.cc:2030
faddeeva_algorithm_916
void faddeeva_algorithm_916(Vector &ls_attenuation, Vector &ls_phase, Vector &xvector, const Numeric f0, const Numeric gamma, const Numeric sigma, ConstVectorView f_grid)
Definition: lineshapes.cc:1699
Array
This can be used to make arrays out of anything.
Definition: array.h:40
lineshape_CO2_lorentz
void lineshape_CO2_lorentz(Vector &ls_attenuation, Vector &ls_phase, Vector &X, const Numeric f0, const Numeric gamma, const Numeric sigma, ConstVectorView f_grid)
Definition: lineshapes.cc:1607
min
#define min(a, b)
Definition: continua.cc:20460
global_data
Definition: agenda_record.cc:33
ConstVectorView
A constant view of a Vector.
Definition: matpackI.h:292
LineshapeRecord
Lineshape related information.
Definition: absorption.h:58
define_lineshape_norm_data
void define_lineshape_norm_data()
Definition: lineshapes.cc:2033
lineshape_CO2_drayson
void lineshape_CO2_drayson(Vector &ls_attenuation, Vector &ls_phase, Vector &X, const Numeric f0, const Numeric gamma, const Numeric sigma, ConstVectorView f_grid)
Definition: lineshapes.cc:1651
_U_
#define _U_
Definition: config.h:167
lineshape_no_shape
void lineshape_no_shape(Vector &ls_attenuation, Vector &ls_phase, Vector &X, const Numeric f0, const Numeric gamma, const Numeric sigma, ConstVectorView f_grid)
Definition: lineshapes.cc:60
lineshape_voigt_kuntz6
void lineshape_voigt_kuntz6(Vector &ls_attenuation, Vector &ls_phase, Vector &x, const Numeric f0, const Numeric gamma, const Numeric sigma, ConstVectorView f_grid)
Definition: lineshapes.cc:236
lineshape_voigt_drayson
void lineshape_voigt_drayson(Vector &ls_attenuation, Vector &ls_phase, Vector &x, const Numeric f0, const Numeric gamma, const Numeric sigma, ConstVectorView f_grid)
Definition: lineshapes.cc:1408
bfun3_
long int bfun3_(Numeric y, Numeric x)
Definition: lineshapes.cc:561
chi_cousin
void chi_cousin(Numeric &chi, const Numeric &df)
Definition: lineshapes.cc:1558