// Implementation file for the object 'mie'. These routines compute
// tables of scattering and extinction cross-sections, albedo
// functions, and entries of the scattering matrix for an average dust
// grain. This average dust grain is given by integration over the
// defined dust composition and grain size distribution. The code for
// these calculations is taken from Bohren & Huffman (1983). It is
// based on solving the Maxwell equations for an electromagnetic wave
// scattered off a spherical grain. We implemented three tables of
// optical data for graphite and silicate grains. The non-isotropic
// optical properties of graphite crystals are considered by using
// two sorts of grains, 'GraphOrt' and 'GraphPar', for electromagnetic
// waves traveling perpendicular and parallel to the crystal's
// symmetry axis. The dielectric functions were taken from Draine &
// Lee (1984). The computed Mie data is stored in a look-up table of
// the class 'mie'.
//
// References:
// Bohren, C. F. & Huffman, D. R. 1983, New York: Wiley
// 1983 Draine, B. T. & Lee, H. M. 1984, ApJ, 285, 89
//
// STOKES, version 1.0, Nov 2004


#include <math.h>
#include <time.h>
#include <iostream.h>
#include <fstream.h>
#include <string.h>
#include <stdlib.h>
#include "mie-stokes-v1.0.h"
#include "specmath-stokes-v1.0.h"


// this routine reads an entry of the dielectric function table for
// the dust component 'subst' at wavelength lambda. It returns the
// complex refraction index n = sqrt(epsilon)
complex GetRef(int subst,double lambda)
{
  int i,j;
  double r,phi;
  complex ref,epsilon;

  i=0;
  j=0;
  switch (subst)
  {
  case 0:
    while(GraphParRe[i].lam<lambda) i++;
    epsilon.re = GraphParRe[i].eps;
    while(GraphParIm[j].lam<lambda) j++;
    epsilon.im = GraphParIm[j].eps;
    break;
  case 1:
    while(GraphOrtRe[i].lam<lambda) i++;
    epsilon.re = GraphOrtRe[i].eps;
    while(GraphOrtIm[j].lam<lambda) j++;
    epsilon.im = GraphOrtIm[j].eps;
    break;
  default:
    while(SilicateRe[i].lam<lambda) i++;
    epsilon.re = SilicateRe[i].eps;
    while(SilicateIm[j].lam<lambda) j++;
    epsilon.im = SilicateIm[j].eps;
    break;
  }
  epsilon = CSqrt(epsilon);
  ref.re = CReal(epsilon);
  ref.im = CImag(epsilon);
  return ref;
}


