GPLIB++: CMTTensor.h Source File

00001 #ifndef CMTTENSOR_H_
00002 #define CMTTENSOR_H_
00003 #include <complex>
00004 #include <gsl/gsl_math.h>
00005 #include "miscfunc.h"
00006 
00007 
00008 //! Stores MT-Tensor components at a single frequency, calculates derived quantities
00009 class CMTTensor
00010 {
00011 private:
00012         //the impedance elements
00013         std::complex<double> Zxx;
00014         std::complex<double> Zxy;
00015         std::complex<double> Zyx;
00016         std::complex<double> Zyy;
00017         //and their errors
00018         double dZxx;
00019         double dZxy;
00020         double dZyx;
00021         double dZyy;
00022         double frequency;
00023         double rotangle;
00024         //! Coherency for the x direction
00025         double Rx;
00026         //! Coherency for the y direction
00027         double Ry;
00028         //! Number of degrees of freedom
00029         double Nu;
00030         //! Calculate apparent resistivity
00031         double CalcRho(const std::complex<double> &Z) const{
00032                 return  mu/(2 * PI * frequency)  * (Z.real()*Z.real()+Z.imag()*Z.imag()) * 1000000.0;}
00033         //! Calculate Phase
00034         double CalcPhase(const std::complex<double> &Z) const{ 
00035                 if (gsl_fcmp(Z.real(),0.0,GSL_DBL_EPSILON) !=0) //if real part significantly different from 0. 
00036                         return atan2(Z.imag(),Z.real())/ PI *180;
00037                 else
00038                         return 0.0;
00039         }
00040         //! Calculate Apparent Resistivity Error
00041         double CalcdRho(const std::complex<double> &Z, const double dZ) const 
00042                 {return mu/(PI * frequency)  * abs(Z) * dZ* 1000000.0;}
00043         //! Calculate Phase Error
00044         double CalcdPhase(const std::complex<double> &Z, const double dZ) const
00045                 {return dZ * abs(Z) * 1./(1.+ pow(Z.imag()/Z.real(),2)) * 1./pow(Z.real(),2)/ PI *180;  }
00046         //! Calculate Schmucker's Z*
00047         double CalcZStar(const std::complex<double> &Z) const
00048                 {return fabs((I/(2.0*PI*frequency)*Z*1000.0).real());}
00049         //! Calculate Schmucker's Rho*
00050         double CalcRhoStar(const std::complex<double> &Z) const
00051                 { double phase = CalcPhase(Z)/180.0 *PI;
00052                         if ( phase > 45.0)
00053                         return 2.0 * pow(cos(phase),2)*CalcRho(Z);
00054                   else
00055                         return 0.5 * 1./pow(sin(phase),2)*CalcRho(Z);}
00056 public:
00057         //! Function for Errors that cannot be calculated analytically when we don't want Jacknife errors
00058         double GetdZero() const {return 0.0;}
00059         //! Set the errors for the impedance elements 
00060         void SetErrors(double dxx, double dxy, double dyx, double dyy)
00061                 {dZxx=dxx; dZxy=dxy; dZyx=dyx; dZyy=dyy;}
00062         //!Rotate by the given angle in radian
00063         void Rotate(double angle);
00064         ///!return the current angle in radian
00065         double GetRotangle() const {return rotangle;}
00066         //! Set the rotation angle, without performing the corresponding rotation
00067         double &SetRotangle() {return rotangle;}
00068         //! Get the frequency for the impedance
00069         double GetFrequency() const{return frequency;}
00070         //! Return tensor elements
00071     std::complex<double> GetZyy() const{return Zyy;}
00072     std::complex<double> GetZxx() const{return Zxx;}
00073     std::complex<double> GetZxy() const{return Zxy;}
00074     std::complex<double> GetZyx() const{return Zyx;}
00075     //! Write access to tensor elements
00076     std::complex<double> &SetZyy() {return Zyy;}
00077     std::complex<double> &SetZxx() {return Zxx;}
00078     std::complex<double> &SetZxy() {return Zxy;}
00079     std::complex<double> &SetZyx() {return Zyx;}
00080     //! Return tensor element errors
00081     double GetdZxx() const{return dZxx;}
00082     double GetdZxy() const{return dZxy;}
00083     double GetdZyx() const{return dZyx;}
00084     double GetdZyy() const{return dZyy;}
00085     //! Write access to errors 
00086     double &SetdZyy() {return dZyy;}
00087     double &SetdZxx() {return dZxx;}
00088     double &SetdZxy() {return dZxy;}
00089     double &SetdZyx() {return dZyx;}
00090         //! Return apparent resistivity
00091         double GetRhoxx() const {return CalcRho(Zxx);}
00092         double GetRhoxy() const {return CalcRho(Zxy);}
00093         double GetRhoyx() const {return CalcRho(Zyx);}
00094         double GetRhoyy() const {return CalcRho(Zyy);}
00095         //! Return phase
00096         double GetPhixx() const {return CalcPhase(Zxx);}
00097         double GetPhixy() const {return CalcPhase(Zxy);}
00098         double GetPhiyx() const {return CalcPhase(Zyx);}
00099          double GetPhiyy() const {return CalcPhase(Zyy);}
00100         //! Return Rho Error for tensor elements
00101         double GetdRhoxx() const {return CalcdRho(Zxx,dZxx);}
00102         double GetdRhoxy() const {return CalcdRho(Zxy,dZxy);}
00103         double GetdRhoyx() const {return CalcdRho(Zyx,dZyx);}
00104         double GetdRhoyy() const {return CalcdRho(Zyy,dZyy);}
00105         //! return phase error for tensor elements
00106         double GetdPhixx() const {return CalcdPhase(Zxx,dZxx);}
00107         double GetdPhixy() const {return CalcdPhase(Zxy,dZxy);}
00108         double GetdPhiyx() const {return CalcdPhase(Zyx,dZyx);}
00109         double GetdPhiyy() const {return CalcdPhase(Zyy,dZyy);}
00110         //! Return Schmucker's rho* for tensor elements 
00111         double GetRhoxxStar() const {return CalcRhoStar(Zxx);}
00112         double GetRhoxyStar() const {return CalcRhoStar(Zxy);}
00113         double GetRhoyxStar() const {return CalcRhoStar(Zyx);}
00114         double GetRhoyyStar() const {return CalcRhoStar(Zyy);}
00115         //! Return Schmucker's z* for tensor elements
00116         double GetZxxStar() const { return CalcZStar(Zxx);}
00117         double GetZxyStar() const { return CalcZStar(Zxy);}
00118         double GetZyxStar() const { return CalcZStar(Zyx);}
00119         double GetZyyStar() const { return CalcZStar(Zyy);}
00120         //! Some invariants and intermediate quantities for strike and skew calculation
00121     std::complex<double> GetS1() const {return Zxx + Zyy;}
00122     std::complex<double> GetS2()const {return Zxy + Zyx;}
00123     std::complex<double> GetD1() const {return Zxx - Zyy;}
00124     std::complex<double> GetD2()  const {return Zxy - Zyx;}
00125     //! The Berdichevskyi invariant
00126     std::complex<double> GetBerd() const {return 0.5*GetD2();}
00127     //! The error of the Berdichevskyi invariant
00128     double GetdBerd() const { return 0.5*sqrt(dZxy*dZxy+dZyx*dZyx);}
00129     //! The determinant
00130     std::complex<double> GetDet() const {return Zxx * Zyy - Zxy*Zyx;}
00131     //! The error of the determinant
00132     double GetdDet() const { return sqrt( abs(Zyy*Zyy)*dZxx*dZxx + abs(Zxx*Zxx)*dZyy*dZyy
00133                 + abs(Zxy*Zxy)*dZyx*dZyx + abs(Zyx*Zyx)*dZxy*dZxy);}
00134     //! The determinant of the real parts of Z
00135     double GetDetreal() const {return Zxx.real() * Zyy.real() - Zxy.real() * Zyx.real();}
00136     //! Get the error of the determinant of the real part
00137     double GetdDetreal() const { return sqrt(gsl_pow_2(Zyy.real()) * gsl_pow_2(dZxx) + 
00138                 gsl_pow_2(Zxx.real()) * gsl_pow_2(dZyy) + gsl_pow_2(Zxy.real()) * gsl_pow_2(dZyx)
00139                 + gsl_pow_2(Zyx.real()) * gsl_pow_2(dZxy));}
00140     //! Rotationally invariant phase difference
00141     double GetMu() const { return sqrt(fabs(Commute(GetD1(),GetS2())+Commute(GetS1(),GetD2())))/abs(GetD2());}
00142     //! Swift's skew
00143     double GetKappa() const {return abs(GetS1())/abs(GetD2());}
00144     double GetSigma() const {return (gsl_pow_2(abs(GetD1()))+gsl_pow_2(abs(GetS2())))/gsl_pow_2(abs(GetD2()));}
00145     //! Bahr's skew
00146     double GetEta() const {
00147                 return sqrt(fabs(Commute(GetD1(),GetS2())-Commute(GetS1(),GetD2())))/abs(GetD2());}
00148     //Bahr's phase sensitive strike angle
00149     double GetAlpha() const{
00150                 double denominator = Commute(GetS1(),GetD1())-Commute(GetS2(),GetD2());
00151                 if (gsl_fcmp(denominator,0.0,GSL_DBL_EPSILON) !=0)
00152                         return 0.5 * atan((Commute(GetS1(),GetS2()) - Commute(GetD1(),GetD2()))/
00153                                         (denominator));
00154                 else
00155                         return 0.0;}
00156         //! Calculate strike angle, so it points to high conductivity direction (Phixy > Phiyx)
00157         double GetAlphaHigh() const{
00158                 double regalpha = GetAlpha(); // calculate strike angle
00159                 CMTTensor tmp(*this); //create local copy
00160                 tmp.Rotate(regalpha); //rotate to new coordinate system
00161                 if (tmp.GetPhiyx() > tmp.GetPhixy())
00162                         regalpha += PI/2.0;
00163                 return regalpha;
00164         }
00165         //! Maximum phase difference
00166         double GetMaxPhiDiff() const{
00167                 CMTTensor tmp(*this); //create local copy
00168                 tmp.Rotate(GetAlpha());
00169                 return fabs(tmp.GetPhixy() - tmp.GetPhiyx());
00170         }
00171         //! All the following quantities are defined in Caldwell GJI 158, 457-469, the phase tensor elements
00172     double GetPhi11() const {
00173                 const double dr = GetDetreal();
00174                 if (gsl_fcmp(dr,0.0,GSL_DBL_EPSILON) !=0)
00175                         return ((Zyy.real() * Zxx.imag() - Zxy.real()*Zyx.imag())/dr);
00176                 else 
00177                         return 0.0;}
00178         double GetPhi12() const {
00179                         const double dr = GetDetreal();
00180                     if (gsl_fcmp(dr,0.0,GSL_DBL_EPSILON) !=0)
00181                         return ((Zyy.real() * Zxy.imag() - Zxy.real()*Zyy.imag())/dr); 
00182                     else 
00183                         return 0.0;}
00184         double GetPhi21() const {
00185                         const double dr = GetDetreal();
00186                     if (gsl_fcmp(dr,0.0,GSL_DBL_EPSILON) !=0)
00187                         return ((Zxx.real() * Zyx.imag() - Zyx.real()*Zxx.imag())/dr); 
00188                     else 
00189                         return 0.0;}
00190         double GetPhi22() const {
00191                         const double dr = GetDetreal();
00192                         if (gsl_fcmp(dr,0.0,GSL_DBL_EPSILON) !=0)
00193                                 return ((Zxx.real() * Zyy.imag() - Zyx.real()*Zxy.imag())/dr);
00194                         else 
00195                             return 0.0;}
00196         double GetAlpha_phi() const {
00197                         const double diff = GetPhi11() - GetPhi22();
00198                         if (gsl_fcmp(diff,0.0,GSL_DBL_EPSILON) !=0)
00199                                 return 0.5* atan2((GetPhi12() + GetPhi21()),diff);
00200                         else
00201                                 return 0.0;}
00202         double GetBeta_phi() const {
00203                         return 0.5* atan2((GetPhi12() - GetPhi21()),(GetPhi11() + GetPhi22()));}
00204         double GetPi1() const { return 0.5 * sqrt( pow(GetPhi11() - GetPhi22(),2) + pow(GetPhi12() + GetPhi21(),2));}
00205         double GetPi2() const {return 0.5 * sqrt( pow(GetPhi11() + GetPhi22(),2) + pow(GetPhi12() - GetPhi21(),2));}
00206         double GetPhiStrike() const {return GetAlpha_phi() - GetBeta_phi();}
00207         double GetPhiMax() const {return GetPi2() + GetPi1();}
00208         double GetPhiMin() const {return GetPi2() - GetPi1();}
00209         double GettrPhi() const {return GetPhi11() + GetPhi22();}
00210         double GetskPhi() const {return GetPhi12() - GetPhi21();}
00211         double GetdetPhi() const {return GetPhi11() * GetPhi22() - GetPhi12()* GetPhi21();}
00212         double GetPhi1() const {return GettrPhi()/2.;}
00213         double GetPhi2() const {return sqrt(fabs(GetdetPhi()));}
00214         double GetPhi2Sq() const {return GetdetPhi();}
00215         double GetPhi3() const {return GetskPhi()/2.;}
00216         double GetPhiEllip() const {
00217                         return (GetPhiMax() - GetPhiMin())/(GetPhiMax()+GetPhiMin()); }
00218         //Here start the Weaver et al. invariants, equations taken from Marti, GJI 2004
00219         double GetXi1() const {return 0.5*(Zxx.real() + Zyy.real());}
00220         double GetXi2() const {return 0.5*(Zxy.real() + Zyx.real());}
00221         double GetXi3() const {return 0.5*(Zxx.real() - Zyy.real());}
00222         double GetXi4() const {return 0.5*(Zxy.real() - Zyx.real());}
00223         double GetEta1() const {return 0.5*(Zxx.imag() + Zyy.imag());}
00224         double GetEta2() const {return 0.5*(Zxy.imag() + Zyx.imag());}
00225         double GetEta3() const {return 0.5*(Zxx.imag() - Zyy.imag());}
00226         double GetEta4() const {return 0.5*(Zxy.imag() - Zyx.imag());}
00227         double GetI1() const {return sqrt(gsl_pow_2(GetXi1()) + gsl_pow_2(GetXi4()));}
00228         double GetI2() const {return sqrt(gsl_pow_2(GetEta1()) + gsl_pow_2(GetEta4()));}
00229         double GetI3() const {return sqrt(gsl_pow_2(GetXi2()) + gsl_pow_2(GetXi3()))/GetI1();}
00230         double GetI4() const {return sqrt(gsl_pow_2(GetEta2()) + gsl_pow_2(GetEta3()))/GetI2();}
00231         double GetI5() const {return (GetXi4()*GetEta1()+GetXi1()*GetEta4())/(GetI1()*GetI2());}
00232         double GetI6() const {return (GetXi4()*GetEta1()-GetXi1()*GetEta4())/(GetI1()*GetI2());}
00233         double Getd13() const {return (GetXi1()*GetEta3()-GetXi3()*GetEta1())/(GetI1()*GetI2());}
00234         double Getd12() const {return (GetXi1()*GetEta2()-GetXi2()*GetEta1())/(GetI1()*GetI2());}
00235         double Getd24() const {return (GetXi2()*GetEta4()-GetXi4()*GetEta2())/(GetI1()*GetI2());}
00236         double Getd34() const {return (GetXi3()*GetEta4()-GetXi4()*GetEta3())/(GetI1()*GetI2());}
00237         double Getd41() const {return (GetXi4()*GetEta1()-GetXi1()*GetEta4())/(GetI1()*GetI2());}
00238         double Getd23() const {return (GetXi2()*GetEta3()-GetXi3()*GetEta2())/(GetI1()*GetI2());}
00239         double GetQ() const {return sqrt(gsl_pow_2(Getd12()-Getd34())+gsl_pow_2(Getd13()+Getd24()));}
00240         double GetI7() const {return (Getd41() + Getd23())/GetQ();}
00241         double Geta() const {return gsl_pow_2(GetI5() - GetI6());}
00242         double Getb() const {return 1.0 - GetI5() * GetI6() + sqrt(1.0+gsl_pow_2(GetI5())*gsl_pow_2(GetI6())
00243                                 - gsl_pow_2(GetI5()) - gsl_pow_2(GetI6()));
00244         }
00245         double Getr() const {return GetI2()/GetI1();}
00246         //const std::complex<double> GetThetaBar(){}
00247         //const std::complex<double> GetdTheta(){}
00248         //const std::complex<double> GetZs(){}
00249         //const std::complex<double> GetZp(){}
00250         //! Return coherency for the x-direction
00251         double GetRx()const{return Rx;}
00252         //! Return coherency for the y-direction
00253         double GetRy()const{return Ry;}
00254         //! The degrees of freedom used for transfer function estimation
00255         double GetNu()const{return Nu;}
00256         friend class CMTStation;
00257     friend class C1DMTSynthData;
00258     friend class JParser;
00259     friend class EDIParser;
00260         CMTTensor();
00261         CMTTensor(const std::complex<double> &xx,const std::complex<double> &xy,const std::complex<double> &yx,
00262                         const std::complex<double> &yy, const double freq = 1., const double angle = 0.0);
00263         CMTTensor& operator= (const CMTTensor& source);
00264         virtual ~CMTTensor();
00265 };
00266 
00267 std::complex<double> inline RhoPhiToZ(const double freq, const double rho, const double phi){
00268                 return sqrt(2 * PI * freq/mu) *0.001* sqrt(rho)*(cos(phi/180*PI)+ I*sin(phi/180*PI));
00269         }
00270 
00271         // This is some remnant code for parallel and serial transforms of the impedance tensor
00272         // only kept so I remember how they are calculated
00273     /*  temp1 = gsl_complex_rect(D1.at(i).real(),D1.at(i).imag());
00274         temp2 = gsl_complex_rect(S2.at(i).real(),S2.at(i).imag()); 
00275         temp3= gsl_complex_mul_real(gsl_complex_arctan(gsl_complex_div(gsl_complex_mul_real(temp1,-1.),temp2)),0.5);
00276         thetabar.at(i) = GSL_REAL(temp3);
00277         thetabar.at(i) += I*GSL_IMAG(temp3);
00278         temp1 = gsl_complex_rect(S1.at(i).real(),S1.at(i).imag());
00279         temp2 = gsl_complex_rect(D2.at(i).real(),D2.at(i).imag());
00280         temp3 = gsl_complex_arctan(gsl_complex_div(temp1,temp2));
00281         dtheta.at(i) = GSL_REAL(temp3);
00282         dtheta.at(i) += I*GSL_IMAG(temp3);
00283         //dtheta.at(i) = atan(S1.at(i)/D2.at(i))
00284         Zs.at(i) = sqrt((pow(DataXX.Z.at(i),2) + pow(DataXY.Z.at(i),2) + pow(DataYX.Z.at(i),2) + pow(DataYY.Z.at(i),2))/2.);
00285         Zp.at(i) = (DataXY.Z.at(i)*DataYX.Z.at(i)-DataXX.Z.at(i)*DataYY.Z.at(i))/Zs.at(i);
00286         rhos.at(i) = 1/(2 * PI * frequency.at(i)) * mu * (pow(Zs.at(i).real(),2)+pow(Zs.at(i).imag(),2)) * pow(1000.0,2);
00287         rhop.at(i) = 1/(2 * PI * frequency.at(i)) * mu * (pow(Zp.at(i).real(),2)+pow(Zp.at(i).imag(),2)) * pow(1000.0,2);
00288         phis.at(i) = atan(Zs.at(i).imag()/Zs.at(i).real());
00289         phip.at(i) = atan(Zp.at(i).imag()/Zp.at(i).real());
00290         dPhips.at(i) = phip.at(i) - phis.at(i);*/
00291         
00292 
00293 
00294 #endif /*CMTTENSOR_H_*/