//******************************************************************************
//** SCATMECH: Polarized Light Scattering C++ Class Library
//** 
//** File: bobvlieg1.cpp
//**
//** Thomas A. Germer
//** Optical Technology Division, National Institute of Standards and Technology
//** 100 Bureau Dr. Stop 8443; Gaithersburg, MD 20899-8443
//** Phone: (301) 975-2876; FAX: (301) 975-6991
//** Email: thomas.germer@nist.gov
//**
//** Version: 6.00 (February 2008)
//**
//******************************************************************************
//#include <float.h>
#include <fstream>
#include <vector>
#include "scatmech.h"
#include "bobvlieg.h"
#include "axipart.h"
#include "dielfunc.h"
#include "askuser.h"

using namespace std;

namespace SCATMECH {
    
    using namespace BobVlieg_Supp;

    // The following code implements the theory described in 
    //   P.A. Bobbert and J. Vlieger, Physica 137A (1986) 209-242.
    // Some details of the theory can be found also in 
    //   P.A. Bobbert, J. Vlieger, and  R. Grief, Physica 137A (1986) 243-257
    // These articles will be referred to as B&V and BV&G, respectively.
    //

    static complex<float> cfI(0.,1.);

    static 
    COMPLEX
    D(int n,int m,int f) 
    {
        int mm = (m==0 ? 1 : 2);

        double d = mm*(2*n+1)*Fact(n-abs(m))/(4*n*(n+1)*Fact(n+abs(m)));

        double e = sqrt((2*n+1)*Fact(n-abs(m))/Fact(n+abs(m)));

        COMPLEX ff = (f==0) ? 1. : -1.;

        return 1./(ff*(double)(e/d));
    }

    int 
    Axisymmetric_Particle_BRDF_Model::
    get_recalc()
    {
        return Local_BRDF_Model::get_recalc()||Shape->get_recalc();
    }

    //
    // Bmatrix() calculates the Mie scattering matrix found 
    // in Eq. (4.2) of B&V. 
    //
    void 
    Axisymmetric_Particle_BRDF_Model::
    Bmatrix(ScatterTMatrix& T)
    {
        //float xdummy;
        //long ldummy;
        complex<float> dummy;

        int m;
        for (m=-MMAX;m<=MMAX;++m) {
         int mm = LMAX+m;
         int beginl = (m==0) ? 1 : abs(m);
         for (int l=beginl;l<=LMAX;++l) {
          for (int f=0;f<=1;++f) {
           int i = index(l,m,f);
           int ll = 2*(LMAX-l)+f;
           for (int l_=beginl;l_<=LMAX;++l_) {
            for (int f_=0;f_<=1;++f_) {
              int i_ = index(l_,m,f_);
              int ll_ = 2*(LMAX-l_)+f_;
              T[mm][ll_][ll]=0;
            }
           }
          }
         }
        }

        Shape->Get_TMatrix(T, LMAX, MMAX, k, N2);

        for (m=-MMAX;m<=MMAX;++m) {
         int mm = LMAX+m;
         int beginl = (m==0) ? 1 : abs(m);
         for (int l=beginl;l<=LMAX;++l) {
          for (int f=0;f<=1;++f) {
           int ll = 2*(LMAX-l)+f;
           for (int l_=beginl;l_<=LMAX;++l_) {
            for (int f_=0;f_<=1;++f_) {
              int ll_ = 2*(LMAX-l_)+f_;
              T[mm][ll_][ll]=T[mm][ll_][ll]*D(l_,m,f_)/D(l,m,f);
            }
           }
          }
         }
        }


    }

    //
    // The complex arccosine is not part of the standard C++ library (!)...
    //
    static COMPLEX arccosine(const COMPLEX& a)
    {
        double x=real(a);
        double y=imag(a);

        return pi/2. - arg(sqrt(1. - sqr(x + cI*y)) + 
                cI*(x + cI*y)) + 
                cI*log(abs(sqrt(1. - sqr(x + cI*y)) + 
                cI*(x + cI*y)));
    }

