### Radiative Transfer Code for a Plane Parallel Multiple Scattering/ Absorbing Atmosphere

```/**************************************************************************************
*   These C/C++ routines calculate the source function produced by multiple scattering  and absorption
*   in a plane-parallel atmosphere, and the intensities on a given line in the medium (in both directions).
*   The former program is based on an iterative solution of the radiative transfer equation, the latter
*   is a straightforward integration based on the known optical depth and source function vectors
*   (see the page http://www.plasmaphysics.org.uk/radiative_transfer.htm for more).
*
*  Input parameters for the routine plansource are:
*    mode:  'c' for continuum- (monochromatic), 'l' for line-scattering
*    N:  number of layers (with the source function assumed constant in each)
*    tau:  vertical optical depth of scatterers (vector of length N+1)
*    t:  vertical optical depth of absorbers (vector of length N+1)
*    S0:  initial source function (vector of length N)
*    xd:  base interval for Gaussian Quadrature (recommended =0.5) (for the kernel function routine)
*    epsi:  required relative accuracy (see also the notes for the kernel function routine)
*
*  Output parameters for plansource are:
*    ST:  total source function (vector of length N)
*    ni:  number of iterations used
*
*
*  Input parameters for the routine intens are:
*    mode:  'c' for continuum- (monochromatic), 'l' for line-scattering
*    N:  number of layers/intervals (with the source function assumed constant in each)
*    tau:  line-of-sight optical depth of scatterers (vector of length N+1)
*    t:  line-of-sight optical depth of absorbers (vector of length N+1)
*    S:  source function (vector of length N)
*    xd:  base interval for Gaussian Quadrature (recommended =0.5) (for the kernel function routine)
*    epsi:  required relative accuracy (see also the notes for the kernel function routine)
*
*  Output parameters for intens are:
*    Iu:  up-looking intensities (vector of length N)    (referred to bottom of each layer/interval)
*    Id:  down-looking intensities (vector of length N)  (referred to top of each layer/interval)
*
*   (due to its design, intens is not restricted to a plane-parallel atmosphere, but can
*   be used for arbitray geometries as long as the optical depth and source function
*   along the line of sight is specified).
*
*   All variables are of double precision floating point type, apart from N and ni (integer) and 'mode'.
*
*   Required here are additonally the routines for the calculation of the kernel functions, which for
*   the case of continuum- (monochromatic) or spectral line transfer with a Gaussian profile
*   can be found under http://www.plasmaphysics.org.uk/programs/kernels.htm .
*   The code should work with different kernel functions (i.e. different line profiles) as well
*   (provided these are properly normalized), but this has not been tested.
*
*
*   In case of any questions or problems, please let me know
*   (for contact details see  http://www.plasmaphysics.org.uk/feedback.htm ).
*
*    Thomas Smid, June 2008
*                 July 2008 (added 'mode' option for line/continuum scattering,
*                            changed Neumann series iteration to direct iteration
*                            because of better convergence if absorption is significant)
*                 August 2008 (split the program into two different routines
*                             for the source functions and intensities
*                             for better flexibility and efficiency)
*****************************************************************************************/

void plansource(char mode, int N, double tau[], double t[], double S0[], double mu, double xd, double epsi,
double S[], int& ni)  {

double S1[N],KS[N][N],al[N],q[N];
double convgc,ddtau,ddt,jktau,jkt,jktau1,jkt1,dev,devmax;
int j,k;

/* calculate absorbing/scattering optical depth ratio, and initialize source function */

for (j=0; j<=N-1; j=j+1)  {
al[j]=(t[j+1]-t[j])/(tau[j+1]-tau[j]) ;
S[j]=S0[j];
}

/* assign the kernel function coefficients  */

for (j=0; j<=N-1; j=j+1)  {  // KS referred to midpoint of layer
for (k=0; k<=N-1; k=k+1)  {
if (k>=j) {ddtau=-(tau[k+1]-tau[k])/2; ddt=-(t[k+1]-t[k])/2; }
else {ddtau=(tau[k+1]-tau[k])/2; ddt=(t[k+1]-t[k])/2; }
jktau=fabs(tau[j]-tau[k])+ddtau; jkt=fabs(t[j]-t[k])+ddt;
if (k==j) {jktau=0; jkt=0; }
jktau1=fabs(tau[j]-tau[k+1])+ddtau; jkt1=fabs(t[j]-t[k+1])+ddt;
KS[j][k]= K2(mode,jktau,jkt,al[k],xd,epsi)- K2(mode,jktau1,jkt1,al[k],xd,epsi) ;
if (KS[j][k]<0) {KS[j][k]=-KS[j][k];}
}
KS[j][j]=2*KS[j][j];
}

/* iteration loop  */

convgc=1; ni=0;

while (convgc >= epsi) {

for (j=0; j<=N-1; j=j+1)  {
S1[j]=S[j];
}

for (j=0; j<=N-1; j=j+1)  {
S[j]=0;
for (k=0; k<=N-1; k=k+1)  {
if (k!=j) { S[j]=S[j]+ S1[k]*KS[j][k]; }
}
S[j]=(S[j]+S0[j])/(1-KS[j][j]);
}

devmax=0;
for (j=0; j<=N-1; j=j+1)  {
if (S1[j]!=0) { dev=fabs(S[j]/S1[j]-1.0); if (dev>devmax) devmax=dev; }
}
ni=ni+1;
convgc=devmax*ni;  // *ni in order to assure actual convergence

}     // end iteration loop

}

void intens(char mode, int N, double tau[], double t[], double S[], double xd, double epsi,
double Iu[], double Id[])  {

double KI[N][N],al[N];
double ddtau,ddt,jktau,jkt,jktau1,jkt1;
int j,k;

/* calculate absorbing/scattering optical depth ratio */

for (j=0; j<=N-1; j=j+1)  {
al[j]=(t[j+1]-t[j])/(tau[j+1]-tau[j]) ;
}

/* calculate transmission function matrix */

for (j=0; j<=N-1; j=j+1)  {  // KI referred to top/bottom of layer
for (k=0; k<=N-1; k=k+1)  {
if (k>=j) {jktau=fabs(tau[k]-tau[j]); jkt=fabs(t[k]-t[j]);
jktau1=fabs(tau[k+1]-tau[j]); jkt1=fabs(t[k+1]-t[j]);
}
else {jktau=fabs(tau[j+1]-tau[k]); jkt=fabs(t[j+1]-t[k]);
jktau1=fabs(tau[j+1]-tau[k+1]); jkt1=fabs(t[j+1]-t[k+1]);
}

KI[j][k]= Z(mode,jktau,jkt,al[k],xd,epsi)- Z(mode,jktau1,jkt1,al[k],xd,epsi) ;
if (KI[j][k]<0) {KI[j][k]=-KI[j][k];}
}
}

/* calculate up- and down looking intensities,
Iu is referred here to the bottom of each layer/interval, Id to the top,
(with 'top to bottom' defined as being in the direction of increasing j  ) */

for (j=0; j<=N-1; j=j+1)  {
Iu[j]=0; Id[j]=0;
for (k=0; k<=N-1; k=k+1)  {
if (k>=j) {Id[j]=Id[j]+S[k]*KI[j][k]; }
if (k<=j) {Iu[j]=Iu[j]+S[k]*KI[j][k]; }
}
}

}

```