sln_faddeeva.c (5226B)
1 /* Copyright (C) 2022, 2026 |Méso|Star> (contact@meso-star.com) 2 * Copyright (C) 2026 Université de Lorraine 3 * Copyright (C) 2022 Centre National de la Recherche Scientifique 4 * Copyright (C) 2022 Université Paul Sabatier 5 * 6 * This program is free software: you can redistribute it and/or modify 7 * it under the terms of the GNU General Public License as published by 8 * the Free Software Foundation, either version 3 of the License, or 9 * (at your option) any later version. 10 * 11 * This program is distributed in the hope that it will be useful, 12 * but WITHOUT ANY WARRANTY; without even the implied warranty of 13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 14 * GNU General Public License for more details. 15 * 16 * You should have received a copy of the GNU General Public License 17 * along with this program. If not, see <http://www.gnu.org/licenses/>. */ 18 19 #include "sln.h" 20 21 /******************************************************************************* 22 * Exported function 23 ******************************************************************************/ 24 double 25 sln_faddeeva(const double x, const double Y) 26 { 27 /* Constants */ 28 const double RRTPI = 0.56418958; 29 30 /* For CPF12 algorithm */ 31 const double Y0 = 1.5; 32 const double Y0PY0 = Y0+Y0; 33 const double Y0Q = Y0*Y0; 34 35 const double C[6] = { 36 1.0117281, 37 -0.75197147, 38 0.012557727, 39 0.010022008, 40 -0.00024206814, 41 0.00000050084806 42 }; 43 const double S[6] = { 44 1.393237, 45 0.23115241, 46 -0.15535147, 47 0.0062183662, 48 0.000091908299, 49 -0.00000062752596 50 }; 51 const double T[6] = { 52 0.31424038, 53 0.94778839, 54 1.5976826, 55 2.2795071, 56 3.0206370, 57 3.8897249 58 }; 59 60 double ABX, XQ, YQ, YRRTPI; 61 double XLIM0, XLIM1, XLIM2, XLIM3, XLIM4; 62 double A0, D0, D2, E0, E2, E4, H0, H2, H4, H6; 63 double P0, P2, P4, P6, P8, Z0, Z2, Z4, Z6, Z8; 64 double XP[6], XM[6], YP[6], YM[6]; 65 double MQ[6], PQ[6], MF[6], PF[6]; 66 double D, YF, YPY0, YPY0Q; 67 68 /* Output */ 69 double k = 0; 70 71 int i; 72 int RG1, RG2, RG3; 73 74 YQ = Y*Y; 75 YRRTPI = Y*RRTPI; 76 77 if(Y >= 70.55) { 78 XQ = x*x; 79 k=YRRTPI/(XQ+YQ); 80 return k; 81 } 82 83 RG1 = RG2 = RG3 = 1; 84 85 XLIM0 = sqrt(15100.0 + Y*(40.0 - Y*3.6)); 86 if(Y >= 8.425){ 87 XLIM1 = 0.0; 88 } else { 89 XLIM1 = sqrt(164.0 - Y*(4.3 + Y*1.8)); 90 } 91 92 XLIM2 = 6.8-Y; 93 XLIM3 = 2.4*Y; 94 XLIM4 = 18.1*Y+1.65; 95 96 if(Y<=1.0e-6) { 97 XLIM1=XLIM0; 98 XLIM2=XLIM0; 99 } 100 101 ABX = sqrt(x*x); 102 XQ = ABX*ABX; 103 if(ABX >= XLIM0) { 104 k = YRRTPI/(XQ + YQ); 105 } else if (ABX >= XLIM1) { 106 if(RG1 != 0) { 107 RG1 = 0; 108 A0 = YQ + 0.5; 109 D0 = A0*A0; 110 D2 = YQ + YQ - 1.0; 111 } 112 D = RRTPI/(D0 + XQ*(D2 + XQ)); 113 k = D*Y*(A0 + XQ); 114 115 } else if (ABX > XLIM2) { 116 if(RG2 != 0) { 117 RG2 = 0; 118 H0 = 0.5625 + YQ*(4.5 + YQ*(10.5 + YQ*(6.0 + YQ))); 119 H2 = -4.5 + YQ*(9.0 + YQ*(6.0 + YQ*4.0)); 120 H4 = 10.5 - YQ*(6.0 - YQ*6.0); 121 H6 = -6.0 + YQ*4.0; 122 E0 = 1.875 + YQ*(8.25 + YQ*(5.5+YQ)); 123 E2 = 5.25 + YQ*(1.0 + YQ*3.0); 124 E4 = 0.75*H6; 125 } 126 D = RRTPI/(H0 + XQ*(H2 + XQ*(H4 + XQ*(H6 + XQ)))); 127 k = D*Y*(E0 + XQ*(E2 + XQ*(E4 + XQ))); 128 129 } else if (ABX<XLIM3) { 130 if(RG3 != 0) { 131 RG3 = 0; 132 Z0 = 272.1014 133 + Y*(1280.829 + Y*(2802.870 + Y*(3764.966 134 + Y*(3447.629 + Y*(2256.981 + Y*(1074.409 135 + Y*(369.1989 + Y*(88.26741 + Y*(13.39880 + Y))))))))); 136 Z2 = 211.678 137 + Y*(902.3066 + Y*(1758.336 + Y*(2037.310 138 + Y*(1549.675 + Y*(793.4273 + Y*(266.2987 139 + Y*(53.59518 + Y*5.0))))))); 140 Z4 = 78.86585 141 + Y*(308.1852 + Y*(497.3014 + Y*(479.2576 142 + Y*(269.2916 + Y*(80.39278 + Y*10.0))))); 143 Z6 = 22.03523 144 + Y*(55.02933 + Y*(92.75679 + Y*(53.59518 + Y*10.0))); 145 Z8 = 1.496460 146 + Y*(13.39880 + Y*5.0); 147 P0 = 153.5168 148 + Y*(549.3954 + Y*(919.4955 + Y*(946.8970 149 + Y*(662.8097 + Y*(328.2151 + Y*(115.3772 150 + Y*(27.93941 + Y*(4.264678 + Y*0.3183291)))))))); 151 P2 = -34.16955 152 + Y*(-1.322256+ Y*(124.5975 + Y*(189.7730 153 + Y*(139.4665 + Y*(56.81652 + Y*(12.79458 + Y*1.2733163)))))); 154 P4 = 2.584042 155 + Y*(10.46332 + Y*(24.01655 + Y*(29.81482 156 + Y*(12.79568 + Y*1.9099744)))); 157 P6 = -0.07272979 158 + Y*(0.9377051 + Y*(4.266322 + Y*1.273316)); 159 P8 = 0.0005480304 160 + Y*0.3183291; 161 } 162 D = 1.7724538/(Z0 + XQ*(Z2 + XQ*(Z4 + XQ*(Z6 + XQ*(Z8 + XQ))))); 163 k = D*(P0 + XQ*(P2 + XQ*(P4 + XQ*(P6 + XQ*P8)))); 164 165 } else { 166 YPY0 = Y+Y0; 167 YPY0Q = YPY0*YPY0; 168 k=0.0; 169 170 FOR_EACH(i, 0, 6) { 171 D = x - T[i]; 172 MQ[i] = D*D; 173 MF[i] = 1.0/(MQ[i] + YPY0Q); 174 XM[i] = MF[i]*D; 175 YM[i] = MF[i]*YPY0; 176 177 D = x + T[i]; 178 PQ[i] = D*D; 179 PF[i] = 1.0/(PQ[i] + YPY0Q); 180 XP[i] = PF[i]*D; 181 YP[i] = PF[i]*YPY0; 182 } 183 184 if(ABX <= XLIM4) { 185 FOR_EACH(i, 0, 6) { 186 k = k + C[i]*(YM[i] + YP[i]) - S[i]*(XM[i] - XP[i]); 187 } 188 } else { 189 YF = Y+Y0PY0; 190 FOR_EACH(i, 0, 6) { 191 k = k 192 + (C[i]*(MQ[i]*MF[i] - Y0*YM[i]) + S[i]*YF*XM[i])/(MQ[i] + Y0Q) 193 + (C[i]*(PQ[i]*PF[i] - Y0*YP[i]) - S[i]*YF*XP[i])/(PQ[i] + Y0Q); 194 } 195 k = Y*k + exp(-XQ); 196 } 197 } 198 return k; 199 }