//******************************************************************************
//** SCATMECH: Polarized Light Scattering C++ Class Library
//** 
//** File: bobvlieg2.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 <assert.h>
#include <vector>
#include "scatmech.h"
#include "axipart.h"
#include "bobvlieg.h"

using namespace std;


namespace SCATMECH {

    using namespace BobVlieg_Supp;

    //
    // The reflection coefficient for p-polarized light (as a function of the 
    // cosine of the complex angle)
    //
    COMPLEX 
    Axisymmetric_Particle_BRDF_Model::
    rp(COMPLEX cosx) const
    {
        COMPLEX cosxf = sqrt(1.-(1.-sqr(cosx))/sqr(N1));
        if (imag(N1*cosxf)<0.) cosxf = -cosxf;
    
        COMPLEX cosxt = sqrt(1.-(1.-sqr(cosx))/sqr(N0));
        if (imag(N0*cosxt)<0.) cosxt = -cosxt;
    
        COMPLEX r12 = (N1*cosx-cosxf)/(N1*cosx+cosxf);
        COMPLEX r23 = (N0*cosxf-N1*cosxt)/(N0*cosxf+N1*cosxt);
        COMPLEX phase23 = exp(2.*cI*tau*N1*cosxf);
        COMPLEX result = (r12+r23*phase23)/(1.0+r12*r23*phase23);
    
        //COMPLEX result = (N0*cosx-cosxt)/(N0*cosx+cosxt);
        return result;
    }

    //
    // The reflection coefficient for s-polarized light (as a function of the 
    // cosine of the complex angle)
    //
    COMPLEX 
    Axisymmetric_Particle_BRDF_Model::
    rs(COMPLEX cosx) const
    {
        COMPLEX cosxf = sqrt(1.-(1.-sqr(cosx))/sqr(N1));
        if (imag(N1*cosxf)<0.) cosxf = -cosxf;
    
        COMPLEX cosxt = sqrt(1.-(1.-sqr(cosx))/sqr(N0)); 
        if (imag(N0*cosxt)<0.) cosxt = -cosxt;

        COMPLEX r12 = (cosx-N1*cosxf)/(cosx+N1*cosxf);
        COMPLEX r23 = (N1*cosxf-N0*cosxt)/(N1*cosxf+N0*cosxt);
        COMPLEX phase23 = exp(2.*cI*tau*N1*cosxf);
        COMPLEX result = (r12+r23*phase23)/(1.0+r12*r23*phase23);

        //COMPLEX result = (cosx-N0*cosxt)/(cosx+N0*cosxt);
        return result;
    }    

    //
    // The element of the incident plane wave vector for p-polarized light
    // from Eq. (2.14) of BV&G...
    //
    COMPLEX 
    Axisymmetric_Particle_BRDF_Model::
    VIp(int l,int m,int f,double thetai) const
    {
        COMPLEX result;
        COMPLEX costheta2=cos(thetai/2.);
        COMPLEX sintheta2=sin(thetai/2.);
        if (f==efield) {
            // Eq. (2.14a) of BV&G
            result =  1./k*ipow(l-1)*lvector[l]*mpow(m-1)*
                            dminus(l,m,costheta2,sintheta2);
        } else {
            // Eq. (2.14b) of BV&G
            result =  -1./k*ipow(l)*lvector[l]*mpow(m-1)*
                            dplus(l,m,costheta2,sintheta2);
        }
        return result;
    }

