//******************************************************************************
//** SCATMECH: Polarized Light Scattering C++ Class Library
//** 
//** File: axifree.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 "axifree.h"
#include "dielfunc.h"
#include "askuser.h"
#include "matrix3d.h"
#include "float.h"

using namespace std;

namespace SCATMECH {

    using namespace BobVlieg_Supp;

    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));
    }


    //
    // Bmatrix() calculates the Mie scattering matrix found 
    // in Eq. (4.2) of B&V. 
    //
    void 
    TMatrix_Axisymmetric_Scatterer::
    Bmatrix(ScatterTMatrix& T)
    {
        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;
            }
           }
          }
         }
        }

        COMPLEX relN2 = N2/N0;
        Shape->Get_TMatrix(T, LMAX, MMAX, k*N0, relN2);

        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)));
    }

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

        Shape->Write("shape.dat");
        double length = Shape->Get_Base_Length();
    
        k = 2.0*pi/lambda;
        double q = k * Shape->Get_MaxRadius();

        // Bohren and Huffman suggested lmax...
        BH_LMAX = (int)(q+4.05*pow(q,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;

        if (medium.k(lambda)!=0) error("Surrounding medium must be non-absorbing.");

        N2 = COMPLEX(particle.index(lambda));
        N0 = medium.n(lambda);

        ScatMatrix.resize(LMAX);
        Bmatrix(ScatMatrix);
        
        // 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.)));
        }

        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);
    
        // 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
    TMatrix_Axisymmetric_Scatterer::
    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>=pi) {
            thetai-=2*pi;
        }
        if (thetai<-pi) {
            thetai+=2*pi;
        }

        if (thetas>=pi) {
            thetas-=2*pi;
        }

        if (thetas<-pi) {
            thetas+=2*pi;
        }

        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;
        }
    }

    JonesMatrix 
    TMatrix_Axisymmetric_Scatterer::
    jones(const Vector& kin,const Vector& kout)
    {
        SETUP();
        double _thetai = acos(-kin.z);
        double _thetas = acos(kout.z);

        double _phis =  ((kout.y==0 && kout.x==0) || (kin.x==0 && kin.y==0)) ?
                            0 : atan2(kout.y,kout.x) - atan2(-kin.y,-kin.x) - pi;

        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);
    
        JonesMatrix scatter = JonesMatrix(pp,ss,ps,sp)*k;

        Vector sr = perpto(kout,Vector(0,0,1));
        Vector si = perpto(kin,Vector(0,0,1));
        Vector pr = perpto(kout,sr);
        Vector pi = perpto(kin,si);

        Vector ai = perpto(kin,kout);
        Vector ar = ai;
        Vector bi = perpto(kin,ai);
        Vector br = perpto(kout,ar);

        JonesMatrix rotatein(ai*pi,bi*si,bi*pi,ai*si);
        JonesMatrix rotateout(ar*pr,br*sr,br*pr,ar*sr);

        return rotateout*scatter*rotatein;
    }

    //
    // Constructor for class TMatrix_Axisymmetric_Scatterer...
    //
    TMatrix_Axisymmetric_Scatterer::
    TMatrix_Axisymmetric_Scatterer()
    {
    }

    //
    // The element of the incident plane wave vector for p-polarized light
    // from Eq. (2.14) of BV&G...
    //
    COMPLEX 
    TMatrix_Axisymmetric_Scatterer::
    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 s-polarized light
    // from Eq. (2.16) of BV&G...
    //
    COMPLEX 
    TMatrix_Axisymmetric_Scatterer::
    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;    
    }

    //
    // Routine to calculate the scattered wave vector for a specific incident
    // angle.  This is evaluation of Eq. (5.3) of B&V...
    //
    void 
    TMatrix_Axisymmetric_Scatterer::
    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);
                    Vs[i]=VIs(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];
                        }
                    }
                }
            }
        }
    }

    //
    // 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
    TMatrix_Axisymmetric_Scatterer::
    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);

        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;
                            } else {
                                Zp[i] = -ipow(l+1)*mpow(l+1)*temp1;
                                Zs[i] = -ipow(l)*mpow(l)*temp2;
                            }
                        }
                    }
                } 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;
                            } else {
                                Zp[i] = -ipow(l+1)*mpow(2*l+1)*temp;
                                Zs[i] = -ipow(l)*mpow(2*l+1)*temp*(double)m;
                            }
                        } 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 
    TMatrix_Axisymmetric_Scatterer::
    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 
    TMatrix_Axisymmetric_Scatterer::
    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;
    }

       
    DEFINE_MODEL(TMatrix_Axisymmetric_Scatterer,Free_Space_Scatterer,
               "TMatrix_Axisymmetric_Scatterer",
               "Axisymmetric particle solved using T-Matrix method");

    DEFINE_PTRPARAMETER(TMatrix_Axisymmetric_Scatterer,
                        Axisymmetric_Shape_Ptr,Shape,
                        "Particle Shape",
                        "Ellipsoid_Axisymmetric_Shape");
    DEFINE_PARAMETER(TMatrix_Axisymmetric_Scatterer,dielectric_function,particle,"Particle","(1.59,0)");
    DEFINE_PARAMETER(TMatrix_Axisymmetric_Scatterer,int,lmax,"Maximum polar order","0");
    DEFINE_PARAMETER(TMatrix_Axisymmetric_Scatterer,int,mmax,"Maximum azimuthal order","0");


} // namespace SCATMECH