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 ).
*   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=tkk.(τ-τ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 ). 
*      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,
   double x2,x3,x4,x5,expx;

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

   expx=exp(-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') {
         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 ;
          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') {
          df=df+qg[i+1]*(df1+df2) ;
        dfm=f/jj ;
        if (dfm>0) {convgc=fabs(df/dfm) ;}
       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);

Programs Home