    //
    // The element of the incident plane wave vector for p-polarized light
    // reflected from the surface from Eq. (2.15) of BV&G...
    //
    COMPLEX 
    Axisymmetric_Particle_BRDF_Model::
    VIRp(int l,int m,int f,double thetai) const
    {
        COMPLEX result;
        COMPLEX cost = cos(thetai);
        COMPLEX phase = exp(2.*cI*qq*cost);
        COMPLEX rp = Axisymmetric_Particle_BRDF_Model::rp(cost);
        COMPLEX costheta2=cos(thetai/2.);
        COMPLEX sintheta2=sin(thetai/2.);
    
        if (f==efield) {
            // Eq. (2.15a) of BV&G
            result = 1./k*phase*rp*
                   ipow(l-1)*lvector[l]*mpow(l)*
                                dminus(l,m,costheta2,sintheta2);
        } else {
            // Eq. (2.15b) of BV&G
            result =  -1./k*phase*rp*
                   ipow(l)*lvector[l]*mpow(l-1)*
                                dplus(l,m,costheta2,sintheta2);
        }
        return result;
    }

    //
    // The element of the incident plane wave vector for s-polarized light
    // from Eq. (2.16) of BV&G...
    //
    COMPLEX 
    Axisymmetric_Particle_BRDF_Model::
    VIs(int l,int m,int f,double thetai) const
    {
        COMPLEX result;
        COMPLEX costheta2=cos(thetai/2.);
        COMPLEX sintheta2=sin(thetai/2.);
    
        if (f==efield) {
            // Eq. (2.16a) of BV&G
            result =  -cI/k*
                    ipow(l-1)*lvector[l]*mpow(m-1)*
                            dplus(l,m,costheta2,sintheta2);
        } else {
            // Eq. (2.16b) of BV&G
            result = cI/k*
                    ipow(l)*lvector[l]*mpow(m-1)*
                            dminus(l,m,costheta2,sintheta2);
        }
        return result;    
    }

    //
    // The element of the incident plane wave vector for s-polarized light
    // reflected from the surface from Eq. (2.17) of BV&G...
    //
    COMPLEX 
    Axisymmetric_Particle_BRDF_Model::
    VIRs(int l,int m,int f,double thetai) const
    {
        COMPLEX result;
        COMPLEX cost = cos(thetai);
        COMPLEX phase = exp(2.*cI*qq*cost);
        COMPLEX rs = Axisymmetric_Particle_BRDF_Model::rs(cost);
        COMPLEX costheta2=cos(thetai/2.);
        COMPLEX sintheta2=sin(thetai/2.);
    
        if (f==efield) {
            // Eq. (2.17a) of BV&G
            result =  -cI/k*phase*rs*
                   ipow(l-1)*lvector[l]*
                    mpow(l-1)*dplus(l,m,costheta2,sintheta2);
        } else {
            // Eq. (2.17b) of BV&G
            result =  cI/k*phase*rs*
                    ipow(l)*lvector[l]*mpow(l)*
                      dminus(l,m,costheta2,sintheta2);
        }
        return result;    
    }

