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

using namespace std;


namespace SCATMECH {

 
    //
    // Conversion from JonesMatrix to MuellerMatrix
    //
    MuellerMatrix::
    MuellerMatrix(const JonesMatrix& jj)
    {
        //
        // from
        // Bohren and Huffman
        // _Absorption_and_Scattering_of_Light_by_Small_Particles_,
        // (Wiley, New York, 1983). 
        //
        //init();
        double norm0 = norm(jj.j[0]);
        double norm1 = norm(jj.j[1]);
        double norm2 = norm(jj.j[2]);
        double norm3 = norm(jj.j[3]);

        m[0][0] = 0.5*(norm0+norm1+norm2+norm3);
        m[0][1] = 0.5*(norm1-norm0+norm3-norm2);
        m[0][2] = real(jj.j[1]*conj(jj.j[2])+jj.j[0]*conj(jj.j[3]));
        m[0][3] = imag(jj.j[1]*conj(jj.j[2])-jj.j[0]*conj(jj.j[3]));
        m[1][0] = 0.5*(norm1-norm0-norm3+norm2);
        m[1][1] = 0.5*(norm1+norm0-norm3-norm2);
        m[1][2] = real(jj.j[1]*conj(jj.j[2])-jj.j[0]*conj(jj.j[3]));
        m[1][3] = imag(jj.j[1]*conj(jj.j[2])+jj.j[0]*conj(jj.j[3]));
        m[2][0] = real(jj.j[1]*conj(jj.j[3])+jj.j[0]*conj(jj.j[2]));
        m[2][1] = real(jj.j[1]*conj(jj.j[3])-jj.j[0]*conj(jj.j[2]));
        m[2][2] = real(jj.j[0]*conj(jj.j[1])+jj.j[2]*conj(jj.j[3]));
        m[2][3] = imag(jj.j[1]*conj(jj.j[0])+jj.j[3]*conj(jj.j[2]));
        m[3][0] = imag(jj.j[3]*conj(jj.j[1])+jj.j[0]*conj(jj.j[2]));
        m[3][1] = imag(jj.j[3]*conj(jj.j[1])-jj.j[0]*conj(jj.j[2]));
        m[3][2] = imag(jj.j[0]*conj(jj.j[1])-jj.j[2]*conj(jj.j[3]));
        m[3][3] = real(jj.j[0]*conj(jj.j[1])-jj.j[2]*conj(jj.j[3]));
    }

    MuellerMatrix::
    MuellerMatrix(const StokesVector& s1,const StokesVector& s2,const StokesVector& s3, const StokesVector& s4)
    {
        for (int i=0;i<4;++i) {
            m[0][i]=s1[i];
            m[1][i]=s2[i];
            m[2][i]=s3[i];
            m[3][i]=s4[i];
        }
    }

    //
    // Routine to change Mueller matrix to one for a parity conversion...
    //
    MuellerMatrix 
    MuellerMatrix::
    parity() const
    {
        MuellerMatrix mm = *this;
        mm[0][2] = -mm[0][2];
        mm[0][3] = -mm[0][3];
        mm[1][2] = -mm[1][2];
        mm[1][3] = -mm[1][3];
        mm[2][0] = -mm[2][0];
        mm[2][1] = -mm[2][1];
        mm[3][0] = -mm[3][0];
        mm[3][1] = -mm[3][1];
        return mm;
    }

    MuellerMatrix& 
    MuellerMatrix::
    operator=(const MuellerMatrix& x) 
    {
        for (int i=0;i<4;++i) {
            for (int j=0;j<4;++j) {
                m[i][j]=x.m[i][j];
            }
        }
        return *this;
    }

    //
    // Define multiplication between two Mueller matrices
    //
    MuellerMatrix 
    MuellerMatrix::
    operator*(const MuellerMatrix& matrix) const 
    {
        MuellerMatrix out;
        int i,j,k;
        for (i=0;i<4;++i) {
            for (j=0;j<4;++j) { 
                out.m[i][j]=0;
                for (k=0;k<4;++k) {
                    out.m[i][j]+=m[i][k]*matrix.m[k][j];
                }
            }
        }
        return out;
    }

    //
    // Define multiplication between a Mueller matrix and a constant...
    //
    MuellerMatrix 
    MuellerMatrix::
    operator*(double d) const
    {
        MuellerMatrix out;

        int i,j;
        for (i=0;i<4;++i) 
            for (j=0;j<4;++j) { 
                out.m[i][j]=m[i][j]*d;
            }
        return out;
    }

    MuellerMatrix 
    MuellerMatrix::
    operator/(double d) const
    {
        MuellerMatrix out;

        int i,j;
        for (i=0;i<4;++i) 
            for (j=0;j<4;++j) { 
                out.m[i][j]=m[i][j]/d;
            }
        return out;
    }