    //
    // Amatrix() calculates the matrix A given by Eq. (8.17) of B&V
    //
    void
    Axisymmetric_Particle_BRDF_Model::
    Amatrix(ScatterTMatrix& A)
    {
        // Choose number of points in Gaussian integration, depending upon
        // LMAX, and round up to nearest 10.
        int npts=((LMAX/10)+1)*10;
        if (npts>100) npts=100;
        int ipt = npts/10-1;

        // Calculate the coefficients in front of the 
        // Legendre polynomials found in Eqs. (8.5a) and (8.5b) of B&V.
        vector<COMPLEX> U1(index(LMAX,LMAX)+1);
        vector<COMPLEX> U2(index(LMAX,LMAX)+1);
        vector<COMPLEX> U3(index(LMAX,LMAX)+1);
        vector<COMPLEX> U4(index(LMAX,LMAX)+1);
        vector<COMPLEX> V1(index(LMAX,LMAX)+1);
        vector<COMPLEX> V2(index(LMAX,LMAX)+1);
        for (int l=1;l<=LMAX;++l) {
            for (int m=0;m<=l;++m) { // We don't need negative m values
                int i=index(l,m);

                double temp1=(2.*l-1.)*(2.*l+1.);
                double temp2=(2.*l+3.)*(2.*l+1.);

                U1[i] = (l-1.)/2.*ipow(-(l-1))*sqrt((l+m-1.)*(l+m)/temp1);
                U2[i] = (l-1.)/2.*ipow(-(l-1))*sqrt((l-m-1.)*(l-m)/temp1);
                U3[i] =-(l+2.)/2.*ipow(-(l+1))*sqrt((l-m+1.)*(l-m+2.)/temp2);
                U4[i] =-(l+2.)/2.*ipow(-(l+1))*sqrt((l+m+1.)*(l+m+2.)/temp2);

                V1[i] = 0.5*ipow(-(l-1))*sqrt((l-m+1.)*(l+m));
                V2[i] =-0.5*ipow(-(l-1))*sqrt((l+m+1.)*(l-m));
            }
        }

        // Normal incidence reflection coefficients...
        COMPLEX rsnormal = rs(COMPLEX(1.,0.));
        COMPLEX rpnormal = rp(COMPLEX(1.,0.));

        vector<COMPLEX> reflect_s(npts);
        vector<COMPLEX> reflect_p(npts);
        vector<vector<COMPLEX> > Umatrix(npts,vector<COMPLEX>(index(LMAX,LMAX)+1));
        vector<vector<COMPLEX> > Vmatrix(npts,vector<COMPLEX>(index(LMAX,LMAX)+1));
        vector<vector<COMPLEX> > dminusmatrix(npts,vector<COMPLEX>(index(LMAX,LMAX)+1));
        vector<vector<COMPLEX> > dplusmatrix(npts,vector<COMPLEX>(index(LMAX,LMAX)+1));

        for (int j=0;j<npts;++j) {
            ostringstream msg;
            msg << "Creating lookup tables: j = " << j << '/' << npts 
                    << ret;
            message(msg.str());
        
            COMPLEX cosa = 1.0 - Gauss_Laguerre_Integration::zeros[ipt][j]/(2.*cI*qq); 
            COMPLEX alpha = arccosine(cosa);              
            COMPLEX cosa2 = cos((pi-alpha)/2.); 
            COMPLEX sina2 = sin((pi-alpha)/2.);

            if (Norm_Inc_Approx) { // Normal incidence approximation
                reflect_s[j] = rsnormal;
                reflect_p[j] = rpnormal;
            } else {               // Exact solution
                reflect_s[j] = rs(cosa);
                reflect_p[j] = rp(cosa);
            }

            // Calculate U, V, d+, and d- for each angle and (l,m) combination...
            for (int l=1;l<=LMAX;++l) {
                int _mmax = (MMAX>l) ? l : MMAX;
                for (int m=0;m<=_mmax;++m) {  // We don't need negative m values
                    int i = index(l,m);
                    Umatrix[j][i] = U1[i]*Ptilde(l-1,m-1,cosa)+
                                    U2[i]*Ptilde(l-1,m+1,cosa)+
                                    U3[i]*Ptilde(l+1,m-1,cosa)+
                                    U4[i]*Ptilde(l+1,m+1,cosa);
                    Vmatrix[j][i] = V1[i]*Ptilde(l,m-1,cosa)+
                                    V2[i]*Ptilde(l,m+1,cosa);

                    dminusmatrix[j][i] = dminus(l,m,cosa2,sina2);
                    dplusmatrix[j][i] = dplus(l,m,cosa2,sina2);
                }
            }
        }
    
        // Constant from change of integration variable (see Sec. 4 of BV&G)...
        COMPLEX expx = exp(2.*cI*qq)/(2.*cI*qq);

        // Construct A matrix...
        for (int m=0;m<=MMAX;++m) {
            ostringstream msg;
            msg << "Constructing A matrix: m = " << m << '/' << MMAX << ret;
            message(msg.str());
            int beginl = (m==0) ? 1 : m;              
            for (int l=beginl;l<=LMAX;++l) {
                for (int l_=beginl;l_<=LMAX;++l_) {
                    if (l>old_LMAX || l_>old_LMAX) {
                        COMPLEX Aee=0,Ahe=0,Aeh=0,Ahh=0;

                        int ii = index(l,m);
                        int ii_ = index(l_,m);

                        // Integration by Gaussian method [see Eq. (4.2) of BV&G]
                        for (int k=0;k<npts;++k) {
                            double weight = Gauss_Laguerre_Integration::weights[ipt][k];
                            Aee += weight*(reflect_p[k]*Vmatrix[k][ii]*dminusmatrix[k][ii_]
                                          +reflect_s[k]*Umatrix[k][ii]*dplusmatrix[k][ii_]);
                            Ahe += weight*(reflect_p[k]*Vmatrix[k][ii]*dplusmatrix[k][ii_]
                                          +reflect_s[k]*Umatrix[k][ii]*dminusmatrix[k][ii_]);
                            Aeh += weight*(reflect_p[k]*Umatrix[k][ii]*dminusmatrix[k][ii_]
                                          +reflect_s[k]*Vmatrix[k][ii]*dplusmatrix[k][ii_]);
                            Ahh += weight*(reflect_p[k]*Umatrix[k][ii]*dplusmatrix[k][ii_]
                                          +reflect_s[k]*Vmatrix[k][ii]*dminusmatrix[k][ii_]);
                        }
                    
                        // Constant value in front of Eq. (8.17) of BV...
                        COMPLEX aa = expx*mpow(m-1)*ipow(l_-1)*lvector[l_];

                        Aee *= aa;
                        Ahe *= -cI*aa;
                        Aeh *= cI*aa;
                        Ahh *= aa;

                        // Put the results into the A matrix...
                    
                        A[LMAX+m][2*(LMAX-l_)][2*(LMAX-l)] = Aee;
                        A[LMAX+m][2*(LMAX-l_)+1][2*(LMAX-l)] = Ahe;
                        A[LMAX+m][2*(LMAX-l_)+1][2*(LMAX-l)+1] = Ahh;
                        A[LMAX+m][2*(LMAX-l_)][2*(LMAX-l)+1] = Aeh;
                    
                        // Use the symmetry relations described in 
                        // Eqs. (4.5) of BV&G...

                        A[LMAX-m][2*(LMAX-l_)][2*(LMAX-l)] = Aee;
                        A[LMAX-m][2*(LMAX-l_)+1][2*(LMAX-l)] = -Ahe;
                        A[LMAX-m][2*(LMAX-l_)+1][2*(LMAX-l)+1] = Ahh;
                        A[LMAX-m][2*(LMAX-l_)][2*(LMAX-l)+1] = -Aeh;
                    }
                }
            }
        }
    }

