Main Lemma Repository
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

modelreadertem1dubc.cpp 11KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389
  1. /* This file is part of Lemma, a geophysical modelling and inversion API */
  2. /* This Source Code Form is subject to the terms of the Mozilla Public
  3. * License, v. 2.0. If a copy of the MPL was not distributed with this
  4. * file, You can obtain one at http://mozilla.org/MPL/2.0/. */
  5. /**
  6. @file
  7. @author M. Andy Kass
  8. @date 02/19/2011
  9. @version $Id: modelreadertem1dubc.cpp 153 2014-03-21 23:54:38Z tirons $
  10. **/
  11. #include "modelreadertem1dubc.h"
  12. namespace Lemma {
  13. std::ostream &operator<<(std::ostream &stream,
  14. const ModelReaderTem1DUBC &ob) {
  15. stream << *(ModelReaderTem1DUBC*)(&ob);
  16. return stream;
  17. }
  18. ModelReaderTem1DUBC::ModelReaderTem1DUBC(const std::string &name) :
  19. ModelReaderTem1D(name), Trans(NULL), Receivers(NULL),
  20. EarthMod(NULL) {
  21. }
  22. ModelReaderTem1DUBC::~ModelReaderTem1DUBC() {
  23. if (NumberOfReferences != 0) {
  24. throw DeleteObjectWithReferences(this);
  25. }
  26. }
  27. ModelReaderTem1DUBC* ModelReaderTem1DUBC::New() {
  28. ModelReaderTem1DUBC* Obj = new
  29. ModelReaderTem1DUBC("ModelReaderTem1DUBC");
  30. Obj->AttachTo(Obj);
  31. return Obj;
  32. }
  33. void ModelReaderTem1DUBC::Delete() {
  34. this->DetachFrom(this);
  35. }
  36. void ModelReaderTem1DUBC::Release() {
  37. delete this;
  38. }
  39. void ModelReaderTem1DUBC::SetEMEarth1D(LayeredEarthEM* earth) {
  40. if (this->EarthMod != NULL) {
  41. this->EarthMod->DetachFrom(this);
  42. }
  43. earth->AttachTo(this);
  44. this->EarthMod=earth;
  45. }
  46. void ModelReaderTem1DUBC::SetTransmitter(WireAntenna* antennae) {
  47. if (this->Trans != NULL) {
  48. this->Trans->DetachFrom(this);
  49. }
  50. antennae->AttachTo(this);
  51. this->Trans=antennae;
  52. }
  53. void ModelReaderTem1DUBC::SetReceiver(ReceiverPoints* receiver) {
  54. if (this->Receivers != NULL) {
  55. this->Receivers->DetachFrom(this);
  56. }
  57. receiver->AttachTo(this);
  58. this->Receivers=receiver;
  59. }
  60. void ModelReaderTem1DUBC::ReadParameters(const std::string &params) {
  61. std::string temp;
  62. std::string temp2;
  63. char* commentBuffer = new char[200];
  64. int bufsize;
  65. bufsize = 200;
  66. std::fstream infile(params.c_str(),std::ios::in);
  67. if (infile.fail()) {
  68. std::cout << "Parameter file I/O error." << std::endl;
  69. }
  70. infile.getline(commentBuffer,bufsize);
  71. temp = commentBuffer;
  72. std::istringstream in(temp);
  73. in >> this->instrumentfile;
  74. in.clear();
  75. //testing only
  76. //std::cout << this->instrumentfile << std::endl;
  77. infile.getline(commentBuffer,bufsize);
  78. temp = " ";
  79. temp = commentBuffer;
  80. std::istringstream in2(temp);
  81. in2 >> this->modelfile;
  82. in2.clear();
  83. //testing only
  84. //std::cout << this->modelfile << std::endl;
  85. infile.getline(commentBuffer,bufsize);
  86. infile.getline(commentBuffer,bufsize);
  87. infile.getline(commentBuffer,bufsize);
  88. temp = " ";
  89. temp = commentBuffer;
  90. std::istringstream in3(temp);
  91. in3 >> temp2;
  92. if (temp2 == "y" || temp2 == "Y" || temp2 == "yes" || temp2 == "Yes" ||
  93. temp2 == "YES" || temp2 == "true" || temp2 == "True" ||
  94. temp2 == "TRUE") {
  95. this->addednoise = true;
  96. }
  97. else {
  98. this->addednoise = false;
  99. }
  100. in3.clear();
  101. //testing only
  102. //std::cout << this->addednoise << std::endl;
  103. if (this->addednoise == true) {
  104. infile.getline(commentBuffer,bufsize);
  105. temp = " ";
  106. temp = commentBuffer;
  107. std::istringstream in4(temp);
  108. in4 >> this->noisepercent;
  109. in4 >> this->noisethresh;
  110. in4 >> this->noiseseed;
  111. in4.clear();
  112. //std::cout << this->noisepercent << " " << this->noisethresh <<
  113. // " " << this->noiseseed << std::endl;
  114. }
  115. infile.close();
  116. delete [] commentBuffer;
  117. // Call the other readers from here.
  118. this->ReadModel();
  119. this->ReadInstrument();
  120. }
  121. void ModelReaderTem1DUBC::ReadModel() {
  122. // Model File cannot have comments!
  123. int numlayers; //includes both halfspace and air layers
  124. VectorXr thicknesses;
  125. VectorXr conductivities;
  126. VectorXcr complexconduct;
  127. std::string temp;
  128. Real temp2;
  129. std::fstream infile(this->modelfile.c_str(),std::ios::in);
  130. if (infile.fail()) {
  131. std::cout << "Model File I/O error." << std::endl;
  132. }
  133. infile >> numlayers;
  134. numlayers = numlayers + 1;
  135. thicknesses.resize(numlayers-2);
  136. conductivities.resize(numlayers);
  137. conductivities(0) = 0.;
  138. //Testing
  139. //std::cout << "Numlayers " << numlayers << std::endl;
  140. //std::cout << conductivities(0) << std::endl;
  141. for (int ii=0;ii<numlayers-2;ii++) {
  142. infile >> thicknesses(ii);
  143. infile >> conductivities(ii+1);
  144. //Testing
  145. //std::cout<<ii<<" "<<thicknesses(ii)<<" "<<conductivities(ii+1)
  146. // << std::endl;
  147. }
  148. infile >> temp2;
  149. //Testing
  150. //std::cout << temp2 << std::endl;
  151. infile >> conductivities(numlayers-1);
  152. //Testing
  153. //std::cout << conductivities(numlayers-1) << std::endl;
  154. //std::cout << "Thicknesses: " << thicknesses << std::endl;
  155. //std::cout << "Conductivities: " << conductivities << std::endl;
  156. //Didn't read the function declaration. Conductivities needs
  157. //to be complex...
  158. complexconduct.resize(conductivities.size());
  159. complexconduct.setZero();
  160. complexconduct.real() = conductivities.array();
  161. //Testing
  162. //std::cout << complexconduct << std::endl;
  163. // Attach all to layered earth class
  164. this->EarthMod->SetNumberOfLayers(numlayers);
  165. this->EarthMod->SetLayerConductivity(complexconduct);
  166. this->EarthMod->SetLayerThickness(thicknesses);
  167. infile.close();
  168. }
  169. void ModelReaderTem1DUBC::ReadInstrument() {
  170. // No comments allowed in this data file either for now.
  171. int tempint;
  172. Real tempreal;
  173. VectorXr transpoints;
  174. MatrixXr receiverpoints;
  175. std::string tempstring;
  176. VectorXr tempvec;
  177. int timeunits;
  178. int ntimes;
  179. int ntimestemp;
  180. int numsounds;
  181. VectorXr orientation;
  182. MatrixXr timestemp;
  183. //Testing
  184. //std::cout << "ModelReaderTem1DUBC::ReadInstrument()" << std::endl;
  185. std::fstream infile(this->instrumentfile.c_str(),std::ios::in);
  186. if (infile.fail()) {
  187. std::cout << "Observations File I/O error." << std::endl;
  188. }
  189. infile >> tempint;
  190. this->Trans->SetNumberOfPoints(tempint+1);
  191. //Testing
  192. //std::cout << "Number of Points: " << tempint << std::endl;
  193. transpoints.resize(tempint*2+2);
  194. //Testing
  195. //std::cout << "Transmitter corners: ";
  196. for (int ii=0;ii<(tempint*2);ii++) {
  197. infile >> transpoints(ii);
  198. //Testing
  199. //std::cout << transpoints(ii) << " ";
  200. }
  201. transpoints(tempint*2) = transpoints(0);
  202. transpoints(tempint*2+1) = transpoints(1);
  203. //Testing
  204. //std::cout << transpoints(tempint*2) << " " <<
  205. // transpoints(tempint*2+1) << std::endl;
  206. infile >> tempreal;
  207. //Testing
  208. //std::cout << "Transmitter height: " << tempreal << std::endl;
  209. int jj;
  210. jj = 0;
  211. for (int ii=0;ii<(tempint+1);ii++) {
  212. this->Trans->SetPoint(ii,Vector3r(transpoints(jj),
  213. transpoints(jj+1),tempreal));
  214. //Testing
  215. //std::cout<< transpoints(jj) << transpoints(jj+1) << std::endl;
  216. jj = jj+2;
  217. }
  218. //Testing
  219. infile >> this->waveformfile;
  220. //Testing
  221. //std::cout << "Waveform file: " << this->waveformfile << std::endl;
  222. infile >> numsounds; // number of soundings
  223. //Testing
  224. //std::cout << "Number of soundings: " << numsounds << std::endl;
  225. this->Trans->SetCurrent(1);
  226. this->Trans->SetNumberOfTurns(1);
  227. infile >> timeunits; // flag for units of time
  228. // 1 = microseconds, 2 = milliseconds, 3 = seconds
  229. //Testing
  230. //std::cout << "Time units: " << timeunits << std::endl;
  231. receiverpoints.resize(numsounds,3);
  232. //Testing
  233. //std::cout << "Receiverpoints resized" << std::endl;
  234. this->receivermoment.resize(numsounds);
  235. //Testing
  236. //std::cout << "Receivermoment resized" << std::endl;
  237. orientation.resize(numsounds);
  238. //Testing
  239. //std::cout << "Orientation resized" << std::endl;
  240. ntimestemp = 30;
  241. timestemp.resize(numsounds,ntimestemp);
  242. //Testing
  243. //std::cout << "timestemp resized" << std::endl;
  244. this->Receivers->SetNumberOfReceivers(numsounds);
  245. //Testing
  246. //std::cout << "Reading receiver points" << std::endl;
  247. for (int ii=0;ii<numsounds;ii++) {
  248. //Testing
  249. //std::cout << "Beginning iteration: " << ii << std::endl;
  250. infile >> this->receivermoment(ii); //receiver moment
  251. //Testing
  252. //std::cout << "Receiver moment: " << this->receivermoment(ii)
  253. // << std::endl;
  254. infile >> receiverpoints(ii,0);
  255. infile >> receiverpoints(ii,1);
  256. infile >> receiverpoints(ii,2);
  257. //Testing
  258. //std::cout << "Receiverpoints: " << receiverpoints.row(ii)
  259. // << std::endl;
  260. this->Receivers->SetLocation(ii,receiverpoints.row(ii));
  261. infile >> tempstring; //receiver orientation
  262. //Testing
  263. //std::cout << "Receiver orientation: " << tempstring << std::endl;
  264. if (tempstring == "z") {
  265. orientation(ii) = 0;
  266. }
  267. if (tempstring == "x") {
  268. orientation(ii) = 1;
  269. }
  270. if (tempstring == "y") {
  271. orientation(ii) = 2;
  272. }
  273. else {
  274. orientation(ii) = 0;
  275. }
  276. //Testing
  277. //std::cout << "Modified orientation: " << orientation(ii)
  278. // << std::endl;
  279. infile >> ntimes;
  280. //Testing
  281. //std::cout << "Number of gates: " << ntimes << std::endl;
  282. if (ntimes != ntimestemp) {
  283. timestemp.resize(Eigen::NoChange,ntimes);
  284. ntimestemp = ntimes;
  285. //Testing
  286. // std::cout << "Better not be here!" << std::endl;
  287. }
  288. //Testing
  289. //std::cout << "Actual size: " << timestemp.size() << std::endl;
  290. infile >> tempint;
  291. //Testing
  292. //std::cout << "Unknown variable: " << tempint << std::endl;
  293. //std::cout << "Times: ";
  294. for (int jj=0;jj<ntimes;jj++) {
  295. infile >> timestemp(ii,jj);
  296. //Testing
  297. std::cout << "time read? " << jj << ". " << timestemp(ii,jj) << "\t " << tempint << std::endl;
  298. infile >> tempint;
  299. }
  300. //Testing
  301. //std::cout << std::endl;
  302. }
  303. infile.close();
  304. //Testing
  305. //std::cout << "Left the loop" << std::endl;
  306. // Make sure times are in seconds
  307. if (timeunits == 1) {
  308. timestemp = timestemp.array()/1e6;
  309. timeunits = 3;
  310. }
  311. if (timeunits == 2) {
  312. timestemp = timestemp.array()/1e3;
  313. timeunits = 3;
  314. }
  315. //Testing
  316. //std::cout << "Scaled times: " << timestemp << std::endl;
  317. this->modtimes.resize(ntimes);
  318. //Testing
  319. //std::cout << "Modtimes size: " << this->modtimes.size() << std::endl;
  320. //std::cout << "timestemp size: " << timestemp.rows() << " x " <<
  321. // timestemp.cols() << std::endl;
  322. tempvec = timestemp.row(0);
  323. this->modtimes = tempvec.array();
  324. //Testing
  325. //std::cout << "Finished with loading modtimes" << std::endl;
  326. // This code has been written so far to allow for the general case
  327. // of soundings with different vectors of times. However, the
  328. // instrument class has not been defined that way yet. So right now
  329. // only one vector of times is returned via an accessor. This will
  330. // be changed later.
  331. }
  332. void ModelReaderTem1DUBC::ReadWaveform() {
  333. }
  334. VectorXr ModelReaderTem1DUBC::GetTimes() {
  335. //Testing
  336. //std::cout << this->modtimes;
  337. return this->modtimes;
  338. }
  339. //Datamembers:
  340. // modelfile, instrumentfile, waveformfile
  341. } // namespace Lemma