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

using namespace std;


namespace SCATMECH {

    //
    // Routine to get the isotropic spatial frequency...
    //
    void
    Roughness_BRDF_Model::
    Bragg_Frequency(double& fx,double& fy)
    {
        double _fx,_fy;
        if (type==REFLECTION) {
            // Calculate magnitude of spatial frequency...
            _fx = sin(thetas)*cos(phis)-sin(thetai);
            _fy = sin(thetas)*sin(phis);
        } else { // type==TRANSMISSION
            double m = substrate.n(lambda);
            _fx = sin(thetas)*cos(phis)*m-sin(thetai);
            _fy = sin(thetas)*sin(phis)*m;
        }
        _fx = _fx/lambda;
        _fy = _fy/lambda;

        double cosrot = cos(rotation);
        double sinrot = sin(rotation);
    
        fx =  cosrot*_fx + sinrot*_fy;
        fy = -sinrot*_fx + cosrot*_fy;
    }

    int
    Roughness_BRDF_Model::
    get_recalc()
    {
        return BRDF_Model::get_recalc()||psd->get_recalc();
    }

    double 
    ABC_PSD_Function::
    psd(double f) 
    {
        SETUP();
        return A/pow(1+sqr(B*f),C/2);
    }

    double 
    Table_PSD_Function::
    psd(double f) 
    {
        SETUP();
        return T.value(f);
    }

    double 
    Fractal_PSD_Function::
    psd(double f) 
    {
        SETUP();
        return A/pow(f,exponent);
    }

    void
    Gaussian_PSD_Function::
    setup()
    {
        PSD_Function::setup();
        A = pi*sqr(sigma*length);
        B = pi*length;
    }

    double 
    Gaussian_PSD_Function::
    psd(double f) 
    {
        SETUP();
        return A*exp(-sqr(B*f));
    }

    static double poly(double x,double* coeff,int n) {
        double accum=0.;
        for (int i=0;i<n;++i) accum = coeff[i]+x*accum;
        return accum;
    }

    double j1pos(double x)
    {
        static double coeff1[] = {1.,376.9991397,99447.43394,18583304.74,2300535178.0,144725228443.0};
        static double coeff2[] = {-30.16036606,15704.48260,-2972611.439,242396853.1,-7895059235.0,72362614232.0};
        static double coeff3[] = {-0.240337019e-6,0.2457520174e-5,-0.3516396496e-4,0.183105e-2,1.0};
        static double coeff4[] = {0.105787412e-6,-0.88228987e-6,0.8449199096e-5,-0.2002690873e-3,0.04687499995};
    
        if (fabs(x) < 8.) {
            return x*poly(x*x,coeff2,6)/poly(x*x,coeff1,6);
        } else {
            return sqrt(0.636619772/x)*(cos(x-2.356194491)*poly(sqr(8./x),coeff3,5)-
                                    8./x*sin(x-2.356194491)*poly(sqr(8./x),coeff4,5));
        }
    }

    double j1(double x) {
        if (x>0) return j1pos(x);
        else return -j1pos(-x);
    }

    double
    Elliptical_Mesa_PSD_Function::
    psd(double fx,double fy)
    {
        SETUP();
        double f = sqrt(sqr(fx*axisx/2.)+sqr(fy*axisy/2.));
        double rx_ry = axisx*axisy/4.;
        //double temp =  (f!=0) ? _j1(2*pi*f)/f : pi;
        double temp = (f!=0) ? j1(2*pi*f)/f : pi;
        return sqr(rx_ry*temp*height)*density;
    }

    static double sinc(double x)
    {
        if (x==0.) return 1;
        else return sin(x)/x;
    }

    double
    Rectangular_Mesa_PSD_Function::
    psd(double fx,double fy)
    {
        SETUP();
        return sqr(lengthx*lengthy*sinc(fx*lengthx*pi)*sinc(fy*lengthy*pi)*height)*density;
    }