    //
    // Matrix inversion
    //
    #define VV(a,b) V[(a)][(b)]

    static
    void 
    cvert(COMPLEX** V, int N)
    {   
        vector<int> W(N);
        
        COMPLEX Y,Z;

        double S,T;
        int   I,J,K,L,M,P;
        //int test;
        if ( N != 1 ) {
            L = 0;
            M = 1;
            Line10: 
            if ( L != N ) {
                K = L;
                L = M;
                M = M + 1;
                //C     ---------------------------------------
                //C     |*** FIND PIVOT AND START ROW SWAP ***|
                //C     ---------------------------------------
                P = L;
                S = abs(VV(L-1,L-1));
                if ( M <= N ) {
                    for (I=M;I<=N;++I) { 
                        T = abs(VV(I-1,L-1));
                        if ( T > S ) {
                            P = I;
                            S = T;
                        }
                    } 
                    W[L-1] = P;
                }
                if ( S == 0. ) throw SCATMECH_exception("Matrix has no inverse (1)");
                Y = VV(P-1,L-1);
                VV(P-1,L-1) = VV(L-1,L-1);
                //C     -----------------------------
                //C     |*** COMPUTE MULTIPLIERS ***|
                //C     -----------------------------
                VV(L-1,L-1) = COMPLEX(-1.,0.);
                Y = COMPLEX(1.,0.)/Y;
                for (I=1;I<=N;++I) { 
                    VV(I-1,L-1) = -Y*VV(I-1,L-1);
                }
                J = L;

                while (1) {
                    do {
                        J = J + 1;
                        if ( J > N ) J = 1;
                        if ( J == L ) goto Line10;
                        Z = VV(P-1,J-1);
                        VV(P-1,J-1) = VV(L-1,J-1);
                        VV(L-1,J-1) = Z;
                    } while ( abs(Z) == 0. ); // goto Line50;
                    //C     ------------------------------
                    //C     |*** ELIMINATE BY COLUMNS ***|
                    //C     ------------------------------
                    if ( K != 0 ) {
                        for (I=1;I<=K;++I) { // DO 60 I = 1,K
                            VV(I-1,J-1) = VV(I-1,J-1) + Z*VV(I-1,L-1);
                        }
                    }
                    VV(L-1,J-1) = Y*Z;
                    if ( M <= N ) {
                        for (I=M;I<=N;++I) {
                            VV(I-1,J-1) = VV(I-1,J-1) + Z*VV(I-1,L-1);
                        }
                    }
                }; // goto Line50;
                //C     -----------------------
                //C     |*** PIVOT COLUMNS ***|
                //C     -----------------------
            }

            do {
                L = W[K-1];
                for (I=1;I<=N;++I) { 
                    Y = VV(I-1,L-1);
                    VV(I-1,L-1) = VV(I-1,K-1);
                    VV(I-1,K-1) = Y;
                }
                K = K - 1;
            } while ( K > 0 );
            return;
        }
        if ( abs(VV(0,0)) != 0. ) {
            VV(0,0) = COMPLEX(1.,0.)/VV(0,0);
            return;
        }

        throw SCATMECH_exception("Matrix has no inverse (2)");
    }