// Calculation of the relevant Mie data for a given size parameter
// x=2*pi*a/lambda, with a being the grain size radius, for a given
// relative refraction index RefRel, and for a given number NAng of
// scattering angles. The results are the entries S1 and S2 of the
// scattering matrix and the extinction and scattering efficiencies
// Qext=Cext/(pi*a^2) and Qsca=Csca/(pi*a^2). All results basically
// depend on two complex functions an and bn that are computed by an
// iterative method. For further details see Bohren and Huffman (1983).
//
// Reference: 
// Bohren, C. F. & Huffman, D. R. 1983, New York: Wiley
void MieCalc(double x,complex RefRel,int NAng,complex S1[],
	     complex S2[],double& Qext,double& Qsca)
{
  int NStop,j,jj,nmx,nn,rn,n,P,T;
  double AMu[NANG],theta[NANG],Pi[NANG],
  	   Pi0[NANG],Pi1[NANG],tau[NANG],
  	   psi,psi0,psi1,APsi,APsi0,APsi1,
  	   chi,chi0,chi1,dn,dx,XStop,YMod,fn,DAng;
  complex D[DMax],y,xi,xi0,xi1,an,bn;

  dx = x;
  y = x*RefRel;
  
  // Series terminated after NStop terms
  XStop = x + 4*pow(x,1.0/3) + 2;
  NStop = int(ceil(XStop));
  YMod = CAbs(y);
  if (XStop >= YMod)
  {
    nmx = int(ceil(XStop) + 15);
  }
  else
  {
    nmx = int(ceil(YMod) + 15);
  }
  if (nmx > DMax-1)
  {
    nmx = DMax-1;
  }

  // Initialization of the array for the scattering angle
  DAng = pi/(2*(NAng-1));
  for(j=0; j<NAng; j++)
  {
    theta[j] = j*DAng;
    AMu[j] = cos(theta[j]);
  }

  // Logarithmic derivative D[j] calculated by downward recurrence
  // beginning with initial value 0+0i at j=nmx
  D[nmx] = C(0,0);
  nn = nmx;
  for(n=0; n<nn; n++)
  {
    rn = nmx-n;
    D[nmx-(n+1)] = rn/y-1/(D[nmx-n]+rn/y);
  }
  for(n=0; n<NAng; n++)
  {
    Pi0[n] = 0;
    Pi1[n] = 1;
  }
  nn = 2*NAng-1;
  for(n=0; n<nn; n++)
  {
    S1[n] = C(0,0);
    S2[n] = C(0,0);
  }
  
  // Riccati-Bessel functions with real argument x 
  // calculated by upward recurrence
  psi0 = cos(dx);
  psi1 = sin(dx);
  chi0 = -sin(dx);
  chi1 = cos(dx);
  APsi0 = psi0;
  APsi1 = psi1;
  xi0 = C(APsi0,-chi0);
  xi1 = C(APsi1,-chi1);
  Qsca = 0;
  n = 1;
  do
  {
    dn = n;
    rn = n;
    fn = double((2*rn+1))/(rn*(rn+1));
    psi = (2*dn-1)*psi1/dx-psi0;
    APsi = psi;
    chi = (2*rn-1)*chi1/x-chi0;
    xi = C(APsi,-chi);
    an = (D[n]/RefRel+rn/x)*APsi-APsi1;
    an = an/((D[n]/RefRel+rn/x)*xi-xi1);
    bn = (D[n]*RefRel+rn/x)*APsi-APsi1;
    bn = bn/((D[n]*RefRel+rn/x)*xi-xi1);

    // (pre-)developement of Qsca 
    Qsca=Qsca+(2*rn+1)*(pow(CAbs(an),2)+pow(CAbs(bn),2));

    // development of S1 and S2
    for(j=0; j<NAng; j++)
    {
      jj = 2*(NAng-1)-j;
      Pi[j]=Pi1[j];
      tau[j]=rn*AMu[j]*Pi[j]-(rn+1)*Pi0[j];
      P = int(pow(-1,n-1));
      S1[j] = S1[j] + fn*(an*Pi[j]+bn*tau[j]);
      T = int(pow(-1,n));
      S2[j] = S2[j] + fn*(an*tau[j]+bn*Pi[j]);
      if (j==jj) break;
      S1[jj] = S1[jj] + fn*(an*Pi[j]*P+bn*tau[j]*T);
      S2[jj] = S2[jj] + fn*(an*tau[j]*T+bn*Pi[j]*P);
    }
    psi0 = psi1;
    psi1 = psi;
    APsi1 = psi1;
    chi0 = chi1;
    chi1 = chi;
    xi1 = C(APsi1,-chi1);
    n++;
    rn = n;
    for(j=0; j<NAng; j++)
    {
      Pi1[j] = double(2*rn-1)/(rn-1)*AMu[j]*Pi[j];
      Pi1[j] = Pi1[j] - rn*Pi0[j]/(rn-1);
      Pi0[j] = Pi[j];
    }
  }
  while(n-1<NStop);

  // Final calculation of Qsca and Qext
  Qsca = 2/(x*x)*Qsca;
  Qext = 4/(x*x)*CReal(S1[0]);
}


// This routine computes or, if already existing, reads in the Mie data for an
// average dust grain. The average grain is obtained by a numerical integration
// over the Mie data for all dust components and grain size radii. The
// Mie data table is stored in a dynamical array to have enough memory
// available. Depending on the model parameters the table may become huge.
void Mie::InitMie(double& RadMin,double& RadMax,int& RadNum,
                  int& LamNum,double Comp[SubNum],double& RadAlpha,
                  char dustfile[80])


