C1DMTSynthData.cpp

Go to the documentation of this file.
00001 #include "C1DMTSynthData.h"
00002 #include <fstream>
00003 #include <algorithm>
00004 #include <cassert>
00005 #include "FatalException.h"
00006 
00007 using namespace std;
00008 
00009 namespace gplib
00010   {
00011     const int frequenzen = 50;
00012     //these are the frequencies we calculate the response at, if no frequencies are set before the calculation
00013     const float T[frequenzen] =
00014       { 0.0025, 0.00308337, 0.00513901, 0.00625, 0.0104167, 0.0125, 0.0208333,
00015           0.025, 0.0416667, 0.05, 0.0833333, 0.1, 0.166667, 0.2, 0.333333, 0.4,
00016           0.667, 0.64, 1.06667, 1.29333, 2.15554, 2.59997, 4.33332, 5.19994,
00017           8.66701, 10.666667, 16.000, 21.333, 32.000, 42.667, 64.0, 85.332,
00018           127.992, 170.678, 256.016, 341.297, 512.033, 682.594, 1023.541,
00019           1366.120, 2049.180, 2732.240, 4098.361, 5464.481, 8196.722,
00020           10869.565, 16393.443, 28800, 43200, 86400 };
00021 
00022     C1DMTSynthData::C1DMTSynthData()
00023       {
00024       }
00025 
00026     C1DMTSynthData::~C1DMTSynthData()
00027       {
00028       }
00029 
00030     C1DMTSynthData::C1DMTSynthData(const C1DMTSynthData &old) :
00031       MTStation(old), Z(old.Z), calc_frequencies(old.calc_frequencies),
00032           resistivity(old.resistivity), thickness(old.thickness), reserrors(
00033               old.reserrors), thickerrors(old.thickerrors)
00034       {
00035 
00036       }
00037     void C1DMTSynthData::WritePlot(std::string filename)
00038       {
00039         ofstream outfile(filename.c_str());
00040         double cumthick = 0;
00041         assert(resistivity.size() == thickness.size());
00042 
00043         for (unsigned int i = 0; i < thickness.size(); ++i)
00044           {
00045             outfile << cumthick << "  " << resistivity.at(i);
00046             if (!thickerrors.empty() || !reserrors.empty())
00047               outfile << " " << thickerrors.at(i) << " " << reserrors.at(i);
00048             outfile << endl;
00049             cumthick += thickness.at(i);
00050             outfile << cumthick << " " << resistivity.at(i);
00051             if (!thickerrors.empty() || !reserrors.empty())
00052               outfile << " " << thickerrors.at(i) << " " << reserrors.at(i);
00053             outfile << endl;
00054           }
00055       }
00056 
00057     void C1DMTSynthData::WriteModel(string filename)
00058       {
00059         assert(resistivity.size() == thickness.size());
00060         ofstream outfile(filename.c_str());
00061 
00062         outfile << thickness.size() << endl;
00063         for (unsigned int i = 0; i < thickness.size(); ++i)
00064           {
00065             outfile << thickness.at(i) << "  " << resistivity.at(i) << endl;
00066           }
00067       }
00068 
00069     gplib::rvec C1DMTSynthData::GetModelVector()
00070       {
00071         assert(resistivity.size() == thickness.size());
00072         const int length = resistivity.size();
00073         gplib::rvec result(length * 2);
00074         for (int i = 0; i < length; ++i)
00075           {
00076             result(i) = log10(resistivity.at(i));
00077             result(i + length) = thickness.at(i);
00078           }
00079         return result;
00080       }
00081 
00082     void C1DMTSynthData::ReadModel(string filename)
00083       {
00084         int layers;
00085         double currthick, curres;
00086         int i = 0;
00087 
00088         ifstream infile(filename.c_str());
00089         infile >> layers;
00090         if (!infile.good() || layers < 0)
00091           throw FatalException("Problem reading file: " + filename);
00092         thickness.reserve(layers);
00093         resistivity.reserve(layers);
00094         while (infile.good())
00095           {
00096             infile >> currthick >> curres;
00097             if (infile.good())
00098               {
00099                 thickness.push_back(currthick);
00100                 resistivity.push_back(curres);
00101                 ++i;
00102               }
00103           }
00104         assert(resistivity.size() == thickness.size());
00105         assert(i==layers);
00106       }
00107     //fill data member variables with data from Cagniard algorithm
00108     void C1DMTSynthData::CalcSynthetic()
00109 
00110       {
00111         assert(resistivity.size() == thickness.size());
00112         calc_frequencies = GetFrequencies();
00113         if (calc_frequencies.empty())
00114           {
00115             for (int i = 0; i < frequenzen; ++i)
00116               calc_frequencies.push_back(1.0 / T[i]);
00117           }
00118         Calc();
00119         Assign(calc_frequencies.size());
00120         for (unsigned int i = 0; i < calc_frequencies.size(); ++i)
00121           {
00122             MTData.at(i).frequency = calc_frequencies.at(i);
00123             MTData.at(i).Zxy = Z.at(i);
00124             MTData.at(i).Zyx = -Z.at(i);
00125           }
00126         Update();
00127       }
00128 
00129     void C1DMTSynthData::Calc()
00130       {
00131         dcomp omegamu;
00132         dcomp kcurr, klow;
00133         dcomp xi;
00134         vector<dcomp> alpha(thickness.size(), 0);
00135         dcomp adm;
00136         double d;
00137         double sigmacurr, sigmalow;
00138         Z.clear();
00139         for (unsigned int i = 0; i < calc_frequencies.size(); ++i)
00140           {
00141             omegamu = -I * 8e-7 * PI * PI * calc_frequencies.at(i);
00142             fill(alpha.begin(), alpha.end(), 0);
00143             sigmacurr = 1.0 / resistivity.back();
00144             kcurr = sqrt(omegamu * sigmacurr);
00145             for (int j = thickness.size() - 2; j >= 0; --j)
00146               {
00147                 sigmalow = 1.0 / resistivity.at(j + 1);
00148                 sigmacurr = 1.0 / resistivity.at(j);
00149                 kcurr = sqrt(omegamu * sigmacurr);
00150                 klow = sqrt(omegamu * sigmalow);
00151                 if (kcurr.real() < 0.0)
00152                   {
00153                     kcurr *= -1.0;
00154                     klow *= -1.0;
00155                   }
00156                 xi = omegamu * (sigmacurr - sigmalow) / std::pow(kcurr + klow, 2);
00157                 alpha.at(j + 1) = (xi + alpha.at(j + 1)) / (1. + xi * alpha.at(
00158                     j + 1));
00159                 d = thickness.at(j) * 1000.0;
00160                 alpha.at(j) = alpha.at(j + 1) * exp(-2. * kcurr * d);
00161               }
00162             adm = kcurr / (-I * 2. * PI * calc_frequencies.at(i)) * ((1.
00163                 - alpha.at(0)) / (1. + alpha.at(0)));
00164             Z.push_back(conj(1. / (1000.0 * adm)));
00165 
00166           }
00167       }
00168   }

Generated on Tue May 4 16:52:14 2010 for GPLIB++ by  doxygen 1.5.8