ARTS: lineshapes.cc Source File

00001 /* Copyright (C) 2000, 2001 Axel von Engeln <engeln@uni-bremen.de>
00002                             Stefan Buehler  <sbuehler@uni-bremen.de>
00003 
00004    This program is free software; you can redistribute it and/or modify it
00005    under the terms of the GNU General Public License as published by the
00006    Free Software Foundation; either version 2, or (at your option) any
00007    later version.
00008 
00009    This program is distributed in the hope that it will be useful,
00010    but WITHOUT ANY WARRANTY; without even the implied warranty of
00011    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012    GNU General Public License for more details.
00013 
00014    You should have received a copy of the GNU General Public License
00015    along with this program; if not, write to the Free Software
00016    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
00017    USA. */
00018 
00031 #include <math.h>
00032 #include "arts.h"
00033 #include "matpackI.h"
00034 #include "array.h"
00035 #include "absorption.h"
00036 
00037 
00057 void lineshape_no_shape(  Vector&       /* ls */,
00058                           Vector&       /* X */,
00059                           Numeric       /* f0 */,
00060                           Numeric       /* gamma */,
00061                           Numeric       /* sigma */,
00062                           VectorView    /* f_mono */,
00063                           const Index   /* nf */)
00064 {
00065   // This function should never be called so throw an error here: 
00066   throw runtime_error("The no_shape lineshape is only a placeholder, but you tried\n"
00067                       "to use it like a real lineshape.");
00068 }
00069 
00070 
00083 void lineshape_lorentz(Vector&       ls,
00084                        Vector&       /* X */,
00085                        Numeric       f0,
00086                        Numeric       gamma,
00087                        Numeric       /* sigma */,
00088                        VectorView f_mono,
00089                        const Index  nf)
00090 {
00091   // FIXME: nf is actually redundant. Could be thrown out in the
00092   // future. For now, let's do an assertion that at least it is
00093   // correct: 
00094   assert( nf==f_mono.nelem() );
00095 
00096   // PI:
00097   extern const Numeric PI;
00098 
00099   //  assert( ls.nelem() == nf );
00100 
00101   Numeric gamma2 = gamma * gamma;
00102   Numeric fac = gamma/PI;
00103 
00104   for ( Index i=0; i<nf; ++i )
00105     {
00106       ls[i] =  fac / ( (f_mono[i]-f0) * (f_mono[i]-f0) + gamma2 );
00107     }
00108 }
00109 
00122 void lineshape_doppler(Vector&       ls,
00123                        Vector&       /* x */,
00124                        Numeric       f0,
00125                        Numeric       /* gamma */,
00126                        Numeric       sigma,
00127                        VectorView f_mono,
00128                        const Index  nf)
00129 {
00130   // FIXME: nf is actually redundant. Could be thrown out in the
00131   // future. For now, let's do an assertion that at least it is
00132   // correct: 
00133   assert( nf==f_mono.nelem() );
00134 
00135   // SQRT(PI):
00136   extern const Numeric PI;
00137   static const Numeric sqrtPI = sqrt(PI);
00138 
00139   //  assert( ls.nelem() == nf );
00140 
00141   Numeric sigma2 = sigma * sigma;
00142   Numeric fac = 1.0 / (sqrtPI * sigma);
00143   
00144   for ( Index i=0; i<nf ; ++i )
00145     {
00146       ls[i] = fac * exp( - pow( f_mono[i]-f0, 2) / sigma2 );
00147     }
00148 }
00149 
00150 
00151 
00152 //------------------------------------------------------------------------
00153 
00154 // help function for lineshape_voigt_kuntz1
00155 long bfun6_(Numeric y, Numeric x)
00156 {
00157   /* System generated locals */
00158   long int ret_val;
00159 
00160   /* Local variables */
00161   static Numeric s;
00162 
00163   /* -------------------------------------------------------------------- */
00164   s = fabs(x) + y;
00165   if (s >= 15.f) {
00166     ret_val = 1;
00167   } else if (s >= 5.5f) {
00168     ret_val = 2;
00169   } else if (y >= fabs(x) * .195f - .176f) {
00170     ret_val = 3;
00171   } else {
00172     ret_val = 4;
00173   }
00174   return ret_val;
00175 } /* bfun6_ */
00176 
00177 
00178 
00239 void lineshape_voigt_kuntz6(Vector&       ls,
00240                             Vector&       x,
00241                             Numeric       f0,
00242                             Numeric       gamma,
00243                             Numeric       sigma,
00244                             VectorView f_mono,
00245                             const Index  nf)
00246 
00247 {
00248   // FIXME: nf is actually redundant. Could be thrown out in the
00249   // future. For now, let's do an assertion that at least it is
00250   // correct: 
00251   assert( nf==f_mono.nelem() );
00252   
00253 
00254   // seems not necessary for Doppler correction
00255   //    extern const Numeric SQRT_NAT_LOG_2;
00256 
00257   // PI
00258   extern const Numeric PI;
00259 
00260   // constant sqrt(1/pi)
00261   const Numeric sqrt_invPI =  sqrt(1/PI);
00262 
00263   // constant normalization factor for voigt
00264   Numeric fac = 1.0 / sigma * sqrt_invPI;
00265 
00266 
00267   /* Initialized data */
00268 
00269   static Numeric yps1 = -1.0;
00270   static Numeric yps2 = -1.0;
00271   static Numeric yps3 = -1.0;
00272   static Numeric yps4 = -1.0;
00273 
00274   /* System generated locals */
00275   long int i__1, i__2;
00276   Numeric r__1;
00277 
00278   /* Local variables */
00279   static long int bmin, lauf[16]        /* was [4][4] */, bmax;
00280   static long int imin, imax, stack[80] /* was [20][4] */;
00281   static Numeric a1, a2, a3, a4, a5, a6, a8, b8, c8, d8, e8, f8, g8, h8, a7, 
00282     b7, c7, d7, e7, f7, o8, p8, q8, r8, s8, t8, g7, h7, o7, p7, q7, 
00283     r7, s7, t7, b6, c6, d6, e6, b5, c5, d5, e5, b4, c4, d4, b3, c3, 
00284     d3, b1, y2;
00285   static long int i2, i1;
00286   static Numeric x2, b2, c1;
00287   static long int stackp, imitte;
00288   static Numeric ym2;
00289 
00290 
00291   // variables needed in original c routine:
00292 
00293   // Ratio of the Lorentz halfwidth to the Doppler halfwidth
00294   Numeric y = gamma / sigma;
00295 
00296   // frequency in units of Doppler 
00297   for (i1=0; i1< (int) nf; i1++)
00298     {
00299       x[i1] = (f_mono[i1] - f0) / sigma;
00300     }
00301 
00302   /* Parameter adjustments */
00303   // this does not work for variables type Vector, corrected by
00304   // adjusting the index to array ls and x
00305   //--ls;
00306   //--x;
00307 
00308   /* Function Body */
00309   y2 = y * y;
00310   if (y >= 15.0 || x[0] >= 15.0 || x[nf-1] <= -15.0) {
00311     lauf[0] = 1;
00312     lauf[4] = nf;
00313     lauf[8] = nf;
00314     lauf[12] = 0;
00315     goto L7;
00316   }
00317   for (i2 = 1; i2 <= 4; ++i2) {
00318     for (i1 = 1; i1 <= 4; ++i1) {
00319       lauf[i1 + (i2 << 2) - 5] = i2 % 2 * (nf + 1);
00320       /* L1: */
00321     }
00322   }
00323   stackp = 1;
00324   stack[stackp - 1] = 1;
00325   stack[stackp + 19] = nf;
00326   stack[stackp + 39] = bfun6_(y, x[0]);
00327   stack[stackp + 59] = bfun6_(y, x[nf-1]);
00328  L2:
00329   imin = stack[stackp - 1];
00330   imax = stack[stackp + 19];
00331   bmin = stack[stackp + 39];
00332   bmax = stack[stackp + 59];
00333   if (bmin == bmax) {
00334     if (x[imax-1] < 0.f) {
00335       /* Computing MIN */
00336       i__1 = imin, i__2 = lauf[bmin - 1];
00337       lauf[bmin - 1] = min(i__1,i__2);
00338       /* Computing MAX */
00339       i__1 = imax, i__2 = lauf[bmax + 3];
00340       lauf[bmax + 3] = max(i__1,i__2);
00341       --stackp;
00342       goto L3;
00343     } else if (x[imin-1] >= 0.f) {
00344       /* Computing MIN */
00345       i__1 = imin, i__2 = lauf[bmin + 7];
00346       lauf[bmin + 7] = min(i__1,i__2);
00347       /* Computing MAX */
00348       i__1 = imax, i__2 = lauf[bmax + 11];
00349       lauf[bmax + 11] = max(i__1,i__2);
00350       --stackp;
00351       goto L3;
00352     }
00353   }
00354   imitte = (imax + imin) / 2;
00355   stack[stackp - 1] = imitte + 1;
00356   stack[stackp + 19] = imax;
00357   stack[stackp + 39] = bfun6_(y, x[imitte]);
00358   stack[stackp + 59] = bmax;
00359   ++stackp;
00360   stack[stackp - 1] = imin;
00361   stack[stackp + 19] = imitte;
00362   stack[stackp + 39] = bmin;
00363   stack[stackp + 59] = bfun6_(y, x[imitte-1]);
00364  L3:
00365   if (stackp > 0) {
00366     goto L2;
00367   }
00368   /* ---- Region 4 */
00369   /* -------------------------------------------------------------------- */
00370   if (lauf[7] >= lauf[3] || lauf[15] >= lauf[11]) {
00371     if ((r__1 = y - yps4, fabs(r__1)) > 1e-8f) {
00372       yps4 = y;
00373       a7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
00374                     y2 * (y2 * (y2 * (2.35944f - y2 * .56419f) - 72.9359f) + 
00375                     571.687f) - 5860.68f) + 40649.2f) - 320772.f) + 1684100.f)
00376                      - 9694630.f) + 40816800.f) - 1.53575e8f) + 4.56662e8f) - 
00377                     9.86604e8f) + 1.16028e9f);
00378       b7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
00379                     y2 * (y2 * (23.0312f - y2 * 7.33447f) - 234.143f) - 
00380                     2269.19f) + 45251.3f) - 234417.f) + 3599150.f) - 
00381                     7723590.f) + 86482900.f) - 2.91876e8f) + 8.06985e8f) - 
00382                     9.85386e8f) - 5.60505e8f);
00383       c7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
00384                     y2 * (97.6203f - y2 * 44.0068f) + 1097.77f) - 25338.3f) + 
00385                     98079.1f) + 576054.f) - 2.3818e7f) + 22930200.f) - 
00386                     2.04467e8f) + 2.94262e8f) + 2.47157e8f) - 6.51523e8f);
00387       d7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
00388                     228.563f - y2 * 161.358f) + 8381.97f) - 66431.2f) - 
00389                     303569.f) + 2240400.f) + 38311200.f) - 41501300.f) - 
00390                     99622400.f) + 2.70167e8f) - 2.63894e8f);
00391       e7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
00392                     296.38f - y2 * 403.396f) + 23507.6f) - 66212.1f) - 
00393                     1.003e6f) + 468142.f) + 24620100.f) + 5569650.f) + 
00394                     1.40677e8f) - 63177100.f);
00395       f7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (125.591f - 
00396                     y2 * 726.113f) + 37544.8f) + 8820.94f) - 934717.f) - 
00397                     1931140.f) - 33289600.f) + 4073820.f) - 16984600.f);
00398       g7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (-260.198f - y2 * 
00399                     968.15f) + 37371.9f) + 79902.5f) - 186682.f) - 900010.f) 
00400                     + 7528830.f) - 1231650.f);
00401       h7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (-571.645f - y2 * 968.15f)
00402                      + 23137.1f) + 72520.9f) + 153468.f) + 86407.6f) - 
00403                     610622.f);
00404       o7 = y * (y2 * (y2 * (y2 * (y2 * (-575.164f - y2 * 726.113f) + 
00405                     8073.15f) + 26538.5f) + 49883.8f) - 23586.5f);
00406       p7 = y * (y2 * (y2 * (y2 * (-352.467f - y2 * 403.396f) + 
00407                     953.655f) + 2198.86f) - 8009.1f);
00408       q7 = y * (y2 * (y2 * (-134.792f - y2 * 161.358f) - 271.202f) - 
00409                     622.056f);
00410       r7 = y * (y2 * (-29.7896f - y2 * 44.0068f) - 77.0535f);
00411       s7 = y * (-2.92264f - y2 * 7.33447f);
00412       t7 = y * -.56419f;
00413       a8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
00414                     y2 * (y2 * (y2 * (y2 - 3.68288f) + 126.532f) - 955.194f) 
00415                     + 9504.65f) - 70946.1f) + 483737.f) - 2857210.f) + 
00416                     14464700.f) - 61114800.f) + 2.11107e8f) - 5.79099e8f) + 
00417                     1.17022e9f) - 1.5599e9f) + 1.02827e9f;
00418       b8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
00419                     y2 * (y2 * (y2 * 14.f - 40.5117f) + 533.254f) + 3058.26f) 
00420                     - 55600.f) + 498334.f) - 2849540.f) + 13946500.f) - 
00421                     70135800.f) + 2.89676e8f) - 7.53828e8f) + 1.66421e9f) - 
00422                     2.28855e9f) + 1.5599e9f;
00423       c8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
00424                     y2 * (y2 * 91 - 198.876f) - 1500.17f) + 48153.3f) - 
00425                     217801.f) - 1063520.f) + 1.4841e7f) - 46039600.f) + 
00426                     63349600.f) - 6.60078e8f) + 1.06002e9f) - 1.66421e9f) + 
00427                     1.17022e9f;
00428       d8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
00429                     y2 * 364 - 567.164f) - 16493.7f) + 161461.f) + 280428.f) 
00430                     - 6890020.f) - 6876560.f) + 1.99846e8f) + 54036700.f) + 
00431                     6.60078e8f) - 7.53828e8f) + 5.79099e8f;
00432       e8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 
00433                     1001 - 1012.79f) - 55582.f) + 240373.f) + 1954700.f) - 
00434                     5257220.f) - 50101700.f) - 1.99846e8f) + 63349600.f) - 
00435                     2.89676e8f) + 2.11107e8f;
00436       f8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 2002 - 
00437                     1093.82f) - 106663.f) + 123052.f) + 3043160.f) + 
00438                     5257220.f) - 6876560.f) + 46039600.f) - 70135800.f) + 
00439                     61114800.f;
00440       g8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 3003 - 
00441                     486.14f) - 131337.f) - 123052.f) + 1954700.f) + 6890020.f)
00442                      + 1.4841e7f) - 13946500.f) + 14464700.f;
00443       h8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 3432 + 486.14f) - 
00444                     106663.f) - 240373.f) + 280428.f) + 1063520.f) - 
00445                     2849540.f) + 2857210.f;
00446       o8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 3003 + 1093.82f) - 
00447                     55582.f) - 161461.f) - 217801.f) - 498334.f) + 483737.f;
00448       p8 = y2 * (y2 * (y2 * (y2 * (y2 * 2002 + 1012.79f) - 16493.7f) - 
00449                     48153.3f) - 55600.f) + 70946.1f;
00450       q8 = y2 * (y2 * (y2 * (y2 * 1001.f + 567.164f) - 1500.17f) - 
00451                     3058.26f) + 9504.65f;
00452       r8 = y2 * (y2 * (y2 * 364 + 198.876f) + 533.254f) + 955.194f;
00453       s8 = y2 * (y2 * 91.f + 40.5117f) + 126.532f;
00454       t8 = y2 * 14.f + 3.68288f;
00455     }
00456     ym2 = y * 2;
00457     for (i2 = 1; i2 <= 3; i2 += 2) {
00458       i__1 = lauf[(i2 + 1 << 2) - 1];
00459       for (i1 = lauf[(i2 << 2) - 1]; i1 <= i__1; ++i1) {
00460         x2 = x[i1-1] * x[i1-1];
00461         ls[i1-1] = fac * (exp(y2 - x2) * cos(x[i1-1] * ym2) - (a7 + x2 *
00462                          (b7 + x2 * (c7 + x2 * (d7 + x2 * (e7 + x2 * (f7 + x2 
00463                         * (g7 + x2 * (h7 + x2 * (o7 + x2 * (p7 + x2 * (q7 + 
00464                         x2 * (r7 + x2 * (s7 + x2 * t7))))))))))))) / (a8 + x2 
00465                         * (b8 + x2 * (c8 + x2 * (d8 + x2 * (e8 + x2 * (f8 + 
00466                         x2 * (g8 + x2 * (h8 + x2 * (o8 + x2 * (p8 + x2 * (q8 
00467                         + x2 * (r8 + x2 * (s8 + x2 * (t8 + x2)))))))))))))));
00468         /* L4: */
00469       }
00470     }
00471   }
00472   /* ---- Region 3 */
00473   /* -------------------------------------------------------------------- */
00474   if (lauf[6] >= lauf[2] || lauf[14] >= lauf[10]) {
00475     if ((r__1 = y - yps3, fabs(r__1)) > 1e-8f) {
00476       yps3 = y;
00477       a5 = y * (y * (y * (y * (y * (y * (y * (y * (y * 
00478                     .564224f + 7.55895f) + 49.5213f) + 204.501f) + 581.746f) 
00479                     + 1174.8f) + 1678.33f) + 1629.76f) + 973.778f) + 272.102f;
00480       b5 = y * (y * (y * (y * (y * (y * (y * 2.25689f + 22.6778f)
00481                      + 100.705f) + 247.198f) + 336.364f) + 220.843f) - 
00482                     2.34403f) - 60.5644f;
00483       c5 = y * (y * (y * (y * (y * 3.38534f + 22.6798f) + 52.8454f)
00484                      + 42.5683f) + 18.546f) + 4.58029f;
00485       d5 = y * (y * (y * 2.25689f + 7.56186f) + 1.66203f) - .128922f;
00486       e5 = y * .564224f + 9.71457e-4f;
00487       a6 = y * (y * (y * (y * (y * (y * (y * (y * (y * (y + 
00488                     13.3988f) + 88.2674f) + 369.199f) + 1074.41f) + 2256.98f) 
00489                     + 3447.63f) + 3764.97f) + 2802.87f) + 1280.83f) + 
00490                     272.102f;
00491       b6 = y * (y * (y * (y * (y * (y * (y * (y * 5.f + 
00492                     53.5952f) + 266.299f) + 793.427f) + 1549.68f) + 2037.31f) 
00493                     + 1758.34f) + 902.306f) + 211.678f;
00494       c6 = y * (y * (y * (y * (y * (y * 10.f + 80.3928f) + 
00495                     269.292f) + 479.258f) + 497.302f) + 308.186f) + 78.866f;
00496       d6 = y * (y * (y * (y * 10.f + 53.5952f) + 92.7568f) + 
00497                     55.0293f) + 22.0353f;
00498       e6 = y * (y * 5.f + 13.3988f) + 1.49645f;
00499     }
00500     for (i2 = 1; i2 <= 3; i2 += 2) {
00501       i__1 = lauf[(i2 + 1 << 2) - 2];
00502       for (i1 = lauf[(i2 << 2) - 2]; i1 <= i__1; ++i1) {
00503         x2 = x[i1-1] * x[i1-1];
00504         ls[i1-1] = fac * (a5 + x2 * (b5 + x2 * (c5 + x2 * (d5 + x2 * 
00505                         e5)))) / (a6 + x2 * (b6 + x2 * (c6 + x2 * (d6 + x2 * (
00506                         e6 + x2)))));
00507         /* L5: */
00508       }
00509     }
00510   }
00511   /* ---- Region 2 */
00512   /* -------------------------------------------------------------------- */
00513   if (lauf[5] >= lauf[1] || lauf[13] >= lauf[9]) {
00514     if ((r__1 = y - yps2, fabs(r__1)) > 1e-8f) {
00515       yps2 = y;
00516       a3 = y * (y2 * (y2 * (y2 * .56419f + 3.10304f) + 4.65456f) + 
00517                     1.05786f);
00518       b3 = y * (y2 * (y2 * 1.69257f + .56419f) + 2.962f);
00519       c3 = y * (y2 * 1.69257f - 2.53885f);
00520       d3 = y * .56419f;
00521       a4 = y2 * (y2 * (y2 * (y2 + 6.f) + 10.5f) + 4.5f) + .5625f;
00522       b4 = y2 * (y2 * (y2 * 4.f + 6.f) + 9.f) - 4.5f;
00523       c4 = y2 * (y2 * 6.f - 6.f) + 10.5f;
00524       d4 = y2 * 4.f - 6.f;
00525     }
00526     for (i2 = 1; i2 <= 3; i2 += 2) {
00527       i__1 = lauf[(i2 + 1 << 2) - 3];
00528       for (i1 = lauf[(i2 << 2) - 3]; i1 <= i__1; ++i1) {
00529         x2 = x[i1-1] * x[i1-1];
00530         ls[i1-1] = fac * (a3 + x2 * (b3 + x2 * (c3 + x2 * d3))) / (a4 
00531                         + x2 * (b4 + x2 * (c4 + x2 * (d4 + x2))));
00532         /* L6: */
00533       }
00534     }
00535   }
00536   /* ---- Region 1 */
00537   /* -------------------------------------------------------------------- */
00538  L7:
00539   if (lauf[4] >= lauf[0] || lauf[12] >= lauf[8]) {
00540     if ((r__1 = y - yps1, fabs(r__1)) > 1e-8f) {
00541       yps1 = y;
00542       a1 = y * .5641896f;
00543       b1 = y2 + .5f;
00544       a2 = y2 * 4;
00545     }
00546 
00547     c1 = fac * a1;
00548     for (i2 = 1; i2 <= 3; i2 += 2) {
00549       i__1 = lauf[(i2 + 1 << 2) - 4];
00550       for (i1 = lauf[(i2 << 2) - 4]; i1 <= i__1; ++i1) {
00551         x2 = x[i1-1] * x[i1-1];
00552         b2 = b1 - x2;
00553         ls[i1-1] = c1 * (b1 + x2) / (b2 * b2 + a2 * x2);
00554         /* L8: */
00555       }
00556     }
00557   }
00558 }
00559 
00560 
00561 //---------------------------------------------------------------
00562 // help function for lineshape_voigt_kuntz3
00563 
00564 long int bfun3_(Numeric y, Numeric x)
00565 {
00566     /* System generated locals */
00567     long int ret_val;
00568 
00569     /* Local variables */
00570     static Numeric x2, y2;
00571 
00572 /* -------------------------------------------------------------------- */
00573     x2 = x * x;
00574     y2 = y * y;
00575     if (x2 * .4081676f + y2 > 21.159543f) {
00576         if (x2 * .7019639f + y2 > 1123.14221f) {
00577             ret_val = 0;
00578         } else {
00579             ret_val = 1;
00580         }
00581     } else {
00582         if (x2 * .20753051f + y2 > 4.20249292f) {
00583             ret_val = 2;
00584         } else if (y >= fabs(x) * .08f - .12f) {
00585             ret_val = 3;
00586         } else {
00587             ret_val = 4;
00588         }
00589     }
00590     return ret_val;
00591 } /* bfun3_ */
00592 
00593 
00594 
00655 void lineshape_voigt_kuntz3(Vector&       ls,
00656                             Vector&       x,
00657                             Numeric       f0,
00658                             Numeric       gamma,
00659                             Numeric       sigma,
00660                             VectorView f_mono,
00661                             const Index  nf)
00662 
00663 {
00664   // FIXME: nf is actually redundant. Could be thrown out in the
00665   // future. For now, let's do an assertion that at least it is
00666   // correct: 
00667   assert( nf==f_mono.nelem() );
00668 
00669   // seems not necessary for Doppler correction
00670   //    extern const Numeric SQRT_NAT_LOG_2;
00671 
00672   // PI
00673   extern const Numeric PI;
00674 
00675   // constant sqrt(1/pi)
00676   const Numeric sqrt_invPI =  sqrt(1/PI);
00677 
00678   // constant normalization factor for voigt
00679   Numeric fac = 1.0 / sigma * sqrt_invPI;
00680 
00681   /* Initialized data */
00682 
00683   static Numeric yps0 = -1.0;
00684   static Numeric yps1 = -1.0;
00685   static Numeric yps2 = -1.0;
00686   static Numeric yps3 = -1.0;
00687   static Numeric yps4 = -1.0;
00688 
00689   /* System generated locals */
00690   long i__1, i__2;
00691   Numeric r__1;
00692 
00693   /* Local variables */
00694   static long bmin, lauf[20]    /* was [5][4] */, bmax, imin, imax;
00695   static long stack[80] /* was [20][4] */;
00696   static Numeric a0, a1, a2, a3, a4, a5, a6, a7, a8, b8, c8, d8, e8, f8, g8, 
00697     h8, b7, c7, d7, e7, f7, g7, o8, p8, q8, r8, s8, t8, h7, o7, p7, 
00698     q7, r7, s7, t7, b6, c6, d6, e6, b5, c5, d5, e5, b4, c4, d4, b3, 
00699     c3, d3, b1, y2;
00700   static long i2, i1;
00701   static Numeric x2, b2, b0, c1;
00702   static long stackp, imitte;
00703   static Numeric ym2;
00704 
00705   // variables needed in original c routine:
00706 
00707   // Ratio of the Lorentz halfwidth to the Doppler halfwidth
00708   Numeric y = gamma / sigma;
00709 
00710   // frequency in units of Doppler 
00711   for (i1=0; i1< (int) nf; i1++)
00712     {
00713       x[i1] = (f_mono[i1] - f0) / sigma;
00714     }
00715 
00716 
00717   /* Parameter adjustments */
00718   // this does not work for variables type Vector, corrected by
00719   // adjusting the index to array ls and x
00720   //--ls;
00721   //--x;
00722 
00723   /* Function Body */
00724   y2 = y * y;
00725   if (y >= 23.0 || x[0] >= 39.0 || x[nf-1] <= -39.0) {
00726     lauf[0] = 1;
00727     lauf[5] = nf;
00728     lauf[10] = nf;
00729     lauf[15] = 0;
00730     goto L8;
00731   }
00732   for (i2 = 1; i2 <= 4; ++i2) {
00733     for (i1 = 0; i1 <= 4; ++i1) {
00734       lauf[i1 + i2 * 5 - 5] = i2 % 2 * (nf + 1);
00735       /* L1: */
00736     }
00737   }
00738   stackp = 1;
00739   stack[stackp - 1] = 1;
00740   stack[stackp + 19] = nf;
00741   stack[stackp + 39] = bfun3_(y, x[0]);
00742   stack[stackp + 59] = bfun3_(y, x[nf-1]);
00743  L2:
00744   imin = stack[stackp - 1];
00745   imax = stack[stackp + 19];
00746   bmin = stack[stackp + 39];
00747   bmax = stack[stackp + 59];
00748   if (bmin == bmax) {
00749     if (x[imax-1] < 0.f) {
00750       /* Computing MIN */
00751       i__1 = imin, i__2 = lauf[bmin];
00752       lauf[bmin] = min(i__1,i__2);
00753       /* Computing MAX */
00754       i__1 = imax, i__2 = lauf[bmax + 5];
00755       lauf[bmax + 5] = max(i__1,i__2);
00756       --stackp;
00757       goto L3;
00758     } else if (x[imin-1] >= 0.f) {
00759       /* Computing MIN */
00760       i__1 = imin, i__2 = lauf[bmin + 10];
00761       lauf[bmin + 10] = min(i__1,i__2);
00762       /* Computing MAX */
00763       i__1 = imax, i__2 = lauf[bmax + 15];
00764       lauf[bmax + 15] = max(i__1,i__2);
00765       --stackp;
00766       goto L3;
00767     }
00768   }
00769   imitte = (imax + imin) / 2;
00770   stack[stackp - 1] = imitte + 1;
00771   stack[stackp + 19] = imax;
00772   stack[stackp + 39] = bfun3_(y, x[imitte]);
00773   stack[stackp + 59] = bmax;
00774   ++stackp;
00775   stack[stackp - 1] = imin;
00776   stack[stackp + 19] = imitte;
00777   stack[stackp + 39] = bmin;
00778   stack[stackp + 59] = bfun3_(y, x[imitte-1]);
00779  L3:
00780   if (stackp > 0) {
00781     goto L2;
00782   }
00783   /* ---- Region 4 */
00784   /* -------------------------------------------------------------------- */
00785   if (lauf[9] >= lauf[4] || lauf[19] >= lauf[14]) {
00786     if ((r__1 = y - yps4, fabs(r__1)) > 1e-8f) {
00787       yps4 = y;
00788       a7 = y * (y2 * (y2 * 4.56662e8f - 9.86604e8f) + 1.16028e9f);
00789       b7 = y * (y2 * (y2 * 8.06985e8f - 9.85386e8f) - 5.60505e8f);
00790       c7 = y * (y2 * (y2 * 2.94262e8f + 2.47157e8f) - 6.51523e8f);
00791       d7 = y * (y2 * (2.70167e8f - y2 * 99622400.f) - 2.63894e8f);
00792       e7 = y * (y2 * (y2 * 5569650.f + 1.40677e8f) - 63177100.f);
00793       f7 = y * (y2 * (4073820.f - y2 * 33289600.f) - 16984600.f);
00794       g7 = y * (y2 * (7528830.f - y2 * 900010) - 1231650.f);
00795       h7 = y * (y2 * (y2 * 153468 + 86407.6f) - 610622.f);
00796       o7 = y * (y2 * (y2 * 26538.5f + 49883.8f) - 23586.5f);
00797       p7 = y * (y2 * (y2 * 953.655f + 2198.86f) - 8009.1f);
00798       q7 = y * (y2 * (-271.202f - y2 * 134.792f) - 622.056f);
00799       r7 = y * (y2 * (-29.7896f - y2 * 44.0068f) - 77.0535f);
00800       s7 = y * (-2.92264f - y2 * 7.33447f);
00801       t7 = y * -.56419f;
00802       a8 = y2 * (y2 * 1.17022e9f - 1.5599e9f) + 1.02827e9f;
00803       b8 = y2 * (y2 * 1.66421e9f - 2.28855e9f) + 1.5599e9f;
00804       c8 = y2 * (y2 * 1.06002e9f - 1.66421e9f) + 1.17022e9f;
00805       d8 = y2 * (y2 * 6.60078e8f - 7.53828e8f) + 5.79099e8f;
00806       e8 = y2 * (y2 * 63349600.f - 2.89676e8f) + 2.11107e8f;
00807       f8 = y2 * (y2 * 46039600.f - 70135800.f) + 61114800.f;
00808       g8 = y2 * (y2 * 1.4841e7f - 13946500.f) + 14464700.f;
00809       h8 = y2 * (y2 * 1063520.f - 2849540.f) + 2857210.f;
00810       o8 = y2 * (-498334.f - y2 * 217801.f) + 483737.f;
00811       p8 = y2 * (-55600.f - y2 * 48153.3f) + 70946.1f;
00812       q8 = y2 * (-3058.26f - y2 * 1500.17f) + 9504.65f;
00813       r8 = y2 * (y2 * 198.876f + 533.254f) + 955.194f;
00814       s8 = y2 * (y2 * 91.f + 40.5117f) + 126.532f;
00815       t8 = y2 * 14.f + 3.68288f;
00816     }
00817     ym2 = y * 2;
00818     for (i2 = 1; i2 <= 3; i2 += 2) {
00819       i__1 = lauf[(i2 + 1) * 5 - 1];
00820       for (i1 = lauf[i2 * 5 - 1]; i1 <= i__1; ++i1) {
00821         x2 = x[i1-1] * x[i1-1];
00822         ls[i1-1] = fac * (exp(y2 - x2) * cos(x[i1-1] * ym2) - (a7 + x2 *
00823                          (b7 + x2 * (c7 + x2 * (d7 + x2 * (e7 + x2 * (f7 + x2 
00824                         * (g7 + x2 * (h7 + x2 * (o7 + x2 * (p7 + x2 * (q7 + 
00825                         x2 * (r7 + x2 * (s7 + x2 * t7))))))))))))) / (a8 + x2 
00826                         * (b8 + x2 * (c8 + x2 * (d8 + x2 * (e8 + x2 * (f8 + 
00827                         x2 * (g8 + x2 * (h8 + x2 * (o8 + x2 * (p8 + x2 * (q8 
00828                         + x2 * (r8 + x2 * (s8 + x2 * (t8 + x2)))))))))))))));
00829         /* L4: */
00830       }
00831     }
00832   }
00833   /* ---- Region 3 */
00834   /* -------------------------------------------------------------------- */
00835   if (lauf[8] >= lauf[3] || lauf[18] >= lauf[13]) {
00836     if ((r__1 = y - yps3, fabs(r__1)) > 1e-8f) {
00837       yps3 = y;
00838       a5 = y * (y * (y * (y * (y * (y * (y * (y * (y * 
00839                     .564224f + 7.55895f) + 49.5213f) + 204.501f) + 581.746f) 
00840                     + 1174.8f) + 1678.33f) + 1629.76f) + 973.778f) + 272.102f;
00841       b5 = y * (y * (y * (y * (y * (y * (y * 2.25689f + 22.6778f)
00842                      + 100.705f) + 247.198f) + 336.364f) + 220.843f) - 
00843                     2.34403f) - 60.5644f;
00844       c5 = y * (y * (y * (y * (y * 3.38534f + 22.6798f) + 52.8454f)
00845                      + 42.5683f) + 18.546f) + 4.58029f;
00846       d5 = y * (y * (y * 2.25689f + 7.56186f) + 1.66203f) - .128922f;
00847       e5 = y * .564224f + 9.71457e-4f;
00848       a6 = y * (y * (y * (y * (y * (y * (y * (y * (y * (y + 
00849                     13.3988f) + 88.2674f) + 369.199f) + 1074.41f) + 2256.98f) 
00850                     + 3447.63f) + 3764.97f) + 2802.87f) + 1280.83f) + 
00851                     272.102f;
00852       b6 = y * (y * (y * (y * (y * (y * (y * (y * 5.f + 
00853                     53.5952f) + 266.299f) + 793.427f) + 1549.68f) + 2037.31f) 
00854                     + 1758.34f) + 902.306f) + 211.678f;
00855       c6 = y * (y * (y * (y * (y * (y * 10.f + 80.3928f) + 
00856                     269.292f) + 479.258f) + 497.302f) + 308.186f) + 78.866f;
00857       d6 = y * (y * (y * (y * 10.f + 53.5952f) + 92.7568f) + 
00858                     55.0293f) + 22.0353f;
00859       e6 = y * (y * 5.f + 13.3988f) + 1.49645f;
00860     }
00861     for (i2 = 1; i2 <= 3; i2 += 2) {
00862       i__1 = lauf[(i2 + 1) * 5 - 2];
00863       for (i1 = lauf[i2 * 5 - 2]; i1 <= i__1; ++i1) {
00864         x2 = x[i1-1] * x[i1-1];
00865         ls[i1-1] = fac * (a5 + x2 * (b5 + x2 * (c5 + x2 * (d5 + x2 * 
00866                         e5)))) / (a6 + x2 * (b6 + x2 * (c6 + x2 * (d6 + x2 * (
00867                         e6 + x2)))));
00868         /* L5: */
00869       }
00870     }
00871   }
00872   /* ---- Region 2 */
00873   /* -------------------------------------------------------------------- */
00874   if (lauf[7] >= lauf[2] || lauf[17] >= lauf[12]) {
00875     if ((r__1 = y - yps2, fabs(r__1)) > 1e-8f) {
00876       yps2 = y;
00877       a3 = y * (y2 * (y2 * (y2 * .56419f + 3.10304f) + 4.65456f) + 
00878                 1.05786f);
00879       b3 = y * (y2 * (y2 * 1.69257f + .56419f) + 2.962f);
00880       c3 = y * (y2 * 1.69257f - 2.53885f);
00881       d3 = y * .56419f;
00882       a4 = y2 * (y2 * (y2 * (y2 + 6.f) + 10.5f) + 4.5f) + .5625f;
00883       b4 = y2 * (y2 * (y2 * 4.f + 6.f) + 9.f) - 4.5f;
00884       c4 = y2 * (y2 * 6.f - 6.f) + 10.5f;
00885       d4 = y2 * 4.f - 6.f;
00886     }
00887     for (i2 = 1; i2 <= 3; i2 += 2) {
00888       i__1 = lauf[(i2 + 1) * 5 - 3];
00889       for (i1 = lauf[i2 * 5 - 3]; i1 <= i__1; ++i1) {
00890         x2 = x[i1-1] * x[i1-1];
00891         ls[i1-1] = fac * (a3 + x2 * (b3 + x2 * (c3 + x2 * d3))) / (a4 
00892                         + x2 * (b4 + x2 * (c4 + x2 * (d4 + x2))));
00893         /* L6: */
00894       }
00895     }
00896   }
00897   /* ---- Region 1 */
00898   /* -------------------------------------------------------------------- */
00899   if (lauf[6] >= lauf[1] || lauf[16] >= lauf[11]) {
00900     if ((r__1 = y - yps1, fabs(r__1)) > 1e-8f) {
00901       yps1 = y;
00902       a1 = y * .5641896f;
00903       b1 = y2 + .5f;
00904       a2 = y2 * 4;
00905     }
00906 
00907     c1 = a1 * fac;
00908     for (i2 = 1; i2 <= 3; i2 += 2) {
00909       i__1 = lauf[(i2 + 1) * 5 - 4];
00910       for (i1 = lauf[i2 * 5 - 4]; i1 <= i__1; ++i1) {
00911         x2 = x[i1-1] * x[i1-1];
00912         b2 = b1 - x2;
00913         ls[i1-1] = c1 * (b1 + x2) / (b2 * b2 + a2 * x2);
00914         /* L7: */
00915       }
00916     }
00917   }
00918   /* ---- Region 0 (Lorentz) */
00919   /* -------------------------------------------------------------------- */
00920 L8:
00921   if (lauf[5] >= lauf[0] || lauf[15] >= lauf[10]) {
00922     if ((r__1 = y - yps0, fabs(r__1)) > 1e-8f) {
00923       yps0 = y;
00924       a0 = y * .5641896f;
00925     }
00926 
00927     b0 = a0 * fac;
00928     for (i2 = 1; i2 <= 3; i2 += 2) {
00929       i__1 = lauf[(i2 + 1) * 5 - 5];
00930       for (i1 = lauf[i2 * 5 - 5]; i1 <= i__1; ++i1) {
00931         ls[i1-1] = b0 / (x[i1-1] * x[i1-1] + y2);
00932         /* L9: */
00933       }
00934     }
00935   }
00936 }
00937 
00938 
00939 //------------------------------------------------------------------
00940 // help function for lineshape_voigt_kuntz4
00941 
00942 long bfun4_(Numeric y, Numeric x)
00943 {
00944     /* System generated locals */
00945     long ret_val;
00946 
00947     /* Local variables */
00948     static Numeric x2, y2;
00949 
00950     x2 = x * x;
00951     y2 = y * y;
00952     if (x2 * .0062f + y2 * .01417f > 1.f) {
00953         if (x2 * 6.2e-5f + y2 * 1.98373e-4f > 1.f) {
00954             ret_val = 0;
00955         } else {
00956             ret_val = 1;
00957         }
00958     } else {
00959         if (x2 * .041649f + y2 * .111111111f > 1.f) {
00960             ret_val = 2;
00961         } else if (y >= fabs(x) * .19487f - .1753846f) {
00962             ret_val = 3;
00963         } else {
00964             ret_val = 4;
00965         }
00966     }
00967     return ret_val;
00968 }
00969 
01029 void lineshape_voigt_kuntz4(Vector&       ls,
01030                             Vector&       x,
01031                             Numeric       f0,
01032                             Numeric       gamma,
01033                             Numeric       sigma,
01034                             VectorView f_mono,
01035                             const Index   nf)
01036 {
01037   // FIXME: nf is actually redundant. Could be thrown out in the
01038   // future. For now, let's do an assertion that at least it is
01039   // correct: 
01040   assert( nf==f_mono.nelem() );
01041 
01042   // seems not necessary for Doppler correction
01043   //    extern const Numeric SQRT_NAT_LOG_2;
01044 
01045   // PI
01046   extern const Numeric PI;
01047 
01048   // constant sqrt(1/pi)
01049   const Numeric sqrt_invPI =  sqrt(1/PI);
01050 
01051   // constant normalization factor for voigt
01052   Numeric fac = 1.0 / sigma * sqrt_invPI;
01053 
01054 
01055     /* Initialized data */
01056 
01057     static float yps0 = -1.f;
01058     static float yps1 = -1.f;
01059     static float yps2 = -1.f;
01060     static float yps3 = -1.f;
01061     static float yps4 = -1.f;
01062 
01063     /* System generated locals */
01064     long i__1, i__2;
01065     float r__1;
01066 
01067     /* Local variables */
01068     static long bmin, lauf[20]  /* was [5][4] */, bmax, imin, imax;
01069     static long stack[80]       /* was [20][4] */;
01070     static Numeric a0, a1, a2, a3, a4, a5, a6, a7, a8, b8, c8, d8, e8, f8, g8, 
01071             h8, b7, c7, d7, e7, f7, g7, o8, p8, q8, r8, s8, t8, h7, o7, p7, 
01072             q7, r7, s7, t7, b6, c6, d6, e6, b5, c5, d5, e5, b4, c4, d4, b3, 
01073             c3, d3, b1, y2;
01074     static long i2, i1;
01075     static Numeric x2, b2, b0, c1;
01076     static long stackp, imitte;
01077     static Numeric ym2;
01078 
01079   // variables needed in original c routine:
01080 
01081   // Ratio of the Lorentz halfwidth to the Doppler halfwidth
01082   Numeric y = gamma / sigma;
01083 
01084   // frequency in units of Doppler 
01085   for (i1=0; i1< (int) nf; i1++)
01086     {
01087       x[i1] = (f_mono[i1] - f0) / sigma;
01088     }
01089 
01090 
01091   /* Parameter adjustments */
01092   // this does not work for variables type Vector, corrected by
01093   // adjusting the index to array ls and x
01094   //--ls;
01095   //--x;
01096 
01097   /* Function Body */
01098   y2 = y * y;
01099   if (y >= 71.f || x[0] >= 123.f || x[nf-1] <= -123.f) {
01100     lauf[0] = 1;
01101     lauf[5] = nf;
01102     lauf[10] = nf;
01103     lauf[15] = 0;
01104     goto L8;
01105   }
01106   for (i2 = 1; i2 <= 4; ++i2) {
01107     for (i1 = 0; i1 <= 4; ++i1) {
01108       lauf[i1 + i2 * 5 - 5] = i2 % 2 * (nf + 1);
01109       /* L1: */
01110     }
01111   }
01112   stackp = 1;
01113   stack[stackp - 1] = 1;
01114   stack[stackp + 19] = nf;
01115   stack[stackp + 39] = bfun4_(y, x[0]);
01116   stack[stackp + 59] = bfun4_(y, x[nf-1]);
01117  L2:
01118   imin = stack[stackp - 1];
01119   imax = stack[stackp + 19];
01120   bmin = stack[stackp + 39];
01121   bmax = stack[stackp + 59];
01122   if (bmin == bmax) {
01123     if (x[imax-1] < 0.f) {
01124       /* Computing MIN */
01125       i__1 = imin, i__2 = lauf[bmin];
01126       lauf[bmin] = min(i__1,i__2);
01127       /* Computing MAX */
01128       i__1 = imax, i__2 = lauf[bmax + 5];
01129       lauf[bmax + 5] = max(i__1,i__2);
01130       --stackp;
01131       goto L3;
01132     } else if (x[imin-1] >= 0.f) {
01133       /* Computing MIN */
01134       i__1 = imin, i__2 = lauf[bmin + 10];
01135       lauf[bmin + 10] = min(i__1,i__2);
01136       /* Computing MAX */
01137       i__1 = imax, i__2 = lauf[bmax + 15];
01138       lauf[bmax + 15] = max(i__1,i__2);
01139       --stackp;
01140       goto L3;
01141     }
01142   }
01143   imitte = (imax + imin) / 2;
01144   stack[stackp - 1] = imitte + 1;
01145   stack[stackp + 19] = imax;
01146   stack[stackp + 39] = bfun4_(y, x[imitte]);
01147   stack[stackp + 59] = bmax;
01148   ++stackp;
01149   stack[stackp - 1] = imin;
01150   stack[stackp + 19] = imitte;
01151   stack[stackp + 39] = bmin;
01152   stack[stackp + 59] = bfun4_(y, x[imitte-1]);
01153  L3:
01154   if (stackp > 0) {
01155     goto L2;
01156   }
01157   /* ---- Region 4 */
01158   /* -------------------------------------------------------------------- */
01159   if (lauf[9] >= lauf[4] || lauf[19] >= lauf[14]) {
01160     if ((r__1 = y - yps4, fabs(r__1)) > 1e-8f) {
01161       yps4 = y;
01162       a7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
01163                     y2 * (y2 * (y2 * (2.35944f - y2 * .56419f) - 72.9359f) + 
01164                     571.687f) - 5860.68f) + 40649.2f) - 320772.f) + 1684100.f)
01165                      - 9694630.f) + 40816800.f) - 1.53575e8f) + 4.56662e8f) - 
01166                     9.86604e8f) + 1.16028e9f);
01167       b7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
01168                     y2 * (y2 * (23.0312f - y2 * 7.33447f) - 234.143f) - 
01169                     2269.19f) + 45251.3f) - 234417.f) + 3599150.f) - 
01170                     7723590.f) + 86482900.f) - 2.91876e8f) + 8.06985e8f) - 
01171                     9.85386e8f) - 5.60505e8f);
01172       c7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
01173                     y2 * (97.6203f - y2 * 44.0068f) + 1097.77f) - 25338.3f) + 
01174                     98079.1f) + 576054.f) - 2.3818e7f) + 22930200.f) - 
01175                     2.04467e8f) + 2.94262e8f) + 2.47157e8f) - 6.51523e8f);
01176       d7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
01177                     228.563f - y2 * 161.358f) + 8381.97f) - 66431.2f) - 
01178                     303569.f) + 2240400.f) + 38311200.f) - 41501300.f) - 
01179                     99622400.f) + 2.70167e8f) - 2.63894e8f);
01180       e7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
01181                     296.38f - y2 * 403.396f) + 23507.6f) - 66212.1f) - 
01182                     1.003e6f) + 468142.f) + 24620100.f) + 5569650.f) + 
01183                     1.40677e8f) - 63177100.f);
01184       f7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (125.591f - 
01185                     y2 * 726.113f) + 37544.8f) + 8820.94f) - 934717.f) - 
01186                     1931140.f) - 33289600.f) + 4073820.f) - 16984600.f);
01187       g7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (-260.198f - y2 * 
01188                     968.15f) + 37371.9f) + 79902.5f) - 186682.f) - 900010.f) 
01189                     + 7528830.f) - 1231650.f);
01190       h7 = y * (y2 * (y2 * (y2 * (y2 * (y2 * (-571.645f - y2 * 968.15f)
01191                      + 23137.1f) + 72520.9f) + 153468.f) + 86407.6f) - 
01192                     610622.f);
01193       o7 = y * (y2 * (y2 * (y2 * (y2 * (-575.164f - y2 * 726.113f) + 
01194                     8073.15f) + 26538.5f) + 49883.8f) - 23586.5f);
01195       p7 = y * (y2 * (y2 * (y2 * (-352.467f - y2 * 403.396f) + 
01196                     953.655f) + 2198.86f) - 8009.1f);
01197       q7 = y * (y2 * (y2 * (-134.792f - y2 * 161.358f) - 271.202f) - 
01198                     622.056f);
01199       r7 = y * (y2 * (-29.7896f - y2 * 44.0068f) - 77.0535f);
01200       s7 = y * (-2.92264f - y2 * 7.33447f);
01201       t7 = y * -.56419f;
01202       a8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
01203                     y2 * (y2 * (y2 * (y2 - 3.68288f) + 126.532f) - 955.194f) 
01204                     + 9504.65f) - 70946.1f) + 483737.f) - 2857210.f) + 
01205                     14464700.f) - 61114800.f) + 2.11107e8f) - 5.79099e8f) + 
01206                     1.17022e9f) - 1.5599e9f) + 1.02827e9f;
01207       b8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
01208                     y2 * (y2 * (y2 * 14.f - 40.5117f) + 533.254f) + 3058.26f) 
01209                     - 55600.f) + 498334.f) - 2849540.f) + 13946500.f) - 
01210                     70135800.f) + 2.89676e8f) - 7.53828e8f) + 1.66421e9f) - 
01211                     2.28855e9f) + 1.5599e9f;
01212       c8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
01213                     y2 * (y2 * 91 - 198.876f) - 1500.17f) + 48153.3f) - 
01214                     217801.f) - 1063520.f) + 1.4841e7f) - 46039600.f) + 
01215                     63349600.f) - 6.60078e8f) + 1.06002e9f) - 1.66421e9f) + 
01216                     1.17022e9f;
01217       d8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (
01218                     y2 * 364 - 567.164f) - 16493.7f) + 161461.f) + 280428.f) 
01219                     - 6890020.f) - 6876560.f) + 1.99846e8f) + 54036700.f) + 
01220                     6.60078e8f) - 7.53828e8f) + 5.79099e8f;
01221       e8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 
01222                     1001 - 1012.79f) - 55582.f) + 240373.f) + 1954700.f) - 
01223                     5257220.f) - 50101700.f) - 1.99846e8f) + 63349600.f) - 
01224                     2.89676e8f) + 2.11107e8f;
01225       f8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 2002 - 
01226                     1093.82f) - 106663.f) + 123052.f) + 3043160.f) + 
01227                     5257220.f) - 6876560.f) + 46039600.f) - 70135800.f) + 
01228                     61114800.f;
01229       g8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 3003 - 
01230                     486.14f) - 131337.f) - 123052.f) + 1954700.f) + 6890020.f)
01231                      + 1.4841e7f) - 13946500.f) + 14464700.f;
01232       h8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 3432 + 486.14f) - 
01233                     106663.f) - 240373.f) + 280428.f) + 1063520.f) - 
01234                     2849540.f) + 2857210.f;
01235       o8 = y2 * (y2 * (y2 * (y2 * (y2 * (y2 * 3003 + 1093.82f) - 
01236                     55582.f) - 161461.f) - 217801.f) - 498334.f) + 483737.f;
01237       p8 = y2 * (y2 * (y2 * (y2 * (y2 * 2002 + 1012.79f) - 16493.7f) - 
01238                     48153.3f) - 55600.f) + 70946.1f;
01239       q8 = y2 * (y2 * (y2 * (y2 * 1001.f + 567.164f) - 1500.17f) - 
01240                     3058.26f) + 9504.65f;
01241       r8 = y2 * (y2 * (y2 * 364 + 198.876f) + 533.254f) + 955.194f;
01242       s8 = y2 * (y2 * 91.f + 40.5117f) + 126.532f;
01243       t8 = y2 * 14.f + 3.68288f;
01244     }
01245     ym2 = y * 2;
01246     for (i2 = 1; i2 <= 3; i2 += 2) {
01247       i__1 = lauf[(i2 + 1) * 5 - 1];
01248       for (i1 = lauf[i2 * 5 - 1]; i1 <= i__1; ++i1) {
01249         x2 = x[i1-1] * x[i1-1];
01250         ls[i1-1] = fac * (exp(y2 - x2) * cos(x[i1-1] * ym2) - (a7 + x2 *
01251                          (b7 + x2 * (c7 + x2 * (d7 + x2 * (e7 + x2 * (f7 + x2 
01252                         * (g7 + x2 * (h7 + x2 * (o7 + x2 * (p7 + x2 * (q7 + 
01253                         x2 * (r7 + x2 * (s7 + x2 * t7))))))))))))) / (a8 + x2 
01254                         * (b8 + x2 * (c8 + x2 * (d8 + x2 * (e8 + x2 * (f8 + 
01255                         x2 * (g8 + x2 * (h8 + x2 * (o8 + x2 * (p8 + x2 * (q8 
01256                         + x2 * (r8 + x2 * (s8 + x2 * (t8 + x2)))))))))))))));
01257         /* L4: */
01258       }
01259     }
01260   }
01261   /* ---- Region 3 */
01262   /* -------------------------------------------------------------------- */
01263   if (lauf[8] >= lauf[3] || lauf[18] >= lauf[13]) {
01264     if ((r__1 = y - yps3, fabs(r__1)) > 1e-8f) {
01265       yps3 = y;
01266       a5 = y * (y * (y * (y * (y * (y * (y * (y * (y * 
01267                     .564224f + 7.55895f) + 49.5213f) + 204.501f) + 581.746f) 
01268                     + 1174.8f) + 1678.33f) + 1629.76f) + 973.778f) + 272.102f;
01269       b5 = y * (y * (y * (y * (y * (y * (y * 2.25689f + 22.6778f)
01270                      + 100.705f) + 247.198f) + 336.364f) + 220.843f) - 
01271                     2.34403f) - 60.5644f;
01272       c5 = y * (y * (y * (y * (y * 3.38534f + 22.6798f) + 52.8454f)
01273                      + 42.5683f) + 18.546f) + 4.58029f;
01274       d5 = y * (y * (y * 2.25689f + 7.56186f) + 1.66203f) - .128922f;
01275       e5 = y * .564224f + 9.71457e-4f;
01276       a6 = y * (y * (y * (y * (y * (y * (y * (y * (y * (y + 
01277                     13.3988f) + 88.2674f) + 369.199f) + 1074.41f) + 2256.98f) 
01278                     + 3447.63f) + 3764.97f) + 2802.87f) + 1280.83f) + 
01279                     272.102f;
01280       b6 = y * (y * (y * (y * (y * (y * (y * (y * 5.f + 
01281                     53.5952f) + 266.299f) + 793.427f) + 1549.68f) + 2037.31f) 
01282                     + 1758.34f) + 902.306f) + 211.678f;
01283       c6 = y * (y * (y * (y * (y * (y * 10.f + 80.3928f) + 
01284                     269.292f) + 479.258f) + 497.302f) + 308.186f) + 78.866f;
01285       d6 = y * (y * (y * (y * 10.f + 53.5952f) + 92.7568f) + 
01286                     55.0293f) + 22.0353f;
01287       e6 = y * (y * 5.f + 13.3988f) + 1.49645f;
01288     }
01289     for (i2 = 1; i2 <= 3; i2 += 2) {
01290       i__1 = lauf[(i2 + 1) * 5 - 2];
01291       for (i1 = lauf[i2 * 5 - 2]; i1 <= i__1; ++i1) {
01292         x2 = x[i1-1] * x[i1-1];
01293         ls[i1-1] = fac * (a5 + x2 * (b5 + x2 * (c5 + x2 * (d5 + x2 * 
01294                         e5)))) / (a6 + x2 * (b6 + x2 * (c6 + x2 * (d6 + x2 * (
01295                         e6 + x2)))));
01296         /* L5: */
01297       }
01298     }
01299   }
01300   /* ---- Region 2 */
01301   /* -------------------------------------------------------------------- */
01302   if (lauf[7] >= lauf[2] || lauf[17] >= lauf[12]) {
01303     if ((r__1 = y - yps2, fabs(r__1)) > 1e-8f) {
01304       yps2 = y;
01305       a3 = y * (y2 * (y2 * (y2 * .56419f + 3.10304f) + 4.65456f) + 
01306                     1.05786f);
01307       b3 = y * (y2 * (y2 * 1.69257f + .56419f) + 2.962f);
01308       c3 = y * (y2 * 1.69257f - 2.53885f);
01309       d3 = y * .56419f;
01310       a4 = y2 * (y2 * (y2 * (y2 + 6.f) + 10.5f) + 4.5f) + .5625f;
01311       b4 = y2 * (y2 * (y2 * 4.f + 6.f) + 9.f) - 4.5f;
01312       c4 = y2 * (y2 * 6.f - 6.f) + 10.5f;
01313       d4 = y2 * 4.f - 6.f;
01314     }
01315     for (i2 = 1; i2 <= 3; i2 += 2) {
01316       i__1 = lauf[(i2 + 1) * 5 - 3];
01317       for (i1 = lauf[i2 * 5 - 3]; i1 <= i__1; ++i1) {
01318         x2 = x[i1-1] * x[i1-1];
01319         ls[i1-1] = fac * (a3 + x2 * (b3 + x2 * (c3 + x2 * d3))) / (a4 
01320                         + x2 * (b4 + x2 * (c4 + x2 * (d4 + x2))));
01321         /* L6: */
01322       }
01323     }
01324   }
01325   /* ---- Region 1 */
01326   /* -------------------------------------------------------------------- */
01327   if (lauf[6] >= lauf[1] || lauf[16] >= lauf[11]) {
01328     if ((r__1 = y - yps1, fabs(r__1)) > 1e-8f) {
01329       yps1 = y;
01330       a1 = y * .5641896f;
01331       b1 = y2 + .5f;
01332       a2 = y2 * 4;
01333     }
01334 
01335     c1 = a1 * fac;
01336     for (i2 = 1; i2 <= 3; i2 += 2) {
01337       i__1 = lauf[(i2 + 1) * 5 - 4];
01338       for (i1 = lauf[i2 * 5 - 4]; i1 <= i__1; ++i1) {
01339         x2 = x[i1-1] * x[i1-1];
01340         b2 = b1 - x2;
01341         ls[i1-1] = c1 * (b1 + x2) / (b2 * b2 + a2 * x2);
01342         /* L7: */
01343       }
01344     }
01345   }
01346   /* ---- Region 0 (Lorentz) */
01347   /* -------------------------------------------------------------------- */
01348  L8:
01349   if (lauf[5] >= lauf[0] || lauf[15] >= lauf[10]) {
01350     if ((r__1 = y - yps0, fabs(r__1)) > 1e-8f) {
01351       yps0 = y;
01352       a0 = y * .5641896f;
01353     }
01354 
01355     b0 = a0 * fac;
01356     for (i2 = 1; i2 <= 3; i2 += 2) {
01357       i__1 = lauf[(i2 + 1) * 5 - 5];
01358       for (i1 = lauf[i2 * 5 - 5]; i1 <= i__1; ++i1) {
01359         ls[i1-1] = b0 / (x[i1-1] * x[i1-1] + y2);
01360         /* L9: */
01361       }
01362     }
01363   }
01364 }
01365 
01366 
01367 
01412 /***  ROUTINE COMPUTES THE VOIGT FUNCTION Y/PI*INTEGRAL FROM ***/
01413 /***   - TO + INFINITY OF EXP(-T*T)/(Y*Y+(X-T)*(X-T)) DT     ***/
01414 void lineshape_voigt_drayson(Vector&       ls,
01415                              Vector&       x,
01416                              Numeric         f0,
01417                              Numeric       gamma,
01418                              Numeric       sigma,
01419                              VectorView f_mono,
01420                              const Index   nf)
01421 
01422 {
01423   // FIXME: nf is actually redundant. Could be thrown out in the
01424   // future. For now, let's do an assertion that at least it is
01425   // correct: 
01426   assert( nf==f_mono.nelem() );
01427 
01428   // seems not necessary for Doppler correction
01429   //    extern const Numeric SQRT_NAT_LOG_2;
01430 
01431   // PI
01432   extern const Numeric PI;
01433 
01434   // constant sqrt(1/pi)
01435   const Numeric sqrt_invPI =  sqrt(1/PI);
01436 
01437   // constant normalization factor for voigt
01438   Numeric fac = 1.0 / sigma * sqrt_invPI;
01439 
01440       static Numeric B[22+1] = {0.,0.,.7093602e-7};
01441       static Numeric RI[15+1];
01442       const  Numeric XN[15+1] = {0.,10.,9.,8.,8.,7.,6.,5.,4.,3.,3.,3.,3.,3.,3.,3.};
01443       const  Numeric YN[15+1] = {0.,.6,.6,.6,.5,.4,.4,.3,.3,.3,.3,1.,.9,.8,.7,.7};
01444       static Numeric D0[25+1], D1[25+1], D2[25+1], D3[25+1], D4[25+1];
01445       static Numeric HN[25+1];
01446       static Numeric H = .201;
01447       const  Numeric XX[3+1] = {0.,.5246476,1.65068,.7071068};
01448       const  Numeric HH[3+1] = {0.,.2562121,.2588268e-1,.2820948};
01449       const  Numeric NBY2[19+1] = {0.,9.5,9.,8.5,8.,7.5,7.,6.5
01450           ,6.,5.5,5.,4.5,4.,3.5,3.,2.5,2.,1.5,1.,.5};
01451       const  Numeric C[21+1] = {0.,.7093602e-7,-.2518434e-6,.856687e-6,
01452           -.2787638e-5,.866074e-5,-.2565551e-4,.7228775e-4,
01453           -.1933631e-3,.4899520e-3,-.1173267e-2,.2648762e-2,
01454           -.5623190e-2,.1119601e-1,-.2084976e-1,.3621573e-1,
01455           -.5851412e-1,.8770816e-1,-.121664,.15584,-.184,.2};
01456       static Numeric CO;
01457       Numeric U, V, UU, VV, Y2, DX;
01458       int I, J, K, MAX, MIN, N, i1;
01459       static int TRU = 0;
01460 
01461 
01462       // variables needed in original c routine:
01463 
01464       // Ratio of the Lorentz halfwidth to the Doppler halfwidth
01465       Numeric Y = gamma / sigma;
01466 
01467       // frequency in units of Doppler, Note: Drayson accepts only
01468       // positive values
01469       for (i1=0; i1< (int) nf; i1++)
01470         {
01471           x[i1] = fabs( (f_mono[i1] - f0) )/ sigma;
01472         }
01473 
01474 
01475       if (TRU) goto L104;
01476 /***  REGION I. COMPUTE DAWSON'S FUNCTION AT MESH POINTS ***/
01477       TRU = 1;
01478       for (I=1; I<=15; I++)
01479         RI[I] = -I/2.;
01480       for (I=1; I<=25; I++)
01481       {
01482         HN[I] = H*(I-.5);
01483         CO = 4.*HN[I]*HN[I]/25.-2.;
01484         for (J=2; J<=21; J++)
01485           B[J+1] = CO*B[J]-B[J-1]+C[J];
01486         D0[I] = HN[I]*(B[22]-B[21])/5.;
01487         D1[I] = 1.-2.*HN[I]*D0[I];
01488         D2[I] = (HN[I]*D1[I]+D0[I])/RI[2];
01489         D3[I] = (HN[I]*D2[I]+D1[I])/RI[3];
01490         D4[I] = (HN[I]*D3[I]+D2[I])/RI[4];
01491       }
01492 L104: for (K=0; K<(int) nf; K++)
01493   {
01494     if ((x[K]-5.) < 0.) goto L105; else goto L112;
01495   L105: if ((Y-1.) <= 0.) goto L110; else goto L106;
01496   L106: if (x[K] > 1.85*(3.6-Y)) goto L112;
01497       /***  REGION II: CONTINUED FRACTION. COMPUTE NUMBER OF TERMS NEEDED. ***/
01498       if (Y < 1.45) goto L107;
01499       I = (int) (Y+Y);
01500       goto L108;
01501   L107: I = (int) (11.*Y);
01502   L108: J = (int) (x[K]+x[K]+1.85);
01503       MAX = (int) (XN[J]*YN[I]+.46);
01504       MIN = (int) ( (16 < 21-2*MAX) ? 16 : 21-2*MAX );
01505       /***  EVALUATE CONTINUED FRACTION ***/
01506       UU = Y;
01507       VV = x[K];
01508       for (J=MIN; J <= 19; J++)
01509         {
01510           U = NBY2[J]/(UU*UU+VV*VV);
01511           UU = Y+U*UU;
01512           VV = x[K]-U*VV;
01513         }
01514       ls[K] = UU/(UU*UU+VV*VV)/1.772454*fac;
01515       continue;
01516   L110: Y2 = Y*Y;
01517       if (x[K]+Y >= 5.) goto L113;
01518       /***  REGION I. COMPUTE DAWSON'S FUNCTION AT X FROM TAYLOR SERIES. ***/
01519       N = (int) (x[K]/H);
01520       DX = x[K]-HN[N+1];
01521       U = (((D4[N+1]*DX+D3[N+1])*DX+D2[N+1])*DX+D1[N+1])*DX+D0[N+1];
01522       V = 1.-2.*x[K]*U;
01523       /***  TAYLOR SERIES EXPANSION ABOUT Y=0.0 ***/
01524       VV = exp(Y2-x[K]*x[K])*cos(2.*x[K]*Y)/1.128379-Y*V;
01525       UU = -Y;
01526       MAX = (int) (5.+(12.5-x[K])*.8*Y);
01527       for (I=2; I<=MAX; I+=2)
01528         {
01529           U = (x[K]*V+U)/RI[I];
01530           V = (x[K]*U+V)/RI[I+1];
01531           UU = -UU*Y2;
01532           VV = VV+V*UU;
01533         }
01534       ls[K] = 1.128379*VV*fac;
01535       continue;
01536   L112: Y2 = Y*Y;
01537       if (Y < 11.-.6875*x[K]) goto L113;
01538       /***  REGION IIIB: 2-POINT GAUSS-HERMITE QUADRATURE. ***/
01539       U = x[K]-XX[3];
01540       V = x[K]+XX[3];
01541       ls[K] = Y*(HH[3]/(Y2+U*U)+HH[3]/(Y2+V*V))*fac;
01542       continue;
01543       /***  REGION IIIA: 4-POINT GAUSS-HERMITE QUADRATURE. ***/
01544   L113: U = x[K]-XX[1];
01545       V = x[K]+XX[1];
01546       UU = x[K]-XX[2];
01547       VV = x[K]+XX[2];
01548       ls[K] = Y*(HH[1]/(Y2+U*U)+HH[1]/(Y2+V*V)+HH[2]/(Y2+UU*UU)+HH[2]/(Y2+
01549                                                                       VV*VV))*fac;
01550       continue;
01551   }
01552 }
01553 
01554 
01555 
01580 void lineshape_rosenkranz_voigt_kuntz6(Vector&       ls,
01581                                        Vector&       x,
01582                                        Numeric       f0,
01583                                        Numeric       gamma,
01584                                        Numeric       sigma,
01585                                        VectorView f_mono,
01586                                        const Index   nf)
01587 {
01588   // FIXME: nf is actually redundant. Could be thrown out in the
01589   // future. For now, let's do an assertion that at least it is
01590   // correct: 
01591   assert( nf==f_mono.nelem() );
01592 
01593   // seems not necessary for Doppler correction
01594   //    extern const Numeric SQRT_NAT_LOG_2;
01595 
01596   extern const Numeric PI;
01597 
01598 
01599   // calculate the required stuff with the parameter of the x array:
01600   Numeric overlap;      // overlap correction
01601   Numeric gamma_o2;     // linewidth 
01602   Numeric gamma_o2_2;   // square of linewidth
01603   {
01604     // x[0] = theta
01605     // x[1] = theta_Nair
01606     // x[2] = total pressure    
01607     // x[3] = O2 VMR
01608     // x[4] = H2O VMR
01609     // x[5] = l_l.Agam()
01610     // x[6] = l_l.Nair()
01611     // x[7] = l_l.Aux()[0] = overlap y
01612     // x[8] = l_l.Aux()[1] = overlap v
01613 
01614     // Pressure Broadening:
01615     // pressure broadening is calculated differently for oxygen, not
01616     // with the partial pressure, but with the dry pressure. the
01617     // passed gamma parameter is calculated with the partial pressure,
01618     // so we correct it here:
01619     gamma *= ( 1.0 - x[4] ) / ( 1 - x[3] );
01620 
01621     // Overlap:
01622     // Get the overlap correction for oxygen, this is always calculated
01623     // even if not needed. FIXME: I assume that theta is to the power of
01624     // Nair (just as in the pressure broadening), otherwise we need
01625     // another parameter to pass with the catalogue.
01626     //cout << "x7, x8: " << x[7] << " " << x[8] << "\n";
01627     overlap = (x[7] + x[8] * ( x[0] - 1.0)) * 
01628       x[2] * x[1];
01629 
01630     // water vapor impact:
01631     // impact of water vapor upon the oxygen pressure broadening. 
01632     gamma_o2 = gamma + 1.1 * x[5] * x[2] * x[4] * x[0];
01633     gamma_o2_2 = gamma_o2 * gamma_o2;
01634   }
01635 
01636   // use voigt lineshape for high altitudes, otherwise lorentz with
01637   // overlap correction
01638   if ( (gamma_o2/sigma - 40) <= 0.0 )
01639     {
01640       /* call voigt function */
01641       lineshape_voigt_kuntz6(ls,
01642                              x,
01643                              f0,
01644                              gamma_o2,
01645                              sigma,
01646                              f_mono,
01647                              nf);
01648     }
01649   else
01650     {
01651       Numeric FD, FS, SF1, SF2;
01652       for (Index J=0; J<(Index) nf; J++)
01653         {
01654           FD = f_mono[J] - f0;
01655           FS = f_mono[J] + f0;
01656           SF1 = (gamma_o2 + FD*overlap) / (FD*FD + gamma_o2_2);
01657           SF2 = (gamma_o2 - FS*overlap) / (FS*FS + gamma_o2_2);
01658           ls[J] = (SF1 + SF2) / PI;
01659         }
01660     }
01661 }
01662 
01663 
01688 void lineshape_rosenkranz_voigt_drayson(Vector&       ls,
01689                                         Vector&       x,
01690                                         Numeric      f0,
01691                                         Numeric       gamma,
01692                                         Numeric       sigma,
01693                                         VectorView f_mono,
01694                                         const Index   nf)
01695 {
01696   // FIXME: nf is actually redundant. Could be thrown out in the
01697   // future. For now, let's do an assertion that at least it is
01698   // correct: 
01699   assert( nf==f_mono.nelem() );
01700 
01701   // seems not necessary for Doppler correction
01702   //    extern const Numeric SQRT_NAT_LOG_2;
01703 
01704   extern const Numeric PI;
01705 
01706 
01707   // calculate the required stuff with the parameter of the x array:
01708   Numeric overlap;      // overlap correction
01709   Numeric gamma_o2;     // linewidth 
01710   Numeric gamma_o2_2;   // square of linewidth
01711   {
01712     // x[0] = theta
01713     // x[1] = theta_Nair
01714     // x[2] = total pressure    
01715     // x[3] = O2 VMR
01716     // x[4] = H2O VMR
01717     // x[5] = l_l.Agam()
01718     // x[6] = l_l.Nair()
01719     // x[7] = l_l.Aux()[0] = overlap y
01720     // x[8] = l_l.Aux()[1] = overlap v
01721 
01722     // Pressure Broadening:
01723     // pressure broadening is calculated differently for oxygen, not
01724     // with the partial pressure, but with the dry pressure. the
01725     // passed gamma parameter is calculated with the partial pressure,
01726     // so we correct it here:
01727     gamma *= ( 1.0 - x[4] ) / ( 1 - x[3] );
01728 
01729     // Overlap:
01730     // Get the overlap correction for oxygen, this is always calculated
01731     // even if not needed. FIXME: I assume that theta is to the power of
01732     // Nair (just as in the pressure broadening), otherwise we need
01733     // another parameter to pass with the catalogue.
01734     overlap = (x[7] + x[8] * ( x[0] - 1.0)) * 
01735       x[2] * x[1];
01736 
01737     // water vapor impact:
01738     // impact of water vapor upon the oxygen pressure broadening. 
01739     gamma_o2 = gamma + 1.1 * x[5] * x[2] * x[4] * x[0];
01740     gamma_o2_2 = gamma_o2 * gamma_o2;
01741   }
01742 
01743   // use voigt lineshape for high altitudes, otherwise lorentz with
01744   // overlap correction
01745   if ( (gamma_o2/sigma - 40) <= 0.0 )
01746     {
01747       /* call voigt function */
01748       lineshape_voigt_drayson(ls,
01749                               x,
01750                               f0,
01751                               gamma_o2,
01752                               sigma,
01753                               f_mono,
01754                               nf);
01755 
01756       // the controlfile uses generally a (f/f0)^2 factor for the
01757       // lineshape function, but we have to remove this normalization
01758       // factor (f/f0)^2 over here, because Rosenkranz does not use it
01759       // for the Voigt part, but for the Lorentz part. don't ask me
01760       // why...
01761       // FIXME: Clarify whether this is correct and if yes use the
01762       // normalization used in controlfile
01763       Numeric f0_2 = f0 * f0;
01764       for (Index J=0; J<(Index) nf; J++)
01765         {
01766           ls[J] *= f0_2 / (f_mono[J] * f_mono[J]);
01767         }
01768     }
01769   else
01770     {
01771       Numeric FD, FS, SF1, SF2;
01772       for (Index J=0; J<(Index) nf; J++)
01773         {
01774           FD = f_mono[J] - f0;
01775           FS = f_mono[J] + f0;
01776           SF1 = (gamma_o2 + FD*overlap) / (FD*FD + gamma_o2_2);
01777           SF2 = (gamma_o2 - FS*overlap) / (FS*FS + gamma_o2_2);
01778           ls[J] = (SF1 + SF2) / PI;
01779         }
01780     }
01781 }
01782 
01783 
01784 
01796 void chi_cousin(
01797             Numeric&   chi,
01798       const Numeric&   df )
01799 {
01800   // Conversion factor from Hz to cm-1
01801   extern const Numeric HZ2CM;
01802 
01803   const Numeric n2 = 0.79;
01804   const Numeric o2 = 0.21;
01805 
01806   const Numeric   df_cm     = df * HZ2CM;
01807   const Numeric   df_cm_abs = fabs( df_cm );
01808 
01809   chi = 0;  
01810 
01811   // N2 term
01812   if( df_cm_abs <= 5 )
01813     { chi += n2; }
01814   else if( df_cm_abs <= 22 )
01815     { chi += n2 * 1.968 * exp( -0.1354 * df_cm_abs ); }
01816   else if( df_cm_abs <= 50 )
01817     { chi += n2 * 0.160 * exp( -0.0214 * df_cm_abs ); }
01818   else
01819     { chi += n2 * 0.162 * exp( -0.0216 * df_cm_abs ); }
01820 
01821   // O2 term
01822   if( df_cm_abs <= 5 )
01823     { chi += o2; }
01824   else if( df_cm_abs <= 22 )
01825     { chi += o2 * 1.968 * exp( -0.1354 * df_cm_abs ); }
01826   else if( df_cm_abs <= 50 )
01827     { chi += o2 * 0.160 * exp( -0.0214 * df_cm_abs ); }
01828   else
01829     { chi += o2 * 0.162 * exp( -0.0216 * df_cm_abs ); }
01830 }
01831 
01832 
01833 
01846 void lineshape_CO2_lorentz(
01847             Vector&   ls,
01848             Vector&   /* X */,
01849             Numeric   f0,
01850             Numeric   gamma,
01851             Numeric   /* sigma */,
01852       VectorView      f_mono,
01853       const Index     nf )
01854 {
01855   assert( f_mono.nelem() == nf );
01856 
01857   // PI:
01858   extern const Numeric PI;
01859 
01860   // Some constant variables
01861   const Numeric gamma2 = gamma * gamma;
01862   const Numeric fac = gamma/PI;
01863 
01864   for ( Index i=0; i<nf; ++i )
01865     {
01866       // f-f0
01867       const Numeric df = f_mono[i] - f0;
01868 
01869       // The chi factor
01870       Numeric chi;
01871       chi_cousin( chi, df );      
01872 
01873       // chi * Lorentz
01874       ls[i] =  chi * fac / ( df * df + gamma2 );
01875     }
01876 }
01877 
01878 
01879 
01892 void lineshape_CO2_drayson(
01893             Vector&   ls,
01894             Vector&   X,
01895             Numeric   f0,
01896             Numeric   gamma,
01897             Numeric   sigma,
01898       VectorView      f_mono,
01899       const Index     nf )
01900 {
01901   lineshape_voigt_drayson(  ls, X, f0, gamma, sigma, f_mono, nf );
01902 
01903   for ( Index i=0; i<nf; ++i )
01904     {
01905       // f-f0
01906       const Numeric df = f_mono[i] - f0;
01907 
01908       // The chi factor
01909       Numeric chi;
01910       chi_cousin( chi, df );      
01911 
01912       // chi * Lorentz
01913       ls[i] *=  chi;
01914     }
01915 }
01916 
01917 
01918 
01919 
01920 
01921 
01922 //------------------------------------------------------------------------
01923 // Normalization Functions 
01924 //------------------------------------------------------------------------
01925 
01926 
01936 void lineshape_norm_no_norm(Vector&       fac,
01937                             Numeric       /* f0 */,
01938                             VectorView    f_mono,
01939                             const Numeric /* T */,
01940                             const Index   nf)
01941 {
01942   // FIXME: nf is actually redundant. Could be thrown out in the
01943   // future. For now, let's do an assertion that at least it is
01944   // correct: 
01945   assert( nf==f_mono.nelem() );
01946 
01947   for ( Index i=0; i<nf; ++i )
01948     {
01949       fac[i] = 1.0;
01950     }
01951 }
01952 
01953 
01954 
01964 void lineshape_norm_linear(Vector&       fac,
01965                            Numeric       f0,
01966                            VectorView    f_mono,
01967                            const Numeric /* T */,
01968                            const Index   nf)
01969 {
01970   // FIXME: nf is actually redundant. Could be thrown out in the
01971   // future. For now, let's do an assertion that at least it is
01972   // correct: 
01973   assert( nf==f_mono.nelem() );
01974 
01975   for ( Index i=0; i<nf; ++i )
01976     {
01977       fac[i] = f_mono[i] / f0;
01978     }
01979 }
01980 
01991 void lineshape_norm_quadratic(Vector&       fac,
01992                               Numeric       f0,
01993                               VectorView    f_mono,
01994                               const Numeric /* T */,
01995                               const Index   nf)
01996 {
01997   // FIXME: nf is actually redundant. Could be thrown out in the
01998   // future. For now, let's do an assertion that at least it is
01999   // correct: 
02000   assert( nf==f_mono.nelem() );
02001 
02002   // don't do this for the whole loop
02003   Numeric f0_2 = f0 * f0;
02004 
02005   for ( Index i=0; i<nf; ++i )
02006     {
02007       fac[i] = (f_mono[i] * f_mono[i]) / f0_2;
02008     }
02009 }
02010 
02027 void lineshape_norm_VVH(Vector&       fac,
02028                         Numeric       f0,
02029                         VectorView    f_mono,
02030                         const Numeric T,
02031                         const Index   nf)
02032 {
02033   extern const Numeric PLANCK_CONST;
02034   extern const Numeric BOLTZMAN_CONST;
02035 
02036   // FIXME: nf is actually redundant. Could be thrown out in the
02037   // future. For now, let's do an assertion that at least it is
02038   // correct: 
02039   assert( nf==f_mono.nelem() );
02040 
02041   // 2kT is constant for the loop
02042   Numeric kT = 2.0 * BOLTZMAN_CONST * T;
02043 
02044   // denominator is constant for the loop
02045   Numeric denom = f0 * tanh( PLANCK_CONST * f0 / kT );
02046 
02047   for ( Index i=0; i<nf; ++i )
02048     {
02049       fac[i] = f_mono[i] * tanh( PLANCK_CONST * f_mono[i] / kT ) /
02050         denom;
02051     }
02052 }
02053 
02054 
02055 
02056 
02057 //------------------------------------------------------------------------
02058 // Available Lineshapes
02059 //------------------------------------------------------------------------
02060 
02062 Array<LineshapeRecord> lineshape_data;
02063 
02064 void define_lineshape_data()
02065 {
02066   // Initialize to empty, just in case.
02067   lineshape_data.resize(0);
02068 
02069   lineshape_data.push_back
02070     (LineshapeRecord
02071      ("no_shape",
02072       "This lineshape does nothing. It only exists, because formally\n"
02073       "you have to specify a lineshape also for continuum tags.", 
02074       lineshape_no_shape));
02075 
02076   lineshape_data.push_back
02077     (LineshapeRecord
02078      ("Lorentz",
02079       "The Lorentz line shape.",
02080       lineshape_lorentz));
02081 
02082   lineshape_data.push_back
02083     (LineshapeRecord
02084      ("Doppler",
02085       "The Doppler line shape.",
02086       lineshape_doppler));
02087 
02088   lineshape_data.push_back
02089     (LineshapeRecord
02090      ("Voigt_Kuntz6",
02091       "The Voigt line shape. Approximation by Kuntz: Accuracy 2*10-6",
02092       lineshape_voigt_kuntz6));
02093 
02094   lineshape_data.push_back
02095     (LineshapeRecord
02096      ("Voigt_Kuntz3",
02097       "The Voigt line shape. Approximation by Kuntz: Accuracy 2*10-3",
02098       lineshape_voigt_kuntz3));
02099 
02100   lineshape_data.push_back
02101     (LineshapeRecord
02102      ("Voigt_Kuntz4",
02103       "The Voigt line shape. Approximation by Kuntz: Accuracy 2*10-4",
02104       lineshape_voigt_kuntz4));
02105 
02106   lineshape_data.push_back
02107     (LineshapeRecord
02108      ("Voigt_Drayson",
02109       "The Voigt line shape. Approximation by Drayson.",
02110       lineshape_voigt_drayson));
02111 
02112   lineshape_data.push_back
02113     (LineshapeRecord
02114      ("Rosenkranz_Voigt_Kuntz6",
02115       "Rosenkranz lineshape for oxygen with overlap correction, "
02116       "at high altitudes Voigt_Kuntz6.",
02117       lineshape_rosenkranz_voigt_kuntz6));
02118 
02119   lineshape_data.push_back
02120     (LineshapeRecord
02121      ("Rosenkranz_Voigt_Drayson",
02122       "Rosenkranz lineshape for oxygen with overlap correction, "
02123       "at high altitudes Drayson.",
02124       lineshape_rosenkranz_voigt_drayson));
02125 
02126   lineshape_data.push_back
02127     (LineshapeRecord
02128      ("CO2_Lorentz",
02129       "Special line shape for CO2 in the infrared, neglecting Doppler\n"
02130       "broadening and details of line mixing. The line shape can be\n"
02131       "expressed as\n"
02132       "   chi(f,f0) * Lorentz(f,f0) \n"
02133       "\n"
02134       "The chi-factor follows Cousin et al. 1985. Broadening by N2 and O2\n"
02135       "is considered, while self-broadening (CO2-CO2) is neglected."
02136       "\n"
02137       "NOTE: Temperature dependency is not yet included. The chi factor is\n"
02138       "valid for 238 K.",
02139       lineshape_CO2_lorentz));
02140 
02141   lineshape_data.push_back
02142     (LineshapeRecord
02143      ("CO2_Drayson",
02144       "Special line shape for CO2 in the infrared, neglecting details of\n"
02145       "line mixing. The line shape can be expressed as\n"
02146       "   chi(f,f0) * Drayson(f,f0) \n"
02147       "\n"
02148       "The chi-factor follows Cousin et al. 1985. Broadening by N2 and O2\n"
02149       "is considered, while self-broadening (CO2-CO2) is neglected.\n"
02150       "\n"
02151       "NOTE: Temperature dependency is not yet included. The chi factor is\n"
02152       "valid for 238 K.",
02153       lineshape_CO2_drayson));
02154 
02155 }
02156 
02159 Array<LineshapeNormRecord> lineshape_norm_data;
02160 
02161 void define_lineshape_norm_data()
02162 {
02163   // Initialize to empty, just in case.
02164   lineshape_norm_data.resize(0);
02165 
02166   lineshape_norm_data.push_back
02167     (LineshapeNormRecord
02168      ("no_norm",
02169       "No normalization of the lineshape.",
02170       lineshape_norm_no_norm));
02171 
02172   lineshape_norm_data.push_back
02173     (LineshapeNormRecord
02174      ("linear",
02175       "Linear normalization of the lineshape with f/f0.",
02176       lineshape_norm_linear));
02177 
02178   lineshape_norm_data.push_back
02179     (LineshapeNormRecord
02180      ("quadratic",
02181       "Quadratic normalization of the lineshape with (f/f0)^2.",
02182       lineshape_norm_quadratic));
02183 
02184   lineshape_norm_data.push_back
02185     (LineshapeNormRecord
02186      ("VVH",
02187       "Van Vleck Huber normalization of the lineshape with\n"
02188       "             (f*tanh(h*f/(2*k*T))) / (f0*tanh(h*f0/(2*k*T))).\n"
02189       "             The denominator is a result of catalogue intensities.",
02190       lineshape_norm_VVH));
02191 }