    //
    // Multiplication with a StokesVector...
    //
    StokesVector 
    MuellerMatrix::
    operator*(const StokesVector& in) const
    {
        StokesVector out;
        int i,j;
        for (i=0;i<4;++i) {
            out.s[i]=0;
            for (j=0;j<4;++j) 
                out.s[i]+=m[i][j]*in.s[j];
        }
        return out;
    }

    //
    // Define addition of two Mueller matrices...
    //
    MuellerMatrix 
    MuellerMatrix::
    operator+(const MuellerMatrix& a) const
    {
        MuellerMatrix c;
    
        int i,j;
        for (i=0;i<4;++i) {
            for (j=0;j<4;++j) {
                c.m[i][j] = m[i][j] + a.m[i][j];
            }
        }
        return c;
    }

    //
    // Define subtraction of two Mueller matrices...
    //
    MuellerMatrix 
    MuellerMatrix::
    operator-(const MuellerMatrix& a) const
    {
        MuellerMatrix c;
    
        int i,j;
        for (i=0;i<4;++i) {
            for (j=0;j<4;++j) {
                c.m[i][j] = m[i][j] - a.m[i][j];
            }
        }
        return c;
    }

    MuellerMatrix 
    MuellerMatrix::
    operator-() const
    {
        MuellerMatrix c;
    
        int i,j;
        for (i=0;i<4;++i) {
            for (j=0;j<4;++j) {
                c.m[i][j] = -m[i][j];
            }
        }
        return c;
    }

    MuellerMatrix 
    MuellerMatrix::
    rotate(const double angle) const
    {
        MuellerMatrix temp(*this);
        return (MuellerMatrix(JonesRotator(angle))*temp)*
                      MuellerMatrix(JonesRotator(-angle));
    }


    MuellerMatrix 
    MuellerMatrix::
    transpose() const
    {
        MuellerMatrix mm;
        int i,j;
        for (i=0;i<4;++i) 
            for (j=0;j<4;++j) 
                mm[i][j]=m[j][i];
        return mm;
    }


    MuellerMatrix 
    MuellerZero()
    {
        MuellerMatrix m;
        for (int i=0;i<4;++i) 
            for (int j=0;j<4;++j) 
                m[i][j]=0.;
        return m;
    }

    MuellerMatrix 
    MuellerUnit(double transmittance)
    {
        MuellerMatrix m=MuellerZero();
        m[0][0]=m[1][1]=m[2][2]=m[3][3]=transmittance;
        return m;
    }

    MuellerMatrix 
    MuellerDepolarizer(double transmittance,double depolarization)
    {
        MuellerMatrix m=MuellerZero();
        m[0][0]=transmittance;
        m[1][1]=m[2][2]=m[3][3]=transmittance*(1-depolarization);
        return m;
    }

    MuellerMatrix 
    MuellerPartialLinearPolarizer(double tmax, double tmin, double angle)
    {
        MuellerMatrix m;
    
        double tpt=0.5*(tmax+tmin);
        double tmt=0.5*(tmax-tmin);
        double tt= sqrt(tmax*tmin);
        double cos2a = cos(2*angle);
        double sin2a = sin(2*angle);

        m[0][0]=tpt;
        m[0][1]=tmt*cos2a;
        m[0][2]=tmt*sin2a;
        m[0][3]=0;
        m[1][0]=m[0][1];
        m[1][1]=tpt*sqr(cos2a)+tt*sqr(sin2a); // Correction: TAG 23 JUN 2004
        m[1][2]=(tpt-tt)*cos2a*sin2a;
        m[1][3]=0;
        m[2][0]=m[0][2];
        m[2][1]=m[1][2]; // Correction: TAG 5 DEC 2003
        m[2][2]=tt*sqr(cos2a)+tpt*sqr(sin2a);
        m[2][3]=0;
        m[3][0]=m[3][1]=m[3][2]=0;
        m[3][3]=tt;

        return m;
    }

    ostream& operator<<(ostream& os,const MuellerMatrix& m) 
    {
        os << '(';
        for (int i=0;i<4;++i) {
            os << '(';
            for (int j=0;j<4;++j) {
                os << m[i][j];
                if (j!=3) os << ',';
            }
            os << ')';
            if (i!=3) os << ',';
        }
        os << ')';
        return os;
    }

    istream& operator>>(istream& is,MuellerMatrix& m)
    {
        StokesVector s1,s2,s3,s4;
        char c;
        is >> c;
        if (c=='(') {
            is >> s1;
            is >> c;
            if (c==',') {
                is >> s2;
                is >> c;
                if (c==',') {
                    is >> s3;
                    is >> c;
                    if (c==',') {
                        if (c==',') {
                            is >> s4;
                            is >> c;
                            if (c==')') {
                                m = MuellerMatrix(s1,s2,s3,s4);
                                return is;
                            }
                        }
                    }
                }
            }
        }
        is.setstate(ios::failbit);
        return is;
    }

} // namespace SCATMECH