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