|
| 1 | +/*************************************************************************** |
| 2 | + * Copyright (C) 2009 by gempa GmbH |
| 3 | + * |
| 4 | + * Author: Jan Becker |
| 5 | + |
| 6 | + * |
| 7 | + * Modifications 2010 - 2011 by Stefan Heimers <[email protected]> |
| 8 | + * Modifications 2014 - 2015 by JeromeSalichon |
| 9 | + * |
| 10 | + * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, |
| 11 | + * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES |
| 12 | + * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND |
| 13 | + * NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT |
| 14 | + * HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, |
| 15 | + * WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING |
| 16 | + * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR |
| 17 | + * OTHER DEALINGS IN THE SOFTWARE. |
| 18 | + ***************************************************************************/ |
| 19 | + |
| 20 | +#define SEISCOMP_COMPONENT MLR |
| 21 | + |
| 22 | +#include "mlr.h" |
| 23 | +#include <seiscomp3/logging/log.h> |
| 24 | +#include <seiscomp3/core/strings.h> |
| 25 | +#include <seiscomp3/processing/magnitudeprocessor.h> |
| 26 | +#include <seiscomp3/math/geo.h> |
| 27 | + |
| 28 | +#include <iostream> |
| 29 | + |
| 30 | +#include <boost/bind.hpp> |
| 31 | +#include <unistd.h> |
| 32 | + |
| 33 | +// Only valid within 0-20 degrees |
| 34 | +#define DELTA_MIN 0. |
| 35 | +#define DELTA_MAX 20 |
| 36 | + |
| 37 | +// Only valid for depths up to 800 km |
| 38 | +#define DEPTH_MAX 800 |
| 39 | + |
| 40 | +#define MAG_TYPE "MLr" |
| 41 | + |
| 42 | +using namespace std; |
| 43 | + |
| 44 | + |
| 45 | +using namespace Seiscomp; |
| 46 | +using namespace Seiscomp::Processing; |
| 47 | + |
| 48 | + |
| 49 | +ADD_SC_PLUGIN( |
| 50 | + "MLr GNS magnitude, J. Ristau method", |
| 51 | + "Derived from MLh magnitude & Sc3 Tutorial methods, J.Salichon, GNS Science New Zealand, J.Becker, Gempa", |
| 52 | + 0, 0, 1 |
| 53 | +); |
| 54 | + |
| 55 | + |
| 56 | +class Magnitude_MLR : public Processing::MagnitudeProcessor { |
| 57 | + public: |
| 58 | + struct param_struct { |
| 59 | + double dist; |
| 60 | + double A; |
| 61 | + bool nomag; |
| 62 | + }; |
| 63 | + |
| 64 | + list<param_struct> list_of_parametersets; |
| 65 | + param_struct selected_parameterset; |
| 66 | + |
| 67 | + Magnitude_MLR() : Processing::MagnitudeProcessor(MAG_TYPE) {} |
| 68 | + |
| 69 | + //Same amplitude measurement as MLv (ML Sc3 amplitude on vertical) |
| 70 | + std::string amplitudeType() const { |
| 71 | + return "MLv" ; |
| 72 | + } |
| 73 | + |
| 74 | + // Do not need to compute another amplitude. |
| 75 | + |
| 76 | + bool setup(const Settings &settings) { |
| 77 | + // Reset the configuration |
| 78 | + list_of_parametersets.clear(); |
| 79 | + try { |
| 80 | + // Read the settings variable MLr.params. This can be found in the applications |
| 81 | + // configuration file at: |
| 82 | + // module.trunk.global.MLr.params |
| 83 | + // or per station (highest priority) |
| 84 | + // module.trunk.NZ.WEL.MLr.params |
| 85 | + if ( !initParameters(list_of_parametersets, settings.getString("MLr.params")) ) |
| 86 | + return false; |
| 87 | + } |
| 88 | + catch ( ... ) {} |
| 89 | + |
| 90 | + return true; |
| 91 | + } |
| 92 | + |
| 93 | + MagnitudeProcessor::Status computeMagnitude( |
| 94 | + double amplitude, // in milimeters |
| 95 | + double period, // in seconds |
| 96 | + double delta, // in degrees |
| 97 | + double depth, // in kilometers |
| 98 | + double &value) |
| 99 | + { |
| 100 | + if ( delta < DELTA_MIN || delta > DELTA_MAX ) |
| 101 | + return DistanceOutOfRange; |
| 102 | + |
| 103 | + if ( depth > DEPTH_MAX ) |
| 104 | + return DepthOutOfRange; |
| 105 | + |
| 106 | + return compute_ML(amplitude, delta, depth, &value); |
| 107 | + } |
| 108 | + |
| 109 | + private: |
| 110 | + // create a C++ list of structs containing all configured |
| 111 | + // parameter sets. |
| 112 | + // run this once at program start |
| 113 | + bool initParameters(list<param_struct> ¶mlist, const string ¶ms) { |
| 114 | + string paramset, range_str, minrange_str; |
| 115 | + string A_str; |
| 116 | + |
| 117 | + // for each parameter set |
| 118 | + istringstream iss_params(params); |
| 119 | + while ( getline(iss_params,paramset,';') ) { |
| 120 | + // extract the parameters from this parameter set |
| 121 | + istringstream iss_paramset(paramset); |
| 122 | + iss_paramset >> range_str; |
| 123 | + iss_paramset >> A_str; |
| 124 | + |
| 125 | + param_struct new_paramset; |
| 126 | + if ( !Core::fromString(new_paramset.dist, range_str) ) { |
| 127 | + SEISCOMP_ERROR("MLr: %s: range is not a numeric value", |
| 128 | + params.c_str()); |
| 129 | + return false; |
| 130 | + } |
| 131 | + |
| 132 | + if ( A_str == "nomag" ) { |
| 133 | + new_paramset.A = 0; |
| 134 | + new_paramset.nomag = true; |
| 135 | + } |
| 136 | + else { |
| 137 | + if ( !Core::fromString(new_paramset.A, A_str) ) { |
| 138 | + SEISCOMP_ERROR("MLr: %s: not a numeric value", |
| 139 | + A_str.c_str()); |
| 140 | + return false; |
| 141 | + } |
| 142 | + |
| 143 | + new_paramset.nomag = false; |
| 144 | + } |
| 145 | + |
| 146 | + paramlist.push_back(new_paramset); |
| 147 | + } |
| 148 | + |
| 149 | + return true; |
| 150 | + } |
| 151 | + |
| 152 | + |
| 153 | + // select the right parameter set for a given distance. init_parameters() must |
| 154 | + // have been called before using this function. |
| 155 | + param_struct selectParameters( double distance, const list<param_struct> ¶mlist) { |
| 156 | + double last_dist = 1000000; // arbitrary number larger than any expected distance; |
| 157 | + param_struct selected_parameters; |
| 158 | + |
| 159 | + // defaults in case the distance is out of the definded ranges |
| 160 | + selected_parameters.nomag = true; |
| 161 | + selected_parameters.dist = 0; |
| 162 | + selected_parameters.A = 0; |
| 163 | + |
| 164 | + list<param_struct>::const_iterator paramit; |
| 165 | + for ( paramit = paramlist.begin(); paramit != paramlist.end(); ++paramit ) { |
| 166 | + if ( (paramit->dist < last_dist) && (paramit->dist >= distance) ) { |
| 167 | + selected_parameters = *paramit; |
| 168 | + last_dist = paramit->dist; |
| 169 | + } |
| 170 | + } |
| 171 | + |
| 172 | + return selected_parameters; |
| 173 | + } |
| 174 | + |
| 175 | + MagnitudeProcessor::Status compute_ML( |
| 176 | + double amplitude, // in micrometers |
| 177 | + double delta, // in degrees |
| 178 | + double depth, // in kilometers |
| 179 | + double *mag) { |
| 180 | + float epdistkm,hypdistkm; |
| 181 | + float LogAmpREF; //MLr parameters |
| 182 | + float magTemp ; //MLr parameters |
| 183 | + |
| 184 | + if ( amplitude <= 0. ) { |
| 185 | + *mag = 0; |
| 186 | + return Error; |
| 187 | + } |
| 188 | + |
| 189 | + // examples: |
| 190 | + // epdistkm <= 60 km => mag=log10(waampl) + 0.018 *epdistkm+1.77 ; |
| 191 | + // 60 < epdistkm <= 700 => mag=log10(waampl) + 0.0038*epdistkm+2.62 ; |
| 192 | + // more generic: mag = log10(waampl) + A * epdistkm |
| 193 | + // |
| 194 | + // a list of distance ranges and corresponding values for A |
| 195 | + // is read from a config file. |
| 196 | + |
| 197 | + // calculate the distance in kilometers from the distance in degrees |
| 198 | + epdistkm = Math::Geo::deg2km(delta); |
| 199 | + hypdistkm = sqrt(epdistkm*epdistkm + depth * depth); |
| 200 | + |
| 201 | + // read the values for A and epdistkm from the config file and |
| 202 | + // select the right set depending on the distance |
| 203 | + selected_parameterset = selectParameters(hypdistkm, list_of_parametersets); |
| 204 | + |
| 205 | + SEISCOMP_DEBUG("epdistkm: %f\n",epdistkm); |
| 206 | + SEISCOMP_DEBUG("hypdistkm: %f\n",hypdistkm); |
| 207 | + |
| 208 | + if ( selected_parameterset.nomag ) { |
| 209 | + SEISCOMP_DEBUG( "epicenter distance out of configured range, no magnitude"); |
| 210 | + return DistanceOutOfRange; |
| 211 | + } |
| 212 | + else { |
| 213 | + SEISCOMP_DEBUG("The selected range is: %f", selected_parameterset.dist); |
| 214 | + SEISCOMP_DEBUG("A: %f", selected_parameterset.A); |
| 215 | + //SEISCOMP_DEBUG("MLr: station %s.%s: %s",settings.networkCode.c_str(), settings.stationCode.c_str(),s.c_str()); |
| 216 | + |
| 217 | + // ORI(sed): *mag = log10(amplitude) + selected_parameterset.A * hypdistkm + selected_parameterset.B; |
| 218 | + // JRistau: *mag = log10(amplitude) - Log(Amp(hypdistkm)) + Selected_parameterset.A ; |
| 219 | + // |
| 220 | + // Computing terms of magnitude |
| 221 | + LogAmpREF = 0.2869 - (1.272 * 1e-3) * hypdistkm - (1.493 * log10(hypdistkm) ) ; |
| 222 | + |
| 223 | + // magTemp = log10(amplitude) - LogAmpREF + StatCorr ; |
| 224 | + // A is defined as STATION dependent ie. module.trunk.NZ.WEL.MLr.params ; |
| 225 | + // Not corrected stations: Default is 0. Use nomag option to disable; |
| 226 | + magTemp = log10(amplitude) - LogAmpREF + selected_parameterset.A ; |
| 227 | + *mag = magTemp ; |
| 228 | + |
| 229 | + SEISCOMP_DEBUG("Mlr: %f", *mag); |
| 230 | + return OK; |
| 231 | + } |
| 232 | + } |
| 233 | +}; |
| 234 | + |
| 235 | +REGISTER_MAGNITUDEPROCESSOR(Magnitude_MLR, MAG_TYPE); |
0 commit comments