{
  int l,r,s,t;
  double a,da,lam,rad,x,Qext,Qsca,Csca,N,Nt;
  double S11,S12,S33,S34;
  double d,dt,Int,ran,Sm1,Sm2;
  complex RefRel,S1[2*NANG-1],S2[2*NANG-1];
  char str[200],filename[80];
  time_t rawtime;
  fstream in_stream,out_stream;

  // try to open an existing dust model
  in_stream.open(dustfile,ios::in);

  // if no dust model exists yet, it is computed
  if (in_stream.fail()!=0)
  {
    dustfile[strlen(dustfile)-4]='\0';
    cout << "Dust model '" << dustfile << "' does not yet exist\n";
    cout << "Computation of dust properties...\n";
    strcat(dustfile,".dst");

    // prepare the grid of wavelengths for the grain model; the
    // wavelength range is given by the available data table of
    // optical constants
    GrainLamMin = 1000;
    GrainLamMax = 10000;
    GrainLamNum = LamNum;
    GrainLamDiff = (GrainLamMax-GrainLamMin)/GrainLamNum;
      
    // define a new Mie data table as a dynamical array
    MD = new Entry[GrainLamNum];
  
    // initialize the table, set everything to zero
    for(l=0; l<GrainLamNum; l++)
    {
      MD[l].Cext = 0;
      MD[l].Albedo = 0;
      for(t=0; t<2*NANG-1; t++)
      {
             MD[l].S11[t] = 0;
             MD[l].S12[t] = 0;
             MD[l].S33[t] = 0;
             MD[l].S34[t] = 0;
      }
      for(t=0; t<ThetaAcc; t++)
      {
             MD[l].Theta[t] = 0;
      }
    }
      
    // loop for the wavelength
    da = (RadMax-RadMin)/RadNum;
    l=0;
    for(lam=GrainLamMin+GrainLamDiff/2; lam<GrainLamMax; lam=lam+GrainLamDiff)
    {
      Csca = 0;
  
      // loop to integrate over the grain radii
      for(a=RadMin+0.5*da; a<RadMax; a=a+da)
      {
  
        // loop to integrate over the different dust components
        for(s=0; s<SubNum; s++)
        {
  
          // calculation of S1, S2, Cext, and Csca
          x = 20000*pi*a/lam;
          RefRel = GetRef(s,lam);
          MieCalc(x,RefRel,NANG,S1,S2,Qext,Qsca);
          MD[l].Cext = MD[l].Cext + Comp[s]*pi*pow(a,RadAlpha+2)*Qext*da;
          Csca = Csca + Comp[s]*pi*pow(a,RadAlpha+2)*Qsca*da;
  	  
       	  // loop over the different scattering angles theta
       	  // and computation of the Mueller matrices
       	  for(t=0; t<2*NANG-1; t++)
       	  {
       	    S11 = 0.5*(pow(CAbs(S2[t]),2) + pow(CAbs(S1[t]),2));
       	    S12 = 0.5*(pow(CAbs(S2[t]),2) - pow(CAbs(S1[t]),2));
       	    S33 = CReal(S1[t]*CConj(S2[t]));
       	    S34 = CImag(S2[t]*CConj(S1[t]));
       	    MD[l].S11[t] = MD[l].S11[t] + Comp[s]*pow(a,RadAlpha)*S11*da;
       	    MD[l].S12[t] = MD[l].S12[t] + Comp[s]*pow(a,RadAlpha)*S12*da;
       	    MD[l].S33[t] = MD[l].S33[t] + Comp[s]*pow(a,RadAlpha)*S33*da;
       	    MD[l].S34[t] = MD[l].S34[t] + Comp[s]*pow(a,RadAlpha)*S34*da;
       	  }
       	  // end loop over theta
  
        }
        // end loop over the components
      }
      // end loop over the grain radii
  
      // normalization of all data 
      N = 0;
      for(a=RadMin+0.5*da; a<RadMax; a=a+da)
      {
        N=N+pow(a,RadAlpha)*da;
      }
      Csca = Csca/N;
      MD[l].Cext = MD[l].Cext/N;
      for(t=0; t<2*NANG-1; t++)
      {
        MD[l].S11[t] = MD[l].S11[t] / N;
        MD[l].S12[t] = MD[l].S12[t] / N;
        MD[l].S33[t] = MD[l].S33[t] / N;
        MD[l].S34[t] = MD[l].S34[t] / N;
      }
  
      // calculation of the albedo function
      MD[l].Albedo = Csca/MD[l].Cext;
      
      // transform Cext from microns^2 to cm^2
      MD[l].Cext = MD[l].Cext*1.0e-8;

      l++;
    }
    // end loop over lambda
  
    // create file for the dust model
    out_stream.open(dustfile,ios::out);
    
    // save global parameters of the dust type
    strcpy(filename,dustfile);
    filename[strlen(dustfile)-4]='\0';
    out_stream << "Dust model " << filename << "\n";
    out_stream << "Stokes version 1.0" << "\n";
    time(&rawtime);
    out_stream << "Created " << ctime(&rawtime) << "\n\n";
    out_stream << "LambdaMin " << GrainLamMin << "\n";
    out_stream << "LambdaMax " << GrainLamMax << "\n";
    out_stream << "GrainRadMin " << RadMin << "\n";
    out_stream << "GrainRadMax " << RadMax << "\n";
    out_stream << "GrainSizeIndex " << RadAlpha << "\n";
    out_stream << "Para-Graphite " << Comp[0]*100 << "\n";
    out_stream << "Ortho-Graphit " << Comp[1]*100 << "\n";
    out_stream << "Silicate " << Comp[2]*100 << "\n";
    out_stream << "GrainRadNum " << RadNum << "\n";
    out_stream << "GrainLambdaNum " << GrainLamNum << "\n\n";
    
    // print the wavelength scale
    out_stream << "***** ";
    for(l=0; l<GrainLamNum; l++) out_stream << (GrainLamMin+l*GrainLamDiff)
                                            << ' ';
    out_stream << "\n";
    
    // save extinction cross sections
    out_stream << "sigma_ext ";
    for(l=0; l<GrainLamNum; l++) out_stream << MD[l].Cext << ' ';
    out_stream << "\n";
    
    // save scattering sections
    out_stream << "sigma_sca ";
    for(l=0; l<GrainLamNum; l++) out_stream << MD[l].Albedo*MD[l].Cext << ' ';
    out_stream << "\n";
    
    // save absorption cross sections
    out_stream << "sigma_abs ";
    for(l=0; l<GrainLamNum; l++) out_stream << (1-MD[l].Albedo)*MD[l].Cext
                                            << ' ';
    out_stream << "\n";
    
    // save albedos
    out_stream << "albedo ";
    for(l=0; l<GrainLamNum; l++) out_stream << MD[l].Albedo << ' ';
    out_stream << "\n\n";
    
    // print a second wavelength scale
    out_stream << "***** ***** ";
    for(l=0; l<GrainLamNum; l++) out_stream << (GrainLamMin+l*GrainLamDiff)
                                            << ' ';
    
    // save the parameters S11
    for(t=0; t<2*NANG-1; t++)
    {
      out_stream << "\n";
      out_stream << "S11 " << t << "  ";
      for(l=0; l<GrainLamNum; l++) out_stream << MD[l].S11[t] << ' ';
    }
    
    // save the parameters S12
    for(t=0; t<2*NANG-1; t++)
    {
      out_stream << "\n";
      out_stream << "S12 " << t << "  ";
      for(l=0; l<GrainLamNum; l++) out_stream << MD[l].S12[t] << ' ';
    }
    
    // save the parameters S33
    for(t=0; t<2*NANG-1; t++)
    {
      out_stream << "\n";
      out_stream << "S33 " << t << "  ";
      for(l=0; l<GrainLamNum; l++) out_stream << MD[l].S33[t] << ' ';
    }
    
    // save the parameters S34
    for(t=0; t<2*NANG-1; t++)
    {
      out_stream << "\n";
      out_stream << "S34 " << t << "  ";
      for(l=0; l<GrainLamNum; l++) out_stream << MD[l].S34[t] << ' ';
    }
    out_stream << "\n";
    
    out_stream.close();
    cout << "...finished and saved\n";
  }

  // if the dust model already exists it is read in from the file
  else
  {
    cout << "Dust model '" << dustfile << "' found...\n";
    
    // jump over title elements
    in_stream >> str;
    in_stream >> str;
    in_stream >> str;
    in_stream >> str;
    in_stream >> str;
    in_stream >> str;
    in_stream >> str;
    in_stream >> str;
    in_stream >> str;
    in_stream >> str;
    in_stream >> str;
    in_stream >> str;

    // read in global parameters of the dust type
    in_stream >> str;
    in_stream >> str;
    GrainLamMin = atoi(str);
    in_stream >> str;
    in_stream >> str;
    GrainLamMax = atoi(str);
    in_stream >> str;
    in_stream >> str;
    RadMin = atof(str);
    in_stream >> str;
    in_stream >> str;
    RadMax = atof(str);
    in_stream >> str;
    in_stream >> str;
    RadAlpha = atof(str);
    in_stream >> str;
    in_stream >> str;
    Comp[0] = atof(str)/100;
    in_stream >> str;
    in_stream >> str;
    Comp[1] = atof(str)/100;
    in_stream >> str;
    in_stream >> str;
    Comp[2] = atof(str)/100;
    in_stream >> str;
    in_stream >> str;
    RadNum = atoi(str);
    in_stream >> str;
    in_stream >> str;
    GrainLamNum = atoi(str);
    GrainLamDiff = (GrainLamMax-GrainLamMin)/GrainLamNum;
    LamNum = GrainLamNum;
 
    // define a new Mie data table as a dynamical array
    MD = new Entry[GrainLamNum];
  
    // jump over wavelengths
    in_stream >> str;
    for (l=0; l<GrainLamNum; l++) in_stream >> str;
    
    // read in extinction cross-sections
    in_stream >> str;
    for (l=0; l<GrainLamNum; l++)
    { 
      in_stream >> str;
      MD[l].Cext = atof(str);
    }
    
    // jump over scattering and absorption cross-sections
    in_stream >> str;
    for (l=0; l<GrainLamNum; l++) in_stream >> str;
    in_stream >> str;
    for (l=0; l<GrainLamNum; l++) in_stream >> str;
    
    // read in extinction cross-sections
    in_stream >> str;
    for (l=0; l<GrainLamNum; l++)
    { 
      in_stream >> str;
      MD[l].Albedo = atof(str);
    }
    
    // jump over wavelengths
    in_stream >> str;
    in_stream >> str;
    for (l=0; l<GrainLamNum; l++) in_stream >> str;
    
    // read in parameters S11
    for (t=0; t<2*NANG-1; t++)
    {
      in_stream >> str;
      in_stream >> str;
      for (l=0; l<GrainLamNum; l++)
      {
        in_stream >> str;
        MD[l].S11[t] = atof(str);
      }
    }
    
    // read in parameters S12
    for (t=0; t<2*NANG-1; t++)
    {
      in_stream >> str;
      in_stream >> str;
      for (l=0; l<GrainLamNum; l++)
      {
        in_stream >> str;
        MD[l].S12[t] = atof(str);
      }
    }
    
    // read in parameters S33
    for (t=0; t<2*NANG-1; t++)
    {
      in_stream >> str;
      in_stream >> str;
      for (l=0; l<GrainLamNum; l++)
      {
        in_stream >> str;
        MD[l].S33[t] = atof(str);
      }
    }
    
    // read in parameters S34
    for (t=0; t<2*NANG-1; t++)
    {
      in_stream >> str;
      in_stream >> str;
      for (l=0; l<GrainLamNum; l++)
      {
        in_stream >> str;
        MD[l].S34[t] = atof(str);
      }
    }
    in_stream.close();

    cout << "...and read in\n";
  }
  
  // Creation of a table to find a scattering angle theta for a given
  // random number. First the normalization constant Nt is computed,
  // then for each random number ran the corresponding theta is
  // derived.
  for (l=0; l<GrainLamNum; l++)
  {
    Nt = 0;
    dt=pi/(2*(NANG-1));
    for(t=0; t<2*NANG-1; t++)
    {
      Nt = Nt+2*MD[l].S11[t]*sin(t*dt)*dt;
    }
    r=0;
    for(ran=1.0/ThetaAcc; ran<1; ran=ran+1.0/ThetaAcc)
    {
      d=1;
      Int=0;
      for(t=0; t<2*NANG-1; t++)
      {
        Int = Int+1/Nt*2*MD[l].S11[t]*sin(t*dt)*dt;
        if(fabs(ran-Int)<d)
        {
          d=fabs(ran-Int);
        }
        else
        {
          break;
        }
      }
      MD[l].Theta[r] = (t-1)*dt;
      r++;
    }
  }
}


