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.

LayeredEarthEM.cpp 14KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368
  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 Trevor Irons
  8. @date 06/24/2009
  9. @version $Id: layeredearthem.cpp 216 2015-03-09 02:26:49Z tirons $
  10. **/
  11. #include "LayeredEarthEM.h"
  12. namespace Lemma {
  13. std::ostream &operator << (std::ostream &stream, const LayeredEarthEM &ob) {
  14. stream << ob.Serialize() << "\n";
  15. return stream;
  16. }
  17. // ==================== LIFECYCLE ===================================
  18. LayeredEarthEM::LayeredEarthEM( const ctor_key& key ) : LayeredEarth( key ) {
  19. }
  20. LayeredEarthEM::LayeredEarthEM( const YAML::Node& node, const ctor_key& key ) : LayeredEarth(node, key) {
  21. SetNumberOfLayers(NumberOfLayers);
  22. LayerThickness = node["LayerThickness"].as<VectorXr>();
  23. if (node["LayerConductivity"]) {
  24. LayerConductivity = node["LayerConductivity"].as<VectorXcr>();
  25. }
  26. if (node["LayerSusceptibility"]) {
  27. LayerSusceptibility = node["LayerSusceptibility"].as<VectorXcr>();
  28. }
  29. if (node["LayerLowFreqSusceptibility"]) {
  30. LayerLowFreqSusceptibility = node["LayerLowFreqSusceptibility"].as<VectorXr>();
  31. }
  32. if (node["LayerHighFreqSusceptibility"]) {
  33. LayerHighFreqSusceptibility = node["LayerHighFreqSusceptibility"].as<VectorXr>();
  34. }
  35. if (node["LayerTauSusceptibility"]) {
  36. LayerTauSusceptibility = node["LayerTauSusceptibility"].as<VectorXr>();
  37. }
  38. if (node["LayerBreathSusceptibility"]) {
  39. LayerBreathSusceptibility = node["LayerBreathSusceptibility"].as<VectorXr>();
  40. }
  41. if (node["LayerPermitivity"]) {
  42. LayerPermitivity = node["LayerPermitivity"].as<VectorXcr>();
  43. }
  44. if (node["LayerLowFreqPermitivity "]) {
  45. LayerLowFreqPermitivity = node["LayerLowFreqPermitivity"].as<VectorXr>();
  46. }
  47. if (node["LayerHighFreqPermitivity "]) {
  48. LayerHighFreqPermitivity = node["LayerHighFreqPermitivity"].as<VectorXr>();
  49. }
  50. if (node["LayerTauPermitivity"]) {
  51. LayerTauPermitivity = node["LayerTauPermitivity"].as<VectorXr>();
  52. }
  53. if (node["LayerBreathPermitivity"]) {
  54. LayerBreathPermitivity = node["LayerBreathPermitivity"].as<VectorXr>();
  55. }
  56. }
  57. LayeredEarthEM::~LayeredEarthEM() {
  58. }
  59. std::shared_ptr<LayeredEarthEM> LayeredEarthEM::NewSP() {
  60. return std::make_shared<LayeredEarthEM> ( ctor_key() );
  61. }
  62. std::shared_ptr<LayeredEarthEM> LayeredEarthEM::DeSerialize( const YAML::Node& node ) {
  63. if (node.Tag() != "LayeredEarthEM") {
  64. throw DeSerializeTypeMismatch( "LayeredEarthEM", node.Tag());
  65. }
  66. return std::make_shared<LayeredEarthEM> ( node, ctor_key() );
  67. }
  68. YAML::Node LayeredEarthEM::Serialize() const {
  69. YAML::Node node = LayeredEarth::Serialize();
  70. node.SetTag( GetName() );
  71. node["LayerConductivity"] = LayerConductivity;
  72. node["LayerSusceptibility"] = LayerSusceptibility;
  73. node["LayerLowFreqSusceptibility"] = LayerLowFreqSusceptibility;
  74. node["LayerHighFreqSusceptibility"] = LayerHighFreqSusceptibility;
  75. node["LayerTauSusceptibility"] = LayerTauSusceptibility;
  76. node["LayerBreathSusceptibility"] = LayerBreathSusceptibility;
  77. node["LayerPermitivity"] = LayerPermitivity;
  78. node["LayerLowFreqPermitivity"] = LayerLowFreqPermitivity;
  79. node["LayerHighFreqPermitivity"] = LayerHighFreqPermitivity;
  80. node["LayerTauPermitivity"] = LayerTauPermitivity;
  81. node["LayerBreathPermitivity"] = LayerBreathPermitivity;
  82. return node;
  83. }
  84. //--------------------------------------------------------------------------------------
  85. // Class: LayeredEarthEM
  86. // Method: GetName
  87. // Description: Class identifier
  88. //--------------------------------------------------------------------------------------
  89. inline std::string LayeredEarthEM::GetName ( ) const {
  90. return CName;
  91. } // ----- end of method LayeredEarthEM::GetName -----
  92. // ==================== OPERATIONS ===================================
  93. void LayeredEarthEM::EvaluateColeColeModel(const Real& omega) {
  94. for (int ilay=0; ilay<GetNumberOfLayers(); ++ilay) {
  95. if ( LayerTauSusceptibility(ilay) > 1e-10) {
  96. LayerSusceptibility(ilay) = LayerHighFreqSusceptibility(ilay) + (LayerLowFreqSusceptibility(ilay) -
  97. LayerHighFreqSusceptibility(ilay)) /
  98. ((Real)(1.) + std::pow(Complex(0, omega*
  99. LayerTauSusceptibility(ilay)),
  100. LayerBreathSusceptibility(ilay)));
  101. } else {
  102. LayerSusceptibility(ilay) = LayerLowFreqSusceptibility(ilay);
  103. }
  104. if ( LayerTauPermitivity(ilay) > 1e-10) {
  105. LayerPermitivity(ilay) = LayerHighFreqPermitivity(ilay) +
  106. (LayerLowFreqPermitivity(ilay) -
  107. LayerHighFreqPermitivity(ilay)) /
  108. ((Real)(1.) + std::pow(Complex(0,
  109. omega*LayerTauPermitivity(ilay)),
  110. LayerBreathPermitivity(ilay)));
  111. } else {
  112. LayerPermitivity(ilay) = LayerLowFreqPermitivity(ilay);
  113. }
  114. }
  115. }
  116. std::shared_ptr<LayeredEarthEM> LayeredEarthEM::Clone() {
  117. auto copy = LayeredEarthEM::NewSP();
  118. copy->LayerConductivity = this->LayerConductivity;
  119. copy->LayerSusceptibility = this->LayerSusceptibility;
  120. copy->LayerLowFreqSusceptibility = this->LayerLowFreqSusceptibility;
  121. copy->LayerHighFreqSusceptibility = this->LayerHighFreqSusceptibility;
  122. copy->LayerTauSusceptibility = this->LayerTauSusceptibility;
  123. copy->LayerBreathSusceptibility = this->LayerBreathSusceptibility;
  124. copy->LayerPermitivity = this->LayerPermitivity;
  125. copy->LayerLowFreqPermitivity = this->LayerLowFreqPermitivity;
  126. copy->LayerHighFreqPermitivity = this->LayerHighFreqPermitivity;
  127. copy->LayerTauPermitivity = this->LayerTauPermitivity;
  128. copy->LayerBreathPermitivity = this->LayerBreathPermitivity;
  129. copy->SetNumberOfLayers( this->GetNumberOfLayers() );
  130. copy->NumberOfInterfaces = this->NumberOfInterfaces;
  131. copy->LayerThickness = this->LayerThickness;
  132. return copy;
  133. }
  134. // ==================== ACCESS ==================================
  135. void LayeredEarthEM::SetLayerConductivity(const VectorXcr &sig) {
  136. if (sig.size() != this->GetNumberOfLayers() )
  137. throw EarthModelParametersDoNotMatchNumberOfLayers( );
  138. LayerConductivity = sig;
  139. }
  140. void LayeredEarthEM::SetLayerConductivity(const int& ilay, const Complex &sig) {
  141. if (ilay > this->GetNumberOfLayers() || ilay < 1 )
  142. throw EarthModelParametersDoNotMatchNumberOfLayers( );
  143. LayerConductivity[ilay] = sig;
  144. }
  145. /* // TODO fix layer 0 problem, 1/0 --> infty, plus layer 0 is ignored anyway
  146. void LayeredEarthEM::SetLayerResistivity(const VectorXcr &rhos) {
  147. if (rhos.size() != this->NumberOfLayers )
  148. throw EarthModelParametersDoNotMatchNumberOfLayers( );
  149. LayerConductivity = 1./rhos.array();
  150. }
  151. */
  152. void LayeredEarthEM::SetLayerHighFreqSusceptibility(const VectorXr &sus) {
  153. if (sus.size() != this->GetNumberOfLayers() )
  154. throw EarthModelParametersDoNotMatchNumberOfLayers( );
  155. LayerHighFreqSusceptibility = sus;
  156. }
  157. void LayeredEarthEM::SetLayerLowFreqSusceptibility(const VectorXr &sus) {
  158. if (sus.size() != this->GetNumberOfLayers() )
  159. throw EarthModelParametersDoNotMatchNumberOfLayers( );
  160. LayerLowFreqSusceptibility = sus;
  161. }
  162. void LayeredEarthEM::SetLayerBreathSusceptibility(const VectorXr &sus) {
  163. if (sus.size() != this->GetNumberOfLayers() )
  164. throw EarthModelParametersDoNotMatchNumberOfLayers( );
  165. LayerBreathSusceptibility = sus;
  166. }
  167. void LayeredEarthEM::SetLayerTauSusceptibility(const VectorXr &sus) {
  168. if (sus.size() != this->GetNumberOfLayers() )
  169. throw EarthModelParametersDoNotMatchNumberOfLayers( );
  170. LayerTauSusceptibility = sus;
  171. }
  172. void LayeredEarthEM::SetLayerHighFreqPermitivity(const VectorXr &per) {
  173. if (per.size() != this->GetNumberOfLayers() )
  174. throw EarthModelParametersDoNotMatchNumberOfLayers( );
  175. LayerHighFreqPermitivity = per;
  176. }
  177. void LayeredEarthEM::SetLayerLowFreqPermitivity(const VectorXr &per) {
  178. if (per.size() != this->GetNumberOfLayers() )
  179. throw EarthModelParametersDoNotMatchNumberOfLayers( );
  180. LayerLowFreqPermitivity = per;
  181. }
  182. void LayeredEarthEM::SetLayerBreathPermitivity(const VectorXr &per) {
  183. if (per.size() != this->GetNumberOfLayers() )
  184. throw EarthModelParametersDoNotMatchNumberOfLayers( );
  185. LayerBreathPermitivity = per;
  186. }
  187. void LayeredEarthEM::SetLayerTauPermitivity(const VectorXr &per) {
  188. if (per.size() != this->GetNumberOfLayers() )
  189. throw EarthModelParametersDoNotMatchNumberOfLayers( );
  190. LayerTauPermitivity = per;
  191. }
  192. void LayeredEarthEM::SetLayerThickness(const VectorXr &thick) {
  193. if (thick.size() != this->GetNumberOfLayers() - 2)
  194. throw EarthModelParametersDoNotMatchNumberOfLayers( );
  195. LayerThickness = thick;
  196. }
  197. void LayeredEarthEM::SetNumberOfLayers(const int&nlay) {
  198. // Check for validity
  199. if (nlay < 2) {
  200. throw EarthModelWithLessThanTwoLayers();
  201. }
  202. //else if (nlay > MAXLAYERS)
  203. // throw EarthModelWithMoreThanMaxLayers();
  204. // Otherwise
  205. this->NumberOfLayers = nlay;
  206. this->NumberOfInterfaces = nlay-1;
  207. // Resize all layers
  208. //////////////////////////////////
  209. // Thicknesses set to zero
  210. LayerThickness = VectorXr::Zero(NumberOfLayers);
  211. ///////////////////////////////////
  212. // Conducitivy set to zero
  213. LayerConductivity = VectorXcr::Zero(NumberOfLayers);
  214. ////////////////////////////////////
  215. // Susceptibility set to One (free space)
  216. LayerSusceptibility = VectorXcr::Ones(NumberOfLayers);
  217. // workaround for Valgrind on AMD which compains about assignment to 1 of VectorXr
  218. LayerLowFreqSusceptibility = LayerSusceptibility.real(); // 1
  219. LayerHighFreqSusceptibility = LayerSusceptibility.real(); // 1
  220. LayerTauSusceptibility = LayerSusceptibility.imag(); // 0
  221. LayerBreathSusceptibility = LayerSusceptibility.imag(); // 0
  222. //////////////////////////////////////
  223. // Permitivity set to 1 (free space)
  224. //this->LayerPermitivity.resize(NumberOfLayers);
  225. //this->LayerPermitivity.setOnes();
  226. LayerPermitivity = VectorXcr::Ones(NumberOfLayers);
  227. // workaround for Valgrind on AMD which compains about assignment to 1 of VectorXr
  228. LayerLowFreqPermitivity = LayerPermitivity.real(); // 1
  229. LayerHighFreqPermitivity = LayerPermitivity.real(); // 1
  230. LayerTauPermitivity = LayerPermitivity.imag(); // 0
  231. LayerBreathPermitivity = LayerPermitivity.imag(); // 0
  232. }
  233. // ==================== INQUIRY ===================================
  234. VectorXcr LayeredEarthEM::GetLayerConductivity() {
  235. return this->LayerConductivity;
  236. }
  237. Complex LayeredEarthEM::GetLayerConductivity(const int &ilay) {
  238. return this->LayerConductivity(ilay);
  239. }
  240. Complex LayeredEarthEM::GetLayerSusceptibility(const int &ilay) {
  241. return this->LayerSusceptibility(ilay);
  242. }
  243. VectorXcr LayeredEarthEM::GetLayerSusceptibility( ) {
  244. return this->LayerSusceptibility;
  245. }
  246. Complex LayeredEarthEM::GetLayerPermitivity(const int &ilay) {
  247. return this->LayerPermitivity(ilay);
  248. }
  249. VectorXcr LayeredEarthEM::GetLayerPermitivity( ) {
  250. return this->LayerPermitivity;
  251. }
  252. Real LayeredEarthEM::GetLayerLowFreqSusceptibility(const int &ilay) {
  253. return this->LayerLowFreqSusceptibility(ilay);
  254. }
  255. VectorXr LayeredEarthEM::GetLayerLowFreqSusceptibility( ) {
  256. return this->LayerLowFreqSusceptibility;
  257. }
  258. Real LayeredEarthEM::GetLayerHighFreqSusceptibility(const int &ilay) {
  259. return this->LayerHighFreqSusceptibility(ilay);
  260. }
  261. VectorXr LayeredEarthEM::GetLayerHighFreqSusceptibility( ) {
  262. return this->LayerHighFreqSusceptibility;
  263. }
  264. Real LayeredEarthEM::GetLayerTauSusceptibility(const int &ilay) {
  265. return this->LayerTauSusceptibility(ilay);
  266. }
  267. VectorXr LayeredEarthEM::GetLayerTauSusceptibility( ) {
  268. return this->LayerTauSusceptibility;
  269. }
  270. Real LayeredEarthEM::GetLayerBreathSusceptibility(const int &ilay) {
  271. return this->LayerBreathSusceptibility(ilay);
  272. }
  273. VectorXr LayeredEarthEM::GetLayerBreathSusceptibility( ) {
  274. return this->LayerBreathSusceptibility;
  275. }
  276. Real LayeredEarthEM::GetLayerLowFreqPermitivity(const int &ilay) {
  277. return this->LayerLowFreqPermitivity(ilay);
  278. }
  279. VectorXr LayeredEarthEM::GetLayerLowFreqPermitivity( ) {
  280. return this->LayerLowFreqPermitivity;
  281. }
  282. Real LayeredEarthEM::GetLayerHighFreqPermitivity(const int &ilay) {
  283. return this->LayerHighFreqPermitivity(ilay);
  284. }
  285. VectorXr LayeredEarthEM::GetLayerHighFreqPermitivity( ) {
  286. return this->LayerHighFreqPermitivity;
  287. }
  288. Real LayeredEarthEM::GetLayerTauPermitivity(const int &ilay) {
  289. return this->LayerTauPermitivity(ilay);
  290. }
  291. VectorXr LayeredEarthEM::GetLayerTauPermitivity( ) {
  292. return this->LayerTauPermitivity;
  293. }
  294. Real LayeredEarthEM::GetLayerBreathPermitivity(const int &ilay) {
  295. return this->LayerBreathPermitivity(ilay);
  296. }
  297. VectorXr LayeredEarthEM::GetLayerBreathPermitivity( ) {
  298. return this->LayerBreathPermitivity;
  299. }
  300. }
  301. /* vim: set tabstop=4 expandtab: */
  302. /* vim: set filetype=cpp: */