ARTS  1.0.222
lineshapes.cc
Go to the documentation of this file.
1 /* Copyright (C) 2000, 2001 Axel von Engeln <engeln@uni-bremen.de>
2  Stefan Buehler <sbuehler@uni-bremen.de>
3 
4  This program is free software; you can redistribute it and/or modify it
5  under the terms of the GNU General Public License as published by the
6  Free Software Foundation; either version 2, or (at your option) any
7  later version.
8 
9  This program is distributed in the hope that it will be useful,
10  but WITHOUT ANY WARRANTY; without even the implied warranty of
11  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12  GNU General Public License for more details.
13 
14  You should have received a copy of the GNU General Public License
15  along with this program; if not, write to the Free Software
16  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
17  USA. */
18 
31 #include <math.h>
32 #include "arts.h"
33 #include "matpackI.h"
34 #include "array.h"
35 #include "absorption.h"
36 
37 
57 void lineshape_no_shape( Vector& /* ls */,
58  Vector& /* X */,
59  Numeric /* f0 */,
60  Numeric /* gamma */,
61  Numeric /* sigma */,
62  VectorView /* f_mono */,
63  const Index /* nf */)
64 {
65  // This function should never be called so throw an error here:
66  throw runtime_error("The no_shape lineshape is only a placeholder, but you tried\n"
67  "to use it like a real lineshape.");
68 }
69 
70 
83 void lineshape_lorentz(Vector& ls,
84  Vector& /* X */,
85  Numeric f0,
86  Numeric gamma,
87  Numeric /* sigma */,
88  VectorView f_mono,
89  const Index nf)
90 {
91  // FIXME: nf is actually redundant. Could be thrown out in the
92  // future. For now, let's do an assertion that at least it is
93  // correct:
94  assert( nf==f_mono.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[i] = fac / ( (f_mono[i]-f0) * (f_mono[i]-f0) + gamma2 );
107  }
108 }
109 
122 void lineshape_doppler(Vector& ls,
123  Vector& /* x */,
124  Numeric f0,
125  Numeric /* gamma */,
126  Numeric sigma,
127  VectorView f_mono,
128  const Index nf)
129 {
130  // FIXME: nf is actually redundant. Could be thrown out in the
131  // future. For now, let's do an assertion that at least it is
132  // correct:
133  assert( nf==f_mono.nelem() );
134 
135  // SQRT(PI):
136  extern const Numeric PI;
137  static const Numeric sqrtPI = sqrt(PI);
138 
139  // assert( ls.nelem() == nf );
140 
141  Numeric sigma2 = sigma * sigma;
142  Numeric fac = 1.0 / (sqrtPI * sigma);
143 
144  for ( Index i=0; i<nf ; ++i )
145  {
146  ls[i] = fac * exp( - pow( f_mono[i]-f0, 2) / sigma2 );
147  }
148 }
149 
150 
151 
152 //------------------------------------------------------------------------
153 
154 // help function for lineshape_voigt_kuntz1
155 long bfun6_(Numeric y, Numeric x)
156 {
157  /* System generated locals */
158  long int ret_val;
159 
160  /* Local variables */
161  static Numeric s;
162 
163  /* -------------------------------------------------------------------- */
164  s = fabs(x) + y;
165  if (s >= 15.f) {
166  ret_val = 1;
167  } else if (s >= 5.5f) {
168  ret_val = 2;
169  } else if (y >= fabs(x) * .195f - .176f) {
170  ret_val = 3;
171  } else {
172  ret_val = 4;
173  }
174  return ret_val;
175 } /* bfun6_ */
176 
177 
178 
239 void lineshape_voigt_kuntz6(Vector& ls,
240  Vector& x,
241  Numeric f0,
242  Numeric gamma,
243  Numeric sigma,
244  VectorView f_mono,
245  const Index nf)
246 
247 {
248  // FIXME: nf is actually redundant. Could be thrown out in the
249  // future. For now, let's do an assertion that at least it is
250  // correct:
251  assert( nf==f_mono.nelem() );
252 
253 
254  // seems not necessary for Doppler correction
255  // extern const Numeric SQRT_NAT_LOG_2;
256 
257  // PI
258  extern const Numeric PI;
259 
260  // constant sqrt(1/pi)
261  const Numeric sqrt_invPI = sqrt(1/PI);
262 
263  // constant normalization factor for voigt
264  Numeric fac = 1.0 / sigma * sqrt_invPI;
265 
266 
267  /* Initialized data */
268 
269  static Numeric yps1 = -1.0;
270  static Numeric yps2 = -1.0;
271  static Numeric yps3 = -1.0;
272  static Numeric yps4 = -1.0;
273 
274  /* System generated locals */
275  long int i__1, i__2;
276  Numeric r__1;
277 
278  /* Local variables */
279  static long int bmin, lauf[16] /* was [4][4] */, bmax;
280  static long int imin, imax, stack[80] /* was [20][4] */;
281  static Numeric a1, a2, a3, a4, a5, a6, a8, b8, c8, d8, e8, f8, g8, h8, a7,
282  b7, c7, d7, e7, f7, o8, p8, q8, r8, s8, t8, g7, h7, o7, p7, q7,
283  r7, s7, t7, b6, c6, d6, e6, b5, c5, d5, e5, b4, c4, d4, b3, c3,
284  d3, b1, y2;
285  static long int i2, i1;
286  static Numeric x2, b2, c1;
287  static long int stackp, imitte;
288  static Numeric ym2;
289 
290 
291  // variables needed in original c routine:
292 
293  // Ratio of the Lorentz halfwidth to the Doppler halfwidth
294  Numeric y = gamma / sigma;
295 
296  // frequency in units of Doppler
297  for (i1=0; i1< (int) nf; i1++)
298  {
299  x[(Index)i1] = (f_mono[(Index)i1] - f0) / sigma;
300  }
301 
302  /* Parameter adjustments */
303  // this does not work for variables type Vector, corrected by
304  // adjusting the index to array ls and x
305  //--ls;
306  //--x;
307 
308  /* Function Body */
309  y2 = y * y;
310  if (y >= 15.0 || x[0] >= 15.0 || x[nf-1] <= -15.0) {
311  lauf[0] = 1;
312  lauf[4] = nf;
313  lauf[8] = nf;
314  lauf[12] = 0;
315  goto L7;
316  }
317  for (i2 = 1; i2 <= 4; ++i2) {
318  for (i1 = 1; i1 <= 4; ++i1) {
319  lauf[i1 + (i2 << 2) - 5] = i2 % 2 * (nf + 1);
320  /* L1: */
321  }
322  }
323  stackp = 1;
324  stack[stackp - 1] = 1;
325  stack[stackp + 19] = nf;
326  stack[stackp + 39] = bfun6_(y, x[0]);
327  stack[stackp + 59] = bfun6_(y, x[nf-1]);
328  L2:
329  imin = stack[stackp - 1];
330  imax = stack[stackp + 19];
331  bmin = stack[stackp + 39];
332  bmax = stack[stackp + 59];
333  if (bmin == bmax) {
334  if (x[(Index)imax-1] < 0.f) {
335  /* Computing MIN */
336  i__1 = imin, i__2 = lauf[bmin - 1];
337  lauf[bmin - 1] = min(i__1,i__2);
338  /* Computing MAX */
339  i__1 = imax, i__2 = lauf[bmax + 3];
340  lauf[bmax + 3] = max(i__1,i__2);
341  --stackp;
342  goto L3;
343  } else if (x[(Index)imin-1] >= 0.f) {
344  /* Computing MIN */
345  i__1 = imin, i__2 = lauf[bmin + 7];
346  lauf[bmin + 7] = min(i__1,i__2);
347  /* Computing MAX */
348  i__1 = imax, i__2 = lauf[bmax + 11];
349  lauf[bmax + 11] = max(i__1,i__2);
350  --stackp;
351  goto L3;
352  }
353  }
354  imitte = (imax + imin) / 2;
355  stack[stackp - 1] = imitte + 1;
356  stack[stackp + 19] = imax;
357  stack[stackp + 39] = bfun6_(y, x[(Index)imitte]);
358  stack[stackp + 59] = bmax;
359  ++stackp;
360  stack[stackp - 1] = imin;
361  stack[stackp + 19] = imitte;
362  stack[stackp + 39] = bmin;
363  stack[stackp + 59] = bfun6_(y, x[(Index)imitte-1]);
364  L3:
365  if (stackp > 0) {
366  goto L2;
367  }
368  /* ---- Region 4 */
369  /* -------------------------------------------------------------------- */
370  if (lauf[7] >= lauf[3] || lauf[15] >= lauf[11]) {
371  if ((r__1 = y - yps4, fabs(r__1)) > 1e-8f) {
372  yps4 = y;
373  a7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
374  y2 * (y2 * (y2 * (2.35944f - y2 * .56419f) - 72.9359f) +
375  571.687f) - 5860.68f) + 40649.2f) - 320772.f) + 1684100.f)
376  - 9694630.f) + 40816800.f) - 1.53575e8f) + 4.56662e8f) -
377  9.86604e8f) + 1.16028e9f);
378  b7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
379  y2 * (y2 * (23.0312f - y2 * 7.33447f) - 234.143f) -
380  2269.19f) + 45251.3f) - 234417.f) + 3599150.f) -
381  7723590.f) + 86482900.f) - 2.91876e8f) + 8.06985e8f) -
382  9.85386e8f) - 5.60505e8f);
383  c7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
384  y2 * (97.6203f - y2 * 44.0068f) + 1097.77f) - 25338.3f) +
385  98079.1f) + 576054.f) - 2.3818e7f) + 22930200.f) -
386  2.04467e8f) + 2.94262e8f) + 2.47157e8f) - 6.51523e8f);
387  d7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
388  228.563f - y2 * 161.358f) + 8381.97f) - 66431.2f) -
389  303569.f) + 2240400.f) + 38311200.f) - 41501300.f) -
390  99622400.f) + 2.70167e8f) - 2.63894e8f);
391  e7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
392  296.38f - y2 * 403.396f) + 23507.6f) - 66212.1f) -
393  1.003e6f) + 468142.f) + 24620100.f) + 5569650.f) +
394  1.40677e8f) - 63177100.f);
395  f7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (125.591f -
396  y2 * 726.113f) + 37544.8f) + 8820.94f) - 934717.f) -
397  1931140.f) - 33289600.f) + 4073820.f) - 16984600.f);
398  g7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (-260.198f - y2 *
399  968.15f) + 37371.9f) + 79902.5f) - 186682.f) - 900010.f)
400  + 7528830.f) - 1231650.f);
401  h7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (-571.645f - y2 * 968.15f)
402  + 23137.1f) + 72520.9f) + 153468.f) + 86407.6f) -
403  610622.f);
404  o7 = y * (y2 * (y2 * (y2 * (y2 * (-575.164f - y2 * 726.113f) +
405  8073.15f) + 26538.5f) + 49883.8f) - 23586.5f);
406  p7 = y * (y2 * (y2 * (y2 * (-352.467f - y2 * 403.396f) +
407  953.655f) + 2198.86f) - 8009.1f);
408  q7 = y * (y2 * (y2 * (-134.792f - y2 * 161.358f) - 271.202f) -
409  622.056f);
410  r7 = y * (y2 * (-29.7896f - y2 * 44.0068f) - 77.0535f);
411  s7 = y * (-2.92264f - y2 * 7.33447f);
412  t7 = y * -.56419f;
413  a8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
414  y2 * (y2 * (y2 * (y2 - 3.68288f) + 126.532f) - 955.194f)
415  + 9504.65f) - 70946.1f) + 483737.f) - 2857210.f) +
416  14464700.f) - 61114800.f) + 2.11107e8f) - 5.79099e8f) +
417  1.17022e9f) - 1.5599e9f) + 1.02827e9f;
418  b8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
419  y2 * (y2 * (y2 * 14.f - 40.5117f) + 533.254f) + 3058.26f)
420  - 55600.f) + 498334.f) - 2849540.f) + 13946500.f) -
421  70135800.f) + 2.89676e8f) - 7.53828e8f) + 1.66421e9f) -
422  2.28855e9f) + 1.5599e9f;
423  c8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
424  y2 * (y2 * 91 - 198.876f) - 1500.17f) + 48153.3f) -
425  217801.f) - 1063520.f) + 1.4841e7f) - 46039600.f) +
426  63349600.f) - 6.60078e8f) + 1.06002e9f) - 1.66421e9f) +
427  1.17022e9f;
428  d8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
429  y2 * 364 - 567.164f) - 16493.7f) + 161461.f) + 280428.f)
430  - 6890020.f) - 6876560.f) + 1.99846e8f) + 54036700.f) +
431  6.60078e8f) - 7.53828e8f) + 5.79099e8f;
432  e8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 *
433  1001 - 1012.79f) - 55582.f) + 240373.f) + 1954700.f) -
434  5257220.f) - 50101700.f) - 1.99846e8f) + 63349600.f) -
435  2.89676e8f) + 2.11107e8f;
436  f8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 2002 -
437  1093.82f) - 106663.f) + 123052.f) + 3043160.f) +
438  5257220.f) - 6876560.f) + 46039600.f) - 70135800.f) +
439  61114800.f;
440  g8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 3003 -
441  486.14f) - 131337.f) - 123052.f) + 1954700.f) + 6890020.f)
442  + 1.4841e7f) - 13946500.f) + 14464700.f;
443  h8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 3432 + 486.14f) -
444  106663.f) - 240373.f) + 280428.f) + 1063520.f) -
445  2849540.f) + 2857210.f;
446  o8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 3003 + 1093.82f) -
447  55582.f) - 161461.f) - 217801.f) - 498334.f) + 483737.f;
448  p8 = y2 * (y2 * (y2 * (y2 * (y2 * 2002 + 1012.79f) - 16493.7f) -
449  48153.3f) - 55600.f) + 70946.1f;
450  q8 = y2 * (y2 * (y2 * (y2 * 1001.f + 567.164f) - 1500.17f) -
451  3058.26f) + 9504.65f;
452  r8 = y2 * (y2 * (y2 * 364 + 198.876f) + 533.254f) + 955.194f;
453  s8 = y2 * (y2 * 91.f + 40.5117f) + 126.532f;
454  t8 = y2 * 14.f + 3.68288f;
455  }
456  ym2 = y * 2;
457  for (i2 = 1; i2 <= 3; i2 += 2) {
458  i__1 = lauf[((i2 + 1) << 2) - 1];
459  for (i1 = lauf[(i2 << 2) - 1]; i1 <= i__1; ++i1) {
460  x2 = x[(Index)i1-1] * x[(Index)i1-1];
461  ls[(Index)i1-1] = fac * (exp(y2 - x2) * cos(x[(Index)i1-1] * ym2) - (a7 + x2 *
462  (b7 + x2 * (c7 + x2 * (d7 + x2 * (e7 + x2 * (f7 + x2
463  * (g7 + x2 * (h7 + x2 * (o7 + x2 * (p7 + x2 * (q7 +
464  x2 * (r7 + x2 * (s7 + x2 * t7))))))))))))) / (a8 + x2
465  * (b8 + x2 * (c8 + x2 * (d8 + x2 * (e8 + x2 * (f8 +
466  x2 * (g8 + x2 * (h8 + x2 * (o8 + x2 * (p8 + x2 * (q8
467  + x2 * (r8 + x2 * (s8 + x2 * (t8 + x2)))))))))))))));
468  /* L4: */
469  }
470  }
471  }
472  /* ---- Region 3 */
473  /* -------------------------------------------------------------------- */
474  if (lauf[6] >= lauf[2] || lauf[14] >= lauf[10]) {
475  if ((r__1 = y - yps3, fabs(r__1)) > 1e-8f) {
476  yps3 = y;
477  a5 = y * (y * (y * (y * (y * (y * (y * (y * (y *
478  .564224f + 7.55895f) + 49.5213f) + 204.501f) + 581.746f)
479  + 1174.8f) + 1678.33f) + 1629.76f) + 973.778f) + 272.102f;
480  b5 = y * (y * (y * (y * (y * (y * (y * 2.25689f + 22.6778f)
481  + 100.705f) + 247.198f) + 336.364f) + 220.843f) -
482  2.34403f) - 60.5644f;
483  c5 = y * (y * (y * (y * (y * 3.38534f + 22.6798f) + 52.8454f)
484  + 42.5683f) + 18.546f) + 4.58029f;
485  d5 = y * (y * (y * 2.25689f + 7.56186f) + 1.66203f) - .128922f;
486  e5 = y * .564224f + 9.71457e-4f;
487  a6 = y * (y * (y * (y * (y * (y * (y * (y * (y * (y +
488  13.3988f) + 88.2674f) + 369.199f) + 1074.41f) + 2256.98f)
489  + 3447.63f) + 3764.97f) + 2802.87f) + 1280.83f) +
490  272.102f;
491  b6 = y * (y * (y * (y * (y * (y * (y * (y * 5.f +
492  53.5952f) + 266.299f) + 793.427f) + 1549.68f) + 2037.31f)
493  + 1758.34f) + 902.306f) + 211.678f;
494  c6 = y * (y * (y * (y * (y * (y * 10.f + 80.3928f) +
495  269.292f) + 479.258f) + 497.302f) + 308.186f) + 78.866f;
496  d6 = y * (y * (y * (y * 10.f + 53.5952f) + 92.7568f) +
497  55.0293f) + 22.0353f;
498  e6 = y * (y * 5.f + 13.3988f) + 1.49645f;
499  }
500  for (i2 = 1; i2 <= 3; i2 += 2) {
501  i__1 = lauf[((i2 + 1) << 2) - 2];
502  for (i1 = lauf[(i2 << 2) - 2]; i1 <= i__1; ++i1) {
503  x2 = x[(Index)i1-1] * x[(Index)i1-1];
504  ls[(Index)i1-1] = fac * (a5 + x2 * (b5 + x2 * (c5 + x2 * (d5 + x2 *
505  e5)))) / (a6 + x2 * (b6 + x2 * (c6 + x2 * (d6 + x2 * (
506  e6 + x2)))));
507  /* L5: */
508  }
509  }
510  }
511  /* ---- Region 2 */
512  /* -------------------------------------------------------------------- */
513  if (lauf[5] >= lauf[1] || lauf[13] >= lauf[9]) {
514  if ((r__1 = y - yps2, fabs(r__1)) > 1e-8f) {
515  yps2 = y;
516  a3 = y * (y2 * (y2 * (y2 * .56419f + 3.10304f) + 4.65456f) +
517  1.05786f);
518  b3 = y * (y2 * (y2 * 1.69257f + .56419f) + 2.962f);
519  c3 = y * (y2 * 1.69257f - 2.53885f);
520  d3 = y * .56419f;
521  a4 = y2 * (y2 * (y2 * (y2 + 6.f) + 10.5f) + 4.5f) + .5625f;
522  b4 = y2 * (y2 * (y2 * 4.f + 6.f) + 9.f) - 4.5f;
523  c4 = y2 * (y2 * 6.f - 6.f) + 10.5f;
524  d4 = y2 * 4.f - 6.f;
525  }
526  for (i2 = 1; i2 <= 3; i2 += 2) {
527  i__1 = lauf[((i2 + 1) << 2) - 3];
528  for (i1 = lauf[(i2 << 2) - 3]; i1 <= i__1; ++i1) {
529  x2 = x[(Index)i1-1] * x[(Index)i1-1];
530  ls[(Index)i1-1] = fac * (a3 + x2 * (b3 + x2 * (c3 + x2 * d3))) / (a4
531  + x2 * (b4 + x2 * (c4 + x2 * (d4 + x2))));
532  /* L6: */
533  }
534  }
535  }
536  /* ---- Region 1 */
537  /* -------------------------------------------------------------------- */
538  L7:
539  if (lauf[4] >= lauf[0] || lauf[12] >= lauf[8]) {
540  if ((r__1 = y - yps1, fabs(r__1)) > 1e-8f) {
541  yps1 = y;
542  a1 = y * .5641896f;
543  b1 = y2 + .5f;
544  a2 = y2 * 4;
545  }
546 
547  c1 = fac * a1;
548  for (i2 = 1; i2 <= 3; i2 += 2) {
549  i__1 = lauf[((i2 + 1) << 2) - 4];
550  for (i1 = lauf[(i2 << 2) - 4]; i1 <= i__1; ++i1) {
551  x2 = x[(Index)i1-1] * x[(Index)i1-1];
552  b2 = b1 - x2;
553  ls[(Index)i1-1] = c1 * (b1 + x2) / (b2 * b2 + a2 * x2);
554  /* L8: */
555  }
556  }
557  }
558 }
559 
560 
561 //---------------------------------------------------------------
562 // help function for lineshape_voigt_kuntz3
563 
564 long int bfun3_(Numeric y, Numeric x)
565 {
566  /* System generated locals */
567  long int ret_val;
568 
569  /* Local variables */
570  static Numeric x2, y2;
571 
572 /* -------------------------------------------------------------------- */
573  x2 = x * x;
574  y2 = y * y;
575  if (x2 * .4081676f + y2 > 21.159543f) {
576  if (x2 * .7019639f + y2 > 1123.14221f) {
577  ret_val = 0;
578  } else {
579  ret_val = 1;
580  }
581  } else {
582  if (x2 * .20753051f + y2 > 4.20249292f) {
583  ret_val = 2;
584  } else if (y >= fabs(x) * .08f - .12f) {
585  ret_val = 3;
586  } else {
587  ret_val = 4;
588  }
589  }
590  return ret_val;
591 } /* bfun3_ */
592 
593 
594 
655 void lineshape_voigt_kuntz3(Vector& ls,
656  Vector& x,
657  Numeric f0,
658  Numeric gamma,
659  Numeric sigma,
660  VectorView f_mono,
661  const Index nf)
662 
663 {
664  // FIXME: nf is actually redundant. Could be thrown out in the
665  // future. For now, let's do an assertion that at least it is
666  // correct:
667  assert( nf==f_mono.nelem() );
668 
669  // seems not necessary for Doppler correction
670  // extern const Numeric SQRT_NAT_LOG_2;
671 
672  // PI
673  extern const Numeric PI;
674 
675  // constant sqrt(1/pi)
676  const Numeric sqrt_invPI = sqrt(1/PI);
677 
678  // constant normalization factor for voigt
679  Numeric fac = 1.0 / sigma * sqrt_invPI;
680 
681  /* Initialized data */
682 
683  static Numeric yps0 = -1.0;
684  static Numeric yps1 = -1.0;
685  static Numeric yps2 = -1.0;
686  static Numeric yps3 = -1.0;
687  static Numeric yps4 = -1.0;
688 
689  /* System generated locals */
690  long i__1, i__2;
691  Numeric r__1;
692 
693  /* Local variables */
694  static long bmin, lauf[20] /* was [5][4] */, bmax, imin, imax;
695  static long stack[80] /* was [20][4] */;
696  static Numeric a0, a1, a2, a3, a4, a5, a6, a7, a8, b8, c8, d8, e8, f8, g8,
697  h8, b7, c7, d7, e7, f7, g7, o8, p8, q8, r8, s8, t8, h7, o7, p7,
698  q7, r7, s7, t7, b6, c6, d6, e6, b5, c5, d5, e5, b4, c4, d4, b3,
699  c3, d3, b1, y2;
700  static long i2, i1;
701  static Numeric x2, b2, b0, c1;
702  static long stackp, imitte;
703  static Numeric ym2;
704 
705  // variables needed in original c routine:
706 
707  // Ratio of the Lorentz halfwidth to the Doppler halfwidth
708  Numeric y = gamma / sigma;
709 
710  // frequency in units of Doppler
711  for (i1=0; i1< (int) nf; i1++)
712  {
713  x[(Index)i1] = (f_mono[(Index)i1] - f0) / sigma;
714  }
715 
716 
717  /* Parameter adjustments */
718  // this does not work for variables type Vector, corrected by
719  // adjusting the index to array ls and x
720  //--ls;
721  //--x;
722 
723  /* Function Body */
724  y2 = y * y;
725  if (y >= 23.0 || x[0] >= 39.0 || x[nf-1] <= -39.0) {
726  lauf[0] = 1;
727  lauf[5] = nf;
728  lauf[10] = nf;
729  lauf[15] = 0;
730  goto L8;
731  }
732  for (i2 = 1; i2 <= 4; ++i2) {
733  for (i1 = 0; i1 <= 4; ++i1) {
734  lauf[i1 + i2 * 5 - 5] = i2 % 2 * (nf + 1);
735  /* L1: */
736  }
737  }
738  stackp = 1;
739  stack[stackp - 1] = 1;
740  stack[stackp + 19] = nf;
741  stack[stackp + 39] = bfun3_(y, x[0]);
742  stack[stackp + 59] = bfun3_(y, x[nf-1]);
743  L2:
744  imin = stack[stackp - 1];
745  imax = stack[stackp + 19];
746  bmin = stack[stackp + 39];
747  bmax = stack[stackp + 59];
748  if (bmin == bmax) {
749  if (x[(Index)imax-1] < 0.f) {
750  /* Computing MIN */
751  i__1 = imin, i__2 = lauf[bmin];
752  lauf[bmin] = min(i__1,i__2);
753  /* Computing MAX */
754  i__1 = imax, i__2 = lauf[bmax + 5];
755  lauf[bmax + 5] = max(i__1,i__2);
756  --stackp;
757  goto L3;
758  } else if (x[(Index)imin-1] >= 0.f) {
759  /* Computing MIN */
760  i__1 = imin, i__2 = lauf[bmin + 10];
761  lauf[bmin + 10] = min(i__1,i__2);
762  /* Computing MAX */
763  i__1 = imax, i__2 = lauf[bmax + 15];
764  lauf[bmax + 15] = max(i__1,i__2);
765  --stackp;
766  goto L3;
767  }
768  }
769  imitte = (imax + imin) / 2;
770  stack[stackp - 1] = imitte + 1;
771  stack[stackp + 19] = imax;
772  stack[stackp + 39] = bfun3_(y, x[(Index)imitte]);
773  stack[stackp + 59] = bmax;
774  ++stackp;
775  stack[stackp - 1] = imin;
776  stack[stackp + 19] = imitte;
777  stack[stackp + 39] = bmin;
778  stack[stackp + 59] = bfun3_(y, x[(Index)imitte-1]);
779  L3:
780  if (stackp > 0) {
781  goto L2;
782  }
783  /* ---- Region 4 */
784  /* -------------------------------------------------------------------- */
785  if (lauf[9] >= lauf[4] || lauf[19] >= lauf[14]) {
786  if ((r__1 = y - yps4, fabs(r__1)) > 1e-8f) {
787  yps4 = y;
788  a7 = y * (y2 * (y2 * 4.56662e8f - 9.86604e8f) + 1.16028e9f);
789  b7 = y * (y2 * (y2 * 8.06985e8f - 9.85386e8f) - 5.60505e8f);
790  c7 = y * (y2 * (y2 * 2.94262e8f + 2.47157e8f) - 6.51523e8f);
791  d7 = y * (y2 * (2.70167e8f - y2 * 99622400.f) - 2.63894e8f);
792  e7 = y * (y2 * (y2 * 5569650.f + 1.40677e8f) - 63177100.f);
793  f7 = y * (y2 * (4073820.f - y2 * 33289600.f) - 16984600.f);
794  g7 = y * (y2 * (7528830.f - y2 * 900010) - 1231650.f);
795  h7 = y * (y2 * (y2 * 153468 + 86407.6f) - 610622.f);
796  o7 = y * (y2 * (y2 * 26538.5f + 49883.8f) - 23586.5f);
797  p7 = y * (y2 * (y2 * 953.655f + 2198.86f) - 8009.1f);
798  q7 = y * (y2 * (-271.202f - y2 * 134.792f) - 622.056f);
799  r7 = y * (y2 * (-29.7896f - y2 * 44.0068f) - 77.0535f);
800  s7 = y * (-2.92264f - y2 * 7.33447f);
801  t7 = y * -.56419f;
802  a8 = y2 * (y2 * 1.17022e9f - 1.5599e9f) + 1.02827e9f;
803  b8 = y2 * (y2 * 1.66421e9f - 2.28855e9f) + 1.5599e9f;
804  c8 = y2 * (y2 * 1.06002e9f - 1.66421e9f) + 1.17022e9f;
805  d8 = y2 * (y2 * 6.60078e8f - 7.53828e8f) + 5.79099e8f;
806  e8 = y2 * (y2 * 63349600.f - 2.89676e8f) + 2.11107e8f;
807  f8 = y2 * (y2 * 46039600.f - 70135800.f) + 61114800.f;
808  g8 = y2 * (y2 * 1.4841e7f - 13946500.f) + 14464700.f;
809  h8 = y2 * (y2 * 1063520.f - 2849540.f) + 2857210.f;
810  o8 = y2 * (-498334.f - y2 * 217801.f) + 483737.f;
811  p8 = y2 * (-55600.f - y2 * 48153.3f) + 70946.1f;
812  q8 = y2 * (-3058.26f - y2 * 1500.17f) + 9504.65f;
813  r8 = y2 * (y2 * 198.876f + 533.254f) + 955.194f;
814  s8 = y2 * (y2 * 91.f + 40.5117f) + 126.532f;
815  t8 = y2 * 14.f + 3.68288f;
816  }
817  ym2 = y * 2;
818  for (i2 = 1; i2 <= 3; i2 += 2) {
819  i__1 = lauf[(i2 + 1) * 5 - 1];
820  for (i1 = lauf[i2 * 5 - 1]; i1 <= i__1; ++i1) {
821  x2 = x[(Index)i1-1] * x[(Index)i1-1];
822  ls[(Index)i1-1] = fac * (exp(y2 - x2) * cos(x[(Index)i1-1] * ym2) - (a7 + x2 *
823  (b7 + x2 * (c7 + x2 * (d7 + x2 * (e7 + x2 * (f7 + x2
824  * (g7 + x2 * (h7 + x2 * (o7 + x2 * (p7 + x2 * (q7 +
825  x2 * (r7 + x2 * (s7 + x2 * t7))))))))))))) / (a8 + x2
826  * (b8 + x2 * (c8 + x2 * (d8 + x2 * (e8 + x2 * (f8 +
827  x2 * (g8 + x2 * (h8 + x2 * (o8 + x2 * (p8 + x2 * (q8
828  + x2 * (r8 + x2 * (s8 + x2 * (t8 + x2)))))))))))))));
829  /* L4: */
830  }
831  }
832  }
833  /* ---- Region 3 */
834  /* -------------------------------------------------------------------- */
835  if (lauf[8] >= lauf[3] || lauf[18] >= lauf[13]) {
836  if ((r__1 = y - yps3, fabs(r__1)) > 1e-8f) {
837  yps3 = y;
838  a5 = y * (y * (y * (y * (y * (y * (y * (y * (y *
839  .564224f + 7.55895f) + 49.5213f) + 204.501f) + 581.746f)
840  + 1174.8f) + 1678.33f) + 1629.76f) + 973.778f) + 272.102f;
841  b5 = y * (y * (y * (y * (y * (y * (y * 2.25689f + 22.6778f)
842  + 100.705f) + 247.198f) + 336.364f) + 220.843f) -
843  2.34403f) - 60.5644f;
844  c5 = y * (y * (y * (y * (y * 3.38534f + 22.6798f) + 52.8454f)
845  + 42.5683f) + 18.546f) + 4.58029f;
846  d5 = y * (y * (y * 2.25689f + 7.56186f) + 1.66203f) - .128922f;
847  e5 = y * .564224f + 9.71457e-4f;
848  a6 = y * (y * (y * (y * (y * (y * (y * (y * (y * (y +
849  13.3988f) + 88.2674f) + 369.199f) + 1074.41f) + 2256.98f)
850  + 3447.63f) + 3764.97f) + 2802.87f) + 1280.83f) +
851  272.102f;
852  b6 = y * (y * (y * (y * (y * (y * (y * (y * 5.f +
853  53.5952f) + 266.299f) + 793.427f) + 1549.68f) + 2037.31f)
854  + 1758.34f) + 902.306f) + 211.678f;
855  c6 = y * (y * (y * (y * (y * (y * 10.f + 80.3928f) +
856  269.292f) + 479.258f) + 497.302f) + 308.186f) + 78.866f;
857  d6 = y * (y * (y * (y * 10.f + 53.5952f) + 92.7568f) +
858  55.0293f) + 22.0353f;
859  e6 = y * (y * 5.f + 13.3988f) + 1.49645f;
860  }
861  for (i2 = 1; i2 <= 3; i2 += 2) {
862  i__1 = lauf[(i2 + 1) * 5 - 2];
863  for (i1 = lauf[i2 * 5 - 2]; i1 <= i__1; ++i1) {
864  x2 = x[(Index)i1-1] * x[(Index)i1-1];
865  ls[(Index)i1-1] = fac * (a5 + x2 * (b5 + x2 * (c5 + x2 * (d5 + x2 *
866  e5)))) / (a6 + x2 * (b6 + x2 * (c6 + x2 * (d6 + x2 * (
867  e6 + x2)))));
868  /* L5: */
869  }
870  }
871  }
872  /* ---- Region 2 */
873  /* -------------------------------------------------------------------- */
874  if (lauf[7] >= lauf[2] || lauf[17] >= lauf[12]) {
875  if ((r__1 = y - yps2, fabs(r__1)) > 1e-8f) {
876  yps2 = y;
877  a3 = y * (y2 * (y2 * (y2 * .56419f + 3.10304f) + 4.65456f) +
878  1.05786f);
879  b3 = y * (y2 * (y2 * 1.69257f + .56419f) + 2.962f);
880  c3 = y * (y2 * 1.69257f - 2.53885f);
881  d3 = y * .56419f;
882  a4 = y2 * (y2 * (y2 * (y2 + 6.f) + 10.5f) + 4.5f) + .5625f;
883  b4 = y2 * (y2 * (y2 * 4.f + 6.f) + 9.f) - 4.5f;
884  c4 = y2 * (y2 * 6.f - 6.f) + 10.5f;
885  d4 = y2 * 4.f - 6.f;
886  }
887  for (i2 = 1; i2 <= 3; i2 += 2) {
888  i__1 = lauf[(i2 + 1) * 5 - 3];
889  for (i1 = lauf[i2 * 5 - 3]; i1 <= i__1; ++i1) {
890  x2 = x[(Index)i1-1] * x[(Index)i1-1];
891  ls[(Index)i1-1] = fac * (a3 + x2 * (b3 + x2 * (c3 + x2 * d3))) / (a4
892  + x2 * (b4 + x2 * (c4 + x2 * (d4 + x2))));
893  /* L6: */
894  }
895  }
896  }
897  /* ---- Region 1 */
898  /* -------------------------------------------------------------------- */
899  if (lauf[6] >= lauf[1] || lauf[16] >= lauf[11]) {
900  if ((r__1 = y - yps1, fabs(r__1)) > 1e-8f) {
901  yps1 = y;
902  a1 = y * .5641896f;
903  b1 = y2 + .5f;
904  a2 = y2 * 4;
905  }
906 
907  c1 = a1 * fac;
908  for (i2 = 1; i2 <= 3; i2 += 2) {
909  i__1 = lauf[(i2 + 1) * 5 - 4];
910  for (i1 = lauf[i2 * 5 - 4]; i1 <= i__1; ++i1) {
911  x2 = x[(Index)i1-1] * x[(Index)i1-1];
912  b2 = b1 - x2;
913  ls[(Index)i1-1] = c1 * (b1 + x2) / (b2 * b2 + a2 * x2);
914  /* L7: */
915  }
916  }
917  }
918  /* ---- Region 0 (Lorentz) */
919  /* -------------------------------------------------------------------- */
920 L8:
921  if (lauf[5] >= lauf[0] || lauf[15] >= lauf[10]) {
922  if ((r__1 = y - yps0, fabs(r__1)) > 1e-8f) {
923  yps0 = y;
924  a0 = y * .5641896f;
925  }
926 
927  b0 = a0 * fac;
928  for (i2 = 1; i2 <= 3; i2 += 2) {
929  i__1 = lauf[(i2 + 1) * 5 - 5];
930  for (i1 = lauf[i2 * 5 - 5]; i1 <= i__1; ++i1) {
931  ls[(Index)i1-1] = b0 / (x[(Index)i1-1] * x[(Index)i1-1] + y2);
932  /* L9: */
933  }
934  }
935  }
936 }
937 
938 
939 //------------------------------------------------------------------
940 // help function for lineshape_voigt_kuntz4
941 
942 long bfun4_(Numeric y, Numeric x)
943 {
944  /* System generated locals */
945  long ret_val;
946 
947  /* Local variables */
948  static Numeric x2, y2;
949 
950  x2 = x * x;
951  y2 = y * y;
952  if (x2 * .0062f + y2 * .01417f > 1.f) {
953  if (x2 * 6.2e-5f + y2 * 1.98373e-4f > 1.f) {
954  ret_val = 0;
955  } else {
956  ret_val = 1;
957  }
958  } else {
959  if (x2 * .041649f + y2 * .111111111f > 1.f) {
960  ret_val = 2;
961  } else if (y >= fabs(x) * .19487f - .1753846f) {
962  ret_val = 3;
963  } else {
964  ret_val = 4;
965  }
966  }
967  return ret_val;
968 }
969 
1029 void lineshape_voigt_kuntz4(Vector& ls,
1030  Vector& x,
1031  Numeric f0,
1032  Numeric gamma,
1033  Numeric sigma,
1034  VectorView f_mono,
1035  const Index nf)
1036 {
1037  // FIXME: nf is actually redundant. Could be thrown out in the
1038  // future. For now, let's do an assertion that at least it is
1039  // correct:
1040  assert( nf==f_mono.nelem() );
1041 
1042  // seems not necessary for Doppler correction
1043  // extern const Numeric SQRT_NAT_LOG_2;
1044 
1045  // PI
1046  extern const Numeric PI;
1047 
1048  // constant sqrt(1/pi)
1049  const Numeric sqrt_invPI = sqrt(1/PI);
1050 
1051  // constant normalization factor for voigt
1052  Numeric fac = 1.0 / sigma * sqrt_invPI;
1053 
1054 
1055  /* Initialized data */
1056 
1057  static float yps0 = -1.f;
1058  static float yps1 = -1.f;
1059  static float yps2 = -1.f;
1060  static float yps3 = -1.f;
1061  static float yps4 = -1.f;
1062 
1063  /* System generated locals */
1064  long i__1, i__2;
1065  float r__1;
1066 
1067  /* Local variables */
1068  static long bmin, lauf[20] /* was [5][4] */, bmax, imin, imax;
1069  static long stack[80] /* was [20][4] */;
1070  static Numeric a0, a1, a2, a3, a4, a5, a6, a7, a8, b8, c8, d8, e8, f8, g8,
1071  h8, b7, c7, d7, e7, f7, g7, o8, p8, q8, r8, s8, t8, h7, o7, p7,
1072  q7, r7, s7, t7, b6, c6, d6, e6, b5, c5, d5, e5, b4, c4, d4, b3,
1073  c3, d3, b1, y2;
1074  static long i2, i1;
1075  static Numeric x2, b2, b0, c1;
1076  static long stackp, imitte;
1077  static Numeric ym2;
1078 
1079  // variables needed in original c routine:
1080 
1081  // Ratio of the Lorentz halfwidth to the Doppler halfwidth
1082  Numeric y = gamma / sigma;
1083 
1084  // frequency in units of Doppler
1085  for (i1=0; i1< (int) nf; i1++)
1086  {
1087  x[(Index)i1] = (f_mono[(Index)i1] - f0) / sigma;
1088  }
1089 
1090 
1091  /* Parameter adjustments */
1092  // this does not work for variables type Vector, corrected by
1093  // adjusting the index to array ls and x
1094  //--ls;
1095  //--x;
1096 
1097  /* Function Body */
1098  y2 = y * y;
1099  if (y >= 71.f || x[0] >= 123.f || x[nf-1] <= -123.f) {
1100  lauf[0] = 1;
1101  lauf[5] = nf;
1102  lauf[10] = nf;
1103  lauf[15] = 0;
1104  goto L8;
1105  }
1106  for (i2 = 1; i2 <= 4; ++i2) {
1107  for (i1 = 0; i1 <= 4; ++i1) {
1108  lauf[i1 + i2 * 5 - 5] = i2 % 2 * (nf + 1);
1109  /* L1: */
1110  }
1111  }
1112  stackp = 1;
1113  stack[stackp - 1] = 1;
1114  stack[stackp + 19] = nf;
1115  stack[stackp + 39] = bfun4_(y, x[0]);
1116  stack[stackp + 59] = bfun4_(y, x[nf-1]);
1117  L2:
1118  imin = stack[stackp - 1];
1119  imax = stack[stackp + 19];
1120  bmin = stack[stackp + 39];
1121  bmax = stack[stackp + 59];
1122  if (bmin == bmax) {
1123  if (x[(Index)imax-1] < 0.f) {
1124  /* Computing MIN */
1125  i__1 = imin, i__2 = lauf[bmin];
1126  lauf[bmin] = min(i__1,i__2);
1127  /* Computing MAX */
1128  i__1 = imax, i__2 = lauf[bmax + 5];
1129  lauf[bmax + 5] = max(i__1,i__2);
1130  --stackp;
1131  goto L3;
1132  } else if (x[(Index)imin-1] >= 0.f) {
1133  /* Computing MIN */
1134  i__1 = imin, i__2 = lauf[bmin + 10];
1135  lauf[bmin + 10] = min(i__1,i__2);
1136  /* Computing MAX */
1137  i__1 = imax, i__2 = lauf[bmax + 15];
1138  lauf[bmax + 15] = max(i__1,i__2);
1139  --stackp;
1140  goto L3;
1141  }
1142  }
1143  imitte = (imax + imin) / 2;
1144  stack[stackp - 1] = imitte + 1;
1145  stack[stackp + 19] = imax;
1146  stack[stackp + 39] = bfun4_(y, x[(Index)imitte]);
1147  stack[stackp + 59] = bmax;
1148  ++stackp;
1149  stack[stackp - 1] = imin;
1150  stack[stackp + 19] = imitte;
1151  stack[stackp + 39] = bmin;
1152  stack[stackp + 59] = bfun4_(y, x[(Index)imitte-1]);
1153  L3:
1154  if (stackp > 0) {
1155  goto L2;
1156  }
1157  /* ---- Region 4 */
1158  /* -------------------------------------------------------------------- */
1159  if (lauf[9] >= lauf[4] || lauf[19] >= lauf[14]) {
1160  if ((r__1 = (float)y - yps4, fabs(r__1)) > 1e-8f) {
1161  yps4 = (float)y;
1162  a7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
1163  y2 * (y2 * (y2 * (2.35944f - y2 * .56419f) - 72.9359f) +
1164  571.687f) - 5860.68f) + 40649.2f) - 320772.f) + 1684100.f)
1165  - 9694630.f) + 40816800.f) - 1.53575e8f) + 4.56662e8f) -
1166  9.86604e8f) + 1.16028e9f);
1167  b7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
1168  y2 * (y2 * (23.0312f - y2 * 7.33447f) - 234.143f) -
1169  2269.19f) + 45251.3f) - 234417.f) + 3599150.f) -
1170  7723590.f) + 86482900.f) - 2.91876e8f) + 8.06985e8f) -
1171  9.85386e8f) - 5.60505e8f);
1172  c7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
1173  y2 * (97.6203f - y2 * 44.0068f) + 1097.77f) - 25338.3f) +
1174  98079.1f) + 576054.f) - 2.3818e7f) + 22930200.f) -
1175  2.04467e8f) + 2.94262e8f) + 2.47157e8f) - 6.51523e8f);
1176  d7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
1177  228.563f - y2 * 161.358f) + 8381.97f) - 66431.2f) -
1178  303569.f) + 2240400.f) + 38311200.f) - 41501300.f) -
1179  99622400.f) + 2.70167e8f) - 2.63894e8f);
1180  e7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
1181  296.38f - y2 * 403.396f) + 23507.6f) - 66212.1f) -
1182  1.003e6f) + 468142.f) + 24620100.f) + 5569650.f) +
1183  1.40677e8f) - 63177100.f);
1184  f7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (125.591f -
1185  y2 * 726.113f) + 37544.8f) + 8820.94f) - 934717.f) -
1186  1931140.f) - 33289600.f) + 4073820.f) - 16984600.f);
1187  g7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (-260.198f - y2 *
1188  968.15f) + 37371.9f) + 79902.5f) - 186682.f) - 900010.f)
1189  + 7528830.f) - 1231650.f);
1190  h7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (-571.645f - y2 * 968.15f)
1191  + 23137.1f) + 72520.9f) + 153468.f) + 86407.6f) -
1192  610622.f);
1193  o7 = y * (y2 * (y2 * (y2 * (y2 * (-575.164f - y2 * 726.113f) +
1194  8073.15f) + 26538.5f) + 49883.8f) - 23586.5f);
1195  p7 = y * (y2 * (y2 * (y2 * (-352.467f - y2 * 403.396f) +
1196  953.655f) + 2198.86f) - 8009.1f);
1197  q7 = y * (y2 * (y2 * (-134.792f - y2 * 161.358f) - 271.202f) -
1198  622.056f);
1199  r7 = y * (y2 * (-29.7896f - y2 * 44.0068f) - 77.0535f);
1200  s7 = y * (-2.92264f - y2 * 7.33447f);
1201  t7 = y * -.56419f;
1202  a8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
1203  y2 * (y2 * (y2 * (y2 - 3.68288f) + 126.532f) - 955.194f)
1204  + 9504.65f) - 70946.1f) + 483737.f) - 2857210.f) +
1205  14464700.f) - 61114800.f) + 2.11107e8f) - 5.79099e8f) +
1206  1.17022e9f) - 1.5599e9f) + 1.02827e9f;
1207  b8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
1208  y2 * (y2 * (y2 * 14.f - 40.5117f) + 533.254f) + 3058.26f)
1209  - 55600.f) + 498334.f) - 2849540.f) + 13946500.f) -
1210  70135800.f) + 2.89676e8f) - 7.53828e8f) + 1.66421e9f) -
1211  2.28855e9f) + 1.5599e9f;
1212  c8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
1213  y2 * (y2 * 91 - 198.876f) - 1500.17f) + 48153.3f) -
1214  217801.f) - 1063520.f) + 1.4841e7f) - 46039600.f) +
1215  63349600.f) - 6.60078e8f) + 1.06002e9f) - 1.66421e9f) +
1216  1.17022e9f;
1217  d8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
1218  y2 * 364 - 567.164f) - 16493.7f) + 161461.f) + 280428.f)
1219  - 6890020.f) - 6876560.f) + 1.99846e8f) + 54036700.f) +
1220  6.60078e8f) - 7.53828e8f) + 5.79099e8f;
1221  e8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 *
1222  1001 - 1012.79f) - 55582.f) + 240373.f) + 1954700.f) -
1223  5257220.f) - 50101700.f) - 1.99846e8f) + 63349600.f) -
1224  2.89676e8f) + 2.11107e8f;
1225  f8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 2002 -
1226  1093.82f) - 106663.f) + 123052.f) + 3043160.f) +
1227  5257220.f) - 6876560.f) + 46039600.f) - 70135800.f) +
1228  61114800.f;
1229  g8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 3003 -
1230  486.14f) - 131337.f) - 123052.f) + 1954700.f) + 6890020.f)
1231  + 1.4841e7f) - 13946500.f) + 14464700.f;
1232  h8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 3432 + 486.14f) -
1233  106663.f) - 240373.f) + 280428.f) + 1063520.f) -
1234  2849540.f) + 2857210.f;
1235  o8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 3003 + 1093.82f) -
1236  55582.f) - 161461.f) - 217801.f) - 498334.f) + 483737.f;
1237  p8 = y2 * (y2 * (y2 * (y2 * (y2 * 2002 + 1012.79f) - 16493.7f) -
1238  48153.3f) - 55600.f) + 70946.1f;
1239  q8 = y2 * (y2 * (y2 * (y2 * 1001.f + 567.164f) - 1500.17f) -
1240  3058.26f) + 9504.65f;
1241  r8 = y2 * (y2 * (y2 * 364 + 198.876f) + 533.254f) + 955.194f;
1242  s8 = y2 * (y2 * 91.f + 40.5117f) + 126.532f;
1243  t8 = y2 * 14.f + 3.68288f;
1244  }
1245  ym2 = y * 2;
1246  for (i2 = 1; i2 <= 3; i2 += 2) {
1247  i__1 = lauf[(i2 + 1) * 5 - 1];
1248  for (i1 = lauf[i2 * 5 - 1]; i1 <= i__1; ++i1) {
1249  x2 = x[(Index)i1-1] * x[(Index)i1-1];
1250  ls[(Index)i1-1] = fac * (exp(y2 - x2) * cos(x[(Index)i1-1] * ym2) - (a7 + x2 *
1251  (b7 + x2 * (c7 + x2 * (d7 + x2 * (e7 + x2 * (f7 + x2
1252  * (g7 + x2 * (h7 + x2 * (o7 + x2 * (p7 + x2 * (q7 +
1253  x2 * (r7 + x2 * (s7 + x2 * t7))))))))))))) / (a8 + x2
1254  * (b8 + x2 * (c8 + x2 * (d8 + x2 * (e8 + x2 * (f8 +
1255  x2 * (g8 + x2 * (h8 + x2 * (o8 + x2 * (p8 + x2 * (q8
1256  + x2 * (r8 + x2 * (s8 + x2 * (t8 + x2)))))))))))))));
1257  /* L4: */
1258  }
1259  }
1260  }
1261  /* ---- Region 3 */
1262  /* -------------------------------------------------------------------- */
1263  if (lauf[8] >= lauf[3] || lauf[18] >= lauf[13]) {
1264  if ((r__1 = (float)y - yps3, fabs(r__1)) > 1e-8f) {
1265  yps3 = (float)y;
1266  a5 = y * (y * (y * (y * (y * (y * (y * (y * (y *
1267  .564224f + 7.55895f) + 49.5213f) + 204.501f) + 581.746f)
1268  + 1174.8f) + 1678.33f) + 1629.76f) + 973.778f) + 272.102f;
1269  b5 = y * (y * (y * (y * (y * (y * (y * 2.25689f + 22.6778f)
1270  + 100.705f) + 247.198f) + 336.364f) + 220.843f) -
1271  2.34403f) - 60.5644f;
1272  c5 = y * (y * (y * (y * (y * 3.38534f + 22.6798f) + 52.8454f)
1273  + 42.5683f) + 18.546f) + 4.58029f;
1274  d5 = y * (y * (y * 2.25689f + 7.56186f) + 1.66203f) - .128922f;
1275  e5 = y * .564224f + 9.71457e-4f;
1276  a6 = y * (y * (y * (y * (y * (y * (y * (y * (y * (y +
1277  13.3988f) + 88.2674f) + 369.199f) + 1074.41f) + 2256.98f)
1278  + 3447.63f) + 3764.97f) + 2802.87f) + 1280.83f) +
1279  272.102f;
1280  b6 = y * (y * (y * (y * (y * (y * (y * (y * 5.f +
1281  53.5952f) + 266.299f) + 793.427f) + 1549.68f) + 2037.31f)
1282  + 1758.34f) + 902.306f) + 211.678f;
1283  c6 = y * (y * (y * (y * (y * (y * 10.f + 80.3928f) +
1284  269.292f) + 479.258f) + 497.302f) + 308.186f) + 78.866f;
1285  d6 = y * (y * (y * (y * 10.f + 53.5952f) + 92.7568f) +
1286  55.0293f) + 22.0353f;
1287  e6 = y * (y * 5.f + 13.3988f) + 1.49645f;
1288  }
1289  for (i2 = 1; i2 <= 3; i2 += 2) {
1290  i__1 = lauf[(i2 + 1) * 5 - 2];
1291  for (i1 = lauf[i2 * 5 - 2]; i1 <= i__1; ++i1) {
1292  x2 = x[(Index)i1-1] * x[(Index)i1-1];
1293  ls[(Index)i1-1] = fac * (a5 + x2 * (b5 + x2 * (c5 + x2 * (d5 + x2 *
1294  e5)))) / (a6 + x2 * (b6 + x2 * (c6 + x2 * (d6 + x2 * (
1295  e6 + x2)))));
1296  /* L5: */
1297  }
1298  }
1299  }
1300  /* ---- Region 2 */
1301  /* -------------------------------------------------------------------- */
1302  if (lauf[7] >= lauf[2] || lauf[17] >= lauf[12]) {
1303  if ((r__1 = (float)y - yps2, fabs(r__1)) > 1e-8f) {
1304  yps2 = (float)y;
1305  a3 = y * (y2 * (y2 * (y2 * .56419f + 3.10304f) + 4.65456f) +
1306  1.05786f);
1307  b3 = y * (y2 * (y2 * 1.69257f + .56419f) + 2.962f);
1308  c3 = y * (y2 * 1.69257f - 2.53885f);
1309  d3 = y * .56419f;
1310  a4 = y2 * (y2 * (y2 * (y2 + 6.f) + 10.5f) + 4.5f) + .5625f;
1311  b4 = y2 * (y2 * (y2 * 4.f + 6.f) + 9.f) - 4.5f;
1312  c4 = y2 * (y2 * 6.f - 6.f) + 10.5f;
1313  d4 = y2 * 4.f - 6.f;
1314  }
1315  for (i2 = 1; i2 <= 3; i2 += 2) {
1316  i__1 = lauf[(i2 + 1) * 5 - 3];
1317  for (i1 = lauf[i2 * 5 - 3]; i1 <= i__1; ++i1) {
1318  x2 = x[(Index)i1-1] * x[(Index)i1-1];
1319  ls[(Index)i1-1] = fac * (a3 + x2 * (b3 + x2 * (c3 + x2 * d3))) / (a4
1320  + x2 * (b4 + x2 * (c4 + x2 * (d4 + x2))));
1321  /* L6: */
1322  }
1323  }
1324  }
1325  /* ---- Region 1 */
1326  /* -------------------------------------------------------------------- */
1327  if (lauf[6] >= lauf[1] || lauf[16] >= lauf[11]) {
1328  if ((r__1 = (float)y - yps1, fabs(r__1)) > 1e-8f) {
1329  yps1 = (float)y;
1330  a1 = y * .5641896f;
1331  b1 = y2 + .5f;
1332  a2 = y2 * 4;
1333  }
1334 
1335  c1 = a1 * fac;
1336  for (i2 = 1; i2 <= 3; i2 += 2) {
1337  i__1 = lauf[(i2 + 1) * 5 - 4];
1338  for (i1 = lauf[i2 * 5 - 4]; i1 <= i__1; ++i1) {
1339  x2 = x[(Index)i1-1] * x[(Index)i1-1];
1340  b2 = b1 - x2;
1341  ls[(Index)i1-1] = c1 * (b1 + x2) / (b2 * b2 + a2 * x2);
1342  /* L7: */
1343  }
1344  }
1345  }
1346  /* ---- Region 0 (Lorentz) */
1347  /* -------------------------------------------------------------------- */
1348  L8:
1349  if (lauf[5] >= lauf[0] || lauf[15] >= lauf[10]) {
1350  if ((r__1 = (float)y - yps0, fabs(r__1)) > 1e-8f) {
1351  yps0 = (float)y;
1352  a0 = y * .5641896f;
1353  }
1354 
1355  b0 = a0 * fac;
1356  for (i2 = 1; i2 <= 3; i2 += 2) {
1357  i__1 = lauf[(i2 + 1) * 5 - 5];
1358  for (i1 = lauf[i2 * 5 - 5]; i1 <= i__1; ++i1) {
1359  ls[(Index)i1-1] = b0 / (x[(Index)i1-1] * x[(Index)i1-1] + y2);
1360  /* L9: */
1361  }
1362  }
1363  }
1364 }
1365 
1366 
1367 
1412 /*** ROUTINE COMPUTES THE VOIGT FUNCTION Y/PI*INTEGRAL FROM ***/
1413 /*** - TO + INFINITY OF EXP(-T*T)/(Y*Y+(X-T)*(X-T)) DT ***/
1414 void lineshape_voigt_drayson(Vector& ls,
1415  Vector& x,
1416  Numeric f0,
1417  Numeric gamma,
1418  Numeric sigma,
1419  VectorView f_mono,
1420  const Index nf)
1421 
1422 {
1423  // FIXME: nf is actually redundant. Could be thrown out in the
1424  // future. For now, let's do an assertion that at least it is
1425  // correct:
1426  assert( nf==f_mono.nelem() );
1427 
1428  // seems not necessary for Doppler correction
1429  // extern const Numeric SQRT_NAT_LOG_2;
1430 
1431  // PI
1432  extern const Numeric PI;
1433 
1434  // constant sqrt(1/pi)
1435  const Numeric sqrt_invPI = sqrt(1/PI);
1436 
1437  // constant normalization factor for voigt
1438  Numeric fac = 1.0 / sigma * sqrt_invPI;
1439 
1440  static Numeric B[22+1] = {0.,0.,.7093602e-7};
1441  static Numeric RI[15+1];
1442  const Numeric XN[15+1] = {0.,10.,9.,8.,8.,7.,6.,5.,4.,3.,3.,3.,3.,3.,3.,3.};
1443  const Numeric YN[15+1] = {0.,.6,.6,.6,.5,.4,.4,.3,.3,.3,.3,1.,.9,.8,.7,.7};
1444  static Numeric D0[25+1], D1[25+1], D2[25+1], D3[25+1], D4[25+1];
1445  static Numeric HN[25+1];
1446  static Numeric H = .201;
1447  const Numeric XX[3+1] = {0.,.5246476,1.65068,.7071068};
1448  const Numeric HH[3+1] = {0.,.2562121,.2588268e-1,.2820948};
1449  const Numeric NBY2[19+1] = {0.,9.5,9.,8.5,8.,7.5,7.,6.5
1450  ,6.,5.5,5.,4.5,4.,3.5,3.,2.5,2.,1.5,1.,.5};
1451  const Numeric C[21+1] = {0.,.7093602e-7,-.2518434e-6,.856687e-6,
1452  -.2787638e-5,.866074e-5,-.2565551e-4,.7228775e-4,
1453  -.1933631e-3,.4899520e-3,-.1173267e-2,.2648762e-2,
1454  -.5623190e-2,.1119601e-1,-.2084976e-1,.3621573e-1,
1455  -.5851412e-1,.8770816e-1,-.121664,.15584,-.184,.2};
1456  static Numeric CO;
1457  Numeric U, V, UU, VV, Y2, DX;
1458  int I, J, K, MAX, MIN, N, i1;
1459  static int TRU = 0;
1460 
1461 
1462  // variables needed in original c routine:
1463 
1464  // Ratio of the Lorentz halfwidth to the Doppler halfwidth
1465  Numeric Y = gamma / sigma;
1466 
1467  // frequency in units of Doppler, Note: Drayson accepts only
1468  // positive values
1469  for (i1=0; i1< (int) nf; i1++)
1470  {
1471  x[i1] = fabs( (f_mono[i1] - f0) )/ sigma;
1472  }
1473 
1474 
1475  if (TRU) goto L104;
1476 /*** REGION I. COMPUTE DAWSON'S FUNCTION AT MESH POINTS ***/
1477  TRU = 1;
1478  for (I=1; I<=15; I++)
1479  RI[I] = -I/2.;
1480  for (I=1; I<=25; I++)
1481  {
1482  HN[I] = H*(I-.5);
1483  CO = 4.*HN[I]*HN[I]/25.-2.;
1484  for (J=2; J<=21; J++)
1485  B[J+1] = CO*B[J]-B[J-1]+C[J];
1486  D0[I] = HN[I]*(B[22]-B[21])/5.;
1487  D1[I] = 1.-2.*HN[I]*D0[I];
1488  D2[I] = (HN[I]*D1[I]+D0[I])/RI[2];
1489  D3[I] = (HN[I]*D2[I]+D1[I])/RI[3];
1490  D4[I] = (HN[I]*D3[I]+D2[I])/RI[4];
1491  }
1492 L104: for (K=0; K<(int) nf; K++)
1493  {
1494  if ((x[K]-5.) < 0.) goto L105; else goto L112;
1495  L105: if ((Y-1.) <= 0.) goto L110; else goto L106;
1496  L106: if (x[K] > 1.85*(3.6-Y)) goto L112;
1497  /*** REGION II: CONTINUED FRACTION. COMPUTE NUMBER OF TERMS NEEDED. ***/
1498  if (Y < 1.45) goto L107;
1499  I = (int) (Y+Y);
1500  goto L108;
1501  L107: I = (int) (11.*Y);
1502  L108: J = (int) (x[K]+x[K]+1.85);
1503  MAX = (int) (XN[J]*YN[I]+.46);
1504  MIN = (int) ( (16 < 21-2*MAX) ? 16 : 21-2*MAX );
1505  /*** EVALUATE CONTINUED FRACTION ***/
1506  UU = Y;
1507  VV = x[K];
1508  for (J=MIN; J <= 19; J++)
1509  {
1510  U = NBY2[J]/(UU*UU+VV*VV);
1511  UU = Y+U*UU;
1512  VV = x[K]-U*VV;
1513  }
1514  ls[K] = UU/(UU*UU+VV*VV)/1.772454*fac;
1515  continue;
1516  L110: Y2 = Y*Y;
1517  if (x[K]+Y >= 5.) goto L113;
1518  /*** REGION I. COMPUTE DAWSON'S FUNCTION AT X FROM TAYLOR SERIES. ***/
1519  N = (int) (x[K]/H);
1520  DX = x[K]-HN[N+1];
1521  U = (((D4[N+1]*DX+D3[N+1])*DX+D2[N+1])*DX+D1[N+1])*DX+D0[N+1];
1522  V = 1.-2.*x[K]*U;
1523  /*** TAYLOR SERIES EXPANSION ABOUT Y=0.0 ***/
1524  VV = exp(Y2-x[K]*x[K])*cos(2.*x[K]*Y)/1.128379-Y*V;
1525  UU = -Y;
1526  MAX = (int) (5.+(12.5-x[K])*.8*Y);
1527  for (I=2; I<=MAX; I+=2)
1528  {
1529  U = (x[K]*V+U)/RI[I];
1530  V = (x[K]*U+V)/RI[I+1];
1531  UU = -UU*Y2;
1532  VV = VV+V*UU;
1533  }
1534  ls[K] = 1.128379*VV*fac;
1535  continue;
1536  L112: Y2 = Y*Y;
1537  if (Y < 11.-.6875*x[K]) goto L113;
1538  /*** REGION IIIB: 2-POINT GAUSS-HERMITE QUADRATURE. ***/
1539  U = x[K]-XX[3];
1540  V = x[K]+XX[3];
1541  ls[K] = Y*(HH[3]/(Y2+U*U)+HH[3]/(Y2+V*V))*fac;
1542  continue;
1543  /*** REGION IIIA: 4-POINT GAUSS-HERMITE QUADRATURE. ***/
1544  L113: U = x[K]-XX[1];
1545  V = x[K]+XX[1];
1546  UU = x[K]-XX[2];
1547  VV = x[K]+XX[2];
1548  ls[K] = Y*(HH[1]/(Y2+U*U)+HH[1]/(Y2+V*V)+HH[2]/(Y2+UU*UU)+HH[2]/(Y2+
1549  VV*VV))*fac;
1550  continue;
1551  }
1552 }
1553 
1554 
1555 
1580 void lineshape_rosenkranz_voigt_kuntz6(Vector& ls,
1581  Vector& x,
1582  Numeric f0,
1583  Numeric gamma,
1584  Numeric sigma,
1585  VectorView f_mono,
1586  const Index nf)
1587 {
1588  // FIXME: nf is actually redundant. Could be thrown out in the
1589  // future. For now, let's do an assertion that at least it is
1590  // correct:
1591  assert( nf==f_mono.nelem() );
1592 
1593  // seems not necessary for Doppler correction
1594  // extern const Numeric SQRT_NAT_LOG_2;
1595 
1596  extern const Numeric PI;
1597 
1598 
1599  // calculate the required stuff with the parameter of the x array:
1600  Numeric overlap; // overlap correction
1601  Numeric gamma_o2; // linewidth
1602  Numeric gamma_o2_2; // square of linewidth
1603  {
1604  // x[0] = theta
1605  // x[1] = theta_Nair
1606  // x[2] = total pressure
1607  // x[3] = O2 VMR
1608  // x[4] = H2O VMR
1609  // x[5] = l_l.Agam()
1610  // x[6] = l_l.Nair()
1611  // x[7] = l_l.Aux()[0] = overlap y
1612  // x[8] = l_l.Aux()[1] = overlap v
1613 
1614  // Pressure Broadening:
1615  // pressure broadening is calculated differently for oxygen, not
1616  // with the partial pressure, but with the dry pressure. the
1617  // passed gamma parameter is calculated with the partial pressure,
1618  // so we correct it here:
1619  gamma *= ( 1.0 - x[4] ) / ( 1 - x[3] );
1620 
1621  // Overlap:
1622  // Get the overlap correction for oxygen, this is always calculated
1623  // even if not needed. FIXME: I assume that theta is to the power of
1624  // Nair (just as in the pressure broadening), otherwise we need
1625  // another parameter to pass with the catalogue.
1626  //cout << "x7, x8: " << x[7] << " " << x[8] << "\n";
1627  overlap = (x[7] + x[8] * ( x[0] - 1.0)) *
1628  x[2] * x[1];
1629 
1630  // water vapor impact:
1631  // impact of water vapor upon the oxygen pressure broadening.
1632  gamma_o2 = gamma + 1.1 * x[5] * x[2] * x[4] * x[0];
1633  gamma_o2_2 = gamma_o2 * gamma_o2;
1634  }
1635 
1636  // use voigt lineshape for high altitudes, otherwise lorentz with
1637  // overlap correction
1638  if ( (gamma_o2/sigma - 40) <= 0.0 )
1639  {
1640  /* call voigt function */
1641  lineshape_voigt_kuntz6(ls,
1642  x,
1643  f0,
1644  gamma_o2,
1645  sigma,
1646  f_mono,
1647  nf);
1648  }
1649  else
1650  {
1651  Numeric FD, FS, SF1, SF2;
1652  for (Index J=0; J<(Index) nf; J++)
1653  {
1654  FD = f_mono[J] - f0;
1655  FS = f_mono[J] + f0;
1656  SF1 = (gamma_o2 + FD*overlap) / (FD*FD + gamma_o2_2);
1657  SF2 = (gamma_o2 - FS*overlap) / (FS*FS + gamma_o2_2);
1658  ls[J] = (SF1 + SF2) / PI;
1659  }
1660  }
1661 }
1662 
1663 
1688 void lineshape_rosenkranz_voigt_drayson(Vector& ls,
1689  Vector& x,
1690  Numeric f0,
1691  Numeric gamma,
1692  Numeric sigma,
1693  VectorView f_mono,
1694  const Index nf)
1695 {
1696  // FIXME: nf is actually redundant. Could be thrown out in the
1697  // future. For now, let's do an assertion that at least it is
1698  // correct:
1699  assert( nf==f_mono.nelem() );
1700 
1701  // seems not necessary for Doppler correction
1702  // extern const Numeric SQRT_NAT_LOG_2;
1703 
1704  extern const Numeric PI;
1705 
1706 
1707  // calculate the required stuff with the parameter of the x array:
1708  Numeric overlap; // overlap correction
1709  Numeric gamma_o2; // linewidth
1710  Numeric gamma_o2_2; // square of linewidth
1711  {
1712  // x[0] = theta
1713  // x[1] = theta_Nair
1714  // x[2] = total pressure
1715  // x[3] = O2 VMR
1716  // x[4] = H2O VMR
1717  // x[5] = l_l.Agam()
1718  // x[6] = l_l.Nair()
1719  // x[7] = l_l.Aux()[0] = overlap y
1720  // x[8] = l_l.Aux()[1] = overlap v
1721 
1722  // Pressure Broadening:
1723  // pressure broadening is calculated differently for oxygen, not
1724  // with the partial pressure, but with the dry pressure. the
1725  // passed gamma parameter is calculated with the partial pressure,
1726  // so we correct it here:
1727  gamma *= ( 1.0 - x[4] ) / ( 1 - x[3] );
1728 
1729  // Overlap:
1730  // Get the overlap correction for oxygen, this is always calculated
1731  // even if not needed. FIXME: I assume that theta is to the power of
1732  // Nair (just as in the pressure broadening), otherwise we need
1733  // another parameter to pass with the catalogue.
1734  overlap = (x[7] + x[8] * ( x[0] - 1.0)) *
1735  x[2] * x[1];
1736 
1737  // water vapor impact:
1738  // impact of water vapor upon the oxygen pressure broadening.
1739  gamma_o2 = gamma + 1.1 * x[5] * x[2] * x[4] * x[0];
1740  gamma_o2_2 = gamma_o2 * gamma_o2;
1741  }
1742 
1743  // use voigt lineshape for high altitudes, otherwise lorentz with
1744  // overlap correction
1745  if ( (gamma_o2/sigma - 40) <= 0.0 )
1746  {
1747  /* call voigt function */
1748  lineshape_voigt_drayson(ls,
1749  x,
1750  f0,
1751  gamma_o2,
1752  sigma,
1753  f_mono,
1754  nf);
1755 
1756  // the controlfile uses generally a (f/f0)^2 factor for the
1757  // lineshape function, but we have to remove this normalization
1758  // factor (f/f0)^2 over here, because Rosenkranz does not use it
1759  // for the Voigt part, but for the Lorentz part. don't ask me
1760  // why...
1761  // FIXME: Clarify whether this is correct and if yes use the
1762  // normalization used in controlfile
1763  Numeric f0_2 = f0 * f0;
1764  for (Index J=0; J<(Index) nf; J++)
1765  {
1766  ls[J] *= f0_2 / (f_mono[J] * f_mono[J]);
1767  }
1768  }
1769  else
1770  {
1771  Numeric FD, FS, SF1, SF2;
1772  for (Index J=0; J<(Index) nf; J++)
1773  {
1774  FD = f_mono[J] - f0;
1775  FS = f_mono[J] + f0;
1776  SF1 = (gamma_o2 + FD*overlap) / (FD*FD + gamma_o2_2);
1777  SF2 = (gamma_o2 - FS*overlap) / (FS*FS + gamma_o2_2);
1778  ls[J] = (SF1 + SF2) / PI;
1779  }
1780  }
1781 }
1782 
1783 
1784 
1796 void chi_cousin(
1797  Numeric& chi,
1798  const Numeric& df )
1799 {
1800  // Conversion factor from Hz to cm-1
1801  extern const Numeric HZ2CM;
1802 
1803  const Numeric n2 = 0.79;
1804  const Numeric o2 = 0.21;
1805 
1806  const Numeric df_cm = df * HZ2CM;
1807  const Numeric df_cm_abs = fabs( df_cm );
1808 
1809  chi = 0;
1810 
1811  // N2 term
1812  if( df_cm_abs <= 5 )
1813  { chi += n2; }
1814  else if( df_cm_abs <= 22 )
1815  { chi += n2 * 1.968 * exp( -0.1354 * df_cm_abs ); }
1816  else if( df_cm_abs <= 50 )
1817  { chi += n2 * 0.160 * exp( -0.0214 * df_cm_abs ); }
1818  else
1819  { chi += n2 * 0.162 * exp( -0.0216 * df_cm_abs ); }
1820 
1821  // O2 term
1822  if( df_cm_abs <= 5 )
1823  { chi += o2; }
1824  else if( df_cm_abs <= 22 )
1825  { chi += o2 * 1.968 * exp( -0.1354 * df_cm_abs ); }
1826  else if( df_cm_abs <= 50 )
1827  { chi += o2 * 0.160 * exp( -0.0214 * df_cm_abs ); }
1828  else
1829  { chi += o2 * 0.162 * exp( -0.0216 * df_cm_abs ); }
1830 }
1831 
1832 
1833 
1846 void lineshape_CO2_lorentz(
1847  Vector& ls,
1848  Vector& /* X */,
1849  Numeric f0,
1850  Numeric gamma,
1851  Numeric /* sigma */,
1852  VectorView f_mono,
1853  const Index nf )
1854 {
1855  assert( f_mono.nelem() == nf );
1856 
1857  // PI:
1858  extern const Numeric PI;
1859 
1860  // Some constant variables
1861  const Numeric gamma2 = gamma * gamma;
1862  const Numeric fac = gamma/PI;
1863 
1864  for ( Index i=0; i<nf; ++i )
1865  {
1866  // f-f0
1867  const Numeric df = f_mono[i] - f0;
1868 
1869  // The chi factor
1870  Numeric chi;
1871  chi_cousin( chi, df );
1872 
1873  // chi * Lorentz
1874  ls[i] = chi * fac / ( df * df + gamma2 );
1875  }
1876 }
1877 
1878 
1879 
1892 void lineshape_CO2_drayson(
1893  Vector& ls,
1894  Vector& X,
1895  Numeric f0,
1896  Numeric gamma,
1897  Numeric sigma,
1898  VectorView f_mono,
1899  const Index nf )
1900 {
1901  lineshape_voigt_drayson( ls, X, f0, gamma, sigma, f_mono, nf );
1902 
1903  for ( Index i=0; i<nf; ++i )
1904  {
1905  // f-f0
1906  const Numeric df = f_mono[i] - f0;
1907 
1908  // The chi factor
1909  Numeric chi;
1910  chi_cousin( chi, df );
1911 
1912  // chi * Lorentz
1913  ls[i] *= chi;
1914  }
1915 }
1916 
1917 
1918 
1919 
1920 
1921 
1922 //------------------------------------------------------------------------
1923 // Normalization Functions
1924 //------------------------------------------------------------------------
1925 
1926 
1936 void lineshape_norm_no_norm(Vector& fac,
1937  Numeric /* f0 */,
1938  VectorView f_mono,
1939  const Numeric /* T */,
1940  const Index nf)
1941 {
1942  // FIXME: nf is actually redundant. Could be thrown out in the
1943  // future. For now, let's do an assertion that at least it is
1944  // correct:
1945  assert( nf==f_mono.nelem() );
1946 
1947  for ( Index i=0; i<nf; ++i )
1948  {
1949  fac[i] = 1.0;
1950  }
1951 }
1952 
1953 
1954 
1964 void lineshape_norm_linear(Vector& fac,
1965  Numeric f0,
1966  VectorView f_mono,
1967  const Numeric /* T */,
1968  const Index nf)
1969 {
1970  // FIXME: nf is actually redundant. Could be thrown out in the
1971  // future. For now, let's do an assertion that at least it is
1972  // correct:
1973  assert( nf==f_mono.nelem() );
1974 
1975  for ( Index i=0; i<nf; ++i )
1976  {
1977  fac[i] = f_mono[i] / f0;
1978  }
1979 }
1980 
1991 void lineshape_norm_quadratic(Vector& fac,
1992  Numeric f0,
1993  VectorView f_mono,
1994  const Numeric /* T */,
1995  const Index nf)
1996 {
1997  // FIXME: nf is actually redundant. Could be thrown out in the
1998  // future. For now, let's do an assertion that at least it is
1999  // correct:
2000  assert( nf==f_mono.nelem() );
2001 
2002  // don't do this for the whole loop
2003  Numeric f0_2 = f0 * f0;
2004 
2005  for ( Index i=0; i<nf; ++i )
2006  {
2007  fac[i] = (f_mono[i] * f_mono[i]) / f0_2;
2008  }
2009 }
2010 
2027 void lineshape_norm_VVH(Vector& fac,
2028  Numeric f0,
2029  VectorView f_mono,
2030  const Numeric T,
2031  const Index nf)
2032 {
2033  extern const Numeric PLANCK_CONST;
2034  extern const Numeric BOLTZMAN_CONST;
2035 
2036  // FIXME: nf is actually redundant. Could be thrown out in the
2037  // future. For now, let's do an assertion that at least it is
2038  // correct:
2039  assert( nf==f_mono.nelem() );
2040 
2041  // 2kT is constant for the loop
2042  Numeric kT = 2.0 * BOLTZMAN_CONST * T;
2043 
2044  // denominator is constant for the loop
2045  Numeric denom = f0 * tanh( PLANCK_CONST * f0 / kT );
2046 
2047  for ( Index i=0; i<nf; ++i )
2048  {
2049  fac[i] = f_mono[i] * tanh( PLANCK_CONST * f_mono[i] / kT ) /
2050  denom;
2051  }
2052 }
2053 
2054 
2055 
2056 
2057 //------------------------------------------------------------------------
2058 // Available Lineshapes
2059 //------------------------------------------------------------------------
2060 
2062 Array<LineshapeRecord> lineshape_data;
2063 
2064 void define_lineshape_data()
2065 {
2066  // Initialize to empty, just in case.
2067  lineshape_data.resize(0);
2068 
2069  lineshape_data.push_back
2070  (LineshapeRecord
2071  ("no_shape",
2072  "This lineshape does nothing. It only exists, because formally\n"
2073  "you have to specify a lineshape also for continuum tags.",
2074  lineshape_no_shape));
2075 
2076  lineshape_data.push_back
2077  (LineshapeRecord
2078  ("Lorentz",
2079  "The Lorentz line shape.",
2080  lineshape_lorentz));
2081 
2082  lineshape_data.push_back
2083  (LineshapeRecord
2084  ("Doppler",
2085  "The Doppler line shape.",
2086  lineshape_doppler));
2087 
2088  lineshape_data.push_back
2089  (LineshapeRecord
2090  ("Voigt_Kuntz6",
2091  "The Voigt line shape. Approximation by Kuntz: Accuracy 2*10-6",
2092  lineshape_voigt_kuntz6));
2093 
2094  lineshape_data.push_back
2095  (LineshapeRecord
2096  ("Voigt_Kuntz3",
2097  "The Voigt line shape. Approximation by Kuntz: Accuracy 2*10-3",
2098  lineshape_voigt_kuntz3));
2099 
2100  lineshape_data.push_back
2101  (LineshapeRecord
2102  ("Voigt_Kuntz4",
2103  "The Voigt line shape. Approximation by Kuntz: Accuracy 2*10-4",
2104  lineshape_voigt_kuntz4));
2105 
2106  lineshape_data.push_back
2107  (LineshapeRecord
2108  ("Voigt_Drayson",
2109  "The Voigt line shape. Approximation by Drayson.",
2110  lineshape_voigt_drayson));
2111 
2112  lineshape_data.push_back
2113  (LineshapeRecord
2114  ("Rosenkranz_Voigt_Kuntz6",
2115  "Rosenkranz lineshape for oxygen with overlap correction, "
2116  "at high altitudes Voigt_Kuntz6.",
2117  lineshape_rosenkranz_voigt_kuntz6));
2118 
2119  lineshape_data.push_back
2120  (LineshapeRecord
2121  ("Rosenkranz_Voigt_Drayson",
2122  "Rosenkranz lineshape for oxygen with overlap correction, "
2123  "at high altitudes Drayson.",
2124  lineshape_rosenkranz_voigt_drayson));
2125 
2126  lineshape_data.push_back
2127  (LineshapeRecord
2128  ("CO2_Lorentz",
2129  "Special line shape for CO2 in the infrared, neglecting Doppler\n"
2130  "broadening and details of line mixing. The line shape can be\n"
2131  "expressed as\n"
2132  " chi(f,f0) * Lorentz(f,f0) \n"
2133  "\n"
2134  "The chi-factor follows Cousin et al. 1985. Broadening by N2 and O2\n"
2135  "is considered, while self-broadening (CO2-CO2) is neglected."
2136  "\n"
2137  "NOTE: Temperature dependency is not yet included. The chi factor is\n"
2138  "valid for 238 K.",
2139  lineshape_CO2_lorentz));
2140 
2141  lineshape_data.push_back
2142  (LineshapeRecord
2143  ("CO2_Drayson",
2144  "Special line shape for CO2 in the infrared, neglecting details of\n"
2145  "line mixing. The line shape can be expressed as\n"
2146  " chi(f,f0) * Drayson(f,f0) \n"
2147  "\n"
2148  "The chi-factor follows Cousin et al. 1985. Broadening by N2 and O2\n"
2149  "is considered, while self-broadening (CO2-CO2) is neglected.\n"
2150  "\n"
2151  "NOTE: Temperature dependency is not yet included. The chi factor is\n"
2152  "valid for 238 K.",
2153  lineshape_CO2_drayson));
2154 
2155 }
2156 
2159 Array<LineshapeNormRecord> lineshape_norm_data;
2160 
2161 void define_lineshape_norm_data()
2162 {
2163  // Initialize to empty, just in case.
2164  lineshape_norm_data.resize(0);
2165 
2166  lineshape_norm_data.push_back
2167  (LineshapeNormRecord
2168  ("no_norm",
2169  "No normalization of the lineshape.",
2170  lineshape_norm_no_norm));
2171 
2172  lineshape_norm_data.push_back
2173  (LineshapeNormRecord
2174  ("linear",
2175  "Linear normalization of the lineshape with f/f0.",
2176  lineshape_norm_linear));
2177 
2178  lineshape_norm_data.push_back
2179  (LineshapeNormRecord
2180  ("quadratic",
2181  "Quadratic normalization of the lineshape with (f/f0)^2.",
2182  lineshape_norm_quadratic));
2183 
2184  lineshape_norm_data.push_back
2185  (LineshapeNormRecord
2186  ("VVH",
2187  "Van Vleck Huber normalization of the lineshape with\n"
2188  " (f*tanh(h*f/(2*k*T))) / (f0*tanh(h*f0/(2*k*T))).\n"
2189  " The denominator is a result of catalogue intensities.",
2190  lineshape_norm_VVH));
2191 }
lineshape_voigt_kuntz6
void lineshape_voigt_kuntz6(Vector &ls, Vector &x, Numeric f0, Numeric gamma, Numeric sigma, VectorView f_mono, const Index nf)
Definition: lineshapes.cc:239
PLANCK_CONST
const Numeric PLANCK_CONST
lineshape_norm_VVH
void lineshape_norm_VVH(Vector &fac, Numeric f0, VectorView f_mono, const Numeric T, const Index nf)
Definition: lineshapes.cc:2027
absorption.h
Declarations required for the calculation of absorption coefficients.
b0
#define b0
Definition: continua.cc:13957
LineshapeRecord
Lineshape related information.
Definition: absorption.h:51
bfun6_
long bfun6_(Numeric y, Numeric x)
Definition: lineshapes.cc:155
lineshape_CO2_lorentz
void lineshape_CO2_lorentz(Vector &ls, Vector &, Numeric f0, Numeric gamma, Numeric, VectorView f_mono, const Index nf)
Definition: lineshapes.cc:1846
array.h
This file contains the definition of Array.
lineshape_norm_linear
void lineshape_norm_linear(Vector &fac, Numeric f0, VectorView f_mono, const Numeric, const Index nf)
Definition: lineshapes.cc:1964
matpackI.h
Array
This can be used to make arrays out of anything.
Definition: array.h:48
lineshape_lorentz
void lineshape_lorentz(Vector &ls, Vector &, Numeric f0, Numeric gamma, Numeric, VectorView f_mono, const Index nf)
Definition: lineshapes.cc:83
lineshape_rosenkranz_voigt_drayson
void lineshape_rosenkranz_voigt_drayson(Vector &ls, Vector &x, Numeric f0, Numeric gamma, Numeric sigma, VectorView f_mono, const Index nf)
Definition: lineshapes.cc:1688
lineshape_norm_quadratic
void lineshape_norm_quadratic(Vector &fac, Numeric f0, VectorView f_mono, const Numeric, const Index nf)
Definition: lineshapes.cc:1991
bfun3_
long int bfun3_(Numeric y, Numeric x)
Definition: lineshapes.cc:564
lineshape_CO2_drayson
void lineshape_CO2_drayson(Vector &ls, Vector &X, Numeric f0, Numeric gamma, Numeric sigma, VectorView f_mono, const Index nf)
Definition: lineshapes.cc:1892
Numeric
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: arts.h:147
ConstVectorView::nelem
Index nelem() const
Returns the number of elements.
Definition: matpackI.h:937
VectorView
The VectorView class.
Definition: matpackI.h:281
define_lineshape_data
void define_lineshape_data()
Definition: lineshapes.cc:2064
bfun4_
long bfun4_(Numeric y, Numeric x)
Definition: lineshapes.cc:942
max
#define max(a, b)
Definition: continua.cc:12949
lineshape_doppler
void lineshape_doppler(Vector &ls, Vector &, Numeric f0, Numeric, Numeric sigma, VectorView f_mono, const Index nf)
Definition: lineshapes.cc:122
lineshape_rosenkranz_voigt_kuntz6
void lineshape_rosenkranz_voigt_kuntz6(Vector &ls, Vector &x, Numeric f0, Numeric gamma, Numeric sigma, VectorView f_mono, const Index nf)
Definition: lineshapes.cc:1580
lineshape_voigt_drayson
void lineshape_voigt_drayson(Vector &ls, Vector &x, Numeric f0, Numeric gamma, Numeric sigma, VectorView f_mono, const Index nf)
Definition: lineshapes.cc:1414
lineshape_voigt_kuntz3
void lineshape_voigt_kuntz3(Vector &ls, Vector &x, Numeric f0, Numeric gamma, Numeric sigma, VectorView f_mono, const Index nf)
Definition: lineshapes.cc:655
Index
INDEX Index
The type to use for all integer numbers and indices.
Definition: arts.h:153
chi_cousin
void chi_cousin(Numeric &chi, const Numeric &df)
Definition: lineshapes.cc:1796
min
#define min(a, b)
Definition: continua.cc:12948
define_lineshape_norm_data
void define_lineshape_norm_data()
Definition: lineshapes.cc:2161
lineshape_data
Array< LineshapeRecord > lineshape_data
Definition: lineshapes.cc:2062
HZ2CM
const Numeric HZ2CM
Global constant, converts Hz to cm-1.
lineshape_no_shape
void lineshape_no_shape(Vector &, Vector &, Numeric, Numeric, Numeric, VectorView, const Index)
Definition: lineshapes.cc:57
lineshape_voigt_kuntz4
void lineshape_voigt_kuntz4(Vector &ls, Vector &x, Numeric f0, Numeric gamma, Numeric sigma, VectorView f_mono, const Index nf)
Definition: lineshapes.cc:1029
lineshape_norm_data
Array< LineshapeNormRecord > lineshape_norm_data
Definition: lineshapes.cc:2159
BOLTZMAN_CONST
const Numeric BOLTZMAN_CONST
Vector
The Vector class.
Definition: matpackI.h:389
lineshape_norm_no_norm
void lineshape_norm_no_norm(Vector &fac, Numeric, VectorView f_mono, const Numeric, const Index nf)
Definition: lineshapes.cc:1936
LineshapeNormRecord
Lineshape related normalization function information.
Definition: absorption.h:92
PI
const Numeric PI
Global constant, pi.
arts.h
The global header file for ARTS.
Generated by doxygen 1.8.20