//******************************************************************************
//** SCATMECH: Polarized Light Scattering C++ Class Library
//** 
//** File: mueller.h
//**
//** 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)
//**
//******************************************************************************
#ifndef SCATMECH_MUELLER_HPP 
#define SCATMECH_MUELLER_HPP

#include "scatmech.h"
#include "optconst.h"

//*****************************************************************************
//**
//** Some attempt was made to include most of the terminology found in 
//**    R. A. Chipman, "Polarimetry", Chapter 22 of _Handbook_of_Optics_,
//**    Volume II, (McGraw-Hill, New York, 1995)
//**
//*****************************************************************************
//
// Class declarations...
//

namespace SCATMECH {

    class JonesMatrix;
    class JonesVector;
    class MuellerMatrix;
    class StokesVector;

    /**
     * Class to handle Jones vectors.
     * Class to handle data storage and the various operations 
     *           associated with Jones vectors
     */
    class JonesVector {
        public:
    
            JonesVector() {}
            JonesVector(const JonesVector& x)
                          {j[0]=x.j[0]; j[1]=x.j[1];}
            /// Constructor taking two COMPLEX values
            JonesVector(const COMPLEX& s,const COMPLEX& p)
                      {j[0]=s; j[1]=p;}
            /// Conversion constructor from StokesVector
            explicit JonesVector(const StokesVector& x);

            JonesVector& operator=(const JonesVector& x)
                      {j[0]=x.j[0]; j[1]=x.j[1]; return *this;}
                
            JonesVector operator*(const COMPLEX& x) const
                                            {return JonesVector(j[0]*x,j[1]*x);}
            friend JonesVector operator*(const COMPLEX& x,const JonesVector& y) 
                                            {return y*x;}
            JonesVector operator/(const COMPLEX& x) const 
                                            {return (1./x)*(*this);}
            JonesVector& operator*=(const COMPLEX& x) 
                                            {return *this = x*(*this);}
            JonesVector& operator/=(const COMPLEX& x) 
                                            {return *this = (1./x)*(*this);}

            // Addition operators...
            JonesVector operator+(const JonesVector& a) const 
                            {return JonesVector(j[0]+a.j[0],j[1]+a.j[1]);}
            JonesVector operator-(const JonesVector& a) const 
                            {return JonesVector(j[0]-a.j[0],j[1]-a.j[1]);}
            JonesVector operator+=(const JonesVector& a) 
                            {return *this = *this + a;}
            JonesVector operator-=(const JonesVector& a) 
                            {return *this = *this - a;}
            JonesVector operator-() const {return JonesVector(-j[0],-j[1]);}
            JonesVector operator+() const {return *this;}

            /// Element indexing
            COMPLEX& operator[](int i) {return j[i];}
            const COMPLEX& operator[](int i) const {return j[i];}

            /// Return intensity
            double intensity() const {return std::norm(j[0])+std::norm(j[1]);}
        
            /// arctan of the component ratio
            double psi() const          
                {return atan(std::abs(j[1])/std::abs(j[0]));}    

            /// phase between two components
            double delta() const 
                {return std::arg(j[1])-std::arg(j[0]);}   

            /// principle angle of polarization
            double eta() const;
            
            /// degree of linear polarization
            double DOLP() const;

            /// degree of polarization (always 1 for Jones!)
            double DOP() const      
                        {return 1.;} 

            /// degree of circular polarization
            double DOCP() const;    

            /// ellipticity (ratio of minor to major axes)
            double e() const;       

            /// eccentricity [sqrt(1-e^2)]
            double epsilon() const  
                {double t = e(); return sqrt(1.-t*t);} 

            /// Vector rotated by angle...
            JonesVector rotate(const double angle) const;
        
            /// S component of wave
            COMPLEX& S() {return j[0];}
            const COMPLEX& S() const {return j[0];}

            /// P component of wave
            COMPLEX& P() {return j[1];}
            const COMPLEX& P() const {return j[1];}

