//******************************************************************************
//** SCATMECH: Polarized Light Scattering C++ Class Library
//** 
//** File: bobvlieg.h
//**
//** 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)
//**
//******************************************************************************

#ifndef SCATMECH_BOBVLIEG_H
#define SCATMECH_BOBVLIEG_H

#include <stdio.h>
#include <vector>
#include "local.h"
#include "scatmatrix.h"
#include "filmtran.h"

namespace SCATMECH {
 
    class Bobbert_Vlieger_BRDF_Model: public Local_BRDF_Model
    {
        public:
    
            COMPLEX Epp(double thetai,double thetas,double phis);
            COMPLEX Eps(double thetai,double thetas,double phis);
            COMPLEX Esp(double thetai,double thetas,double phis);
            COMPLEX Ess(double thetai,double thetas,double phis);
        
            DECLARE_MODEL();
        
            // Index of particle...
            DECLARE_PARAMETER(dielectric_function,n2);   
            // Index of particle film... 
            DECLARE_PARAMETER(dielectric_function,n3);  

            // Radius of particle...
            DECLARE_PARAMETER(double,r);
            // Particle film thickness...
            DECLARE_PARAMETER(double,d);
            // Substrate coatings... 
            DECLARE_PARAMETER(dielectric_stack,stack);
            // Distance between particle and substrate...
            DECLARE_PARAMETER(double,delta);
            // Scattering order (-1 means exact)...
            DECLARE_PARAMETER(int,order);
            // Maximum chosen l value (negative value signifies add to nominal)...
            DECLARE_PARAMETER(int,lmax); 
            // Norm_Inc_Approx =0 for full solution, =1 to assume that r(theta)=r(0)
            DECLARE_PARAMETER(int,Norm_Inc_Approx); 
            // Number of iterative improvement steps
            DECLARE_PARAMETER(int,improve);

        public:
            Bobbert_Vlieger_BRDF_Model();
            
        protected:
            void setup();

            JonesMatrix jonesDSC();

        private:

            COMPLEX N0;          // Index of substrate
            COMPLEX N2;          // Index of particle
            COMPLEX N3;          // Index of particle film

            //Private bookkeeping...
            double  old_thetai;  // Last value of incident angle
            double  old_thetas;  // Last value of scattering angle
            double  old_phis;    // Last value of scattering azimuth angle
            int     LMAX;        // Maximum order of spherical harmonic
            int     BH_LMAX;     // Bohren and Huffman recommended LMAX (used for output)
            double  q;           // 2*PI/lambda*particle radius
            double  qq;          // 2*PI/lambda*particle-surface separation
            double  k;           // 2*PI/lambda
            int     sqrsize;     // Dimension of matrices and vectors
            //int     old_LMAX;
        
            //Data arrays...
            //The following are kept after setup()...
        
            ScatterTMatrix ScatMatrix;          // The scattering matrix
            ScatterTMatrix ScatMatrixInverse;   // The inverse of the scattering matrix
            std::vector<COMPLEX> Wp,Ws;     // The vectors W (depends upon thetai)
            std::vector<COMPLEX> Zp,Zs;     // The vectors Z (depends upon thetas)
            std::vector<COMPLEX> Vp,Vs;     // The vectors V (scattering vector)
            std::vector<COMPLEX> eIP;       // The vector Phi (depends upon phis)
        
            std::vector<double> lvector;    // sqrt((2.*ll+1.)/(ll*(ll+1.)))

            // Function which returns a unique number for the set (l,m,f)
            int index(int l,int m, int f) const {return 2*(l*l+l-1+m)+f;}
            int index(int l,int m) const {return l*l+l-1+m;}
        
            //Private Functions...
            COMPLEX VIp(int l,int m,int f,double thetai) const;
            COMPLEX VIRp(int l,int m,int f,double thetai) const;
            COMPLEX VIs(int l,int m,int f,double thetai) const;
            COMPLEX VIRs(int l,int m,int f,double thetai) const;
        
            COMPLEX Bmatrix(int l_,int m_,int f_,int l,int m,int f) const;
            void    Amatrix(ScatterTMatrix& A);
            void    invert_block_diagonal(ScatterTMatrix& AA);
            void    set_geometry(double thetai,double thetas,double phis);
        
            COMPLEX E(std::vector<COMPLEX>& W,std::vector<COMPLEX>& Z);
            void    calculate_W(double thetai);
            void    calculate_Z(double thetas);
            void    calculate_eIP(double phis);
            COMPLEX rp(COMPLEX cost) const;
            COMPLEX rs(COMPLEX cost) const;
         
            void    iterative_improvement(ScatterTMatrix& Ainv, ScatterTMatrix& A, 
                                          std::vector<COMPLEX>& b,std::vector<COMPLEX>& x);
        
            static COMPLEX arccosine(const COMPLEX& a);
    };

    class Gauss_Laguerre_Integration {
      public:
        static double* weights[];
        static double* zeros[];
    };


    namespace BobVlieg_Supp {
        
        static const int efield=0;
        static const int hfield=1;
        
        double  Fact(int j);
        double  sqrtFact(int j);
        double  mpow(int m);  // (-1)**m
        COMPLEX ipow(int m);  // (i)**m
        COMPLEX LegendreP(int l,int m,COMPLEX x);
        COMPLEX Legendre(int l,int m,COMPLEX x);
        COMPLEX Ptilde(int l,int m,COMPLEX x);
        COMPLEX Jhalf(int l,COMPLEX rho);
        COMPLEX Yhalf(int l,COMPLEX rho);
        COMPLEX j(int l,COMPLEX rho);
        COMPLEX y(int l,COMPLEX rho);
        COMPLEX h(int l,COMPLEX rho);
        COMPLEX psi(int l,COMPLEX rho);
        COMPLEX zeta(int l,COMPLEX rho);
        COMPLEX psi_(int l,COMPLEX rho);
        COMPLEX zeta_(int l,COMPLEX rho);
        COMPLEX chi(int l,COMPLEX rho);
        COMPLEX chi_(int l,COMPLEX rho);  
        COMPLEX dLegendreCosx_dx(int l,int m,COMPLEX cosx);
        COMPLEX P_tilde(int l,int m,COMPLEX x);
        COMPLEX d(int l,int m,int m_,COMPLEX cosalpha2,COMPLEX sinalpha2);
        COMPLEX dplus(int l,int m,COMPLEX cosalpha2,COMPLEX sinalpha2);
        COMPLEX dminus(int l,int m,COMPLEX cosalpha2,COMPLEX sinalpha2);

        static const COMPLEX cI(0.,1.);
    }

} // namespace SCATMECH 



#endif