    //
    // Matrix 1-BA in Eq. 5.3 of B&V is block diagonal. 
    // The following routine inverts that matrix.
    //
    void 
    Axisymmetric_Particle_BRDF_Model::
    invert_block_diagonal(ScatterTMatrix& AA) 
    {
        int m;
        for (m=-MMAX;m<=MMAX;++m) {
            ostringstream msg;
            msg << "Inverting matrix (m = " << m << ") " << ret;
            message(msg.str());

            int n= 2*((m==0) ? LMAX : LMAX-abs(m)+1);
      
            cvert(AA[LMAX+m],n);
        }
        message("Done Inverting Matrices");
    }

    //
    // setup() is called whenever the model needs to be "recalculated"
    //
    void 
    Axisymmetric_Particle_BRDF_Model::
    setup()
    {
        Local_BRDF_Model::setup();   

        Shape->Write("shape.dat");
        double length = Shape->Get_Base_Length();
    
        ostringstream msg;
        msg << "Distance of origin from surface: " << length << endl;
        message(msg.str());

        k   = 2.0*pi/lambda;
        double q = k * length;
        qq  = q+k*delta;
        double qqq = Shape->Get_MaxRadius();
        // Bohren and Huffman recommended lmax...
        BH_LMAX = (int)(qqq+4.05*pow(qqq,0.3333)+2.);
    
        if (lmax==0) LMAX = BH_LMAX;
        if (lmax>0) LMAX = lmax;
        if (lmax<0) LMAX = BH_LMAX+abs(lmax);

        if (mmax==0) MMAX = LMAX;
        if (mmax>0) MMAX = mmax;
        if (mmax<0) MMAX = LMAX+abs(mmax);

        if (MMAX>LMAX) MMAX = LMAX;

        tau = k*t;
        N0  = COMPLEX(substrate.index(lambda));
        N2 = COMPLEX(n2.index(lambda));

        ScatterTMatrix BMatrix(LMAX);

        Bmatrix(BMatrix);

        N1 = n1.index(lambda);
        N2 = n2.index(lambda);
        N3 = n3.index(lambda);
    
        int m,l,l_,f,f_,l__,f__;
    
        if (type==TRANSMISSION) 
             error("No code for TRANSMISSION mode"); 

        if (delta<0.) 
             error("delta cannot be less than zero.");

    
        // The Bohren and Huffman recommended LMAX...
        // BH_LMAX = (int)(q+4.*pow(q,0.3333)+2.);
        BH_LMAX = LMAX;

        sqrsize = index(LMAX,LMAX,1)+1;

        lvector.resize(LMAX+1);
        for (int ll=1;ll<=LMAX;++ll) {
            lvector[ll]=sqrt((2.*ll+1.)/(ll*(ll+1.)));
        }

        message("Allocating Memory");



        //Non-temporary-arrays...
        ScatMatrix.resize(LMAX);
        ScatMatrixInverse.resize(LMAX);

        Zp.resize(sqrsize);
        Zs.resize(sqrsize);
        eIP.resize(sqrsize);
        Wp.resize(sqrsize);
        Ws.resize(sqrsize);
        Vp.resize(sqrsize);
        Vs.resize(sqrsize);

        vector<COMPLEX> B(sqrsize);

        //
        // Calculate A matrix or read it from file and 
        // write it to a file...
        //
        if (order!=0) {        
            int retvalue=1;
            old_LMAX=0;
                
            if (old_LMAX<LMAX) {
                message("Calculating A Matrix");
                Amatrix(ScatMatrix);
            }
        }
    
        //
        // Matrix A becomes B^(-1)*(1-B*A)...
        //
    
        message("Calculating matrix BA");
        for (m=-mmax;m<=mmax;++m) {
         int mm = LMAX+m;
         int beginl = (m==0) ? 1 : abs(m);
         for (l=beginl;l<=LMAX;++l) {
          for (f=0;f<=1;++f) {
           int i = index(l,m,f);
           int ll = 2*(LMAX-l)+f;
           for (l_=beginl;l_<=LMAX;++l_) {
            for (f_=0;f_<=1;++f_) {
             int i_ = index(l_,m,f_);
             int ll_ = 2*(LMAX-l_)+f_;
                 ScatMatrixInverse[mm][ll_][ll]=0.;
             if (order!=0) {
              for (l__=beginl;l__<=LMAX;++l__) {
                   for (f__=0;f__<=1;++f__) {
                      int i__ = index(l__,m,f__);
                      int ll__ = 2*(LMAX-l__)+f__;

                ScatMatrixInverse[mm][ll_][ll] += BMatrix[mm][ll_][ll__]*ScatMatrix[mm][ll__][ll];
               }
              }
             }
            }
           }
          }
         }
        }
    
        message("Calculating matrix 1-BA");

        for (m=-mmax;m<=mmax;++m) {
         int mm = LMAX+m;
         int beginl = (m==0) ? 1 : abs(m);
         for (l=beginl;l<=LMAX;++l) {
          for (f=0;f<=1;++f) {
           int i = index(l,m,f);
           int ll = 2*(LMAX-l)+f;
           for (l_=beginl;l_<=LMAX;++l_) {
            for (f_=0;f_<=1;++f_) {
              int i_ = index(l_,m,f_);
              int ll_ = 2*(LMAX-l_)+f_;
                   ScatMatrixInverse[mm][ll_][ll]=(double)(ll_==ll)-ScatMatrixInverse[mm][ll_][ll];
            }
           }
          }
         }
        }

        message("Inverting B");
        invert_block_diagonal(BMatrix);

        message("Creating B^(-1)(1-BA)");
        for (m=-mmax;m<=mmax;++m) {
         int mm = LMAX+m;
         int beginl = (m==0) ? 1 : abs(m);
         for (l=beginl;l<=LMAX;++l) {
          for (f=0;f<=1;++f) {
           int i = index(l,m,f);
           int ll = 2*(LMAX-l)+f;
           for (l_=beginl;l_<=LMAX;++l_) {
            for (f_=0;f_<=1;++f_) {
              int i_ = index(l_,m,f_);
              int ll_ = 2*(LMAX-l_)+f_;
              ScatMatrix[mm][ll_][ll] =0.;
              for (l__=beginl;l__<=LMAX;++l__) {
                for (f__=0;f__<=1;++f__) {
                  int i__ = index(l__,m,f__);
                  int ll__ = 2*(LMAX-l__)+f__;

                  ScatMatrix[mm][ll_][ll] += BMatrix[mm][ll_][ll__]*ScatMatrixInverse[mm][ll__][ll];
                }
              }
            }
           }
          }
         }
        }

        for (m=-mmax;m<=mmax;++m) {
         int mm = LMAX+m;
         int beginl = (m==0) ? 1 : abs(m);
         for (l=beginl;l<=LMAX;++l) {
          for (f=0;f<=1;++f) {
           int i = index(l,m,f);
           int ll = 2*(LMAX-l)+f;
           for (l_=beginl;l_<=LMAX;++l_) {
            for (f_=0;f_<=1;++f_) {
              int i_ = index(l_,m,f_);
              int ll_ = 2*(LMAX-l_)+f_;
              ScatMatrixInverse[mm][ll_][ll]=ScatMatrix[mm][ll_][ll];
            }
           }
          }
         }
        }

        //
        // Invert the matrix B^(-1)*(1-B*A)
        //
        invert_block_diagonal(ScatMatrix);
    
        // Reset previous geometry...
        old_thetai=-1000;
        old_thetas=-1000;
        old_phis=-1000;
    }

