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

using namespace std;


namespace SCATMECH {


    //
    // The function which returns the Jones matrix for scattering...
    //
    // The scattering matrix is described by Thomas A. Germer in "Polarized light
    // scattering by microroughness and small defects in dielectric layers," 
    // J. Opt. Soc. Am. A 18(6), 1279-1288 (2001).
    //
    JonesMatrix 
    Two_Face_BRDF_Model::
    jones()
    {
        // All BRDF_Model::jones(,,) should call setup() if any parameters 
        // have changed...
        SETUP();
    
        if (type==TRANSMISSION) {
             error("No code for TRANSMISSION mode"); 
        }
    
        COMPLEX i(0,1);   
        COMPLEX eps1=substrate.epsilon(lambda);
        COMPLEX eps2=film.epsilon(lambda);
        double  tau=thickness;
        double  k=2.*pi/lambda;
        double  ki=k*sin(thetai);
        double  kr=k*sin(thetas);
        COMPLEX qi1=k*sqrt(eps1-sqr(sin(thetai)));
        COMPLEX qi2=k*sqrt(eps2-sqr(sin(thetai)));
        double  qi3=k*cos(thetai);
        COMPLEX qr1=k*sqrt(eps1-sqr(sin(thetas)));
        COMPLEX qr2=k*sqrt(eps2-sqr(sin(thetas)));
        double  qr3=k*cos(thetas);
    
        COMPLEX Ki_plus = exp(2.*i*qi2*tau)+1.;
        COMPLEX Kr_plus = exp(2.*i*qr2*tau)+1.;
        COMPLEX Ki_minus = Ki_plus-2.;
        COMPLEX Kr_minus = Kr_plus-2.;
    
        COMPLEX F_pi_plus  = eps2*Ki_minus*qi1- eps1*Ki_plus*qi2;
        COMPLEX F_pr_plus  = eps2*Kr_minus*qr1- eps1*Kr_plus*qr2;
        COMPLEX F_pi_minus = eps2*Ki_plus*qi1 - eps1*Ki_minus*qi2;
        COMPLEX F_pr_minus = eps2*Kr_plus*qr1 - eps1*Kr_minus*qr2;

        COMPLEX F_si_plus  = Ki_minus*qi1 - Ki_plus*qi2;
        COMPLEX F_sr_plus  = Kr_minus*qr1 - Kr_plus*qr2;
        COMPLEX F_si_minus = Ki_plus*qi1  - Ki_minus*qi2;
        COMPLEX F_sr_minus = Kr_plus*qr1  - Kr_minus*qr2;    
    
        COMPLEX Gpi = eps2*F_pi_plus*qi3 - F_pi_minus*qi2;
        COMPLEX Gpr = eps2*F_pr_plus*qr3 - F_pr_minus*qr2;
        COMPLEX Gsi =      F_si_plus*qi3 - F_si_minus*qi2;
        COMPLEX Gsr =      F_sr_plus*qr3 - F_sr_minus*qr2;
    
        double fx,fy;
        Bragg_Frequency(fx,fy);
    
        JonesMatrix j1,j2;
     
        // Buried interface...
        if (layer==1) {
            j1.PP() = qi2*qr2*eps2*(-eps1*ki*kr+eps2*qi1*qr1*cos(phis))/(Gpi*Gpr);
            j1.PS() = eps2*k*qi2*qr2*qi1*sin(phis)/(Gpi*Gsr);
            j1.SP() = eps2*k*qi2*qr2*qr1*sin(phis)/(Gsi*Gpr);
            j1.SS() = -sqr(k)*qi2*qr2*cos(phis)/(Gsi*Gsr);
        
            COMPLEX temp = 4.*(eps1-eps2)*exp(i*(qi2+qr2-qr3-qi3)*tau)*qr3*qi3/(COMPLEX)(pi);
            temp = temp*sqrt(psd->psd(fx,fy)/(cos(thetai)*cos(thetas)));
            j1 = j1*temp;
        }

        if (layer==1) return j1;
    
        // Exposed interface...
        if (layer==2) {
            j2.PP() = (-eps2*F_pi_plus*F_pr_plus*ki*kr+
                        F_pi_minus*F_pr_minus*qi2*qr2*cos(phis))/(Gpi*Gpr);
            j2.PS() = -F_pi_minus*F_sr_plus*k*qi2*sin(phis)/(Gpi*Gsr);
            j2.SP() = -F_pr_minus*F_si_plus*k*qr2*sin(phis)/(Gsi*Gpr);
            j2.SS() = -sqr(k)*F_si_plus*F_sr_plus*cos(phis)/(Gsi*Gsr);
        
            COMPLEX temp = (eps2-1.)*qr3*qi3*exp(-i*(qi3+qr3)*tau)/((COMPLEX)(pi));
            temp = temp*sqrt(psd->psd(fx,fy)/(cos(thetai)*cos(thetas)));
            j2 = j2*temp;
        }
    
        if (layer==2) return j2;
    
        error("Invalid interface number in jones(,,)");
        return JonesMatrix(0.,0.,0.,0.);
    }

    DEFINE_MODEL(Two_Face_BRDF_Model,Roughness_BRDF_Model,
                "Two_Face_BRDF_Model",
                "Scattering by a single rough interface of a single dielectric film.");


    DEFINE_PARAMETER(Two_Face_BRDF_Model,dielectric_function,film,"Layer","(1.46,0)");

    DEFINE_PARAMETER(Two_Face_BRDF_Model,double,thickness,"Thickness [um]","0.05");

    DEFINE_PARAMETER(Two_Face_BRDF_Model,int,layer,"Interface number (buried:1, exposed:2)","1");




} // namespace SCATMECH