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.

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254
  1. // ===========================================================================
  2. //
  3. // Filename: hantenna.cpp
  4. //
  5. // Created: 10/07/2010 08:57:04 AM
  6. // Modified: 11 April 2018
  7. // Compiler: Tested with g++, icpc, and MSVC 2017
  8. //
  9. // Author: Trevor Irons (ti)
  10. //
  11. // Copyright (C) 2012,2018 Trevor Irons
  12. //
  13. // Organisation: Lemma Software
  14. //
  15. // Email: Trevor.Irons@lemmasoftware.org
  16. //
  17. // ===========================================================================
  18. /**
  19. @file
  20. @author Trevor Irons
  21. @date 10/07/2010
  22. $Id$
  23. **/
  24. #include "LemmaCore"
  25. #include "FDEM1D"
  26. #include "timer.h"
  27. #if defined(__clang__)
  28. /* Clang/LLVM. ---------------------------------------------- */
  29. const char* compiler = "clang";
  30. const char* ver = __VERSION__;
  31. #elif defined(__ICC) || defined(__INTEL_COMPILER)
  32. /* Intel ICC/ICPC. ------------------------------------------ */
  33. const char* compiler = "icpc";
  34. #elif defined(__GNUC__) || defined(__GNUG__)
  35. /* GNU GCC/G++. --------------------------------------------- */
  36. const char* compiler = "gcc (GCC) ";// __VERSION__;
  37. const char* ver = __VERSION__;
  38. #elif defined(_MSC_VER)
  39. /* Microsoft Visual Studio. --------------------------------- */
  40. const char* compiler = "msvc ";
  41. const int ver = _MSC_FULL_VER;
  42. #elif defined(__PGI)
  43. /* Portland Group PGCC/PGCPP. ------------------------------- */
  44. const char* compiler = "pgc";
  45. #endif
  46. using namespace Lemma;
  47. std::vector<Real> readinpfile(const std::string& fname);
  48. std::vector<std::string> readinpfile2(const std::string& fname);
  49. int main(int argc, char** argv) {
  50. const char *buildString = __DATE__ ", " __TIME__;
  51. std::cout
  52. << "===========================================================================\n"
  53. << "Lemma " << LEMMA_VERSION << "\n"
  54. << "[" << compiler << " " << ver << " " << buildString << "]\n"
  55. << "This program is part of Lemma, a geophysical modelling and inversion API. \n"
  56. << " This Source Code Form is subject to the terms of the Mozilla Public\n"
  57. << " License, v. 2.0. If a copy of the MPL was not distributed with this\n"
  58. << " file, You can obtain one at http://mozilla.org/MPL/2.0/. \n"
  59. << "Copyright (C) 2018 Lemma Software \n"
  60. << "More information may be found at: https://lemmasoftware.org\n"
  61. << " info@lemmasoftware.org\n"
  62. << "===========================================================================\n\n"
  63. << "Hantenna calculates the harmonic H field from polygonal wire loop sources\n";
  64. if (argc < 5) {
  65. std::cout << "usage: hantenna.exe trans.inp cond.inp points.inp config.inp \n";
  66. exit(0);
  67. }
  68. #ifdef LEMMAUSEOMP
  69. std::cout << "OpenMP is using " << omp_get_max_threads() << " threads" << std::endl;
  70. #endif
  71. std::vector<Real> Trans = readinpfile(std::string(argv[1]));
  72. std::vector<Real> CondMod = readinpfile(std::string(argv[2]));
  73. std::vector<Real> Points = readinpfile(std::string(argv[3]));
  74. std::vector<std::string> config = readinpfile2(std::string(argv[4]));
  75. //////////////////////////////////////
  76. // Define transmitter
  77. auto trans = PolygonalWireAntenna::NewSP();
  78. trans->SetNumberOfPoints((int)(Trans[0]));
  79. int ip=1;
  80. for ( ; ip<=(int)(Trans[0])*2; ip+=2) {
  81. trans->SetPoint(ip/2, Vector3r (Trans[ip], Trans[ip+1], -1e-3));
  82. //trans->SetPoint(ip/2, Vector3r (Trans[ip], Trans[ip+1], 50.));
  83. }
  84. trans->SetNumberOfFrequencies(1);
  85. trans->SetFrequency(0, Trans[ip]);
  86. trans->SetCurrent(Trans[ip+1]);
  87. trans->SetMinDipoleRatio(atof(config[1].c_str()));
  88. trans->SetMinDipoleMoment(atof(config[2].c_str()));
  89. trans->SetMaxDipoleMoment(atof(config[3].c_str()));
  90. /////////////////////////////////////
  91. // Field calculations
  92. auto receivers = FieldPoints::NewSP();
  93. int nx = (int)Points[0];
  94. int ny = (int)Points[1];
  95. int nz = (int)Points[2];
  96. Real ox = Points[3];
  97. Real oy = Points[4];
  98. Real oz = Points[5];
  99. Vector3r loc;
  100. VectorXr dx(nx-1);
  101. VectorXr dy(ny-1);
  102. VectorXr dz(nz-1);
  103. ip = 6;
  104. int ir = 0;
  105. for ( ; ip <6+nx-1; ++ip) {
  106. dx[ir] = Points[ip];
  107. ++ir;
  108. }
  109. ir = 0;
  110. for ( ; ip <6+ny-1+nx-1; ++ip) {
  111. dy[ir] = Points[ip];
  112. ++ir;
  113. }
  114. ir = 0;
  115. for ( ; ip <6+nz-1+ny-1+nx-1; ++ip) {
  116. dz[ir] = Points[ip];
  117. ++ir;
  118. }
  119. receivers->SetNumberOfPoints(nx*ny*nz);
  120. ir = 0;
  121. Real pz = oz;
  122. for (int iz=0; iz<nz; ++iz) {
  123. Real py = oy;
  124. for (int iy=0; iy<ny; ++iy) {
  125. Real px = ox;
  126. for (int ix=0; ix<nx; ++ix) {
  127. loc << px, py, pz;
  128. receivers->SetLocation(ir, loc);
  129. if (ix < nx-1) px += dx[ix];
  130. ++ ir;
  131. }
  132. if (iy<ny-1) py += dy[iy];
  133. }
  134. if (iz<nz-1) pz += dz[iz];
  135. }
  136. ////////////////////////////////////
  137. // Define model
  138. auto earth = LayeredEarthEM::NewSP();
  139. VectorXcr sigma;
  140. VectorXr thick;
  141. earth->SetNumberOfLayers(static_cast<int>(CondMod[0])+1);
  142. sigma.resize(static_cast<int>(CondMod[0])+1); sigma(0) = 0; // airlayer
  143. thick.resize(static_cast<int>(CondMod[0])-1);
  144. int ilay=1;
  145. for ( ; ilay/2<CondMod[0]-1; ilay+=2) {
  146. sigma(ilay/2+1) = 1./CondMod[ilay];
  147. thick(ilay/2) = CondMod[ilay+1];
  148. }
  149. sigma(ilay/2+1) = 1./ CondMod[ilay];
  150. earth->SetLayerConductivity(sigma);
  151. if (thick.size() > 0) earth->SetLayerThickness(thick);
  152. auto EmEarth = EMEarth1D::NewSP();
  153. EmEarth->AttachWireAntenna(trans);
  154. EmEarth->AttachLayeredEarthEM(earth);
  155. EmEarth->AttachFieldPoints(receivers);
  156. EmEarth->SetFieldsToCalculate(H);
  157. EmEarth->SetHankelTransformMethod(string2Enum<HANKELTRANSFORMTYPE>(config[0]));
  158. ///////////////////////////////////////////////
  159. // Keep track of time
  160. jsw_timer timer;
  161. timer.begin();
  162. clock_t launch = clock();
  163. EmEarth->CalculateWireAntennaFields(true); // true=status bar
  164. Real paTime = timer.end();
  165. std::cout << "\n\n===========================================\ncalc. real time: " << paTime/60. << "\t[m]\n";
  166. std::cout << "calc. user time: " << (clock()-launch)/CLOCKS_PER_SEC/60. << "\t[CPU m]"
  167. << std::endl;
  168. ////////////////////////////////////
  169. // Report
  170. std::fstream hrep("hfield.yaml", std::ios::out);
  171. std::fstream hreal("hfield.dat", std::ios::out);
  172. hrep << *EmEarth << std::endl;
  173. hrep.close();
  174. //hreal << *trans << std::endl;
  175. //hreal << *earth << std::endl;
  176. hreal << "// Right hand coordinate system, z is positive down\n";
  177. hreal << "// x[m]\ty[m]\tz[m]\tHx[A/m]\tHy[A/m]\tHz[A/m]\n";
  178. hreal.precision(8);
  179. int i=0;
  180. for (int iz=0; iz<nz; ++iz) {
  181. for (int iy=0; iy<ny; ++iy) {
  182. for (int ix=0; ix<nx; ++ix) {
  183. hreal << receivers->GetLocation(i).transpose() << "\t";
  184. //hreal << receivers->GetHfield(0, i).transpose() << "\n"; ( complex, notation )
  185. hreal << receivers->GetHfield(0, i).transpose().real() << "\t";
  186. hreal << receivers->GetHfield(0, i).transpose().imag() << "\n";
  187. ++i;
  188. }
  189. }
  190. }
  191. hreal.close();
  192. // Clean up
  193. }
  194. std::vector<Real> readinpfile(const std::string& fname) {
  195. std::string buf;
  196. char dump[255];
  197. std::vector<Real> vals;
  198. std::fstream input(fname.c_str(), std::ios::in);
  199. if (input.fail()) {
  200. std::cerr << "Input file " << fname << " failed to open\n";
  201. exit(EXIT_FAILURE);
  202. }
  203. while (input >> buf) {
  204. if (buf.substr(0,2) == "//") {
  205. input.getline(dump, 255);
  206. } else {
  207. vals.push_back( atof(buf.c_str() ));
  208. }
  209. }
  210. return vals;
  211. }
  212. std::vector<std::string> readinpfile2(const std::string& fname) {
  213. std::string buf;
  214. char dump[255];
  215. std::vector<std::string> vals;
  216. std::fstream input(fname.c_str(), std::ios::in);
  217. if (input.fail()) {
  218. std::cerr << "Input file " << fname << " failed to open\n";
  219. exit(EXIT_FAILURE);
  220. }
  221. while (input >> buf) {
  222. if (buf.substr(0,2) == "//") {
  223. input.getline(dump, 255);
  224. } else {
  225. vals.push_back( std::string(buf.c_str() ));
  226. }
  227. }
  228. return vals;
  229. }