c kspectrum (http://www.meso-star.com/en_Products.html) - This file is part of kspectrum c Copyright (C) 2008-2015 - Méso-Star - Vincent Eymet c c This file must be used under the terms of the CeCILL license. c This source file is licensed as described in the file COPYING, which c you should have received as part of this distribution. The terms c are also available at c http://www.cecill.info/licences/Licence_CeCILL_V2.1-en.txt c subroutine Lorentz(nu,nu0,P,gamma,delta_air,f) implicit none c c Purpose: to compute lineshape at a given frequency, according to Lorentz's profile c c Inputs: c + nu: frequency at which Lorentz's profile has to be computed [cm¯¹] c + nu0: vacuum wavenumber of the given line [cm¯¹] c + P: total pressure [atm] c + gamma: half-width [cm¯¹] c + delta_air: air pressure-induced line shift [cm¯¹.atm¯¹] c c Ouput: c + f: value of Lorentz's profile c include 'parameters.inc' double precision gamma,nu,nu0,delta_air,P double precision f f=(1.0D+0/pi) & *gamma/(gamma**2+(nu-(nu0+delta_air*P))**2) return end subroutine Voigt(x,y,k) implicit none c c Purpose: to calculate the Faddeeva (Voigt) function with relative error less than 10^(-4). c c Inputs: x and y, parameters for the Voigt function c - x is defined as x=(nu-nu_c)/gamma_D*sqrt(ln(2)) with nu the current wavenumber, c nu_c the wavenumber at line center, gamma_D the Doppler linewidth. c - y is defined as y=gamma_L/gamma_D*sqrt(ln(2)) with gamma_L the Lorentz linewith c and gamma_D the Doppler linewidth c c Output: k, the Voigt function; it has to be multiplied by sqrt(ln(2)/pi)*1/gamma_D c so that the result may be interpretable in terms of line profile. c double precision x,y double precision k c Constants double precision RRTPI ! 1/SQRT(pi) parameter(RRTPI=0.56418958D+0) double precision Y0,Y0PY0,Y0Q ! for CPF12 algorithm parameter ( Y0 = 1.5D+0, Y0PY0 = Y0+Y0, Y0Q = Y0*Y0 ) double precision C(0:5),S(0:5),T(0:5) save C,S,T c save preserves values of C, S and T (static) arrays between procedure calls DATA C / 1.0117281, -0.75197147, 0.012557727, & 0.010022008, -0.00024206814, 0.00000050084806 / DATA S / 1.393237, 0.23115241, -0.15535147, & 0.0062183662, 0.000091908299, -0.00000062752596 / DATA T / 0.31424038, 0.94778839, 1.5976826, & 2.2795071, 3.0206370, 3.8897249 / c Local variables c integer I integer j ! Loop variables integer RG1, RG2, RG3 ! y polynomial flags double precision ABX, XQ, YQ, YRRTPI ! |x|, x^2, y^2, y/SQRT(pi) double precision XLIM0, XLIM1, XLIM2, XLIM3, XLIM4 ! |x| on region boundaries double precision A0, D0, D2, E0, E2, E4, H0, H2, H4, H6 ! W4 temporary variables double precision P0, P2, P4, P6, P8, Z0, Z2, Z4, Z6, Z8 double precision XP(0:5), XM(0:5), YP(0:5), YM(0:5) ! CPF12 temporary values double precision MQ(0:5), PQ(0:5), MF(0:5), PF(0:5) double precision D, YF, YPY0, YPY0Q ***** Start of executable code ***************************************** YQ=Y**2 ! y^2 YRRTPI=Y*RRTPI ! y/SQRT(pi) if (Y.ge.70.55D+0) then ! All points XQ=x**2 k=YRRTPI/(XQ+YQ) RETURN endif RG1=1 ! Set flags RG2=1 RG3=1 XLIM0=SQRT ( 15100.0D+0+Y*(40.0D+0-Y*3.6D+0) ) ! y<70.55 if (Y.ge.8.425D+0) then XLIM1=0.0D+0 else XLIM1=SQRT ( 164.0D+0-Y*(4.3D+0+Y*1.8D+0) ) endif XLIM2=6.8D+0-Y XLIM3=2.4D+0*Y XLIM4=18.1D+0*Y+1.65D+0 if (Y.le.1.0D-6) then ! When y<10^-6 XLIM1=XLIM0 ! avoid W4 algorithm XLIM2=XLIM0 endif *..... ABX=dabs(x) ! |x| XQ =ABX*ABX ! x^2 if (ABX .ge. XLIM0) then ! Region 0 algorithm k=YRRTPI/(XQ+YQ) elseif (ABX .ge. XLIM1) then ! Humlicek W4 Region 1 if (RG1.ne.0) then ! First point in Region 1 RG1=0 A0=YQ+0.5D+0 ! Region 1 y-dependents D0=A0*A0 D2=YQ+YQ-1.0D+0 endif D=RRTPI/(D0+XQ*(D2+XQ)) k=D*Y*(A0+XQ) elseif (ABX.gt.XLIM2) then ! Humlicek W4 Region 2 if (RG2.ne.0) then ! First point in Region 2 RG2=0 H0=0.5625D+0+YQ*(4.5D+0+YQ*(10.5D+0+YQ*(6.0D+0+YQ))) ! Region 2 y-dependents H2=-4.5D+0+YQ*(9.0D+0+YQ*(6.0D+0+YQ* 4.0D+0)) H4=10.5D+0-YQ*(6.0D+0-YQ*6.0D+0) H6=-6.0D+0+YQ* 4.0D+0 E0=1.875D+0+YQ*(8.25D+0+YQ*(5.5D+0+YQ)) E2=5.25D+0+YQ*(1.0D+0+YQ* 3.0D+0) E4=0.75D+0*H6 endif D=RRTPI/(H0+XQ*(H2+XQ*(H4+XQ*(H6+XQ)))) k=D*Y*(E0+XQ*(E2+XQ*(E4+XQ))) elseif (ABX.lt.XLIM3) then ! Humlicek W4 Region 3 if (RG3.ne.0) then ! First point in Region 3 RG3=0 Z0=272.1014D+0+Y*(1280.829D+0+Y*(2802.870D+0+Y*(3764.966D+0 ! Region 3 y-dependents & +Y*(3447.629D+0+Y*(2256.981D+0+Y*(1074.409D+0 & +Y*(369.1989D+0+Y*(88.26741D+0 & +Y*(13.39880D+0+Y))))))))) Z2=211.678D+0+Y*(902.3066D+0+Y*(1758.336D+0+Y*(2037.310D+0 & +Y*(1549.675D+0+Y*(793.4273D+0+Y*(266.2987D+0 & +Y*(53.59518D+0+Y*5.0D+0))))))) Z4=78.86585D+0+Y*(308.1852D+0+Y*(497.3014D+0+Y*(479.2576D+0 & +Y*(269.2916D+0+Y*(80.39278D+0+Y*10.0D+0))))) Z6=22.03523D+0+Y*(55.02933D+0+Y*(92.75679D+0+Y*(53.59518D+0 & +Y*10.0D+0))) Z8=1.496460D+0+Y*(13.39880D+0+Y*5.0D+0) P0=153.5168D+0+Y*(549.3954D+0+Y*(919.4955D+0+Y*(946.8970D+0 & +Y*(662.8097D+0+Y*(328.2151D+0+Y*(115.3772D+0 & +Y*(27.93941D+0 & +Y*(4.264678D+0+Y*0.3183291D+0)))))))) P2=-34.16955D+0+Y*(-1.322256D+0+Y*(124.5975D+0 & +Y*(189.7730D+0 & +Y*(139.4665D+0+Y*(56.81652D+0+Y*(12.79458D+0 & +Y*1.2733163D+0)))))) P4=2.584042D+0+Y*(10.46332D+0+Y*(24.01655D+0+Y*(29.81482D+0 & +Y*(12.79568D+0+Y*1.9099744D+0)))) P6=-0.07272979D+0+Y*(0.9377051D+0+Y*(4.266322D+0 & +Y*1.273316D+0)) P8=0.0005480304D+0+Y*0.3183291D+0 endif D=1.7724538D+0/(Z0+XQ*(Z2+XQ*(Z4+XQ*(Z6+XQ*(Z8+XQ))))) k=D*(P0+XQ*(P2+XQ*(P4+XQ*(P6+XQ*P8)))) else ! Humlicek CPF12 algorithm YPY0=Y+Y0 YPY0Q=YPY0*YPY0 k=0.0D+0 do j=0,5 D=x-T(j) MQ(j)=D*D MF(j)=1.0D+0/(MQ(j)+YPY0Q) XM(j)=MF(j)*D YM(j)=MF(j)*YPY0 D=x+T(j) PQ(j)=D*D PF(j)=1.0D+0/(PQ(j)+YPY0Q) XP(j)=PF(j)*D YP(j)=PF(j)*YPY0 enddo if (ABX.le.XLIM4) then ! Humlicek CPF12 Region I do j=0,5 k=k+C(j)*(YM(j)+YP(j))-S(j)*(XM(j)-XP(j)) enddo else ! Humlicek CPF12 Region II YF=Y+Y0PY0 do j=0,5 k=k & +(C(j)*(MQ(j)*MF(j)-Y0*YM(j))+S(j)*YF*XM(j)) & /(MQ(j)+Y0Q) & +(C(j)*(PQ(j)*PF(j)-Y0*YP(j))-S(j)*YF*XP(j)) & /(PQ(j)+Y0Q) enddo k=Y*k+dexp(-XQ) endif endif *..... return end