/**************************************************************************************
*   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 https://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 https://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  https://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 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]; }
     }
   }


}