// This function returns the extinction cross-section for a given
// wavelength lambda. Checks are included to make sure that the
// wavelength index l does not leave the table limits. Since l is
// computed from a float number this might happen due to improper
// conversion.
double Mie::GetCext(double lambda)
{
  int l;

  l = int(floor((lambda-GrainLamMin)/GrainLamDiff));
  if (l<0)
  {
    l = 0;
  }
  if (l>=GrainLamNum)
  {
    l = GrainLamNum-1;
  }
  return MD[l].Cext;
}


// This function returns the albedo for a given wavelength
// lambda. Checks are included to make sure that the wavelength index
// l does not leave the table limits. Since l is computed from a float
// number this might happen due to improper conversion.
double Mie::GetAlbedo(double lambda)
{
  int l;

  l = int(floor((lambda-GrainLamMin)/GrainLamDiff));
  if (l<0)
  {
    l = 0;
  }
  if (l>=GrainLamNum)
  {
    l = GrainLamNum-1;
  }
  return MD[l].Albedo;
}


// This function returns the scattering angle for a given wavelength
// lambda and a random number r. Checks are included to make sure that
// the wavelength index l and r do not leave the table limits. Since
// they are computed from float numbers this might happen due to
// improper conversion.
double Mie::GetTheta(double lambda,double ran)
{
  int l,r;

  l = int(floor((lambda-GrainLamMin)/GrainLamDiff));
  if (l<0)
  {
    l = 0;
  }
  if (l>=GrainLamNum)
  {
    l = GrainLamNum-1;
  }
  r = int(floor(ran*ThetaAcc));
  if (r<0)
  {
    r = 0;
  }
  if (r>=ThetaAcc)
  {
    r = ThetaAcc-1;
  }
  return MD[l].Theta[r];
}