    //
    // The following function must be called before any attempt to evaluate
    // the scattering field.  It sets the geometry for the measurement...
    //
    void
    Axisymmetric_Particle_BRDF_Model::
    set_geometry(double thetai,double thetas,double phis)
    { 
        SETUP();
    
        if (thetai<0) {
            thetai=-thetai;
            phis = pi+phis;
        }
    
        if (thetas<0) {
            thetas=-thetas;
            phis=pi+phis;
        }
    
        if (thetai!=old_thetai) {
            calculate_W(thetai);
            old_thetai=thetai;
        }
    
        if (thetas!=old_thetas) {
            calculate_Z(pi-thetas);
            old_thetas=thetas;
        }
    
        if (phis!=old_phis) {
            calculate_eIP(phis);
            old_phis=phis;
        }
    }

    // The electric field version of the p-->p 
    // differential scattering cross-section...
    COMPLEX
    Axisymmetric_Particle_BRDF_Model::
    Epp(double thetai,double thetas,double phis) 
    {
        set_geometry(thetai,thetas,phis);   
        return E(Wp,Zp);
    }

    // The electric field version of the p-->s 
    // differential scattering cross-section...
    COMPLEX
    Axisymmetric_Particle_BRDF_Model::
    Eps(double thetai,double thetas,double phis) 
    {
        set_geometry(thetai,thetas,phis);   
        return E(Wp,Zs);
    }

