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


namespace SCATMECH {

    // Parameters for area integration over a circle obtained from
    // Section 25.4.61 (Multidimensional Integration) of _Handbook_
    // of_Mathematical_Functions_, by M. Abramowitz and I.A. Stegun,
    // (Dover, New York, 1965).
    double Circle_Integral::xx[]=
    {0.,0.48206214668854613,0.18413135534528705,-0.18413135534528705,
       -0.48206214668854613,-0.5958615826865181,-0.48206214668854613,
       -0.18413135534528705,0.18413135534528705,0.48206214668854613,
       0.5958615826865181,0.7436573695963749,0.28405183920147253,
       -0.28405183920147253,-0.7436573695963749,-0.9192110607898046,
       -0.7436573695963749,-0.28405183920147253,0.28405183920147253,
       0.7436573695963749,0.9192110607898046};

    double Circle_Integral::yy[]=
    {0.,0.35023865071078736,0.5666980410239565,0.5666980410239565,
       0.35023865071078736,0.,-0.35023865071078736,-0.5666980410239565,
       -0.5666980410239565,-0.35023865071078736,0.,0.5402987052763671,
       0.8742216692147242,0.8742216692147242,0.5402987052763671,0.,
       -0.5402987052763671,-0.8742216692147242,-0.8742216692147242,
       -0.5402987052763671,0.};

    double Circle_Integral::ww[]=
    {0.1111111111111111,0.051248582618842166,0.051248582618842166,
       0.051248582618842166,0.051248582618842166,0.051248582618842166,
       0.051248582618842166,0.051248582618842166,0.051248582618842166,
       0.051248582618842166,0.051248582618842166,0.037640306270046726,
       0.037640306270046726,0.037640306270046726,0.037640306270046726,
       0.037640306270046726,0.037640306270046726,0.037640306270046726,
       0.037640306270046726,0.037640306270046726,0.037640306270046726};


    //
    // Constructor...
    //
    Focussed_Beam_Instrument_BRDF_Model::
    Focussed_Beam_Instrument_BRDF_Model()
    {
        model_cs = xyxy;
    }

    MuellerMatrix 
    Focussed_Beam_Instrument_BRDF_Model::
    mueller()
    {
        synchronize(model);
        SETUP();

        if (alpha==0.) return model->Mueller(thetai,thetas,phis,rotation,BRDF_Model::xyxy);

        MuellerMatrix m=MuellerZero();

        for (int i=0;i<Circle_Integral::n();++i) {
            double theta = Circle_Integral::theta(i);
            double r = Circle_Integral::r(i);
            double w = Circle_Integral::w(i);

            double ri = r*alpha*deg;
            double rs = ri*focal_point;

            Vector vi = four_angles(thetai,pi,ri,theta);
            Vector vs = four_angles(thetas,phis,rs,theta);
            Vector z(0,0,1);

            double thetas_ = acos(vs.z);
            double thetai_ = acos(vi.z);

            double weight = w*cos(thetai_)*cos(thetas_);
            m += model->Mueller(vi,vs,z,Vector(1,0,0),BRDF_Model::xyxy)*weight;
        }
        return (m/cos(thetai)/cos(thetas));
    }

    //
    // The following returns a unit vector, defined by the angles theta, phi, alpha, and beta:
    // Adds planar vector (alpha, beta) on unit sphere to vector (theta,phi).
    //
    Vector 
    Focussed_Beam_Instrument_BRDF_Model::
    four_angles(double theta,double phi,double alpha,double beta) 
    {
        return Vector(
                cos(beta)*sin(alpha)*sin(phi) + 
                    cos(phi)*(cos(theta)*sin(alpha)*sin(beta) + cos(alpha)*sin(theta)),
                -(cos(beta)*cos(phi)*sin(alpha)) + 
                    sin(phi)*(cos(theta)*sin(alpha)*sin(beta) + cos(alpha)*sin(theta)),
                cos(alpha)*cos(theta) - sin(alpha)*sin(beta)*sin(theta));
    }

    DEFINE_MODEL(Focussed_Beam_Instrument_BRDF_Model,Instrument_BRDF_Model, 
                "Focussed_Beam_Instrument_BRDF_Model",
                "A BRDF model measured with a focussed beam.");

    DEFINE_PTRPARAMETER(Focussed_Beam_Instrument_BRDF_Model,BRDF_Model_Ptr,model,"Model to be integrated","Microroughness_BRDF_Model");

    DEFINE_PARAMETER(Focussed_Beam_Instrument_BRDF_Model,double,alpha,"Incident beam half-angle [deg]","0");

    DEFINE_PARAMETER(Focussed_Beam_Instrument_BRDF_Model,double,focal_point,"Focal point to detector ratio","1");



}