GPLIB++: CMTStation.cpp Source File

00001 #include "CMTStation.h"
00002 #include "miscfunc.h"
00003 #include "EDILexer.hpp"
00004 #include "EDIParser.hpp"
00005 #include "JLexer.hpp"
00006 #include "JParser.hpp"
00007 #include <iterator>
00008 #include <gsl/gsl_complex.h>
00009 #include <gsl/gsl_complex_math.h>
00010 #include <string>
00011 #include "CFatalException.h"
00012 #include <boost/filesystem/operations.hpp>
00013 #include <boost/bind.hpp>
00014 #include <functional>
00015 #include <fstream>
00016 #include <cassert>
00017 using namespace std;
00018 
00019 void TrimFilename(string &name)
00020 {
00021         unsigned int slashpos = name.find('/',0);
00022     if (slashpos != string::npos)
00023                 name = name.substr(slashpos+1);
00024     unsigned int dotpos = name.find('.',0);
00025     if (dotpos != string::npos)
00026                 name.erase(dotpos);
00027 }
00028 
00029 void CMTStation::InitialSetup()
00030 {
00031         latitude = 0;
00032         longitude = 0;
00033         elevation = 0;
00034         azimuth = 0;
00035 }
00036 
00037 CMTStation::CMTStation()
00038 {
00039         InitialSetup();
00040 }
00041 
00042 CMTStation::~CMTStation()
00043 {}
00044 
00045 CMTStation::CMTStation(const std::string filename)
00046 {
00047         InitialSetup();
00048         GetData(filename);
00049 }
00050 
00051 CMTStation::CMTStation(const int size):
00052 MTData(size)
00053 {
00054         InitialSetup();
00055         Assign(size);
00056 }
00057 void CMTStation::AssignAll(const int nfreq)
00058 {
00059         Assign(nfreq);
00060 }
00061 void CMTStation::Assign(const int nfreq)
00062 {
00063         MTData.assign(nfreq,CMTTensor());
00064         TFData.assign(nfreq,CMagneticTF());
00065 }
00066 
00067 /*! The Update() method first updates all components individually and then 
00068  * recalculates all derived quantities that depend on more than one component
00069  */
00070  
00071 void CMTStation::Update()
00072 {
00073         
00074 }
00075 
00076 /*! Rotates the component indicated by i by angle rotangle. Does not Update. */
00077 void CMTStation::DoRotate(const int i, const double angle)
00078 {
00079          dcomp newxx, newxy, newyx, newyy;
00080          newxx = MTData.at(i).Zxx * pow(cos(angle),2) + (MTData.at(i).Zxy + MTData.at(i).Zyx)*  sin(angle) * cos(angle) 
00081             + MTData.at(i).Zyy * pow(sin(angle),2);
00082         newxy = MTData.at(i).Zxy * pow(cos(angle),2) - (MTData.at(i).Zxx - MTData.at(i).Zyy)*  sin(angle) * cos(angle) 
00083             - MTData.at(i).Zyx * pow(sin(angle),2);
00084         newyx = MTData.at(i).Zyx * pow(cos(angle),2) - (MTData.at(i).Zxx - MTData.at(i).Zyy)*  sin(angle) * cos(angle) 
00085             - MTData.at(i).Zxy * pow(sin(angle),2);
00086         newyy = MTData.at(i).Zyy * pow(cos(angle),2) - (MTData.at(i).Zxy + MTData.at(i).Zyx)*  sin(angle) * cos(angle) 
00087             + MTData.at(i).Zxx * pow(sin(angle),2);
00088         MTData.at(i).Zxx = newxx;
00089         MTData.at(i).Zxy = newxy;
00090         MTData.at(i).Zyx = newyx;
00091         MTData.at(i).Zyy = newyy;
00092         MTData.at(i).rotangle += angle;
00093 }
00094 
00095 /*! Rotate(const double rotangle) rotates the impedance tensor by the angle rotangle in radian
00096  * and updates all derived quantities
00097  */
00098 void CMTStation::Rotate(const double rotangle)
00099 {
00100     for (unsigned int i = 0; i < MTData.size() ; ++i)
00101     {
00102            DoRotate(i,rotangle);
00103     }
00104     Update();
00105 }
00106 
00107 /*! Rotate(void) rotates the impedance tensor by the negative of the angle given in the rotangles field
00108  */
00109 void CMTStation::Rotate(void)
00110 {
00111     for (unsigned int i = 0; i < MTData.size() ; ++i)
00112     {
00113                 DoRotate(i,-MTData.at(i).rotangle);
00114     }
00115     Update();
00116 }
00117 /*! Return a vector containing the available frequencies; This should be rewritten to be more efficient
00118  * */
00119 trealdata CMTStation::GetFrequencies() const
00120 {
00121         trealdata temp(MTData.size());
00122         transform(MTData.begin(),MTData.end(),temp.begin(),mem_fun_ref(&CMTTensor::GetFrequency));
00123         return temp;
00124 }
00125 
00126 void CMTStation::SetFrequencies(const trealdata &freqs)
00127 {
00128         //assert(MTData.size() == freqs.size());
00129         const unsigned int nfreqs = freqs.size();
00130         if (MTData.size() != nfreqs) //this invalidates the data stored there before
00131         {
00132                 MTData.resize(nfreqs);
00133                 TFData.resize(nfreqs);
00134         }
00135         for (unsigned int i = 0; i < nfreqs; ++i)
00136                 MTData.at(i).frequency = freqs.at(i);
00137         assert(MTData.size() == TFData.size());
00138         
00139         
00140 }
00141 
00142 void CMTStation::ReadEdi(const std::string filename)
00143 {
00144         using namespace antlr;
00145         ifstream infile(filename.c_str());
00146         const double invalid = 1e30;
00147         
00148         if (infile)
00149         {
00150                 try
00151                 {
00152                         EDILexer lexer(infile);
00153                         EDIParser parser(lexer);
00154                         parser.edi_file();
00155                         latitude = parser.latitude;
00156                         longitude = parser.longitude;
00157                         elevation = parser.elevation;
00158                         azimuth = parser.azimuth;
00159                         name = filename;
00160                         TrimFilename(name);
00161                     unsigned int valfreqs = 0;
00162                         for (unsigned int i= 0; i < parser.frequency.size(); ++i)
00163                         {
00164                                 if (abs(parser.DataXX.Z.at(i)) > invalid)
00165                                 {
00166                                         parser.frequency.at(i) = 0;
00167                                 }
00168                                 else
00169                                 {
00170                                         valfreqs++;
00171                                 }
00172                         }
00173                         //frequency.assign(parser.frequency.size(),0);
00174                         //copy(parser.frequency.begin(),parser.frequency.end(),frequency.begin());
00175                         //rotangles.assign(parser.frequency.size(),0);
00176                         azimuth += parser.rotangles.at(0)*180./PI;
00177                         Assign(valfreqs);
00178                         for (unsigned int i = 0; i < parser.rotangles.size(); ++i)
00179                         {
00180                                 if (abs(parser.DataXX.Z.at(i)) < invalid)
00181                                 {
00182                                         MTData.at(i).frequency = parser.frequency.at(i);
00183                                         MTData.at(i).rotangle = parser.rotangles.at(i) * PI/180.;
00184                                         MTData.at(i).Zxx = parser.DataXX.Z.at(i);
00185                                         MTData.at(i).Zxy = parser.DataXY.Z.at(i);
00186                                         MTData.at(i).Zyx = parser.DataYX.Z.at(i);
00187                                         MTData.at(i).Zyy = parser.DataYY.Z.at(i);
00188                                         MTData.at(i).dZxx = parser.DataXX.dZ.at(i);
00189                                         MTData.at(i).dZxy = parser.DataXY.dZ.at(i);
00190                                         MTData.at(i).dZyx = parser.DataYX.dZ.at(i);
00191                                         MTData.at(i).dZyy = parser.DataYY.dZ.at(i);
00192                                 }
00193                         }
00194                         
00195                         if (parser.DataZY.Z.empty() == false)
00196                         {
00197                                 for (unsigned int i = 0; i < parser.DataZY.Z.size(); ++i)
00198                                 {
00199                                         if (abs(parser.DataZX.Z.at(i)) < invalid)
00200                                         {
00201                                                 TFData.at(i).frequency = parser.frequency.at(i);
00202                                                 TFData.at(i).Tx = parser.DataZX.Z.at(i);
00203                                                 TFData.at(i).Ty = parser.DataZY.Z.at(i);
00204                                                 TFData.at(i).dTx = parser.DataZX.dZ.at(i);
00205                                                 TFData.at(i).dTy = parser.DataZY.dZ.at(i);
00206                                         }
00207                                 }
00208                         }
00209                 }
00210                 catch(ANTLRException& e)
00211                 {
00212                         cerr << "Parse exception: " << e.toString() << endl;
00213                 }
00214             Update();
00215         }
00216         else
00217     {
00218         throw CFatalException("File not found: "+filename); 
00219     }
00220 }
00221 
00222 void CMTStation::ReadPek1D(const std::string filename)
00223 {
00224         ifstream infile(filename.c_str());
00225         if (infile)
00226         {
00227                 const double convfactor = 1./(1000. * mu);
00228                 double period, rxx, imxx, rxy, imxy, ryx, imyx, ryy, imyy;
00229                 while (infile.good())
00230                 {
00231                         infile >> period >> rxx >> imxx >> rxy >> imxy >> ryx >> imyx >> ryy >> imyy;
00232                         
00233                         if (infile.good())
00234                         {
00235                                 CMTTensor CurrentImp;
00236                                 CurrentImp.frequency = 1./period;
00237                                 CurrentImp.Zxx = (rxx + I * imxx) * convfactor;
00238                                 CurrentImp.Zxy = (rxy + I * imxy) * convfactor;
00239                                 CurrentImp.Zyx = (ryx + I * imyx) * convfactor;
00240                                 CurrentImp.Zyy = (ryy + I * imyy) * convfactor;
00241                                 MTData.push_back(CurrentImp);
00242                                 TFData.push_back(CMagneticTF());
00243                         }
00244                 }
00245                 name = filename;
00246                 TrimFilename(name);
00247         }
00248         else
00249     {
00250         throw CFatalException("File not found: "+filename); 
00251     }
00252 }
00253 
00254 void CMTStation::ReadJ(const std::string filename)
00255 {
00256                 
00257         using namespace antlr;
00258         ifstream infile(filename.c_str());
00259         double factor;
00260         
00261         if (infile)
00262         {
00263                 try
00264            {
00265                         JLexer lexer(infile);
00266                         JParser parser(lexer);
00267                         parser.jfile();
00268                         latitude = parser.latitude;
00269                         longitude = parser.longitude;
00270                         elevation = parser.elevation;
00271                         azimuth = parser.azimuth;
00272                         name = parser.name;
00273                         TrimFilename(name);
00274                         Assign( parser.frequency.size());
00275                         
00276                         //frequency = parser.frequency;
00277                         if ( parser.zassigned == false)
00278                         {
00279                                 if (parser.rassigned == false)
00280                                 {
00281                                         cerr << "No MT data in file !" << endl;
00282                                 }
00283                                 for (unsigned int i =0; i < parser.frequency.size(); ++i)
00284                                 {               
00285                                         MTData.at(i).frequency = parser.frequency.at(i);
00286                                         factor = sqrt(2 * PI * parser.frequency.at(i)/mu);
00287                                         if ((parser.DataXX.rhoa.size() > 0) )
00288                                                 MTData.at(i).Zxx = factor * sqrt(parser.DataXX.rhoa.at(i))
00289                                                         *(cos(parser.DataXX.phi.at(i)/180*PI)+ I*sin(parser.DataXX.phi.at(i)/180*PI));
00290                                         if ((parser.DataXY.rhoa.size() > 0) )
00291                                                 MTData.at(i).Zxy = factor * sqrt(parser.DataXY.rhoa.at(i))
00292                                                         *(cos(parser.DataXY.phi.at(i)/180*PI)+ I*sin(parser.DataXY.phi.at(i)/180*PI));
00293                                         if ((parser.DataYX.rhoa.size() > 0) )
00294                                                 MTData.at(i).Zyx = factor * sqrt(parser.DataYX.rhoa.at(i))
00295                                                         *(cos(parser.DataYX.phi.at(i)/180*PI)+ I*sin(parser.DataYX.phi.at(i)/180*PI));
00296                                         if ((parser.DataYY.rhoa.size() > 0) )
00297                                                 MTData.at(i).Zyy = factor * sqrt(parser.DataYY.rhoa.at(i))
00298                                                         *(cos(parser.DataYY.phi.at(i)/180*PI)+ I*sin(parser.DataYY.phi.at(i)/180*PI));
00299                                 }       
00300                         }
00301                         else
00302                         {
00303                                 for (unsigned int i =0; i < parser.frequency.size(); ++i)
00304                                 {
00305                                         MTData.at(i).frequency = parser.frequency.at(i);
00306                                         MTData.at(i).rotangle = parser.azimuth;
00307                                         MTData.at(i).Zxx = parser.DataXX.Z.at(i);
00308                                         MTData.at(i).Zxy = parser.DataXY.Z.at(i);
00309                                         MTData.at(i).Zyx = parser.DataYX.Z.at(i);
00310                                         MTData.at(i).Zyy = parser.DataYY.Z.at(i);
00311                                         MTData.at(i).dZxx = parser.DataXX.dZ.at(i);
00312                                         MTData.at(i).dZxy = parser.DataXY.dZ.at(i);
00313                                         MTData.at(i).dZyx = parser.DataYX.dZ.at(i);
00314                                         MTData.at(i).dZyy = parser.DataYY.dZ.at(i);
00315                                         MTData.at(i).Rx = parser.Rx.at(i);
00316                                         MTData.at(i).Ry = parser.Ry.at(i);
00317                                 }
00318                         }
00319                         if (parser.tassigned)
00320                                 {
00321                                         for (unsigned int i = 0; i < parser.DataZY.Z.size(); ++i)
00322                                         {
00323                                                 TFData.at(i).frequency = parser.frequency.at(i);
00324                                                 TFData.at(i).Tx = parser.DataZX.Z.at(i);
00325                                                 TFData.at(i).Ty = parser.DataZY.Z.at(i);
00326                                                 TFData.at(i).dTx = parser.DataZX.dZ.at(i);
00327                                                 TFData.at(i).dTy = parser.DataZY.dZ.at(i);
00328                                                 TFData.at(i).Rz = parser.Rz.at(i);
00329                                         }
00330                                 }
00331                 }
00332                 catch(ANTLRException& e)
00333                 {
00334                         cerr << "Parse exception: " << e.toString() << endl;    
00335                 }
00336             Update();
00337         }
00338         else
00339     {
00340         throw CFatalException("File not found: "+filename); 
00341     }
00342         
00343 }
00344 
00345 void CMTStation::ReadMtt(const std::string filename)
00346 {
00347         ifstream infile;
00348     double currentreal, currentimag;
00349         int nentries = 0;
00350         infile.open(filename.c_str());
00351         while ( infile.good())
00352         {
00353                 infile >> currentreal;
00354                 ++nentries;
00355         }
00356         infile.close();
00357         if ( ((nentries-1) % 23) != 0)
00358                 throw CFatalException("Number of records does not match expected: "+filename); 
00359         const int nrecords = (nentries-1) / 23;
00360         Assign(nrecords);
00361         infile.open(filename.c_str());
00362         int currentrecord = 0;
00363     if (infile)
00364         {
00365         while ( infile.good()) // read in inputfile
00366                 {       
00367                         infile >> currentreal;
00368                         if (infile.good())
00369                         {
00370                         MTData.at(currentrecord).frequency = currentreal;
00371                         TFData.at(currentrecord).frequency = currentreal;
00372                         infile >> currentreal;
00373                         MTData.at(currentrecord).Nu = currentreal;
00374                         infile >> currentreal >> currentimag;
00375                         MTData.at(currentrecord).Zxx = (currentreal + I * currentimag); 
00376                         infile >> currentreal >> currentimag;
00377                         MTData.at(currentrecord).Zxy = (currentreal + I * currentimag); 
00378                         infile >> currentreal >> currentimag;
00379                                 MTData.at(currentrecord).Zyx = (currentreal + I * currentimag); 
00380                             infile >> currentreal >> currentimag;
00381                                 MTData.at(currentrecord).Zyy = (currentreal + I * currentimag); 
00382                             infile >> currentreal;
00383                             MTData.at(currentrecord).dZxx = currentreal; 
00384                             infile >> currentreal;
00385                             MTData.at(currentrecord).dZxy = currentreal; 
00386                             infile >> currentreal;
00387                             MTData.at(currentrecord).dZyx = currentreal; 
00388                             infile >> currentreal;
00389                             MTData.at(currentrecord).dZyy = currentreal; 
00390                             infile >> currentreal >> currentimag;
00391                             TFData.at(currentrecord).Tx = currentreal + I * currentimag;
00392                             infile >> currentreal >> currentimag;
00393                             TFData.at(currentrecord).Ty = currentreal + I * currentimag;
00394                             infile >> currentreal;
00395                             TFData.at(currentrecord).dTx = currentreal; 
00396                             infile >> currentreal;
00397                             TFData.at(currentrecord).dTy = currentreal; 
00398                             infile >> currentreal;
00399                             MTData.at(currentrecord).Rx = currentreal;
00400                             infile >> currentreal;
00401                             MTData.at(currentrecord).Ry = currentreal;
00402                             infile >> currentreal;
00403                             TFData.at(currentrecord).Rz = currentreal;   
00404                             ++currentrecord;
00405                         }  
00406                 }
00407         infile.close();
00408         name = filename;
00409         TrimFilename(name);
00410     }
00411     else
00412     {
00413         throw CFatalException("File not found: "+filename); 
00414     }
00415     Update();
00416 }
00417 
00418 
00419 void CMTStation::GetData(const string filename)
00420 {
00421         if (!boost::filesystem::exists(filename))
00422                 throw CFatalException("File does not exist : "+filename);
00423         string ending;
00424         unsigned int dotpos = filename.find('.',0);
00425         if (dotpos != string::npos)
00426                 ending = filename.substr(dotpos);
00427         else
00428         {
00429                 cerr << "File has no extension ! " << filename << endl;
00430                 dataformat = unknown;
00431                 throw CFatalException("File has no extension ! " + filename);
00432         }
00433         if ( ending == ".mtt")
00434                                         {
00435                                                 ReadMtt(filename);
00436                                                 dataformat = mtt;
00437                                         }
00438                                         else
00439                                                 if      (ending == ".dat" || ending == ".j")
00440                                                 {
00441                                                         ReadJ(filename);
00442                                                         dataformat = j;
00443                                                 }
00444                                                 else
00445                                                         if      (filename.substr(filename.size()-4) == ".edi")
00446                                                         {
00447                                                                 ReadEdi(filename);
00448                                                                 dataformat = edi;
00449                                                         }
00450                                                         else
00451                                                                 if (ending == ".pek")
00452                                                                 {
00453                                                                         ReadPek1D(filename);
00454                                                                         dataformat = pek;
00455                                                                 }
00456                                                                 else
00457                                                                 {
00458                                                                         dataformat = unknown;
00459                                                                         throw CFatalException("File not recognised: "+filename); 
00460                                                                 }
00461         assert(MTData.size() == TFData.size());
00462         Rotate(); //rotate data to 0 degree 
00463                 
00464 }
00465 
00466 void CMTStation::WriteData(const std::string filename)
00467 {
00468         name = filename;
00469         WriteBack();
00470 }
00471 
00472 void CMTStation::WriteBack()
00473 {
00474         switch (dataformat)
00475         {
00476                 case mtt:
00477                         WriteAsMtt(name.c_str());
00478                         break;
00479                 case j:
00480                         WriteAsJ(name.c_str());
00481                         break;
00482                 case edi:
00483                         WriteAsEdi(name.c_str());
00484                         break;
00485                 default:
00486                         throw CFatalException("Unknown File Format for Writing ! This should not happen");
00487         }
00488         
00489 }
00490 
00491 void CMTStation::WriteAsMtt(const string filename)
00492 {
00493         Update();
00494     WriteMtt((filename+".mtt").c_str());
00495 }
00496 
00497 void CMTStation::WriteAsEdi(const string filename)
00498 {
00499         throw CFatalException("WriteAsEdi not implemented yet !");
00500 }
00501 
00502 
00503 
00504 void CMTStation::WriteJBlock(boost::function<complex<double> (const CMTTensor*)> Comp, 
00505         boost::function<double (const CMTTensor*)> Err, ofstream &outfile, const double convfactor)
00506 {
00507         for (unsigned int i = 0; i < MTData.size(); ++i)
00508         {
00509                 if (MTData.at(i).frequency != 0)
00510                 {
00511                         outfile <<  setfill(' ') << setw(15) <<resetiosflags(ios::fixed)<<  1./MTData.at(i).frequency << " ";
00512                         outfile <<  setfill(' ') << setw(15) <<  convfactor * boost::bind(Comp,&MTData.at(i))().real() << " ";
00513                         outfile <<  setfill(' ') << setw(15) <<  convfactor * boost::bind(Comp,&MTData.at(i))().imag() << " ";
00514                         outfile <<  setfill(' ') << setw(15) <<  convfactor * boost::bind(Err,&MTData.at(i))() << " ";
00515                         outfile << setfill(' ')  << setw(15)  <<setiosflags(ios::fixed)<<  1.000 << " " << endl;
00516                 }
00517         }
00518 }
00519 
00520 void CMTStation::WriteAsJ(const string filename)
00521 {
00522         ofstream outfile;
00523         outfile.open((filename+".j").c_str());
00524         int actfreqs = 0;
00525         const double convfactor = 4. * acos(-1.)*1e-4;
00526         
00527         outfile << "# J-File Produced by CJData.cpp" << endl;
00528         outfile << ">LATITUDE  = " << latitude  << endl; 
00529         outfile << ">LONGITUDE = " << longitude << endl; 
00530         outfile << ">ELEVATION = " << elevation << endl;
00531         outfile << ">AZIMUTH   = " << azimuth << endl;
00532         outfile << name << endl;
00533         for (unsigned int i = 0; i < MTData.size(); ++i)
00534                 if (MTData.at(i).frequency != 0)
00535                         actfreqs++;
00536         outfile << "ZXX S.I." << endl;
00537         outfile << actfreqs << endl;
00538         //WriteJBlock(DataXX,outfile,convfactor);
00539         WriteJBlock(&CMTTensor::GetZxx,&CMTTensor::GetdZxx,outfile,convfactor);
00540         outfile << "ZXY S.I." << endl;
00541         outfile << actfreqs << endl;
00542         WriteJBlock(&CMTTensor::GetZxy,&CMTTensor::GetdZxy,outfile,convfactor);
00543         
00544         outfile << "ZYX S.I." << endl;
00545         outfile << actfreqs << endl;
00546         WriteJBlock(&CMTTensor::GetZyx,&CMTTensor::GetdZyx,outfile,convfactor);
00547         
00548         outfile << "ZYY S.I." << endl;
00549         outfile << actfreqs << endl;
00550         WriteJBlock(&CMTTensor::GetZyy,&CMTTensor::GetdZyy,outfile,convfactor);
00551         
00552 /*      if (norm(DataZY.Z.at(0)) !=0)
00553         {
00554                 outfile << "TZX " << endl;
00555                 outfile << actfreqs << endl;
00556                 WriteJBlock(DataZX,outfile,1.);
00557         }
00558         
00559         if (norm(DataZY.Z.at(0)) !=0)
00560         {
00561                 outfile << "TZY " << endl;
00562                 outfile << actfreqs << endl;
00563                 WriteJBlock(DataZY,outfile,1.);
00564         }*/
00565         outfile.close();
00566 }
00567 
00568 void CMTStation::WriteMtt(const std::string filename)
00569 {
00570     ofstream outfile;
00571 
00572     outfile.open(filename.c_str());
00573     vector<int> mttindex(MTData.size());
00574    // index(frequency, mttindex);
00575     assert(MTData.size() == TFData.size());
00576     for (unsigned int i = 0; i < MTData.size(); ++i) //write mtt-file
00577     {
00578         //k = mttindex.at(i);
00579         if (MTData.at(i).frequency > 0)
00580         {
00581                 outfile << setw(9) << setfill(' ') << setprecision(4) << setiosflags(ios::scientific)  << MTData.at(i).frequency;
00582                 outfile  << "  " << resetiosflags(ios::scientific) << setprecision(5) << MTData.at(i).Nu << "\n";
00583                 
00584             outfile << setw(9) << setfill(' ') << setprecision(4) <<setiosflags(ios::fixed)  <<  MTData.at(i).Zxx.real() << " ";
00585                 outfile << setw(9) << setfill(' ') << setprecision(4) <<setiosflags(ios::fixed)  <<  MTData.at(i).Zxx.imag() << " ";
00586                 outfile << setw(9) << setfill(' ') << setprecision(4) <<setiosflags(ios::fixed)  <<  MTData.at(i).Zxy.real() << " ";
00587                 outfile << setw(9) << setfill(' ') << setprecision(4) <<setiosflags(ios::fixed)  <<  MTData.at(i).Zxy.imag() << " ";
00588                 outfile << setw(9) << setfill(' ') << setprecision(4) <<setiosflags(ios::fixed)  <<  MTData.at(i).Zyx.real() << " ";
00589                 outfile << setw(9) << setfill(' ') << setprecision(4) <<setiosflags(ios::fixed)  <<  MTData.at(i).Zyx.imag() << " ";
00590                 outfile << setw(9) << setfill(' ') << setprecision(4) <<setiosflags(ios::fixed)  <<  MTData.at(i).Zyy.real() << " ";
00591                 outfile << setw(9) << setfill(' ') << setprecision(4) <<setiosflags(ios::fixed)  <<  MTData.at(i).Zyy.imag() << " ";
00592                 outfile << "\n";
00593                 outfile << setw(9) << setfill(' ') << setprecision(4) <<setiosflags(ios::fixed)  <<  MTData.at(i).dZxx << " ";
00594                 outfile << setw(9) << setfill(' ') << setprecision(4) <<setiosflags(ios::fixed)  <<  MTData.at(i).dZxy<< " ";
00595                 outfile << setw(9) << setfill(' ') << setprecision(4) <<setiosflags(ios::fixed)  <<  MTData.at(i).dZyx << " ";
00596                 outfile << setw(9) << setfill(' ') << setprecision(4) <<setiosflags(ios::fixed)  <<  MTData.at(i).dZyy << " ";
00597                 outfile << setw(9) << setfill(' ') << setprecision(4) <<setiosflags(ios::fixed)  <<  TFData.at(i).Tx.real() << " ";
00598                 outfile << setw(9) << setfill(' ') << setprecision(4) <<setiosflags(ios::fixed)  <<  TFData.at(i).Tx.imag() << " ";
00599                 outfile << setw(9) << setfill(' ') << setprecision(4) <<setiosflags(ios::fixed)  <<  TFData.at(i).Ty.real() << " ";
00600                 outfile << setw(9) << setfill(' ') << setprecision(4) <<setiosflags(ios::fixed)  <<  TFData.at(i).Ty.imag() << " ";
00601                 outfile << "\n";
00602                 outfile << setw(9) << setfill(' ') << setprecision(4) <<setiosflags(ios::fixed)  <<  TFData.at(i).dTx << " ";
00603                 outfile << setw(9) << setfill(' ') << setprecision(4) <<setiosflags(ios::fixed)  <<  TFData.at(i).dTy << " ";
00604                 outfile << setw(9) << setfill(' ') << setprecision(4) <<setiosflags(ios::fixed)  <<  MTData.at(i).Rx << " ";
00605                 outfile << setw(9) << setfill(' ') << setprecision(4) <<setiosflags(ios::fixed)  <<  MTData.at(i).Ry << " ";
00606                 outfile << setw(9) << setfill(' ') << setprecision(4) <<setiosflags(ios::fixed)  <<  TFData.at(i).Rz << " ";
00607                 outfile << "\n" << resetiosflags(ios::fixed);
00608                 //write out mtt file entries
00609         }
00610         
00611     } 
00612     outfile.close();
00613 }
00614 CMTStation::CMTStation(const CMTStation &old):
00615         MTData(old.MTData),
00616         TFData(old.TFData),
00617         latitude(old.latitude),
00618         longitude(old.longitude),
00619         elevation(old.elevation),
00620         azimuth(old.azimuth),
00621         name(old.name),
00622         dataformat(old.dataformat)
00623 {
00624         Update();
00625 }
00626 CMTStation& CMTStation::operator= (const CMTStation& source)
00627 {
00628         if (this == &source) return *this;
00629         
00630         copy(source.MTData.begin(),source.MTData.end(),back_inserter(this->MTData));
00631         copy(source.TFData.begin(),source.TFData.end(),back_inserter(this->TFData));
00632 
00633         
00634         this->latitude = source.latitude;
00635         this->longitude = source.longitude;
00636         this->elevation = source.elevation;
00637         this->azimuth = source.azimuth;
00638         this->name = source.name;
00639         this->dataformat = source.dataformat;
00640         this->Update();
00641         return *this;
00642 }