            friend class JonesMatrix;
            friend class MuellerMatrix;
            friend class StokesVector; 

        private:
    
            // The elements of the vector...
            COMPLEX j[2];
    };

    /**
     * Class to handle Jones Matrices...
     * Class to store the elements and handle operations of a Jones matrix
     */
    class JonesMatrix {
        public:
            JonesMatrix() {}
            JonesMatrix(const JonesMatrix& x)
                { j[0]=x.j[0]; j[1]=x.j[1]; j[2]=x.j[2]; j[3]=x.j[3];}

            JonesMatrix(const COMPLEX& pp,const COMPLEX& ss,
                        const COMPLEX& ps,const COMPLEX& sp)
                        { j[0]=pp; j[1]=ss; j[2]=ps, j[3]=sp;}
            explicit JonesMatrix(const MuellerMatrix& x);

            JonesMatrix& operator=(const JonesMatrix& x) 
                { j[0]=x.j[0]; j[1]=x.j[1]; j[2]=x.j[2]; j[3]=x.j[3]; return *this;}
        
            // Multiplication operators...
            JonesMatrix operator*(const JonesMatrix& matrix) const
                { return JonesMatrix(j[3]*matrix.j[2]+j[0]*matrix.j[0],
                                     j[1]*matrix.j[1]+j[2]*matrix.j[3],
                                     j[1]*matrix.j[2]+j[2]*matrix.j[0],
                                     j[3]*matrix.j[1]+j[0]*matrix.j[3]); }

            JonesMatrix operator*(const COMPLEX& x) const
                {return JonesMatrix(j[0]*x,j[1]*x,j[2]*x,j[3]*x);}

            friend JonesMatrix operator*(const COMPLEX& x,const JonesMatrix& y) 
                                        {return y*x;}
            JonesMatrix operator/(const COMPLEX& x) const 
                        {return JonesMatrix(j[0]/x,j[1]/x,j[2]/x,j[3]/x);}
            JonesMatrix& operator*=(const JonesMatrix& a) 
                        {return *this = *this * a;}
            JonesMatrix& operator*=(const COMPLEX& a) {return *this = *this * a;}
            JonesMatrix& operator/=(const COMPLEX& a) {return *this = *this / a;}
            JonesVector operator*(const JonesVector& a) const 
                { return JonesVector(j[1]*a.j[0]+j[2]*a.j[1],j[3]*a.j[0]+j[0]*a.j[1]); }
        
            // Addition operators...
            JonesMatrix operator+(const JonesMatrix& a) const
                { return JonesMatrix(j[0]+a.j[0],j[1]+a.j[1],j[2]+a.j[2],j[3]+a.j[3]);}

            JonesMatrix operator-(const JonesMatrix& a) const
                { return JonesMatrix(j[0]-a.j[0],j[1]-a.j[1],j[2]-a.j[2],j[3]-a.j[3]);}

            JonesMatrix& operator+=(const JonesMatrix& a) 
                                        {return *this = *this + a;}
            JonesMatrix& operator-=(const JonesMatrix& a) 
                                        {return *this = *this - a;}
            JonesMatrix operator-() const {return JonesMatrix(-j[0],-j[1],-j[2],-j[3]);}

            JonesMatrix operator+() const {return *this;}
        
            /// Element indexing
            COMPLEX& operator[](int i) {return j[i];};

            /// Constant element indexing
            const COMPLEX& operator[](int i) const {return j[i];};

            /// Return matrix rotated by angle
            JonesMatrix rotate(double angle) const;
        
            /// Return matrix transpose
            JonesMatrix transpose() const {return JonesMatrix(j[0],j[1],j[3],j[2]);}
        
            /// p->p matrix element
            COMPLEX& PP() {return j[0];}
            const COMPLEX& PP() const {return j[0];}

            /// s->s matrix element
            COMPLEX& SS() {return j[1];}
            const COMPLEX& SS() const {return j[1];}
        
            /// p->s matrix element
            COMPLEX& PS() {return j[2];}
            const COMPLEX& PS() const {return j[2];}