    double
    Triangular_Mesa_PSD_Function::
    psd(double fx,double fy)
    {
        SETUP();

        static const double sqrt3=sqrt(3.);
        static const COMPLEX i(0,1);
        COMPLEX s;
        double eps = 1./side*0.000001;
        if (fabs(fy)<eps) {
            s = (sqrt3 - 3.*i*fx*pi*side - sqrt3*cos(sqrt3*fx*pi*side) + 
                i*sqrt3*sin(sqrt3*fx*pi*side))/(6.*sqr(fx*pi));
            if (fabs(fx) < eps) {
                s = sqrt3*sqr(side)/4.;
            }
        } else if (fabs(sqrt3*fx - fy)<eps) {
            (-i/4.*sqrt3*(exp(i*fy*pi*side)*fy*pi*side - sin(fy*pi*side)))/sqr(fy*pi);
        } else {
            s = (sqrt3*fy*(cos(sqrt3*fx*side*pi)-cos(fy*side*pi)-i*sin(sqrt3*fx*side*pi)) +
            3.*i*fx*sin(fy*side*pi))/(2.*(-3.*sqr(fx)*fy+cube(fy))*sqr(pi));
        }
        return density * sqr(height) * norm(s);
    }

    double
    Rectangular_Pyramid_PSD_Function::
    psd(double fx,double fy)
    {
        SETUP();
        double s = (fy*height*lengthy*cos(fy*lengthy*pi)*sin(fx*lengthx*pi) - 
                    fx*height*lengthx*cos(fx*lengthx*pi)*sin(fy*lengthy*pi))/
                    (cube(fx*pi)*fy*sqr(lengthx) - fx*cube(fy*pi)*sqr(lengthy));

        if (fx==0) {
            if (fy==0) {
                s = (height*lengthx*lengthy)/3.;
            } else {
                s = (height*lengthx*(-(fy*lengthy*pi*cos(fy*lengthy*pi)) + sin(fy*lengthy*pi)))/(cube(fy*pi)*sqr(lengthy));
            }
        } else if (fy==0) {
            s = (height*lengthy*(-(fx*lengthx*pi*cos(fx*lengthx*pi)) + sin(fx*lengthx*pi)))/(cube(fx*pi)*sqr(lengthx));
        } else {
            s = (fy*height*lengthy*cos(fy*lengthy*pi)*sin(fx*lengthx*pi) - 
                    fx*height*lengthx*cos(fx*lengthx*pi)*sin(fy*lengthy*pi))/
                    (cube(fx*pi)*fy*sqr(lengthx) - fx*cube(fy*pi)*sqr(lengthy));
        }
        return density * sqr(s);
    }


    double
    Double_PSD_Function::
    psd(double fx,double fy)
    {
        SETUP();
        return psd1->psd(fx,fy) + psd2->psd(fx,fy);
    }

    void Register(const PSD_Function* x)
    {
        static bool Models_Registered = false;
        if (!Models_Registered) {
            Models_Registered=true;

            Register_Model(PSD_Function);
            Register_Model(Unit_PSD_Function);
            Register_Model(ABC_PSD_Function);
            Register_Model(Table_PSD_Function);
            Register_Model(Fractal_PSD_Function);
            Register_Model(Gaussian_PSD_Function);
            Register_Model(Elliptical_Mesa_PSD_Function);
            Register_Model(Rectangular_Mesa_PSD_Function);
            Register_Model(Triangular_Mesa_PSD_Function);
            Register_Model(Rectangular_Pyramid_PSD_Function);
            Register_Model(Double_PSD_Function);
        }
    }

    DEFINE_VIRTUAL_MODEL(Roughness_BRDF_Model,BRDF_Model,
                        "Roughness_BRDF_Model","Scattering by a rough surface in the smooth surface limit.");
    DEFINE_PTRPARAMETER(Roughness_BRDF_Model,PSD_Function_Ptr,psd,"Power spectral density function","ABC_PSD_Function");


    DEFINE_VIRTUAL_MODEL(PSD_Function,Model,
                        "PSD_Function","Virtual class for Power Spectral Density functions");

    DEFINE_MODEL(Unit_PSD_Function,PSD_Function,
                "Unit_PSD_Function","Unit power spectrum");

    DEFINE_MODEL(ABC_PSD_Function,PSD_Function,  
                "ABC_PSD_Function","PSD = A/pow(1+sqr(B*f),C/2)");
    DEFINE_PARAMETER(ABC_PSD_Function,double,A,"Program assumes power spectrum S2(f) = "
                                                "A/(1+(B*f)^2)^(C/2)\n"
                                                "Power Spectrum A Parameter [um^4]","0.01");
    DEFINE_PARAMETER(ABC_PSD_Function,double,B,"Power Spectrum B Parameter [um]","362");
    DEFINE_PARAMETER(ABC_PSD_Function,double,C,"Power Spectrum C Parameter","2.5");

