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

using namespace std;


namespace SCATMECH {


    //
    // Jones matrix for scattering...
    //
    JonesMatrix 
    Correlated_Roughness_BRDF_Model::
    jones()
    {
        SETUP();

        if (type == REFLECTION) {
            // From J. M. Elson, 
            //   "Light scattering from surfaces with a single dielectric
            //    overlayer," J. Opt. Soc. Am. 66 (7) 682--694 (1976)
            // ---Solves problem of correlated roughness in a single film
            // 
            // Comments:
            //
            //  * Error found in manuscript(?) Eq. (2.26b): Should K_0 by K_0minus? 
            //    (Correction is assumed in code below.)
            //
            //  * Since the manuscript has modulus-squared on the p vs. s 
            //    scattering efficiency, a factor of -i has been added to the ps 
            //    and ss terms in order to match the results for microroughness 
            //    in the limit of tau->0 or epsilon0->1 or epsilon0->epsilon.
            //
            //  * Note that phis is phi in paper,
            //    and phi is phi(tau) in paper.
            // 
        
            // All routines should call setup()...
            SETUP();
        
            double tau = 2*pi/lambda*thickness;
        
            COMPLEX I(0,1);
            double theta0=thetai;
            double theta =thetas;
        
            COMPLEX eps0 = (COMPLEX)(film.epsilon(lambda));
            COMPLEX eps = (COMPLEX)(substrate.epsilon(lambda));
        
            COMPLEX q  =cos(theta0);
            COMPLEX k0 =sin(theta0);
            COMPLEX q0 =sqrt(eps0-sqr(k0));
            // Error found in SCATMECH v. 1.0 (Fixed TAG 30 JUNE 2000)
            // Change the following
            // COMPLEX nu0=sqrt(sqr(k0)-eps);
            // to
            COMPLEX nu0=-I*sqrt(eps-sqr(k0));
        
            COMPLEX k=sin(theta);
            COMPLEX eta=-I*cos(theta);
            COMPLEX eta0=-I*sqrt(eps0-sqr(k));
            // Error found in SCATMECH v. 1.0 (Fixed TAG 30 JUNE 2000)
            // Change the following
            // COMPLEX nu=sqrt(sqr(k)-eps);
            // to
            COMPLEX nu=-I*sqrt(eps-sqr(k));
        
            COMPLEX Gamma_0plus =(nu0*eps0-I*q0*eps)*exp(-I*q0*tau)+
                                 (nu0*eps0+I*q0*eps)*exp(I*q0*tau);
            COMPLEX Gamma_0minus=(nu0*eps0-I*q0*eps)*exp(-I*q0*tau)-
                                 (nu0*eps0+I*q0*eps)*exp(I*q0*tau);
            COMPLEX Gamma_plus= (nu*eps0+eta0*eps)*exp(eta0*tau)+
                                (nu*eps0-eta0*eps)*exp(-eta0*tau);
            COMPLEX Gamma_minus=(nu*eps0+eta0*eps)*exp(eta0*tau)-
                                (nu*eps0-eta0*eps)*exp(-eta0*tau);

        
            COMPLEX K_0plus = (nu0-I*q0)*exp(-I*q0*tau)+(nu0+I*q0)*exp(I*q0*tau);
            COMPLEX K_0minus= (nu0-I*q0)*exp(-I*q0*tau)-(nu0+I*q0)*exp(I*q0*tau);
            COMPLEX K_plus  = (nu+eta0)*exp(eta0*tau)+(nu-eta0)*exp(-eta0*tau);
            COMPLEX K_minus = (nu+eta0)*exp(eta0*tau)-(nu-eta0)*exp(-eta0*tau);
        
            COMPLEX phi = (nu*eps0+eta0*eps)*
                          (eta0+eta*eps0)*exp(eta0*tau)+
                          (nu*eps0-eta0*eps)*
                          (eta0-eta*eps0)*exp(-eta0*tau);
            COMPLEX beta = (eta+eta0)*(nu+eta0)*exp(eta0*tau)+
                           (eta0-eta)*(nu-eta0)*exp(-eta0*tau);
        
            COMPLEX phi0 =  (nu0*eps0-I*q0*eps)*
                            (q0+q*eps0)*exp(-I*q0*tau)+
                            (nu0*eps0+I*q0*eps)*
                            (q0-q*eps0)*exp(I*q0*tau);
            COMPLEX beta0 = (q+q0)*(nu0-I*q0)*exp(-I*q0*tau)+
                            (q0-q)*(nu0+I*q0)*exp(I*q0*tau);
           
            COMPLEX B_theta=q0*(4.*eta0*nu0*eps0*(eps0-eps)+
                            (1.-eps0)*Gamma_0plus*K_minus);
            COMPLEX B_phi=4.*eta0*q0*(eps0-eps)+
                          I*K_minus*K_0minus*(1.-eps0);
            COMPLEX P_theta=4.*eta0*q0*eps0*(eps0-eps)*
                            (k*k0*eps+nu*nu0*eps0*cos(phis))
                            +(I*k*k0*eps0*Gamma_minus*Gamma_0minus+
                            eta0*q0*Gamma_plus*Gamma_0plus*cos(phis))*(1.-eps0);
            COMPLEX P_phi=-eta0*(4.*q0*eps0*nu*(eps0-eps)+
                          I*(1.-eps0)*Gamma_plus*K_0minus);
        
            COMPLEX pp = P_theta/phi0/phi;
            COMPLEX sp = I*P_phi/beta0/phi*sin(phis);
            COMPLEX ps = (-I)* B_theta/phi0/beta*sin(phis);
            COMPLEX ss = (-I)* I*B_phi/beta0/beta*cos(phis);

            // Calculate spatial frequency...
            double fx,fy;
            Bragg_Frequency(fx,fy);

            // The fourier transform of the surface roughness...
            double ft=sqrt(psd->psd(fx,fy));

            COMPLEX other_factors=sqr(2.*pi/lambda)/
                                    pi*sqrt(cos(thetai)*cos(thetas));

            return JonesMatrix(pp,ss,ps,sp)*(ft*other_factors);
        
        } else { //TRANSMISSION...
            error("No code for TRANSMISSION mode"); 
            return JonesMatrix(0.,0.,0.,0.);
        }
    }

    DEFINE_MODEL(Correlated_Roughness_BRDF_Model,Roughness_BRDF_Model,
                "Correlated_Roughness_BRDF_Model",
                "Scattering from a single dielectric layer on a substrate with equal and correlated roughness.");

    DEFINE_PARAMETER(Correlated_Roughness_BRDF_Model,double,thickness,"Film thickness [um]","0.05");

    DEFINE_PARAMETER(Correlated_Roughness_BRDF_Model,dielectric_function,film,"Overlayer","(1.46,0)");


} // namespace SCATMECH