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

using namespace std;


namespace SCATMECH {


    //
    // Jones Matrix for scattering...
    //
    JonesMatrix 
    Microroughness_BRDF_Model::
    jones()
    {
        // All routines should call setup()...
        SETUP();

        if (type==REFLECTION) {
            // This routine follows the Rayleigh-Rice first-order vector 
            // perturbation theory. For details, see Thomas A. Germer, 
            // "Angular dependence and polarization of out-of-plane optical 
            // scattering from particulate contamination, subsurface defects, 
            // and surface microroughness," 
            // Applied Optics vol. 36 no. 33, pp 8798-8805 (1997).
        
            // An abreviation...
            COMPLEX e = (COMPLEX)(substrate.epsilon(lambda));
        
            // p-->p scattering
            COMPLEX pp = (e-1.)*
                        (
                        e*sin(thetai)*sin(thetas)-
                        sqrt(e-sqr(sin(thetas)))*sqrt(e-sqr(sin(thetai)))*cos(phis)
                        )/
                        (e*cos(thetai)+sqrt(e-sqr(sin(thetai))))/
                        (e*cos(thetas)+sqrt(e-sqr(sin(thetas))));
        
            // s-->s scattering
            COMPLEX ss = (e-1.)*cos(phis)/
                        (cos(thetai)+sqrt(e-sqr(sin(thetai))))/
                        (cos(thetas)+sqrt(e-sqr(sin(thetas))));
        
            // p-->s scattering
            COMPLEX ps = -(e-1.)*sqrt(e-sqr(sin(thetai)))*sin(phis)/
                        (e*cos(thetai)+sqrt(e-sqr(sin(thetai))))/
                        (cos(thetas)+sqrt(e-sqr(sin(thetas))));  
            // s-->p scattering
            COMPLEX sp = -(e-1.)*sqrt(e-sqr(sin(thetas)))*sin(phis)/
                        (cos(thetai)+sqrt(e-sqr(sin(thetai))))/
                        (e*cos(thetas)+sqrt(e-sqr(sin(thetas))));
        
            // Calculate spatial frequency...
            double fx,fy;
            Bragg_Frequency(fx,fy);
        
            // Various factors common to all the scattering elements...
            COMPLEX factor=sqrt(
                                16.*sqr(pi)/sqr(sqr(lambda))*
                                cos(thetai)*cos(thetas)*
                                psd->psd(fx,fy)
                                );
        
            // Return the Jones matrix...
            return JonesMatrix(pp,ss,ps,sp)*factor;
        
        } else /* TRANSMISSION */ { 

            // The solution for transmissive scatter was obtained from 
            //      J. M. Elson, "Multilayer-coated optics: guided-wave coupling 
            //                    and scattering by means of interface random 
            //                    roughness," J. Opt. Soc. Am. A 12(4) 729 (1995).
            //
            // The case for L=0 simplifies considerably to the following...
            //
            COMPLEX q0[2];
            COMPLEX q1[2];
        
            COMPLEX eps=(COMPLEX)(substrate.epsilon(lambda));    
            COMPLEX index = sqrt(eps);
            COMPLEX k = sin(thetas)*index;
            double k0 = sin(thetai);
 
            q0[0]=sqrt(eps-k0*k0);
            q1[0]=sqrt(eps-k*k);
            q0[1]=sqrt(1.-k0*k0);
            q1[1]=sqrt(1.-k*k);
        
            double fx,fy;
            Bragg_Frequency(fx,fy);

            double PSD = psd->psd(fx,fy);

            JonesMatrix J;
            J.PP() =  -index*(k*k0+cos(phis)*q0[0]*q1[1])/(q0[0]+eps*q0[1])/(q1[0]+eps*q1[1]);
            J.SS() =  cos(phis)/(q0[0]+q0[1])/(q1[0]+q1[1]);
            J.SP() =  -index*q1[1]*sin(phis)/(q0[0]+q0[1])/(q1[0]+eps*q1[1]);
            J.PS() =  q0[0]*sin(phis)/(q0[0]+eps*q0[1])/(q1[0]+q1[1]);

            COMPLEX common = 4.*pi*(1.-eps)*index*sqrt(q1[0]*q0[1]*PSD)/sqr(lambda);
            return J*common;
        }
    }

    DEFINE_MODEL(Microroughness_BRDF_Model,  
                Roughness_BRDF_Model,       
                "Microroughness_BRDF_Model",
                "Scattering by a single rough surface in the smooth-surface limit (no films).");


} // namespace SCATMECH