    //
    // The following should iteratively improve the solution to A.x=b
    //
    void
    Axisymmetric_Particle_BRDF_Model::
    iterative_improvement(ScatterTMatrix& Ainv, ScatterTMatrix& A, vector<COMPLEX>& b, vector<COMPLEX>& x)
    {
        vector<COMPLEX> temp1(sqrsize);
        int m,l,f,l_,f_;
        
        for (m=-MMAX;m<=MMAX;++m) {
            COMPLEX temp2;
            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_);
                    temp1[i_]=0;
                    for (l=beginl;l<=LMAX;++l) {
                        for (f=0;f<=1;++f) {
                            int i = index(l,m,f);
                            int mm = LMAX+m;
                            int ll = 2*(LMAX-l)+f;
                            int ll_ = 2*(LMAX-l_)+f_;
                        
                            temp1[i_]+=A[mm][ll_][ll]*x[i];
                        }
                    }
                    temp1[i_]-=b[i_];
                }
            }
            for (l_=beginl;l_<=LMAX;++l_) {
                for (f_=0;f_<=1;++f_) {
                    int i_ = index(l_,m,f_);
                    temp2=0;
                    for (l=beginl;l<=LMAX;++l) {
                        for (f=0;f<=1;++f) {
                            int i = index(l,m,f);
                            int mm = LMAX+m;
                            int ll = 2*(LMAX-l)+f;
                            int ll_ = 2*(LMAX-l_)+f_;
                        
                            temp2+=Ainv[mm][ll_][ll]*temp1[i];
                        }
                    }
                    x[i_]-=temp2;
                }
            }
        }
    }

    //
    // Routine to calculate the scattered wave vector for a specific incident
    // angle.  This is evaluation of Eq. (5.3) of B&V...
    //
    void 
    Axisymmetric_Particle_BRDF_Model::
    calculate_W(double thetai)
    {
        int i,m,l,f,l_,f_;
    
        // First calculate the V vector of the incident wave...
        for (l=1;l<=LMAX;++l) {
            for (m=-l;m<=l;++m) {
                for (f=0;f<=1;++f) {
                    i=index(l,m,f);
                    Vp[i]=VIp(l,m,f,thetai)+VIRp(l,m,f,thetai);
                    Vs[i]=VIs(l,m,f,thetai)+VIRs(l,m,f,thetai);
                }
            }
        }
    
        // Then multiply the incident vector V by the scattering matrix 
        for (m=-MMAX;m<=MMAX;++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_);
                    Wp[i_]=0;
                    Ws[i_]=0;
                    for (l=beginl;l<=LMAX;++l) {
                        for (f=0;f<=1;++f) {
                            int i = index(l,m,f);
                            int mm = LMAX+m;
                            int ll = 2*(LMAX-l)+f;
                            int ll_ = 2*(LMAX-l_)+f_;
                        
                            Wp[i_]+=ScatMatrix[mm][ll_][ll]*Vp[i];
                            Ws[i_]+=ScatMatrix[mm][ll_][ll]*Vs[i];
                        }
                    }
                }
            }
        }
    
        // Iteratively improve the solution, since the matrix 
        // inversion isn't perfect...
        if (order<0) {
            for (i=0;i<improve;++i) {
                iterative_improvement(ScatMatrix,ScatMatrixInverse,Vp,Wp);
                iterative_improvement(ScatMatrix,ScatMatrixInverse,Vs,Ws);
            }
        }
    }

    //
    // The vector Z is that part of the scattering function which depends
    // upon the scattering angle thetas.  It must be evaluated anytime 
    // that thetas changes...
    //
    void
    Axisymmetric_Particle_BRDF_Model::
    calculate_Z(double thetas) 
    {
        double d_cost = cos(thetas);
        COMPLEX cost= d_cost;
        double sint = real(sqrt(1.-sqr(cost)));
       
        COMPLEX _cost=cos(pi-thetas);
        COMPLEX rp = Axisymmetric_Particle_BRDF_Model::rp(_cost);
        COMPLEX rs = Axisymmetric_Particle_BRDF_Model::rs(_cost);    
        COMPLEX phase = exp(2.*cI*qq*_cost);
        COMPLEX rpphase = rp*phase;
        COMPLEX rsphase = rs*phase;

        if (d_cost > -0.999999) {
            for (int i=0,l=1;l<=LMAX;++l) {
                if (l<=BH_LMAX) {
                    int _mmax = (MMAX>l) ? l : MMAX;
                    for (int m=-_mmax;m<=_mmax;++m) {
                        COMPLEX temp1 = Ptilde(l,m,cost)*(double)m/sint;
                        COMPLEX temp2 = P_tilde(l,m,cost);
                        for (int f=0;f<=1;++f,++i) {
                            if (f==efield) {
                                Zp[i] = ipow(l)*mpow(l)*temp2; 
                                Zs[i] = -ipow(l+1)*mpow(l+1)*temp1;
                                Zp[i] = Zp[i]*(1.+mpow(l-m+1)*rpphase);
                                Zs[i] = Zs[i]*(1.+mpow(l-m)*rsphase);
                            } else {
                                Zp[i] = -ipow(l+1)*mpow(l+1)*temp1;
                                Zs[i] = -ipow(l)*mpow(l)*temp2;
                                Zp[i] = Zp[i]*(1.+mpow(l-m)*rpphase);
                                Zs[i] = Zs[i]*(1.+mpow(l-m+1)*rsphase);
                            }
                        }
                    }
                } else {
                    int _mmax = (MMAX>l) ? l : MMAX;
                    for (int m=-_mmax;m<=_mmax;++m) {
                        for (int f=0;f<=1;++f,++i) {
                                Zp[i] = 0.; 
                                Zs[i] = 0.;
                                Zp[i] = 0.;
                                Zs[i] = 0.;
                        }
                    }
                }
            }
        } else {  // Take limits as cos(thetas) -> -1...
            for (int i=0,l=1;l<=LMAX;++l) {
                double temp = sqrt((double)((2*l+1)*(l+1)*l))/2.;
                  int _mmax = (MMAX>l) ? l : MMAX;
                  for (int m=-_mmax;m<=_mmax;++m) {
                    for (int f=0;f<=1;++f,++i) {
                        if (l<=BH_LMAX && (m==-1 || m==1)) {
                            if (f==efield) {
                                Zp[i] = ipow(l)*mpow(2*l+1)*temp*(double)m; 
                                Zs[i] = -ipow(l+1)*mpow(2*l+1)*temp;
                                Zp[i] = Zp[i]*(1.+mpow(l-m+1)*rpphase);
                                Zs[i] = Zs[i]*(1.+mpow(l-m)*rsphase);
                            } else {
                                Zp[i] = -ipow(l+1)*mpow(2*l+1)*temp;
                                Zs[i] = -ipow(l)*mpow(2*l+1)*temp*(double)m;
                                Zp[i] = Zp[i]*(1.+mpow(l-m)*rpphase);
                                Zs[i] = Zs[i]*(1.+mpow(l-m+1)*rsphase);
                            }
                        } else {
                                Zp[i] = 0.; 
                                Zs[i] = 0.;
                                Zp[i] = 0.;
                                Zs[i] = 0.;
                        }
                    }
                }
            }
        }
    
    }

    //
    // The vector eIP is that part of the scattering function which depends
    // upon the azimuthal scattering angle phis.  It must be evaluated anytime 
    // that phis changes...
    //
    // This function basically calculates exp(i m phis)...
    //
    void 
    Axisymmetric_Particle_BRDF_Model::
    calculate_eIP(double phis) 
    {
        vector<COMPLEX> expmphi(2*LMAX+2);
    
        COMPLEX expphi=exp(cI*phis);
        expmphi[LMAX]=1.;
    
        for (int m=1;m<=LMAX;++m) {
            expmphi[LMAX+m]=expmphi[LMAX+m-1]*expphi;
            expmphi[LMAX-m]=expmphi[LMAX-m+1]/expphi;
        }
    
        for (int i=0,l=1;l<=BH_LMAX;++l) {
            int _mmax = (MMAX>l) ? l : MMAX;
            for (int m=-_mmax;m<=_mmax;++m) {
                for (int f=0;f<=1;++f,++i) {
                    eIP[i] = expmphi[m+LMAX];
                }
            }
        }
    }

    //
    // The total scattered electric field depends upon the incident vector W and
    // the scattered vector Z. 
    //
    COMPLEX 
    Axisymmetric_Particle_BRDF_Model::
    E(vector<COMPLEX>& W,vector<COMPLEX>& Z) 
    {
        COMPLEX E(0.);
    
        for (int i=0,l=1;l<=BH_LMAX;++l) {
            int _mmax = (MMAX>l) ? l : MMAX;
            for (int m=-_mmax;m<=_mmax;++m) {
                for (int f=0;f<=1;++f,++i) {
                    E += W[i]*Z[i]*eIP[i];          
                }
            }
        }
    
        return E;
    }


} // namespace SCATMECH