    // The electric field version of the s-->p 
    // differential scattering cross-section...
    COMPLEX
    Axisymmetric_Particle_BRDF_Model::
    Esp(double thetai,double thetas,double phis) 
    {
        set_geometry(thetai,thetas,phis);   
        return E(Ws,Zp);
    }

    // The electric field version of the s-->s 
    // differential scattering cross-section...
    COMPLEX
    Axisymmetric_Particle_BRDF_Model::
    Ess(double thetai,double thetas,double phis) 
    {
        set_geometry(thetai,thetas,phis);   
        return E(Ws,Zs);
    }

    //
    // Jones matrix differential scattering cross section...
    //
    JonesMatrix 
    Axisymmetric_Particle_BRDF_Model::
    jonesDSC() 
    {
        set_geometry(thetai,thetas,phis);
       
        COMPLEX pp=E(Wp,Zp);
        COMPLEX ps=E(Wp,Zs);
        COMPLEX sp=E(Ws,Zp);
        COMPLEX ss=E(Ws,Zs);
    
        return JonesMatrix(pp,ss,ps,sp);
    }

    //
    // Constructor for class Axisymmetric_Particle_BRDF_Model...
    //
    Axisymmetric_Particle_BRDF_Model::
    Axisymmetric_Particle_BRDF_Model()
    {
        n3=vacuum_function;
        // Initialize array addresses to zero ...
        ScatMatrix=0;
        ScatMatrixInverse=0;

        // Initialize old angles ...
        old_thetai=-1000.;
        old_thetas=-1000.;
        old_phis=-1000;
        LMAX = 0;
        MMAX = 0;
    }

        
    DEFINE_MODEL(Axisymmetric_Particle_BRDF_Model,Local_BRDF_Model,
               "Axisymmetric_Particle_BRDF_Model",
               "Axisymmetric particle on a substrate.");

