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
* complete frequency redistribution (where an adaptive 16-point Gaussian quadrature is
* 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:
* (see http://www.plasmaphysics.org.uk/radiative_transfer.htm ).
*
*
* 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[14]= {-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[0]+a[1]*x+a[2]*x2+a[3]*x3+a[4]*x4+a[5]*x5 -log(x) ; }
else { e1=expx*(x4+a[6]*x3+a[7]*x2+a[8]*x+a[9])/
(x5+a[10]*x4+a[11]*x3+a[12]*x2+a[13]*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[16]=
{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);
}
|