### Kernel Functions for Radiative Transfer in Spectral Lines

```/*********************************************************************************
*   These C/C++ functions compute the Kernel functions required for the radiative transfer of
*   spectral lines, applicable either to line- scattering with a Doppler (Gaussian) profile and
*   used for the frequency intergration) or to continuum- (monochromatic) scattering.
*
*   The fundamental kernel in the integral equation for the source function as a function of
*   the vertical optical depths of scatterers τ and absorbers t in a plane-parallel
*   atmosphere is for the case of spectral line transfer:
*
*
*   K1(τ,t)= 1/2/√π.-∞∫∞dx e-2x2.E1(τ.e-x2+t)
*
*   with E1 the first exponential integral.
*   This function is itself rarely needed however, but only its successive integrals over τ,
*   i.e.
*
*   K2(τ,t,α)= 1/2/√π.-∞∫∞dx e-2x2/(e-x2+α).E2(τ.e-x2+t)
*
*   and
*
*   K3(τ,t,α)= 1/2/√π.-∞∫∞dx e-2x2/(e-x2+α)2.[E3(t-ατ) -E3(τ.e-x2+t)]
*
*
*   where E2 and E3 are the second and third exponential integrals respectively.
*   The additional parameter α arises from a linear parameterization of the absorber optical
*   depth t in terms of the scattering optical depth τ between 2 discretization points k and k+1,
*   i.e.  t=tk+αk.(τ-τk) .
*   (note that K3 is only required for a linear interpolation of the source function
*   between two discretization points, not for a 'step-function' approximation).
*
*
*   The calculation of the resultant intensity from the source function on the other hand
*   involves the transmission function
*
*   T(τ,t)= 1/√π.-∞∫∞dx e-x2.e-(τ.e-x2+t)
*
*   and its integral, the 'curve-of-growth' function
*
*   Z(τ,t,α)= 1/√π.-∞∫∞dx e-x2/(e-x2+α).[e-(t-ατ) -e-(τ.e-x2+t)]   ,
*
*   where here however τ and t are line-of-sight optical depths, not vertical optical
*   depths.
*   (similarly to K3 above, one can here also define a further function for
*   a higher order approximation, but as the program for this requires extended
*   numerical accuracy, it has been omitted here).
*
*
*   The corresponding functions for continuum- (monochromatic) scattering are simply
*   obtained from the above by leaving out the frequency integration and setting e-x2=1
*   (resulting thus in straight exponential functions and exponential integrals).
*
*
*   The above functions are called as dummy functions K1,K2,K3,T,Z with the arguments
*   (all double accuracy variables apart from 'mode'):
*
*    mode:  'c' for continuum- (monochromatic), 'l' for line-scattering
*    tau:  optical depth of scatterers to point k
*    tauab:  optical depth of absorbers to point k
*    al:  the ratio Δt/Δτ in the interval k-1,k (not required for K1 and T).
*    xd:  base interval for Gaussian Quadrature (recommended =0.5)
*    epsi:  required relative accuracy (note that this is limited by the
*           accuracy of the approximation for the exponential integral in
*           case of the functions K1,K2,K3).
*
*
*   In case of any questions or problems, please let me know
*   (for contact details see  http://www.plasmaphysics.org.uk/feedback.htm ).
*
*      Thomas Smid, April 2008
*                   July 2008 (added option for continuum- (monochromatic) scattering)
*********************************************************************************/

// this function computes the exponential integrals required for K1- K3;
// the approximation is taken from Abramowitz and Stegun's
// Handbook of Mathematical Functions (p.231);
// note that its accuracy is of the order of 10-7 for x≤1
// and of the order of 10-4 for x>1.

void expin(double x,double& e1,double& e2,double& e3,double& e4) {

double a= {-0.57721566,0.99999193,-0.24991055,0.05519968,
-0.00976004,0.00107857,8.5733287401,18.0590169730,
8.6347608925,0.2677737343,9.5733223454,25.6329561486,
21.0996530827,3.9584969228};
double x2,x3,x4,x5,expx;

if (x<=0) { e2=1; e3=0.5; e4=1/3.0 ; return; }

expx=exp(-x) ;
x2=x*x;
x3=x2*x;
x4=x3*x;
x5=x4*x;

if (x<=1)  {e1=a+a*x+a*x2+a*x3+a*x4+a*x5 -log(x) ; }
else { e1=expx*(x4+a*x3+a*x2+a*x+a)/
(x5+a*x4+a*x3+a*x2+a*x) ;
}

e2=expx-x*e1 ;
e3=(expx-x*e2)/2 ;
e4=(expx-x*e3)/3 ;

}

// this is the actual Gaussian Quadrature routine called from the 'dummy functions' below

double kernels(char func, char mode, double tau, double tauab, double al, double xd, double epsi)    {

double xl,xu,df,dfm,f,convgc,pi=3.141592653589793e+00,taual,eal;
double a,b,c,ly,x1,x2,expx1,expx2,taux1,taux2,mf1,mf2,df1,df2,e1,e2,e3,e4;
double qg=
{4.947004674958250e-01, 1.357622970587705e-02,
4.722875115366163e-01, 3.112676196932395e-02,
4.328156011939159e-01, 4.757925584124639e-02,
3.777022041775015e-01, 6.231448562776694e-02,
3.089381222013219e-01, 7.479799440828837e-02,
2.290083888286137e-01, 8.457825969750127e-02,
1.408017753896295e-01, 9.130170752246179e-02,
4.750625491881872e-02, 9.472530522753425e-02};

int i,jj;

if (mode=='c') {   // continuum (monochromatic) scattering

if (func=='1') { expin(tau+tauab,e1,e2,e3,e4); return e1/2; }
else if (func=='2')  { expin(tau+tauab,e1,e2,e3,e4); return e2/(1+al)/2; }
else if (func=='3')  { expin(tau+tauab,e1,e2,e3,e4); return (1-e3)/(1+al)/(1+al)/2; }
else if (func=='T')  { return exp(-(tau+tauab)); }
else if (func=='Z')  { return (1-exp(-(tau+tauab)))/(1+al); }

}

else if (mode=='l') { // line scattering (with Doppler profile)

if (tau==0 && tauab==0) {
if (func=='1') { return -100; }   // some nonsense value is returned as K1 diverges at 0
else if (func=='2'&& al==0) {return 0.5;}
else if (func=='T') {return 1.0;}
else if (func=='3' || func=='Z') {return 0.0;}
}

if (func=='Z' || func=='3') {
taual=tauab-al*tau;
if (taual<=0) { if (func=='Z') {eal=1;}
else {eal=0.5;}
}
else   { if (func=='Z') {eal=exp(-taual);}
else {expin(taual,e1,e2,eal,e4);}
}
}

f=0 ;
xu=0 ;
convgc=1 ; jj=1 ;

while(convgc>=epsi)  {

xl=xu ;
xu=xl+xd ;   // note: xd should be made variable for other (e.g. Voigt-) line profiles
a=0.5*(xu+xl) ;
b=xd ;
df=0 ;
for (i=0; i<=14; i=i+2)	{
c=qg[i]*b ;
x1=a+c;
x2=a-c;
expx1=exp(-x1*x1);
expx2=exp(-x2*x2);
taux1=tau*expx1+tauab;
taux2=tau*expx2+tauab;

if (func=='1') {
expin(taux1,e1,e2,e3,e4); df1=expx1*expx1*e1;
expin(taux2,e1,e2,e3,e4); df2=expx2*expx2*e1;
}

else if (func=='2') {
expin(taux1,e1,e2,e3,e4); df1=expx1*expx1/(expx1+al)*e2;
expin(taux2,e1,e2,e3,e4); df2=expx2*expx2/(expx2+al)*e2;
}

else if (func=='3') {
mf1=expx1/(expx1+al);  mf2=expx2/(expx2+al);
expin(taux1,e1,e2,e3,e4); df1=mf1*mf1*(eal-e3);
expin(taux2,e1,e2,e3,e4); df2=mf2*mf2*(eal-e3);
}

else if (func=='T') {df1=expx1*exp(-taux1); df2=expx2*exp(-taux2);}

else if (func=='Z') {
df1=expx1/(expx1+al)*(eal-exp(-taux1));
df2=expx2/(expx2+al)*(eal-exp(-taux2));
}

df=df+qg[i+1]*(df1+df2) ;
}
df=df*b;
f=f+df;
dfm=f/jj ;
if (dfm>0) {convgc=fabs(df/dfm) ;}
jj=jj+1;

}

if (func=='T' || func=='Z') { return 2*f/sqrt(pi) ; }
else  { return f/sqrt(pi) ; }

}

}

// these are just the dummy functions according to the above definitions

double K1(char mode, double tau, double tauab, double xd, double epsi)  {
return kernels('1', mode, tau, tauab, 0, xd, epsi);
}

double K2(char mode, double tau, double tauab, double al, double xd, double epsi)  {
return kernels('2', mode, tau, tauab, al, xd, epsi);
}

double K3(char mode, double tau, double tauab, double al, double xd, double epsi)  {
return kernels('3', mode, tau, tauab, al, xd, epsi);
}

double T(char mode, double tau, double tauab, double xd, double epsi)  {
return kernels('T', mode, tau, tauab, 0, xd, epsi);
}

double Z(char mode, double tau, double tauab, double al, double xd, double epsi)  {
return kernels('Z', mode, tau, tauab, al, xd, epsi);
}

```