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