star-line

sln_faddeeva.c (5264B)

---

 /* Copyright (C) 2022, 2026 |Méso|Star> (contact@meso-star.com)
  * Copyright (C) 2026 Université de Lorraine
  * Copyright (C) 2022 Centre National de la Recherche Scientifique
  * Copyright (C) 2022 Université Paul Sabatier
  *
  * This file is part of Star-Line.
  *
  * This program is free software: you can redistribute it and/or modify
  * it under the terms of the GNU General Public License as published by
  * the Free Software Foundation, either version 3 of the License, or
  * (at your option) any later version.
  *
  * This program is distributed in the hope that it will be useful,
  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  * GNU General Public License for more details.
  *
  * You should have received a copy of the GNU General Public License
  * along with this program. If not, see <http://www.gnu.org/licenses/>. */
 
 #include "sln.h"
 
 /*******************************************************************************
  * Exported function
  ******************************************************************************/
 double
 sln_faddeeva(const double x, const double Y)
 {
   /* Constants */
   const double RRTPI = 0.56418958;
 
   /* For CPF12 algorithm */
   const double Y0 = 1.5;
   const double Y0PY0 = Y0+Y0;
   const double Y0Q = Y0*Y0;
 
   const double C[6] = {
     1.0117281,
     -0.75197147,
     0.012557727,
     0.010022008,
     -0.00024206814,
     0.00000050084806
   };
   const double S[6] = {
     1.393237,
     0.23115241,
     -0.15535147,
     0.0062183662,
     0.000091908299,
     -0.00000062752596
   };
   const double T[6] = {
     0.31424038,
     0.94778839,
     1.5976826,
     2.2795071,
     3.0206370,
     3.8897249
   };
 
   double ABX, XQ, YQ, YRRTPI;
   double XLIM0, XLIM1, XLIM2, XLIM3, XLIM4;
   double A0, D0, D2, E0, E2, E4, H0, H2, H4, H6;
   double P0, P2, P4, P6, P8, Z0, Z2, Z4, Z6, Z8;
   double XP[6], XM[6], YP[6], YM[6];
   double MQ[6], PQ[6], MF[6], PF[6];
   double D, YF, YPY0, YPY0Q;
 
   /* Output */
   double k = 0;
 
   int i;
   int RG1, RG2, RG3;
 
   YQ = Y*Y;
   YRRTPI = Y*RRTPI;
 
   if(Y >= 70.55) {
     XQ = x*x;
     k=YRRTPI/(XQ+YQ);
     return k;
   }
 
   RG1 = RG2 = RG3 = 1;
 
   XLIM0 = sqrt(15100.0 + Y*(40.0 - Y*3.6));
   if(Y >= 8.425){
     XLIM1 = 0.0;
   } else {
     XLIM1 = sqrt(164.0 - Y*(4.3 + Y*1.8));
   }
 
   XLIM2 = 6.8-Y;
   XLIM3 = 2.4*Y;
   XLIM4 = 18.1*Y+1.65;
 
   if(Y<=1.0e-6) {
     XLIM1=XLIM0;
     XLIM2=XLIM0;
   }
 
   ABX = sqrt(x*x);
   XQ = ABX*ABX;
   if(ABX >= XLIM0) {
     k = YRRTPI/(XQ + YQ);
   } else if (ABX >= XLIM1) {
     if(RG1 != 0) {
       RG1 = 0;
       A0 = YQ + 0.5;
       D0 = A0*A0;
       D2 = YQ + YQ - 1.0;
     }
     D = RRTPI/(D0 + XQ*(D2 + XQ));
     k = D*Y*(A0 + XQ);
 
   } else if (ABX > XLIM2) {
     if(RG2 != 0) {
       RG2 = 0;
       H0 = 0.5625 + YQ*(4.5 + YQ*(10.5 + YQ*(6.0 + YQ)));
       H2 = -4.5 + YQ*(9.0 + YQ*(6.0 + YQ*4.0));
       H4 = 10.5 - YQ*(6.0 - YQ*6.0);
       H6 = -6.0 + YQ*4.0;
       E0 = 1.875 + YQ*(8.25 + YQ*(5.5+YQ));
       E2 = 5.25 + YQ*(1.0 + YQ*3.0);
       E4 = 0.75*H6;
     }
     D = RRTPI/(H0 + XQ*(H2 + XQ*(H4 + XQ*(H6 + XQ))));
     k = D*Y*(E0 + XQ*(E2 + XQ*(E4 + XQ)));
 
   } else if (ABX<XLIM3) {
     if(RG3 != 0) {
       RG3 = 0;
       Z0 = 272.1014
          + Y*(1280.829 + Y*(2802.870 + Y*(3764.966
          + Y*(3447.629 + Y*(2256.981 + Y*(1074.409
          + Y*(369.1989 + Y*(88.26741 + Y*(13.39880 + Y)))))))));
       Z2 = 211.678
          + Y*(902.3066 + Y*(1758.336 + Y*(2037.310
          + Y*(1549.675 + Y*(793.4273 + Y*(266.2987
          + Y*(53.59518 + Y*5.0)))))));
       Z4 = 78.86585
          + Y*(308.1852 + Y*(497.3014 + Y*(479.2576
          + Y*(269.2916 + Y*(80.39278 + Y*10.0)))));
       Z6 = 22.03523
          + Y*(55.02933 + Y*(92.75679 + Y*(53.59518 + Y*10.0)));
       Z8 = 1.496460
          + Y*(13.39880 + Y*5.0);
       P0 = 153.5168
          + Y*(549.3954 + Y*(919.4955 + Y*(946.8970
          + Y*(662.8097 + Y*(328.2151 + Y*(115.3772
          + Y*(27.93941 + Y*(4.264678 + Y*0.3183291))))))));
       P2 = -34.16955
          + Y*(-1.322256+ Y*(124.5975 + Y*(189.7730
          + Y*(139.4665 + Y*(56.81652 + Y*(12.79458 + Y*1.2733163))))));
       P4 = 2.584042
          + Y*(10.46332 + Y*(24.01655 + Y*(29.81482
          + Y*(12.79568 + Y*1.9099744))));
       P6 = -0.07272979
          + Y*(0.9377051 + Y*(4.266322 + Y*1.273316));
       P8 = 0.0005480304
          + Y*0.3183291;
     }
     D = 1.7724538/(Z0 + XQ*(Z2 + XQ*(Z4 + XQ*(Z6 + XQ*(Z8 + XQ)))));
     k = D*(P0 + XQ*(P2 + XQ*(P4 + XQ*(P6 + XQ*P8))));
 
   } else {
     YPY0 = Y+Y0;
     YPY0Q = YPY0*YPY0;
     k=0.0;
 
     FOR_EACH(i, 0, 6) {
       D = x - T[i];
       MQ[i] = D*D;
       MF[i] = 1.0/(MQ[i] + YPY0Q);
       XM[i] = MF[i]*D;
       YM[i] = MF[i]*YPY0;
 
       D = x + T[i];
       PQ[i] = D*D;
       PF[i] = 1.0/(PQ[i] + YPY0Q);
       XP[i] = PF[i]*D;
       YP[i] = PF[i]*YPY0;
     }
 
     if(ABX <= XLIM4) {
       FOR_EACH(i, 0, 6) {
         k = k + C[i]*(YM[i] + YP[i]) - S[i]*(XM[i] - XP[i]);
       }
     } else {
       YF = Y+Y0PY0;
       FOR_EACH(i, 0, 6) {
         k = k
           + (C[i]*(MQ[i]*MF[i] - Y0*YM[i]) + S[i]*YF*XM[i])/(MQ[i] + Y0Q)
           + (C[i]*(PQ[i]*PF[i] - Y0*YP[i]) - S[i]*YF*XP[i])/(PQ[i] + Y0Q);
       }
       k = Y*k + exp(-XQ);
     }
   }
   return k;
 }