// This function returns the entry S11 of the Mueller matrix for a
// given wavelength lambda and a scattering angle theta. Checks are
// included to make sure that the wavelength index l and the angle
// index t do not leave the table limits. Since they are computed from
// float numbers this might happen due to improper conversion.
double Mie::GetS11(double lambda,double theta)
{
  int l,t;
  double DAng;

  l = int(floor((lambda-GrainLamMin)/GrainLamDiff));
  if (l<0)
  {
    l = 0;
  }
  if (l>=GrainLamNum)
  {
    l = GrainLamNum-1;
  }
  DAng = pi/(2*(NANG-1));
  t = int(floor((theta/DAng+0.5)));
  if (t<0)
  {
    t = 0;
  }
  if (t>=2*NANG)
  {
          t = 2*NANG-1;
  }
  return MD[l].S11[t];
}


// This function returns the entry S12 of the Mueller matrix for a
// given wavelength lambda and a scattering angle theta. Checks are
// included to make sure that the wavelength index l and the angle
// index t do not leave the table limits. Since they are computed from
// float numbers this might happen due to improper conversion.
double Mie::GetS12(double lambda,double theta)
{
  int l,t;
  double DAng;

  l = int(floor((lambda-GrainLamMin)/GrainLamDiff));
  if (l<0)
  {
    l = 0;
  }
  if (l>=GrainLamNum)
  {
    l = GrainLamNum-1;
  }
  DAng = pi/(2*(NANG-1));
  t = int(floor((theta/DAng+0.5)));
  if (t<0)
  {
    t = 0;
  }
  if (t>=2*NANG)
  {
          t = 2*NANG-1;
  }
  return MD[l].S12[t];
}


