Main Lemma Repository
Du kannst nicht mehr als 25 Themen auswählen Themen müssen mit entweder einem Buchstaben oder einer Ziffer beginnen. Sie können Bindestriche („-“) enthalten und bis zu 35 Zeichen lang sein.

edipole.cpp 7.7KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231
  1. // ===========================================================================
  2. //
  3. // Filename: utsnmrinversion1d.cpp
  4. //
  5. // Created: 10/07/2010 08:57:04 AM
  6. // Compiler: Tested with g++, icpc, and MSVC 2010
  7. //
  8. // Author: Trevor Irons (ti)
  9. //
  10. // Organisation: Colorado School of Mines (CSM)
  11. // United States Geological Survey (USGS)
  12. //
  13. // Email: tirons@mines.edu, tirons@usgs.gov
  14. //
  15. // This program is free software: you can redistribute it and/or modify
  16. // it under the terms of the GNU General Public License as published by
  17. // the Free Software Foundation, either version 3 of the License, or
  18. // (at your option) any later version.
  19. //
  20. // This program is distributed in the hope that it will be useful,
  21. // but WITHOUT ANY WARRANTY; without even the implied warranty of
  22. // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  23. // GNU General Public License for more details.
  24. //
  25. // You should have received a copy of the GNU General Public License
  26. // along with this program. If not, see <http://www.gnu.org/licenses/>.
  27. //
  28. // ===========================================================================
  29. /**
  30. @file
  31. @author Trevor Irons
  32. @date 10/07/2010
  33. @version 0.0
  34. **/
  35. #include "receiverpoints.h"
  36. #include "emearth1d.h"
  37. #include "dipolesource.h"
  38. using namespace Lemma;
  39. std::vector<Real> readinpfile(const std::string& fname);
  40. int main(int argc, char** argv) {
  41. std::cout <<
  42. "\n"
  43. << "ediple.exe Copyright (C) 2011 Broken Spoke Development, LLC\n\n"
  44. << "This program is free software: you can redistribute it and/or modify\n"
  45. << "it under the terms of the GNU General Public License as published by\n"
  46. << "the Free Software Foundation, either version 3 of the License, or\n"
  47. << "(at your option) any later version.\n\n"
  48. << "This program is distributed in the hope that it will be useful,\n"
  49. << "but WITHOUT ANY WARRANTY; without even the implied warranty of\n"
  50. << "MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the\n"
  51. << "GNU General Public License for more details.\n\n"
  52. << "You should have received a copy of the GNU General Public License\n"
  53. << "along with this program. If not, see <http://www.gnu.org/licenses/>.\n\n";
  54. if (argc < 4) {
  55. std::cout << "usage: edipole.exe dipole.inp cond.inp points.inp \n";
  56. exit(0);
  57. }
  58. std::vector<Real> Trans = readinpfile(std::string(argv[1]));
  59. std::vector<Real> CondMod = readinpfile(std::string(argv[2]));
  60. std::vector<Real> Points = readinpfile(std::string(argv[3]));
  61. //////////////////////////////////////
  62. // Define dipole source
  63. DipoleSource* dipole = DipoleSource::New();
  64. switch ( (int)(Trans[0]) ) {
  65. case 0:
  66. dipole->SetType(MAGNETICDIPOLE);
  67. break;
  68. case 1:
  69. dipole->SetType(GROUNDEDELECTRICDIPOLE);
  70. break;
  71. default:
  72. std::cerr << "Undefined dipole type, must be:\n"
  73. << "\t0 -> Magnetic\n"
  74. << "\t1 -> Grounded Electric\n";
  75. exit(EXIT_FAILURE);
  76. }
  77. switch ( (int)(Trans[1]) ) {
  78. case 0:
  79. dipole->SetPolarisation(XPOLARISATION);
  80. break;
  81. case 1:
  82. dipole->SetPolarisation(YPOLARISATION);
  83. break;
  84. case 2:
  85. dipole->SetPolarisation(ZPOLARISATION);
  86. break;
  87. default:
  88. std::cerr << "Undefined orientation, must be:\n"
  89. << "\t0 -> x\n"
  90. << "\t1 -> y\n"
  91. << "\t2 -> z\n";
  92. exit(EXIT_FAILURE);
  93. }
  94. dipole->SetNumberOfFrequencies(1);
  95. dipole->SetFrequency(0, Trans[2]);
  96. dipole->SetLocation(Trans[3], Trans[4], Trans[5]);
  97. dipole->SetMoment(Trans[6]);
  98. // Receivers
  99. ReceiverPoints *receivers = ReceiverPoints::New();
  100. int nx = (int)Points[0];
  101. int ny = (int)Points[1];
  102. int nz = (int)Points[2];
  103. Real ox = Points[3];
  104. Real oy = Points[4];
  105. Real oz = Points[5];
  106. Vector3r loc;
  107. VectorXr dx(nx-1); // TODO map the dx, dy, dz vectors
  108. VectorXr dy(ny-1);
  109. VectorXr dz(nz-1);
  110. int ip = 6;
  111. int ir = 0;
  112. for ( ; ip <6+nx-1; ++ip) {
  113. dx[ir] = Points[ip];
  114. ++ir;
  115. }
  116. ir = 0;
  117. for ( ; ip <6+ny-1+nx-1; ++ip) {
  118. dy[ir] = Points[ip];
  119. ++ir;
  120. }
  121. ir = 0;
  122. for ( ; ip <6+nz-1+ny-1+nx-1; ++ip) {
  123. dz[ir] = Points[ip];
  124. ++ir;
  125. }
  126. receivers->SetNumberOfReceivers(nx*ny*nz);
  127. ir = 0;
  128. Real pz = oz;
  129. for (int iz=0; iz<nz; ++iz) {
  130. Real py = oy;
  131. for (int iy=0; iy<ny; ++iy) {
  132. Real px = ox;
  133. for (int ix=0; ix<nx; ++ix) {
  134. loc << px, py, pz;
  135. receivers->SetLocation(ir, loc);
  136. if (ix < nx-1) px += dx[ix];
  137. ++ ir;
  138. }
  139. if (iy<ny-1) py += dy[iy];
  140. }
  141. if (iz<nz-1) pz += dz[iz];
  142. }
  143. ////////////////////////////////////
  144. // Define model
  145. LayeredEarthEM *earth = LayeredEarthEM::New();
  146. VectorXcr sigma;
  147. VectorXr thick;
  148. earth->SetNumberOfLayers(CondMod[0]+1);
  149. sigma.resize(CondMod[0]+1); sigma(0) = 0; // airlayer
  150. thick.resize(CondMod[0]-1);
  151. int ilay=1;
  152. for ( ; ilay/2<CondMod[0]-1; ilay+=2) {
  153. sigma(ilay/2+1) = 1./CondMod[ilay];
  154. thick(ilay/2) = CondMod[ilay+1];
  155. }
  156. sigma(ilay/2+1) = 1./ CondMod[ilay];
  157. earth->SetLayerConductivity(sigma);
  158. if (thick.size() > 0) earth->SetLayerThickness(thick);
  159. EMEarth1D *EmEarth = EMEarth1D::New();
  160. EmEarth->AttachDipoleSource(dipole);
  161. EmEarth->AttachLayeredEarthEM(earth);
  162. EmEarth->AttachReceiverPoints(receivers);
  163. EmEarth->SetFieldsToCalculate(E);
  164. EmEarth->MakeCalc3();
  165. ////////////////////////////////////
  166. // Report
  167. std::fstream hreal("efield.dat", std::ios::out);
  168. hreal << *dipole << std::endl;
  169. hreal << *earth << std::endl;
  170. hreal << "Right hand coordinate system, z is positive down\n\n";
  171. hreal << "x[m]\ty[m]\tz[m]\tRe(Ex) [V/m]\t Im(Ex) [V/m]\tRe(Ey) [V/m]\tIm(Ey) [V/m]\tRe(Ez) [V/m] Im(Ez) [V/m]\n";
  172. hreal << "\\\\ ====================================================================================================\n";
  173. hreal.precision(8);
  174. int i=0;
  175. for (int iz=0; iz<nz; ++iz) {
  176. for (int iy=0; iy<ny; ++iy) {
  177. for (int ix=0; ix<nx; ++ix) {
  178. hreal << receivers->GetLocation(i).transpose() << "\t";
  179. hreal << std::real(receivers->GetEfield(0, i).transpose()[0]) << "\t";
  180. hreal << std::imag(receivers->GetEfield(0, i).transpose()[0]) << "\t";
  181. hreal << std::real(receivers->GetEfield(0, i).transpose()[1]) << "\t";
  182. hreal << std::imag(receivers->GetEfield(0, i).transpose()[1]) << "\t";
  183. hreal << std::real(receivers->GetEfield(0, i).transpose()[2]) << "\t";
  184. hreal << std::imag(receivers->GetEfield(0, i).transpose()[2]) << "\n";
  185. ++i;
  186. }
  187. }
  188. }
  189. hreal.close();
  190. // Clean up
  191. EmEarth->Delete();
  192. earth->Delete();
  193. receivers->Delete();
  194. dipole->Delete();
  195. }
  196. std::vector<Real> readinpfile(const std::string& fname) {
  197. std::string buf;
  198. char dump[255];
  199. std::vector<Real> vals;
  200. std::fstream input(fname.c_str(), std::ios::in);
  201. if (input.fail()) {
  202. std::cerr << "Input file " << fname << " failed to open\n";
  203. exit(EXIT_FAILURE);
  204. }
  205. while (input >> buf) {
  206. if (buf.substr(0,2) == "//") {
  207. input.getline(dump, 255);
  208. } else {
  209. vals.push_back( atof(buf.c_str() ));
  210. }
  211. }
  212. return vals;
  213. }