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

using namespace std;


namespace SCATMECH {


    static
    COMPLEX 
    ArcSin(const COMPLEX& a)
    {    
        COMPLEX temp = sqrt(1. - sqr(a)) + COMPLEX(0,1)*a;
        return COMPLEX(arg(temp),-log(abs(temp)));
    }

    void
    Subsurface_Particle_BRDF_Model::
    setup() 
    {
        Local_BRDF_Model::setup();
        if (substrate.k(lambda)!=0) 
            error("Substrate cannot be absorbing"); 
        if (scatterer->get_lambda()!=lambda) scatterer->set_lambda(lambda);
        if (scatterer->get_medium().n(lambda)!=substrate.n(lambda)) scatterer->set_medium(substrate);

        Euler = Matrix(cos(alpha),      -(cos(beta)*sin(alpha)),          sin(alpha)*sin(beta),
                          sin(alpha),      cos(alpha)*cos(beta),           -(cos(alpha)*sin(beta)),
                          0,                      sin(beta),                      cos(beta));
    }

    Vector 
    Subsurface_Particle_BRDF_Model::
    rotatez(const Vector& v) const
    {
        Vector result;
        result.x = cos(rotation)*v.x+sin(rotation)*v.y;
        result.y = -sin(rotation)*v.x+cos(rotation)*v.y;
        result.z = v.z;
        return result;
    }

    JonesMatrix
    Subsurface_Particle_BRDF_Model::
    jonesDSC()
    {
        SETUP();

        double n = substrate.n(lambda);

        double k=2.*pi/lambda;

        // Angles of refraction internal to the material...
        double sin_thetai_internal = sin(thetai)/n; 
        double cos_thetai_internal = sqrt(1.-sqr(sin_thetai_internal));
    
        Vector inKi(sin_thetai_internal,0.,-cos_thetai_internal);

        // The surface normal...    
        Vector normal(0.,0.,1.);
        Vector inSi =  perpto(inKi,normal);
        Vector inPi =  perpto(inKi,inSi);
    
    
        if (type==REFLECTION) {
              
            double sin_thetas_internal = sin(thetas)/n; 
            double cos_thetas_internal = sqrt(1.-sqr(sin_thetas_internal));

            Vector outKi(sin_thetas_internal*cos(phis),
                          sin_thetas_internal*sin(phis),
                          cos_thetas_internal);
        
            // The following definitions make {s,p,k} right handed...
            Vector outSi = perpto(outKi,normal);
            Vector outPi = perpto(outKi,outSi);
                
            JonesMatrix ti = stack.t12(thetai,lambda,vacuum,substrate)*(COMPLEX)sqrt(n*cos_thetai_internal/cos(thetai));
            JonesMatrix ts = stack.t12(thetas,lambda,vacuum,substrate)*(COMPLEX)sqrt(n*cos_thetas_internal/cos(thetas));

            Vector perp = perpto(inKi,outKi);  // Polarization unit vectors in local 
            Vector pari = cross(inKi,perp);   // coordinates.
            Vector pars = cross(outKi,perp);

            JonesMatrix matrixin = JonesMatrix(perp*inPi,pari*inSi,pari*inPi,perp*inSi);
            JonesMatrix matrixout = JonesMatrix(outPi*perp,outSi*pars,outSi*perp,outPi*pars);

            //double scatterangle = acos(inKi*outKi);
            JonesMatrix scatter = scatterer->jones(Euler*rotatez(inKi),Euler*rotatez(outKi));
            scatter=matrixout*scatter*matrixin;

            double j1 = cos(thetas)/cos_thetas_internal/sqr(n);
            double j2 = cos(thetai)/cos_thetai_internal;

            // The scatter matrix globally includes transmission to/from defect...
            scatter = ts*scatter*ti;
        
            COMPLEX common=j1*j2/sqr(k*n);
            // Return the whole value with factors common to all elements...
            return scatter*sqrt(common);
      
        } else { // TRANSMISSION...

            double sin_thetas = sin(thetas); 
            double cos_thetas = sqrt(1.-sqr(sin_thetas));

            double sin_thetas_external = sin_thetas*n;
            COMPLEX thetas_external = ArcSin((COMPLEX)sin_thetas_external);

            JonesMatrix scatter_direct;
            {
                Vector outKi(sin_thetas*cos(phis),sin_thetas*sin(phis),-cos_thetas);
                Vector outSi = perpto(outKi,normal);
                Vector outPi = cross(outKi,outSi);

                Vector perp = perpto(inKi,outKi);  // Polarization unit vectors in local 
                Vector pari = cross(inKi,perp);   // coordinates.
                Vector pars = cross(outKi,perp);

                JonesMatrix matrixin = JonesMatrix(perp*inPi,pari*inSi,pari*inPi,perp*inSi);
                JonesMatrix matrixout = JonesMatrix(outPi*perp,outSi*pars,outSi*perp,outPi*pars);

                //double scatterangle = acos(inKi*outKi);
                JonesMatrix scatter = scatterer->jones(Euler*rotatez(inKi),Euler*rotatez(outKi));

                scatter_direct=matrixout*scatter*matrixin;
            }

            JonesMatrix scatter_indirect; 
            {
                Vector outKr(sin_thetas*cos(phis),sin_thetas*sin(phis),cos_thetas);
                Vector outSr = perpto(outKr,normal);
                Vector outPr = cross(outKr,outSr);

                Vector perp = perpto(inKi,outKr);  // Polarization unit vectors in local 
                Vector pari = cross(inKi,perp);   // coordinates.
                Vector pars = cross(outKr,perp);

                JonesMatrix matrixin = JonesMatrix(perp*inPi,pari*inSi,pari*inPi,perp*inSi);
                JonesMatrix matrixout = JonesMatrix(outPr*perp,outSr*pars,outSr*perp,outPr*pars);

                //double scatterangle = acos(inKi*outKr);
                JonesMatrix scatter = scatterer->jones(Euler*rotatez(inKi),Euler*rotatez(outKr));

                JonesMatrix r = stack.r21(thetas_external,lambda,substrate,vacuum);
                COMPLEX phase = exp(COMPLEX(0,2)*cos_thetas*k*depth*n);

                scatter_indirect=matrixout*scatter*matrixin*r*phase;
            }

            JonesMatrix ti = stack.t12(thetai,lambda,vacuum,substrate)*(COMPLEX)sqrt(n);
        
            JonesMatrix scatter = (scatter_direct+scatter_indirect)*ti/k/n;
        
            return scatter;
        }
    }

    DEFINE_MODEL(Subsurface_Particle_BRDF_Model,Local_BRDF_Model,
                "Subsurface_Particle_BRDF_Model",
                "Scattering by a particle beneath an interface of a nonabsorbing media");

    DEFINE_PARAMETER(Subsurface_Particle_BRDF_Model,dielectric_stack,stack,"Film stack on substrate","");
    DEFINE_PARAMETER(Subsurface_Particle_BRDF_Model,double,depth,"Depth of center of particle [um]","0");
    DEFINE_PTRPARAMETER(Subsurface_Particle_BRDF_Model,Free_Space_Scatterer_Ptr,scatterer,"The scattering function","MieScatterer");
    DEFINE_PARAMETER(Subsurface_Particle_BRDF_Model,double,alpha,"Rotation of particle about z-axis [rad]","0");
    DEFINE_PARAMETER(Subsurface_Particle_BRDF_Model,double,beta,"Rotation of particle about y-axis [rad]","0");



} //namespace SCATMECH