    DEFINE_MODEL(Table_PSD_Function,PSD_Function,
                "Table_PSD_Function","PSD given by a table of values.");
    DEFINE_PARAMETER(Table_PSD_Function,Table,T,"Table of PSD versus spatial frequency","1");

    DEFINE_MODEL(Fractal_PSD_Function,
                PSD_Function,
                "Fractal_PSD_Function",
                "A PSD with a fractal behavior");
    DEFINE_PARAMETER(Fractal_PSD_Function,double,A,"Power spectrum amplitude parameter","0.01");
    DEFINE_PARAMETER(Fractal_PSD_Function,double,exponent,"Power spectrum exponent","2.5");

    DEFINE_MODEL(Gaussian_PSD_Function,
                PSD_Function,
                "Gaussian_PSD_Function",
                "A PSD with gaussian statistics");
    DEFINE_PARAMETER(Gaussian_PSD_Function,double,sigma,"Standard deviation of height [um]","0.05");
    DEFINE_PARAMETER(Gaussian_PSD_Function,double,length,"Correlation length [um]","5");

    DEFINE_MODEL(Elliptical_Mesa_PSD_Function,
                 PSD_Function,
                "Elliptical_Mesa_PSD_Function",
                "PSD for elliptical mesas or depressions");
    DEFINE_PARAMETER(Elliptical_Mesa_PSD_Function,double,axisx,"Axis of mesas along x-direction [um]","1.5");
    DEFINE_PARAMETER(Elliptical_Mesa_PSD_Function,double,axisy,"Axis of mesas along y-direction [um]","1.5");
    DEFINE_PARAMETER(Elliptical_Mesa_PSD_Function,double,height,"Height or depth  of mesas [um]","0.01");
    DEFINE_PARAMETER(Elliptical_Mesa_PSD_Function,double,density,"Surface number density of mesas [um^-2]","0.01");

    DEFINE_MODEL(Rectangular_Mesa_PSD_Function,
                 PSD_Function,
                "Rectangular_Mesa_PSD_Function",
                "PSD for rectangular mesas or depressions");
    DEFINE_PARAMETER(Rectangular_Mesa_PSD_Function,double,lengthx,"Length of mesas along x-direction [um]","1.5");
    DEFINE_PARAMETER(Rectangular_Mesa_PSD_Function,double,lengthy,"Length of mesas along y-direction [um]","1.5");
    DEFINE_PARAMETER(Rectangular_Mesa_PSD_Function,double,height,"Height or depth  of mesas [um]","0.01");
    DEFINE_PARAMETER(Rectangular_Mesa_PSD_Function,double,density,"Surface number density of mesas [um^-2]","0.01");

    DEFINE_MODEL(Triangular_Mesa_PSD_Function,
                 PSD_Function,
                "Triangular_Mesa_PSD_Function",
                "PSD for isosceles mesas or depressions");
    DEFINE_PARAMETER(Triangular_Mesa_PSD_Function,double,side,"Length of side of triangle [um]","1.5");
    DEFINE_PARAMETER(Triangular_Mesa_PSD_Function,double,height,"Height or depth of mesas [um]","0.01");
    DEFINE_PARAMETER(Triangular_Mesa_PSD_Function,double,density,"Surface number density of mesas [um^-2]","0.01");

    DEFINE_MODEL(Rectangular_Pyramid_PSD_Function,
                 PSD_Function,
                "Rectangular_Pyramid_PSD_Function",
                "PSD for rectangular pyramids.");
    DEFINE_PARAMETER(Rectangular_Pyramid_PSD_Function,double,lengthx,"Length of side of pyramid in x-direction [um]","1.5");
    DEFINE_PARAMETER(Rectangular_Pyramid_PSD_Function,double,lengthy,"Length of side of pyramid in y-direction [um]","1.5");
    DEFINE_PARAMETER(Rectangular_Pyramid_PSD_Function,double,height,"Height or depth of pyramid [um]","0.01");
    DEFINE_PARAMETER(Rectangular_Pyramid_PSD_Function,double,density,"Surface number density of pyramids [um^-2]","0.01");

    DEFINE_MODEL(Double_PSD_Function,
                 PSD_Function,
                "Double_PSD_Function",
                "Sum of two PSD functions");
    DEFINE_PTRPARAMETER(Double_PSD_Function,PSD_Function_Ptr,psd1,"First PSD Function","ABC_PSD_Function");
    DEFINE_PTRPARAMETER(Double_PSD_Function,PSD_Function_Ptr,psd2,"Second PSD Function","ABC_PSD_Function");


} // namespace SCATMECH