    DEFINE_PTRPARAMETER(Axisymmetric_Particle_BRDF_Model,
                        Axisymmetric_Shape_Ptr,
                        Shape,
                        "Particle Shape",
                        "Ellipsoid_Axisymmetric_Shape");
    DEFINE_PARAMETER(Axisymmetric_Particle_BRDF_Model,dielectric_function,n2,"Particle","(1.46,0)");
    DEFINE_PARAMETER(Axisymmetric_Particle_BRDF_Model,dielectric_function,n1,"Substrate Film","(1.59,0)");
    DEFINE_PARAMETER(Axisymmetric_Particle_BRDF_Model,double,t,"Film thickness [um]","0");
    DEFINE_PARAMETER(Axisymmetric_Particle_BRDF_Model,double,delta,"Particle-substrate separation [um]","0");
    DEFINE_PARAMETER(Axisymmetric_Particle_BRDF_Model,int,lmax,"Maximum polar order (LMAX)","0");
    DEFINE_PARAMETER(Axisymmetric_Particle_BRDF_Model,int,mmax,"Maximum azimuthal order (mmax)","0");
    DEFINE_PARAMETER(Axisymmetric_Particle_BRDF_Model,int,order,"Surface interaction order (-1=exact)","-1");
    DEFINE_PARAMETER(Axisymmetric_Particle_BRDF_Model,int,Norm_Inc_Approx,"Normal incidence approximation","0");
    DEFINE_PARAMETER(Axisymmetric_Particle_BRDF_Model,int,improve,"Iterative improvement steps","0");


} // namespace SCATMECH