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.

LayeredEarth.cpp 4.1KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134
  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: layeredearth.cpp 210 2015-02-25 02:57:03Z tirons $
  10. **/
  11. #include "LayeredEarth.h"
  12. namespace Lemma {
  13. // ==================== FRIENDS ======================
  14. std::ostream &operator << (std::ostream &stream, const LayeredEarth &ob) {
  15. stream << ob.Serialize() << "\n";
  16. return stream;
  17. }
  18. // ==================== LIFECYCLE ===================================
  19. LayeredEarth::LayeredEarth( const ctor_key& key ) : EarthModel( key ),
  20. NumberOfLayers(0), NumberOfInterfaces(0) {
  21. }
  22. LayeredEarth::~LayeredEarth() {
  23. }
  24. LayeredEarth::LayeredEarth(const YAML::Node& node, const ctor_key& key) : EarthModel(node, key) {
  25. NumberOfLayers = node["NumberOfLayers"].as<int>();
  26. NumberOfInterfaces = node["NumberOfInterfaces"].as<int>();
  27. LayerThickness = node["LayerThickness"].as<VectorXr>();
  28. }
  29. YAML::Node LayeredEarth::Serialize() const {
  30. YAML::Node node = EarthModel::Serialize();
  31. node["NumberOfLayers"] = NumberOfLayers;
  32. node["NumberOfInterfaces"] = NumberOfInterfaces;
  33. node["LayerThickness"] = LayerThickness;
  34. node.SetTag( this->GetName() );
  35. return node;
  36. }
  37. // ==================== OPERATIONS ===================================
  38. // ==================== ACCESS ===================================
  39. //--------------------------------------------------------------------------------------
  40. // Class: LayeredEarth
  41. // Method: GetName
  42. // Description: Class identifier
  43. //--------------------------------------------------------------------------------------
  44. inline std::string LayeredEarth::GetName ( ) const {
  45. return CName;
  46. } // ----- end of method LayeredEarth::get_GetName -----
  47. void LayeredEarth::SetLayerThickness(const VectorXr &thick) {
  48. if (thick.size() != this->NumberOfLayers - 2)
  49. throw EarthModelParametersDoNotMatchNumberOfLayers( );
  50. LayerThickness = thick;
  51. }
  52. // ==================== INQUIRY ===================================
  53. int LayeredEarth::GetNumberOfLayers () {
  54. return this->NumberOfLayers;
  55. }
  56. int LayeredEarth::GetNumberOfNonAirLayers () {
  57. return this->NumberOfLayers - 1;
  58. }
  59. VectorXr LayeredEarth::GetLayerThickness( ) {
  60. return LayerThickness;
  61. }
  62. Real LayeredEarth::GetLayerThickness(const int &ilay) {
  63. // Take into account infinite top and bottom layers
  64. // estimate infinity by 1000 m
  65. if (ilay < 0 || ilay > NumberOfLayers - 1) {
  66. throw RequestForNonValidEarthModelParameter( );
  67. } else if (ilay == 0) {
  68. return 1000.;
  69. } else if (ilay == NumberOfLayers - 1) {
  70. return 1000.;
  71. } else {
  72. return this->LayerThickness(ilay-1);
  73. }
  74. }
  75. Real LayeredEarth::GetLayerDepth(const int &ilay) {
  76. Real depth = 0;
  77. if (ilay == 0) {
  78. return depth;
  79. } else {
  80. for (int i=1; i<=ilay; ++i) {
  81. depth += GetLayerThickness(i);
  82. }
  83. }
  84. return depth;
  85. }
  86. int LayeredEarth::GetLayerAtThisDepth(const Real& depth) {
  87. if (depth <= 0 || NumberOfLayers < 2) {
  88. return 0;
  89. }
  90. Real laydep = 0;
  91. for (int ilay=0; ilay<NumberOfLayers-2; ++ilay) {
  92. laydep += LayerThickness[ilay];
  93. if (laydep >= depth) { return ilay+1; }
  94. }
  95. return NumberOfLayers-1;
  96. }
  97. EarthModelWithLessThanTwoLayers::EarthModelWithLessThanTwoLayers() :
  98. runtime_error( "EARTH MODEL WITH LESS THAN TWO LAYERS") { }
  99. EarthModelWithMoreThanMaxLayers::EarthModelWithMoreThanMaxLayers() :
  100. runtime_error( "EARTH MODEL WITH MORE THAN MAX LAYERS") { }
  101. EarthModelParametersDoNotMatchNumberOfLayers::
  102. EarthModelParametersDoNotMatchNumberOfLayers( ) :
  103. runtime_error( "EARTH MODEL PARAMETERS DO NOT MATCH NUMBER OF LAYERS")
  104. {}
  105. RequestForNonValidEarthModelParameter::
  106. RequestForNonValidEarthModelParameter() :
  107. runtime_error( "REQUEST FOR NON VALID EARTH MODEL PARAMETER") {}
  108. }