            /// s->p matrix element
            COMPLEX& SP() {return j[3];}
            const COMPLEX& SP() const {return j[3];}

            friend class JonesVector;
            friend class MuellerMatrix;
            friend class StokesVector;
        private:
    
            // Array of matrix elements...
            // (Note that the elements are p->p, s->s, p->s, and s->p, in order)
            COMPLEX j[4]; 
    };

    /**
     * Class to handle Mueller matrices
     * Class to handle the storage of the elements and operations of a Mueller matrix
     */
    class MuellerMatrix {
    
        public:
    
            // Constructors and Assignments...
            MuellerMatrix() {}
            MuellerMatrix(const JonesMatrix& j);
            MuellerMatrix(const MuellerMatrix& x) {*this=x;}
            MuellerMatrix(const StokesVector& s1,const StokesVector& s2,const StokesVector& s3, const StokesVector& s4);
            MuellerMatrix& operator=(const MuellerMatrix& x);
    
            // Element indexing...
            double* operator[](int i) {return m[i];}
            const double* operator[](int i) const {return m[i];}
        
            /// Maximum "transmission"...
            double Tmax() const 
            {
             return m[0][0]+sqrt(m[0][1]*m[0][1]+m[0][2]*m[0][2]+m[0][3]*m[0][3]);
            }
            /// Minimum "transmission"...
            double Tmin() const 
            {
             return m[0][0]-sqrt(m[0][1]*m[0][1]+m[0][2]*m[0][2]+m[0][3]*m[0][3]);
            }
            /// Diattenuation...
            double diattenuation() const 
            {
             return sqrt(m[0][1]*m[0][1]+m[0][2]*m[0][2]+m[0][3]*m[0][3])/m[0][0];
            }
            /// Linear diattenuation...
            double linear_diattenuation() const 
            {
             return sqrt(m[0][1]*m[0][1]+m[0][2]*m[0][2])/m[0][0];
            }
            /// Polarization dependent loss...
            double polarization_dependent_loss() const 
            {
             return 10.*log(Tmax()/Tmin())/log(10.);
            }
            /// Polarizance...
            double polarizance() const 
            {
             return m[0][0]+sqrt(m[1][0]*m[1][0]+m[2][0]*m[2][0]+m[3][0]*m[3][0]);
            }
            /// Extinction ratio...
            double extinction_ratio() const {return Tmax()/Tmin();}
        
            /// Return matrix rotated by angle...
            MuellerMatrix rotate(const double angle) const;
            /// Return scattering matrix for a parity conversion
            MuellerMatrix parity() const;
            /// Return transpose of matrix...
            MuellerMatrix transpose() const;
        
            // Binary Arithmetic Operations...
            MuellerMatrix operator*(const MuellerMatrix& matrix) const;
            MuellerMatrix& operator*=(const MuellerMatrix& matrix) 
                            {return *this = *this * matrix;}
        
            MuellerMatrix operator*(double d) const;
            friend MuellerMatrix operator*(double d,const MuellerMatrix &v) 
                            {return v*d;}
            StokesVector operator*(const StokesVector &s) const;
            MuellerMatrix& operator*=(const double d) {return *this = *this * d;}
            MuellerMatrix operator/(double d) const;
            MuellerMatrix& operator/=(double d) {return *this = *this * (1./d);}
            MuellerMatrix operator+(const MuellerMatrix& a) const;
            MuellerMatrix operator-(const MuellerMatrix& a) const;
            MuellerMatrix operator-() const;
            MuellerMatrix operator+() const {return *this;}
            MuellerMatrix& operator+=(const MuellerMatrix& matrix) 
                                    {return *this = *this + matrix;}
            MuellerMatrix& operator-=(const MuellerMatrix& matrix) 
                                    {return *this = *this - matrix;}
        
            friend class JonesVector;
            friend class JonesMatrix;
            friend class StokesVector;
            
        private:
    
            // The elements of the Mueller matrix...
            double m[4][4];
    };

