//****************************************************************************** //** 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