Lemma is an Electromagnetics API
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.

Hantenna.cpp 7.6KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234
  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. using namespace Lemma;
  28. std::vector<Real> readinpfile(const std::string& fname);
  29. std::vector<std::string> readinpfile2(const std::string& fname);
  30. int main(int argc, char** argv) {
  31. std::cout <<
  32. "\n"
  33. << "Hantenna \n\n
  34. << "Hantenna is a programme for computing the H field from polygonal wire\n"
  35. << "loop sources \n\n"
  36. << "Hantenna was built using Lemma (Lemma is an Electromagnetics Modelling API)\n"
  37. << "Lemma is Free and Open Source Software (FOSS) and is released under\n"
  38. << "the MPL, it is covered by the following copyrights:\n"
  39. << "Copyright (C) 2009, 2010, 2011, 2012, 218 Trevor P. Irons\n"
  40. << "Copyright (C) 2011, 2012 M. Andy Kass\n\n"
  41. << "More information may be found at: https://lemmasoftware.org\n"
  42. << " info@lemmasoftware.org\n\n"
  43. << "=====================================================================\n"
  44. << "This programme is part of Lemma, a geophysical modelling and inversion API \n"
  45. << "This Source Code Form is subject to the terms of the Mozilla Public\n"
  46. << "License, v. 2.0. If a copy of the MPL was not distributed with this\n"
  47. << "file, You can obtain one at http://mozilla.org/MPL/2.0/. \n"
  48. << "=====================================================================\n\n\n";
  49. if (argc < 5) {
  50. std::cout << "usage: hantenna.exe trans.inp cond.inp points.inp config.inp \n";
  51. exit(0);
  52. }
  53. #ifdef LEMMAUSEOMP
  54. std::cout << "OpenMP is using " << omp_get_max_threads() << " threads" << std::endl;
  55. #endif
  56. std::vector<Real> Trans = readinpfile(std::string(argv[1]));
  57. std::vector<Real> CondMod = readinpfile(std::string(argv[2]));
  58. std::vector<Real> Points = readinpfile(std::string(argv[3]));
  59. std::vector<std::string> config = readinpfile2(std::string(argv[4]));
  60. //////////////////////////////////////
  61. // Define transmitter
  62. auto trans = PolygonalWireAntenna::NewSP();
  63. trans->SetNumberOfPoints((int)(Trans[0]));
  64. int ip=1;
  65. for ( ; ip<=(int)(Trans[0])*2; ip+=2) {
  66. trans->SetPoint(ip/2, Vector3r (Trans[ip], Trans[ip+1], -1e-3));
  67. }
  68. trans->SetNumberOfFrequencies(1);
  69. trans->SetFrequency(0, Trans[ip]);
  70. trans->SetCurrent(Trans[ip+1]);
  71. trans->SetMinDipoleRatio(atof(config[1].c_str()));
  72. trans->SetMinDipoleMoment(atof(config[2].c_str()));
  73. trans->SetMaxDipoleMoment(atof(config[3].c_str()));
  74. /////////////////////////////////////
  75. // Field calculations
  76. auto receivers = FieldPoints::NewSP();
  77. int nx = (int)Points[0];
  78. int ny = (int)Points[1];
  79. int nz = (int)Points[2];
  80. Real ox = Points[3];
  81. Real oy = Points[4];
  82. Real oz = Points[5];
  83. Vector3r loc;
  84. VectorXr dx(nx-1);
  85. VectorXr dy(ny-1);
  86. VectorXr dz(nz-1);
  87. ip = 6;
  88. int ir = 0;
  89. for ( ; ip <6+nx-1; ++ip) {
  90. dx[ir] = Points[ip];
  91. ++ir;
  92. }
  93. ir = 0;
  94. for ( ; ip <6+ny-1+nx-1; ++ip) {
  95. dy[ir] = Points[ip];
  96. ++ir;
  97. }
  98. ir = 0;
  99. for ( ; ip <6+nz-1+ny-1+nx-1; ++ip) {
  100. dz[ir] = Points[ip];
  101. ++ir;
  102. }
  103. receivers->SetNumberOfPoints(nx*ny*nz);
  104. ir = 0;
  105. Real pz = oz;
  106. for (int iz=0; iz<nz; ++iz) {
  107. Real py = oy;
  108. for (int iy=0; iy<ny; ++iy) {
  109. Real px = ox;
  110. for (int ix=0; ix<nx; ++ix) {
  111. loc << px, py, pz;
  112. receivers->SetLocation(ir, loc);
  113. if (ix < nx-1) px += dx[ix];
  114. ++ ir;
  115. }
  116. if (iy<ny-1) py += dy[iy];
  117. }
  118. if (iz<nz-1) pz += dz[iz];
  119. }
  120. ////////////////////////////////////
  121. // Define model
  122. auto earth = LayeredEarthEM::NewSP();
  123. VectorXcr sigma;
  124. VectorXr thick;
  125. earth->SetNumberOfLayers(CondMod[0]+1);
  126. sigma.resize(CondMod[0]+1); sigma(0) = 0; // airlayer
  127. thick.resize(CondMod[0]-1);
  128. int ilay=1;
  129. for ( ; ilay/2<CondMod[0]-1; ilay+=2) {
  130. sigma(ilay/2+1) = 1./CondMod[ilay];
  131. thick(ilay/2) = CondMod[ilay+1];
  132. }
  133. sigma(ilay/2+1) = 1./ CondMod[ilay];
  134. earth->SetLayerConductivity(sigma);
  135. if (thick.size() > 0) earth->SetLayerThickness(thick);
  136. auto EmEarth = EMEarth1D::NewSP();
  137. EmEarth->AttachWireAntenna(trans);
  138. EmEarth->AttachLayeredEarthEM(earth);
  139. EmEarth->AttachFieldPoints(receivers);
  140. EmEarth->SetFieldsToCalculate(H);
  141. EmEarth->SetHankelTransformMethod(string2Enum<HANKELTRANSFORMTYPE>(config[0]));
  142. EmEarth->SetHankelTransformMethod(FHTKEY201);
  143. ///////////////////////////////////////////////
  144. // Keep track of time
  145. jsw_timer timer;
  146. timer.begin();
  147. clock_t launch = clock();
  148. EmEarth->CalculateWireAntennaFields(true); // true=status bar
  149. Real paTime = timer.end();
  150. std::cout << "\n\n===========================================\ncalc. real time: " << paTime/60. << "\t[m]\n";
  151. std::cout << "calc. user time: " << (clock()-launch)/CLOCKS_PER_SEC/60. << "\t[CPU m]"
  152. << std::endl;
  153. ////////////////////////////////////
  154. // Report
  155. std::fstream hrep("hfield.yaml", std::ios::out);
  156. std::fstream hreal("hfield.dat", std::ios::out);
  157. hrep << *EmEarth << std::endl;
  158. hrep.close();
  159. //hreal << *trans << std::endl;
  160. //hreal << *earth << std::endl;
  161. hreal << "// Right hand coordinate system, z is positive down\n";
  162. hreal << "// x[m]\ty[m]\tz[m]\tHx[A/m]\tHy[A/m]\tHz[A/m]\n";
  163. hreal.precision(8);
  164. int i=0;
  165. for (int iz=0; iz<nz; ++iz) {
  166. for (int iy=0; iy<ny; ++iy) {
  167. for (int ix=0; ix<nx; ++ix) {
  168. hreal << receivers->GetLocation(i).transpose() << "\t";
  169. hreal << receivers->GetHfield(0, i).transpose() << "\n";
  170. ++i;
  171. }
  172. }
  173. }
  174. hreal.close();
  175. // Clean up
  176. }
  177. std::vector<Real> readinpfile(const std::string& fname) {
  178. std::string buf;
  179. char dump[255];
  180. std::vector<Real> vals;
  181. std::fstream input(fname.c_str(), std::ios::in);
  182. if (input.fail()) {
  183. std::cerr << "Input file " << fname << " failed to open\n";
  184. exit(EXIT_FAILURE);
  185. }
  186. while (input >> buf) {
  187. if (buf.substr(0,2) == "//") {
  188. input.getline(dump, 255);
  189. } else {
  190. vals.push_back( atof(buf.c_str() ));
  191. }
  192. }
  193. return vals;
  194. }
  195. std::vector<std::string> readinpfile2(const std::string& fname) {
  196. std::string buf;
  197. char dump[255];
  198. std::vector<std::string> vals;
  199. std::fstream input(fname.c_str(), std::ios::in);
  200. if (input.fail()) {
  201. std::cerr << "Input file " << fname << " failed to open\n";
  202. exit(EXIT_FAILURE);
  203. }
  204. while (input >> buf) {
  205. if (buf.substr(0,2) == "//") {
  206. input.getline(dump, 255);
  207. } else {
  208. vals.push_back( std::string(buf.c_str() ));
  209. }
  210. }
  211. return vals;
  212. }