GPLIB++
MTFuncs.h
Go to the documentation of this file.
1 #ifndef MTFUNCS_H_
2 #define MTFUNCS_H_
3 #include <complex>
4 #include "miscfunc.h"
5 
6 namespace gplib
7  {
8  //double CalcRhoa(const MTTensor &Z){
9  //return 1./(2 * PI * Z.frequency) * mu * (pow(Z.real(),2)+pow(Z.at(k).imag(),2)) * pow(1000.0,2);
10  //}
11  class MTTensor;
12  double CalcPhi(const std::complex<double> &Z)
13  {
14  return atan(Z.imag() / Z.real()) / PI * 180;
15  }
16 
17  std::complex<double> CalcS1(const MTTensor &Z)
18  {
19  return Z.GetZxx() + Z.GetZyy();
20  }
21  /*
22 
23 
24  S1.at(i) = DataXX.Z.at(i) + DataYY.Z.at(i);
25  S2.at(i) = DataXY.Z.at(i) + DataYX.Z.at(i);
26  D1.at(i) = DataXX.Z.at(i) - DataYY.Z.at(i);
27  D2.at(i) = DataXY.Z.at(i) - DataYX.Z.at(i);
28  Berd.at(i) = 0.5 * D2.at(i);
29  dBerd.at(i) = sqrt(0.5*pow(DataXY.dZ.at(i),2)+pow(DataYX.dZ.at(i),2));
30  alpha.at(i) = 0.5 * atan((Commute(S1.at(i),S2.at(i))-Commute(D1.at(i),D2.at(i)))/(Commute(S1.at(i),D1.at(i))-Commute(S2.at(i),D2.at(i))));
31  eta.at(i) = sqrt(abs(Commute(D1.at(i),S2.at(i))-Commute(S1.at(i),D2.at(i))))/abs(D2.at(i));
32  det = DataXX.Z.at(i).real() * DataYY.Z.at(i).real() - DataXY.Z.at(i).real() * DataYX.Z.at(i).real();
33  detreal.at(i) = det;
34  Phi11.at(i) = (DataYY.Z.at(i).real() * DataXX.Z.at(i).imag() - DataXY.Z.at(i).real() * DataYX.Z.at(i).imag())/det;
35  Phi12.at(i) = (DataYY.Z.at(i).real() * DataXY.Z.at(i).imag() - DataXY.Z.at(i).real() * DataYY.Z.at(i).imag())/det;
36  Phi21.at(i) = (DataXX.Z.at(i).real() * DataYX.Z.at(i).imag() - DataYX.Z.at(i).real() * DataXX.Z.at(i).imag())/det;
37  Phi22.at(i) = (DataXX.Z.at(i).real() * DataYY.Z.at(i).imag() - DataYX.Z.at(i).real() * DataXY.Z.at(i).imag())/det;
38  alpha_phi.at(i) = 0.5 * atan((Phi12.at(i) + Phi21.at(i))/(Phi11.at(i) - Phi22.at(i)));
39  beta_phi.at(i) = 0.5 * atan((Phi12.at(i) - Phi21.at(i))/(Phi11.at(i) + Phi22.at(i)));
40  trPhi.at(i) = (Phi11.at(i) + Phi22.at(i));
41  skPhi.at(i) = Phi12.at(i) - Phi21.at(i);
42  detPhi.at(i) = Phi11.at(i) * Phi22.at(i) - Phi12.at(i) * Phi21.at(i);
43  Phi1.at(i) = trPhi.at(i) / 2.;
44  Phi2.at(i) = sqrt(detPhi.at(i));
45  Phi3.at(i) = skPhi.at(i) / 2.;
46  Phimax.at(i) = sqrt( pow(Phi1.at(i),2) + pow(Phi3.at(i),2)) + sqrt(pow(Phi1.at(i),2)+pow(Phi3.at(i),2) - pow(Phi2.at(i),2));
47  Phimin.at(i) = sqrt( pow(Phi1.at(i),2) + pow(Phi3.at(i),2)) - sqrt(pow(Phi1.at(i),2)+pow(Phi3.at(i),2) - pow(Phi2.at(i),2));
48  PhiEllip.at(i) = (Phimax.at(i)- Phimin.at(i))/ (Phimax.at(i)+Phimin.at(i));
49  temp1 = gsl_complex_rect(D1.at(i).real(),D1.at(i).imag());
50  temp2 = gsl_complex_rect(S2.at(i).real(),S2.at(i).imag());
51  temp3= gsl_complex_mul_real(gsl_complex_arctan(gsl_complex_div(gsl_complex_mul_real(temp1,-1.),temp2)),0.5);
52  thetabar.at(i) = GSL_REAL(temp3);
53  thetabar.at(i) += I*GSL_IMAG(temp3);
54  temp1 = gsl_complex_rect(S1.at(i).real(),S1.at(i).imag());
55  temp2 = gsl_complex_rect(D2.at(i).real(),D2.at(i).imag());
56  temp3 = gsl_complex_arctan(gsl_complex_div(temp1,temp2));
57  dtheta.at(i) = GSL_REAL(temp3);
58  dtheta.at(i) += I*GSL_IMAG(temp3);
59  //dtheta.at(i) = atan(S1.at(i)/D2.at(i))
60  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.);
61  Zp.at(i) = (DataXY.Z.at(i)*DataYX.Z.at(i)-DataXX.Z.at(i)*DataYY.Z.at(i))/Zs.at(i);
62  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);
63  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);
64  phis.at(i) = atan(Zs.at(i).imag()/Zs.at(i).real());
65  phip.at(i) = atan(Zp.at(i).imag()/Zp.at(i).real());
66  dPhips.at(i) = phip.at(i) - phis.at(i);
67  */
68  }
69 #endif /*MTFUNCS_H_*/
Stores MT-Tensor components at a single frequency, calculates derived quantities. ...
Definition: MTTensor.h:16
std::complex< double > CalcS1(const MTTensor &Z)
Definition: MTFuncs.h:17
std::complex< double > GetZxx() const
Definition: MTTensor.h:122
double CalcPhi(const std::complex< double > &Z)
Definition: MTFuncs.h:12
std::complex< double > GetZyy() const
Return tensor elements.
Definition: MTTensor.h:118