// This function returns the entry S33 of the Mueller matrix for a
// given wavelength lambda and a scattering angle theta. Checks are
// included to make sure that the wavelength index l and the angle
// index t do not leave the table limits. Since they are computed from
// float numbers this might happen due to improper conversion.
double Mie::GetS33(double lambda,double theta)
{
  int l,t;
  double DAng;

  l = int(floor((lambda-GrainLamMin)/GrainLamDiff));
  if (l<0)
  {
    l = 0;
  }
  if (l>=GrainLamNum)
  {
    l = GrainLamNum-1;
  }
  DAng = pi/(2*(NANG-1));
  t = int(floor((theta/DAng+0.5)));
  if (t<0)
  {
    t = 0;
  }
  if (t>=2*NANG)
  {
    t = 2*NANG-1;
  }
  return MD[l].S33[t];
}


// This function returns the entry S34 of the Mueller matrix for a
// given wavelength lambda and a scattering angle theta. Checks are
// included to make sure that the wavelength index l and the angle
// index t do not leave the table limits. Since they are computed from
// float numbers this might happen due to improper conversion.
double Mie::GetS34(double lambda,double theta)
{
  int l,t;
  double DAng;

  l = int(floor((lambda-GrainLamMin)/GrainLamDiff));
  if (l<0)
  {
    l = 0;
  }
  if (l>=GrainLamNum)
  {
    l = GrainLamNum-1;
  }
  DAng = pi/(2*(NANG-1));
  t = int(floor((theta/DAng+0.5)));
  if (t<0)
  {
    t = 0;
  }
  if (t>=2*NANG)
  {
          t = 2*NANG-1;
  }
  return MD[l].S34[t];
}


// This function sets S11, S12, S33, and S34 to the entries of the
// Mueller matrix at the wavelength lambda and for scattering angle theta
void Mie::GetMueller(double lambda,double theta,double& S11,double& S12,
		     double& S33,double& S34)
{
  S11 = GetS11(lambda,theta);
  S12 = GetS12(lambda,theta);
  S12 = S12/S11;
  S33 = GetS33(lambda,theta);
  S33 = S33/S11;
  S34 = GetS34(lambda,theta);
  S34 = S34/S11;
  S11 = S11/S11;
}