    /**
     * Class to handle Stokes vectors
     * Class to store the elements and handle the operations of a Stokes vector
     */
    class StokesVector {
  
        public:
        
            StokesVector() {}
            StokesVector(double I,double Q,double U,double V) 
                            {s[0]=I;s[1]=Q; s[2]=U; s[3]=V;}
            StokesVector(const StokesVector& x)
                {s[0]=x.s[0]; s[1]=x.s[1]; s[2]=x.s[2]; s[3]=x.s[3];}
            StokesVector(const JonesVector& j)
                {
                    using namespace std;
                    s[0] = norm(j.j[0])+norm(j.j[1]);
                    s[1] = norm(j.j[0])-norm(j.j[1]);
                    s[2] = 2*real(conj(j.j[0])*j.j[1]);  // (s+p) minus (s-p)
                    s[3] = 2*imag(conj(j.j[0])*j.j[1]);  // Left minus Right
                }  

            StokesVector& operator=(const StokesVector& x)
                {s[0]=x.s[0]; s[1]=x.s[1]; s[2]=x.s[2]; s[3]=x.s[3]; return *this;}
        
            // Aliases for each element...

            /// The first elements (intensity) of the Stokes vector
            double& I() {return s[0];}
            double I() const {return s[0];}

            /// The second element (Is - Ip) of the Stokes vector
            double& Q() {return s[1];}
            double Q() const {return s[1];}

            /// The third element [I(45) - I(-45)] of the Stokes vector
            double& U() {return s[2];}
            double U() const {return s[2];}

            /// The fourth element [I(lcp) - I(rcp)] of the Stokes vector
            double& V() {return s[3];}
            double V() const {return s[3];}
        
            /// The i-th element (zero indexing)
            double& operator[](int i) {return s[i];}
            double operator[](int i) const {return s[i];}

            // Binary arithmetic operations...
            StokesVector operator*(const MuellerMatrix& matrix) const;
            double operator*(const StokesVector& a) const
                {return a.s[0]*s[0] + a.s[1]*s[1] + a.s[2]*s[2] + a.s[3]*s[3];}

            StokesVector operator+(const StokesVector& a) const
                {return StokesVector(s[0]+a.s[0],s[1]+a.s[1],s[2]+a.s[2],s[3]+a.s[3]);}
            StokesVector operator-(const StokesVector& a) const 
                {return StokesVector(s[0]-a.s[0],s[1]-a.s[1],s[2]-a.s[2],s[3]-a.s[3]);}
            StokesVector operator-() const 
                {return StokesVector(-s[0],-s[1],-s[2],-s[3]);}
            StokesVector operator+() const {return *this;}

            StokesVector operator*(double d) const 
                {return StokesVector(s[0]*d,s[1]*d,s[2]*d,s[3]*d);}
            friend StokesVector operator*(double d,const StokesVector& s) 
                {return StokesVector(s[0]*d,s[1]*d,s[2]*d,s[3]*d);}
            StokesVector& operator*=(double d) {return *this = *this *  d;}
            StokesVector operator/(double d) const 
                {return StokesVector(s[0]/d,s[1]/d,s[2]/d,s[3]/d);}
            StokesVector& operator/=(double d) {return *this = *this /  d;}
        
            /// Return vector rotated by angle...
            StokesVector rotate(double angle) const;
        
            /// Principle angle of polarization
            double eta() const         
                {return atan2(s[2],s[1])/2.; }

            /// Intensity
            double intensity() const    
                {return s[0];}
            
            /// Degree of linear polarization
            double DOLP() const     
                {return sqrt(sqr(s[2])+sqr(s[1]))/s[0];}
        
            /// Degree of polarization 
            double DOP() const      
                {return sqrt(sqr(s[2])+sqr(s[1])+sqr(s[3]))/s[0];}

            /// Degree of circular polarization
            double DOCP() const     
                {return s[3]/s[0];}

            /// Ellipticity (ratio of minor to major axes)
            double e() const        
                {return s[3]/(s[0]+sqrt(sqr(s[1])+sqr(s[2])));}

            /// Phase between two components
            double delta() const;    

            /// Arctangent of the component ratio
            double psi() const;
            
            /// eccentricity [sqrt(1-e^2)]
            double eccentricity() const 
                {return sqrt(1.-sqr(this->e()));}  
            
            /// Polarized part of Stokes Vector
            StokesVector pol_part() const;  

            /// Unpolarized part of Stokes Vector
            StokesVector unpol_part() const; 
        
            friend class JonesVector;
            friend class MuellerMatrix;
            friend class JonesMatrix;

        private:
    
            // The elements of the Stokes vector...
            double s[4];
    };

    std::ostream& operator<<(std::ostream& os,const JonesVector& j);
    std::istream& operator>>(std::istream& is,JonesVector& j);
    std::ostream& operator<<(std::ostream& os,const StokesVector& j);
    std::istream& operator>>(std::istream& is, StokesVector& j);
    std::ostream& operator<<(std::ostream& os,const MuellerMatrix& mm);
    std::istream& operator>>(std::istream& is, MuellerMatrix& mm);
    std::ostream& operator<<(std::ostream& os,const JonesMatrix& j);
    std::istream& operator>>(std::istream& is, JonesMatrix& j);

    // Some common elements...
    inline JonesMatrix JonesZero() {return JonesMatrix(0.,0.,0.,0.);}

    JonesMatrix JonesRotator(double angle);
    JonesMatrix JonesLinearRetarder(double phase, double angle=0);
    JonesMatrix JonesCircularRetarder(double phase);
    JonesMatrix JonesLinearPolarizer(double angle=0,double diattenuation=1);
    JonesMatrix JonesCircularPolarizer(double diattenuation=1);                               

    // Returns Jones matrix given eigenvectors and eigenvalues...
    JonesMatrix JonesGeneralized(const JonesVector& a,const JonesVector& b, 
                        const COMPLEX& ma, const COMPLEX& mb);
 
    MuellerMatrix MuellerZero();
    MuellerMatrix MuellerUnit(double attenuation=1.);
    MuellerMatrix MuellerDepolarizer(double attenuation=1.,double depolarization=1.); 
    MuellerMatrix MuellerPartialLinearPolarizer(double tmax,double tmin,double angle);

    inline StokesVector StokesVectorUnitUnpol() {return StokesVector(1,0,0,0);}
    inline StokesVector StokesVectorUnitS() {return StokesVector(1,1,0,0);}
    inline StokesVector StokesVectorUnitP() {return StokesVector(1,-1,0,0);}
    inline StokesVector StokesVectorUnitLinear(double eta) {return StokesVector(1,cos(2*eta),sin(2*eta),0);}
    inline StokesVector StokesVectorUnitRCP() {return StokesVector(1,0,0,-1);}
    inline StokesVector StokesVectorUnitLCP() {return StokesVector(1,0,0,1);}
    inline StokesVector StokesVectorUnitGeneral(double eta,double DOCP=0.,double DOP=1.)
    {
        double temp=sqrt(1.-sqr(DOCP));
        return StokesVector(1,
                            DOP*temp*cos(2*eta),
                            DOP*temp*sin(2*eta),
                            DOP*DOCP);
    }

    inline JonesVector JonesVectorUnitS() {return JonesVector(1.,0.);}
    inline JonesVector JonesVectorUnitP() {return JonesVector(0.,1.);}
    inline JonesVector JonesVectorUnitLinear(double eta) {return JonesVector(cos(eta),sin(eta));}
    inline JonesVector JonesVectorUnitRCP() {return JonesVector(1.,COMPLEX(0,-1));}
    inline JonesVector JonesVectorUnitLCP() {return JonesVector(1.,COMPLEX(0,1));}
    inline JonesVector JonesVectorUnitGeneral(double eta,double DOCP) 
    {
        return JonesVector(StokesVectorUnitGeneral(eta,DOCP));
    }

} // namespace SCATMECH 

#endif