diff --git a/libs/seiscomp/processing/amplitudes/CMakeLists.txt b/libs/seiscomp/processing/amplitudes/CMakeLists.txt index 31564fff..5d9813c3 100644 --- a/libs/seiscomp/processing/amplitudes/CMakeLists.txt +++ b/libs/seiscomp/processing/amplitudes/CMakeLists.txt @@ -9,7 +9,10 @@ SET(AMPS_SOURCES ML.cpp MLc.cpp MLh.cpp - MLv.cpp + MLv.cpp + mb_Lg.cpp + MLcan.cpp + mb_VC.cpp msbb.cpp Ms20.cpp Mwp.cpp @@ -27,6 +30,9 @@ SET(AMPS_HEADERS MLc.h MLh.h MLv.h + mb_Lg.h + MLcan.h + mb_VC.h msbb.h Ms20.h Mwp.h diff --git a/libs/seiscomp/processing/amplitudes/MLcan.cpp b/libs/seiscomp/processing/amplitudes/MLcan.cpp new file mode 100644 index 00000000..415ec1c1 --- /dev/null +++ b/libs/seiscomp/processing/amplitudes/MLcan.cpp @@ -0,0 +1,603 @@ +/* ################################################################################ +* # Copyright (C) 2024 by IGN Spain # +* # # +* # author: J. Barco, E. Suarez # +* # email: jbarco@transportes.gob.es , eadiaz@transportes.gob.es # +* # last modified: 2024-03-20 # +* # # +* # This program is free software; you can redistribute it and/or modify # +* # it under the terms of the GNU General Public License as published by # +* # the Free Software Foundation; either version 2 of the License, or # +* # (at your option) any later version. # +* # # +* # This program is distributed in the hope that it will be useful, # +* # but WITHOUT ANY WARRANTY; without even the implied warranty of # +* # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # +* # GNU General Public License for more details. # +* # # +* # You should have received a copy of the GNU General Public License # +* # along with this program; if not, write to the # +* # Free Software Foundation, Inc., # +* # 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. # +* ################################################################################ */ + + +#define SEISCOMP_COMPONENT MLcan + +#include +#include + +#include +#include + + +using namespace std; + + +namespace Seiscomp { +namespace Processing { + + +namespace { +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +AmplitudeProcessor::AmplitudeValue average( + const AmplitudeProcessor::AmplitudeValue &v0, + const AmplitudeProcessor::AmplitudeValue &v1) +{ + AmplitudeProcessor::AmplitudeValue v; + // Average both values + v.value = (v0.value + v1.value) * 0.5; + + // Compute lower and upper uncertainty + double v0l = v0.value; + double v0u = v0.value; + double v1l = v1.value; + double v1u = v1.value; + + if ( v0.lowerUncertainty ) v0l -= *v0.lowerUncertainty; + if ( v0.upperUncertainty ) v0u += *v0.upperUncertainty; + if ( v1.lowerUncertainty ) v1l -= *v1.lowerUncertainty; + if ( v1.upperUncertainty ) v1u += *v1.upperUncertainty; + + double l = 0, u = 0; + + l = max(l, v.value - v0l); + l = max(l, v.value - v1l); + + u = max(u, v0u - v.value); + u = max(u, v1u - v.value); + + v.lowerUncertainty = l; + v.upperUncertainty = u; + + return v; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +AmplitudeProcessor::AmplitudeTime average( + const AmplitudeProcessor::AmplitudeTime &t0, + const AmplitudeProcessor::AmplitudeTime &t1) +{ + AmplitudeProcessor::AmplitudeTime t; + t.reference = Core::Time((double(t0.reference) + double(t1.reference)) * 0.5); + + // Compute lower and upper uncertainty + Core::Time t0b = t0.reference + Core::TimeSpan(t0.begin); + Core::Time t0e = t0.reference + Core::TimeSpan(t0.end); + Core::Time t1b = t1.reference + Core::TimeSpan(t1.begin); + Core::Time t1e = t1.reference + Core::TimeSpan(t1.end); + + Core::Time minTime = t.reference; + Core::Time maxTime = t.reference; + + minTime = min(minTime, t0b); + minTime = min(minTime, t0e); + minTime = min(minTime, t1b); + minTime = min(minTime, t1e); + + maxTime = max(maxTime, t0b); + maxTime = max(maxTime, t0e); + maxTime = max(maxTime, t1b); + maxTime = max(maxTime, t1e); + + t.begin = (double)(minTime - t.reference); + t.end = (double)(maxTime - t.reference); + + return t; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +AmplitudeProcessor::AmplitudeValue gmean( + const AmplitudeProcessor::AmplitudeValue &v0, + const AmplitudeProcessor::AmplitudeValue &v1) +{ + AmplitudeProcessor::AmplitudeValue v; + // Average both values + v.value = sqrt(v0.value * v1.value); + + // Compute lower and upper uncertainty + double v0l = v0.value; + double v0u = v0.value; + double v1l = v1.value; + double v1u = v1.value; + + if ( v0.lowerUncertainty ) v0l -= *v0.lowerUncertainty; + if ( v0.upperUncertainty ) v0u += *v0.upperUncertainty; + if ( v1.lowerUncertainty ) v1l -= *v1.lowerUncertainty; + if ( v1.upperUncertainty ) v1u += *v1.upperUncertainty; + + double l = 0, u = 0; + + l = max(l, v.value - v0l); + l = max(l, v.value - v1l); + + u = max(u, v0u - v.value); + u = max(u, v1u - v.value); + + v.lowerUncertainty = l; + v.upperUncertainty = u; + + return v; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +REGISTER_AMPLITUDEPROCESSOR(AmplitudeProcessor_MLcan2h, "MLcan"); +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +AmplitudeProcessor_MLcan::AmplitudeProcessor_MLcan() +: AbstractAmplitudeProcessor_ML("MLcan") {} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +AmplitudeProcessor_MLcan2h::AmplitudeProcessor_MLcan2h() +: Processing::AmplitudeProcessor("MLcan") { + setSignalEnd("min(R / 4 , 150)"); + setMinSNR(0); + // Maximum distance is 8 degrees + setMaxDist(8); + // Maximum depth is 80 km + setMaxDepth(80); + + setUsedComponent(Horizontal); + + _combiner = TakeAverage; + + _ampN.setUsedComponent(FirstHorizontal); + _ampE.setUsedComponent(SecondHorizontal); + + _ampE.setPublishFunction(bind(&AmplitudeProcessor_MLcan2h::newAmplitude, this, placeholders::_1, placeholders::_2)); + _ampN.setPublishFunction(bind(&AmplitudeProcessor_MLcan2h::newAmplitude, this, placeholders::_1, placeholders::_2)); + + // Propagate configuration to single processors + _ampN.setConfig(config()); + _ampE.setConfig(config()); +} + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +const AmplitudeProcessor *AmplitudeProcessor_MLcan2h::componentProcessor(Component comp) const { + switch ( comp ) { + case FirstHorizontalComponent: + return &_ampN; + case SecondHorizontalComponent: + return &_ampE; + default: + break; + } + + return nullptr; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +const DoubleArray *AmplitudeProcessor_MLcan2h::processedData(Component comp) const { + switch ( comp ) { + case FirstHorizontalComponent: + return _ampN.processedData(comp); + case SecondHorizontalComponent: + return _ampE.processedData(comp); + default: + break; + } + + return nullptr; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +void AmplitudeProcessor_MLcan2h::reprocess(OPT(double) searchBegin, OPT(double) searchEnd) { + setStatus(WaitingForData, 0); + _ampN.setConfig(config()); + _ampE.setConfig(config()); + + _results[0] = _results[1] = Core::None; + + _ampN.reprocess(searchBegin, searchEnd); + _ampE.reprocess(searchBegin, searchEnd); + + if ( !isFinished() ) { + if ( _ampN.status() > Finished ) + setStatus(_ampN.status(), _ampN.statusValue()); + else + setStatus(_ampE.status(), _ampE.statusValue()); + } +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +int AmplitudeProcessor_MLcan2h::capabilities() const { + return _ampN.capabilities() | Combiner; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +AmplitudeProcessor::IDList +AmplitudeProcessor_MLcan2h::capabilityParameters(Capability cap) const { + if ( cap == Combiner ) { + IDList params; + params.push_back("Average"); + params.push_back("Max"); + params.push_back("Min"); + params.push_back("Geometric mean"); + return params; + } + + return _ampN.capabilityParameters(cap); +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +bool AmplitudeProcessor_MLcan2h::setParameter(Capability cap, const string &value) { + if ( cap == Combiner ) { + if ( value == "Min" ) { + _combiner = TakeMin; + return true; + } + else if ( value == "Max" ) { + _combiner = TakeMax; + return true; + } + else if ( value == "Average" ) { + _combiner = TakeAverage; + return true; + } + else if ( value == "Geometric mean" ) { + _combiner = TakeGeometricMean; + return true; + } + + return false; + } + + _ampN.setParameter(cap, value); + return _ampE.setParameter(cap, value); +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +std::string AmplitudeProcessor_MLcan2h::parameter(Capability cap) const { + if ( cap == Combiner ) { + switch ( _combiner ) { + case TakeMin: + return "Min"; + case TakeMax: + return "Max"; + case TakeAverage: + return "Average"; + case TakeGeometricMean: + return "Geometric mean"; + default: + break; + } + } + + return _ampN.parameter(cap); +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +bool AmplitudeProcessor_MLcan2h::setup(const Settings &settings) { + // Copy the stream configurations (gain, orientation, responses, ...) to + // the horizontal processors + _ampN.streamConfig(FirstHorizontalComponent) = streamConfig(FirstHorizontalComponent); + _ampE.streamConfig(SecondHorizontalComponent) = streamConfig(SecondHorizontalComponent); + + _combiner = TakeAverage; + + try { + string s = settings.getString("amplitudes." + _type + ".combiner"); + if ( s == "average" ) + _combiner = TakeAverage; + else if ( s == "max" ) + _combiner = TakeMax; + else if ( s == "min" ) + _combiner = TakeMin; + else if ( s == "geometric_mean" ) + _combiner = TakeGeometricMean; + else { + SEISCOMP_ERROR("%s: invalid combiner type for station %s.%s: %s", + _type.c_str(), + settings.networkCode.c_str(), settings.stationCode.c_str(), + s.c_str()); + return false; + } + } + catch ( ... ) {} + + if ( !AmplitudeProcessor::setup(settings) ) return false; + + // Setup each component + if ( !_ampN.setup(settings) || !_ampE.setup(settings) ) return false; + + return true; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +void AmplitudeProcessor_MLcan2h::setTrigger(const Core::Time& trigger) { + AmplitudeProcessor::setTrigger(trigger); + _ampE.setTrigger(trigger); + _ampN.setTrigger(trigger); +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +void AmplitudeProcessor_MLcan2h::computeTimeWindow() { + // Copy configuration to each component + _ampN.setConfig(config()); + _ampE.setConfig(config()); + + _ampN.computeTimeWindow(); + _ampE.computeTimeWindow(); + + // computeTimeWindow evaluates the signal times. This copies back the + // evaluated times. + setConfig(_ampE.config()); + + setTimeWindow(_ampE.timeWindow() | _ampN.timeWindow()); +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +void AmplitudeProcessor_MLcan2h::setEnvironment(const DataModel::Origin *hypocenter, + const DataModel::SensorLocation *receiver, + const DataModel::Pick *pick) { + _ampN.setEnvironment(hypocenter, receiver, pick); + _ampE.setEnvironment(hypocenter, receiver, pick); +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +void AmplitudeProcessor_MLcan2h::reset() { + AmplitudeProcessor::reset(); + + _results[0] = _results[1] = Core::None; + + _ampE.reset(); + _ampN.reset(); +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +void AmplitudeProcessor_MLcan2h::close() const { + // TODO: Check for best available amplitude here +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +bool AmplitudeProcessor_MLcan2h::feed(const Record *record) { + // Both processors finished already? + if ( _ampE.isFinished() && _ampN.isFinished() ) return false; + + // Did an error occur? + if ( status() > WaveformProcessor::Finished ) return false; + + if ( record->channelCode() == _streamConfig[FirstHorizontalComponent].code() ) { + if ( !_ampN.isFinished() ) { + _ampN.feed(record); + if ( _ampN.status() == InProgress ) + setStatus(WaveformProcessor::InProgress, _ampN.statusValue()); + else if ( _ampN.isFinished() && _ampE.isFinished() ) { + if ( !isFinished() ) { + if ( _ampE.status() == Finished ) + setStatus(_ampN.status(), _ampN.statusValue()); + else + setStatus(_ampE.status(), _ampE.statusValue()); + } + } + } + } + else if ( record->channelCode() == _streamConfig[SecondHorizontalComponent].code() ) { + if ( !_ampE.isFinished() ) { + _ampE.feed(record); + if ( _ampE.status() == InProgress ) + setStatus(WaveformProcessor::InProgress, _ampE.statusValue()); + else if ( _ampE.isFinished() && _ampN.isFinished() ) { + if ( !isFinished() ) { + if ( _ampN.status() == Finished ) + setStatus(_ampE.status(), _ampE.statusValue()); + else + setStatus(_ampN.status(), _ampN.statusValue()); + } + } + } + } + + return true; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +bool AmplitudeProcessor_MLcan2h::computeAmplitude(const DoubleArray &data, + size_t i1, size_t i2, + size_t si1, size_t si2, + double offset, + AmplitudeIndex *dt, + AmplitudeValue *amplitude, + double *period, double *snr) { + return false; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +void AmplitudeProcessor_MLcan2h::newAmplitude(const AmplitudeProcessor *proc, + const AmplitudeProcessor::Result &res) { + if ( isFinished() ) return; + + int idx = 0; + + if ( proc == &_ampE ) { + idx = 0; + } + else if ( proc == &_ampN ) { + idx = 1; + } + + _results[idx] = ComponentResult(); + _results[idx]->value = res.amplitude; + _results[idx]->time = res.time; + _results[idx]->snr = res.snr; + _results[idx]->period = res.period; + + if ( _results[0] && _results[1] ) { + setStatus(Finished, 100.); + Result newRes; + newRes.record = res.record; + + switch ( _combiner ) { + case TakeAverage: + newRes.amplitude = average(_results[0]->value, _results[1]->value); + newRes.time = average(_results[0]->time, _results[1]->time); + newRes.snr = (_results[0]->snr + _results[1]->snr) * 0.5; + newRes.period = (_results[0]->period + _results[1]->period) * 0.5; + newRes.component = Horizontal; + break; + case TakeGeometricMean: + newRes.amplitude = gmean(_results[0]->value, _results[1]->value); + newRes.time = average(_results[0]->time, _results[1]->time); + newRes.snr = (_results[0]->snr + _results[1]->snr) * 0.5; + newRes.period = (_results[0]->period + _results[1]->period) * 0.5; + newRes.component = Horizontal; + break; + case TakeMin: + if ( _results[0]->value.value <= _results[1]->value.value ) { + newRes.amplitude = _results[0]->value; + newRes.time = _results[0]->time; + newRes.snr = _results[0]->snr; + newRes.period = _results[0]->period; + newRes.component = _ampE.usedComponent(); + } + else { + newRes.amplitude = _results[1]->value; + newRes.time = _results[1]->time; + newRes.snr = _results[1]->snr; + newRes.period = _results[1]->period; + newRes.component = _ampN.usedComponent(); + } + break; + case TakeMax: + if ( _results[0]->value.value >= _results[1]->value.value ) { + newRes.amplitude = _results[0]->value; + newRes.time = _results[0]->time; + newRes.snr = _results[0]->snr; + newRes.period = _results[0]->period; + newRes.component = _ampE.usedComponent(); + } + else { + newRes.amplitude = _results[1]->value; + newRes.time = _results[1]->time; + newRes.snr = _results[1]->snr; + newRes.period = _results[1]->period; + newRes.component = _ampN.usedComponent(); + } + break; + }; + + emitAmplitude(newRes); + } +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +} +} diff --git a/libs/seiscomp/processing/amplitudes/MLcan.h b/libs/seiscomp/processing/amplitudes/MLcan.h new file mode 100644 index 00000000..d4ebbfcd --- /dev/null +++ b/libs/seiscomp/processing/amplitudes/MLcan.h @@ -0,0 +1,118 @@ +/* ################################################################################ +* # Copyright (C) 2024 by IGN Spain # +* # # +* # author: J. Barco, E. Suarez # +* # email: jbarco@transportes.gob.es , eadiaz@transportes.gob.es # +* # last modified: 2024-03-20 # +* # # +* # This program is free software; you can redistribute it and/or modify # +* # it under the terms of the GNU General Public License as published by # +* # the Free Software Foundation; either version 2 of the License, or # +* # (at your option) any later version. # +* # # +* # This program is distributed in the hope that it will be useful, # +* # but WITHOUT ANY WARRANTY; without even the implied warranty of # +* # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # +* # GNU General Public License for more details. # +* # # +* # You should have received a copy of the GNU General Public License # +* # along with this program; if not, write to the # +* # Free Software Foundation, Inc., # +* # 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. # +* ################################################################################ */ + + +#ifndef SEISCOMP_PROCESSING_AMPLITUDEPROCESSOR_MLcan_H +#define SEISCOMP_PROCESSING_AMPLITUDEPROCESSOR_MLcan_H + + +#include + + +namespace Seiscomp { +namespace Processing { + + +//! Wrapper class that allows access to protected methods for +//! the proxy (see below). +class SC_SYSTEM_CLIENT_API AmplitudeProcessor_MLcan : public AbstractAmplitudeProcessor_ML { + public: + AmplitudeProcessor_MLcan(); + + friend class AmplitudeProcessor_MLcan2h; +}; + + +//! Proxy amplitude processor that holds two MLh processors to calculate +//! the amplitudes on both horizontals and then averages the result. +//! This class does not handle waveforms itself. It directs them to the +//! corresponding amplitude processors instead. +class SC_SYSTEM_CLIENT_API AmplitudeProcessor_MLcan2h : public AmplitudeProcessor { + public: + AmplitudeProcessor_MLcan2h(); + + public: + int capabilities() const override; + IDList capabilityParameters(Capability cap) const override; + bool setParameter(Capability cap, const std::string &value) override; + std::string parameter(Capability cap) const override; + + bool setup(const Settings &settings) override; + + void setTrigger(const Core::Time& trigger) override; + + void setEnvironment(const DataModel::Origin *hypocenter, + const DataModel::SensorLocation *receiver, + const DataModel::Pick *pick) override; + + void computeTimeWindow() override; + + void reset() override; + void close() const override; + + bool feed(const Record *record) override; + + bool computeAmplitude(const DoubleArray &data, + size_t i1, size_t i2, + size_t si1, size_t si2, + double offset, + AmplitudeIndex *dt, AmplitudeValue *amplitude, + double *period, double *snr) override; + + const AmplitudeProcessor *componentProcessor(Component comp) const override; + const DoubleArray *processedData(Component comp) const override; + + void reprocess(OPT(double) searchBegin, OPT(double) searchEnd) override;; + + + private: + void newAmplitude(const AmplitudeProcessor *proc, + const AmplitudeProcessor::Result &res); + + + private: + struct ComponentResult { + AmplitudeValue value; + AmplitudeTime time; + double snr; + double period; + }; + + enum CombinerProc { + TakeMin, + TakeMax, + TakeAverage, + TakeGeometricMean + }; + + mutable AmplitudeProcessor_MLcan _ampE, _ampN; + CombinerProc _combiner; + OPT(ComponentResult) _results[2]; +}; + + +} +} + + +#endif diff --git a/libs/seiscomp/processing/amplitudes/mb_Lg.cpp b/libs/seiscomp/processing/amplitudes/mb_Lg.cpp new file mode 100644 index 00000000..41479486 --- /dev/null +++ b/libs/seiscomp/processing/amplitudes/mb_Lg.cpp @@ -0,0 +1,625 @@ +/* ################################################################################ +* # Copyright (C) 2024 by IGN Spain # +* # # +* # author: J. Barco, E. Suarez # +* # email: jbarco@transportes.gob.es , eadiaz@transportes.gob.es # +* # last modified: 2024-03-20 # +* # # +* # This program is free software; you can redistribute it and/or modify # +* # it under the terms of the GNU General Public License as published by # +* # the Free Software Foundation; either version 2 of the License, or # +* # (at your option) any later version. # +* # # +* # This program is distributed in the hope that it will be useful, # +* # but WITHOUT ANY WARRANTY; without even the implied warranty of # +* # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # +* # GNU General Public License for more details. # +* # # +* # You should have received a copy of the GNU General Public License # +* # along with this program; if not, write to the # +* # Free Software Foundation, Inc., # +* # 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. # +* ################################################################################ */ + + +#define SEISCOMP_COMPONENT Amplitudemb_Lg + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + +using namespace std; +using namespace Seiscomp::Math; + + +namespace Seiscomp { +namespace Processing { + + +REGISTER_AMPLITUDEPROCESSOR(AmplitudeProcessor_mb_Lg2h, "mb_Lg"); + + +namespace { + + +AmplitudeProcessor::AmplitudeValue average( + const AmplitudeProcessor::AmplitudeValue &v0, + const AmplitudeProcessor::AmplitudeValue &v1) +{ + AmplitudeProcessor::AmplitudeValue v; + // Average both values + v.value = (v0.value + v1.value) * 0.5; + + // Compute lower and upper uncertainty + double v0l = v0.value; + double v0u = v0.value; + double v1l = v1.value; + double v1u = v1.value; + + if ( v0.lowerUncertainty ) v0l -= *v0.lowerUncertainty; + if ( v0.upperUncertainty ) v0u += *v0.upperUncertainty; + if ( v1.lowerUncertainty ) v1l -= *v1.lowerUncertainty; + if ( v1.upperUncertainty ) v1u += *v1.upperUncertainty; + + double l = 0, u = 0; + + l = max(l, v.value - v0l); + l = max(l, v.value - v1l); + + u = max(u, v0u - v.value); + u = max(u, v1u - v.value); + + v.lowerUncertainty = l; + v.upperUncertainty = u; + + return v; +} + + +AmplitudeProcessor::AmplitudeTime average( + const AmplitudeProcessor::AmplitudeTime &t0, + const AmplitudeProcessor::AmplitudeTime &t1) +{ + AmplitudeProcessor::AmplitudeTime t; + t.reference = Core::Time((double(t0.reference) + double(t1.reference)) * 0.5); + + // Compute lower and upper uncertainty + Core::Time t0b = t0.reference + Core::TimeSpan(t0.begin); + Core::Time t0e = t0.reference + Core::TimeSpan(t0.end); + Core::Time t1b = t1.reference + Core::TimeSpan(t1.begin); + Core::Time t1e = t1.reference + Core::TimeSpan(t1.end); + + Core::Time minTime = t.reference; + Core::Time maxTime = t.reference; + + minTime = min(minTime, t0b); + minTime = min(minTime, t0e); + minTime = min(minTime, t1b); + minTime = min(minTime, t1e); + + maxTime = max(maxTime, t0b); + maxTime = max(maxTime, t0e); + maxTime = max(maxTime, t1b); + maxTime = max(maxTime, t1e); + + t.begin = double(minTime - t.reference); + t.end = double(maxTime - t.reference); + + return t; +} + + +AmplitudeProcessor::AmplitudeValue gmean( + const AmplitudeProcessor::AmplitudeValue &v0, + const AmplitudeProcessor::AmplitudeValue &v1) +{ + AmplitudeProcessor::AmplitudeValue v; + // Average both values + v.value = sqrt(v0.value * v1.value); + + // Compute lower and upper uncertainty + double v0l = v0.value; + double v0u = v0.value; + double v1l = v1.value; + double v1u = v1.value; + + if ( v0.lowerUncertainty ) v0l -= *v0.lowerUncertainty; + if ( v0.upperUncertainty ) v0u += *v0.upperUncertainty; + if ( v1.lowerUncertainty ) v1l -= *v1.lowerUncertainty; + if ( v1.upperUncertainty ) v1u += *v1.upperUncertainty; + + double l = 0, u = 0; + + l = max(l, v.value - v0l); + l = max(l, v.value - v1l); + + u = max(u, v0u - v.value); + u = max(u, v1u - v.value); + + v.lowerUncertainty = l; + v.upperUncertainty = u; + + return v; +} + + +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +AmplitudeProcessor_mb_Lg::AmplitudeProcessor_mb_Lg() +: AbstractAmplitudeProcessor_ML("mb_Lg") {} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +AmplitudeProcessor_mb_Lg2h::AmplitudeProcessor_mb_Lg2h() +: AmplitudeProcessor("mb_Lg") { + setDefaultConfiguration(); + AmplitudeProcessor_mb_Lg2h::reset(); + + setUsedComponent(Horizontal); + + _ampN.setUsedComponent(FirstHorizontal); + _ampE.setUsedComponent(SecondHorizontal); + + _ampE.setPublishFunction(bind(&AmplitudeProcessor_mb_Lg2h::newAmplitude, this, placeholders::_1, placeholders::_2)); + _ampN.setPublishFunction(bind(&AmplitudeProcessor_mb_Lg2h::newAmplitude, this, placeholders::_1, placeholders::_2)); +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +void AmplitudeProcessor_mb_Lg2h::setDefaultConfiguration() { + setSignalEnd("min(R / 3 + 30, 150)"); + setMinSNR(0); + setMinDist(0); + setMaxDist(8); + setMinDepth(-10); + setMaxDepth(80); + + _ampE._preFilter = _ampN._preFilter = "BW(3,0.5,12)"; + + // Propagate configuration to single processors + _ampE.setConfig(config()); + _ampN.setConfig(config()); +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +int AmplitudeProcessor_mb_Lg2h::capabilities() const { + return _ampN.capabilities() | Combiner; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +AmplitudeProcessor::IDList +AmplitudeProcessor_mb_Lg2h::capabilityParameters(Capability cap) const { + if ( cap == Combiner ) { + IDList params; + params.push_back("Max"); + params.push_back("Average"); + params.push_back("Min"); + params.push_back("Geometric mean"); + return params; + } + + return _ampN.capabilityParameters(cap); +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +bool AmplitudeProcessor_mb_Lg2h::setParameter(Capability cap, const std::string &value) { + if ( cap == Combiner ) { + if ( value == "Min" ) { + _combiner = TakeMin; + return true; + } + else if ( value == "Max" ) { + _combiner = TakeMax; + return true; + } + else if ( value == "Average" ) { + _combiner = TakeAverage; + return true; + } + else if ( value == "Geometric mean" ) { + _combiner = TakeGeometricMean; + return true; + } + + return false; + } + + _ampN.setParameter(cap, value); + return _ampE.setParameter(cap, value); +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +std::string AmplitudeProcessor_mb_Lg2h::parameter(Capability cap) const { + if ( cap == Combiner ) { + switch ( _combiner ) { + case TakeMin: + return "Min"; + case TakeMax: + return "Max"; + case TakeAverage: + return "Average"; + case TakeGeometricMean: + return "Geometric mean"; + default: + break; + } + } + + return _ampN.parameter(cap); +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +void AmplitudeProcessor_mb_Lg2h::reset() { + AmplitudeProcessor::reset(); + + _results[0] = _results[1] = Core::None; + + _ampE.reset(); + _ampN.reset(); +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +bool AmplitudeProcessor_mb_Lg2h::setup(const Settings &settings) { + setDefaultConfiguration(); + + // Copy the stream configurations (gain, orientation, responses, ...) to + // the horizontal processors + _ampN.streamConfig(FirstHorizontalComponent) = streamConfig(FirstHorizontalComponent); + _ampE.streamConfig(SecondHorizontalComponent) = streamConfig(SecondHorizontalComponent); + + _combiner = TakeMax; + + try { + string s = settings.getString("amplitudes." + _type + ".combiner"); + if ( s == "average" ) + _combiner = TakeAverage; + else if ( s == "max" ) + _combiner = TakeMax; + else if ( s == "min" ) + _combiner = TakeMin; + else if ( s == "geometric_mean" ) + _combiner = TakeGeometricMean; + else { + SEISCOMP_ERROR("%s: invalid combiner type for station %s.%s: %s", + _type.c_str(), + settings.networkCode.c_str(), settings.stationCode.c_str(), + s.c_str()); + return false; + } + } + catch ( ... ) {} + + try { + _amplitudeScale = settings.getDouble("amplitudes." + _type + ".amplitudeScale"); + } + catch ( ... ) {} + + if ( !AmplitudeProcessor::setup(settings) ) { + return false; + } + + // Setup each component + return _ampN.setup(settings) && _ampE.setup(settings); +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +void AmplitudeProcessor_mb_Lg2h::setTrigger(const Core::Time &trigger) { + AmplitudeProcessor::setTrigger(trigger); + _ampE.setTrigger(trigger); + _ampN.setTrigger(trigger); +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +void AmplitudeProcessor_mb_Lg2h::setEnvironment(const DataModel::Origin *hypocenter, + const DataModel::SensorLocation *receiver, + const DataModel::Pick *pick) { + _ampN.setEnvironment(hypocenter, receiver, pick); + _ampE.setEnvironment(hypocenter, receiver, pick); +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +void AmplitudeProcessor_mb_Lg2h::computeTimeWindow() { + // Copy configuration to each component + _ampN.setConfig(config()); + _ampE.setConfig(config()); + + _ampE.computeTimeWindow(); + _ampN.computeTimeWindow(); + + // computeTimeWindow evaluates the signal times. This copies back the + // evaluated times. + setConfig(_ampE.config()); + + setTimeWindow(_ampE.timeWindow() | _ampN.timeWindow()); +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +bool AmplitudeProcessor_mb_Lg2h::computeAmplitude( + const DoubleArray &, + size_t, size_t, + size_t, size_t, + double, + AmplitudeIndex *, AmplitudeValue *, + double *, double *) { + return false; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +void AmplitudeProcessor_mb_Lg2h::close() const {} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +bool AmplitudeProcessor_mb_Lg2h::feed(const Record *record) { + // Both processors finished already? + if ( _ampE.isFinished() && _ampN.isFinished() ) return false; + + // Did an error occur? + if ( status() > WaveformProcessor::Finished ) return false; + + if ( record->channelCode() == _streamConfig[FirstHorizontalComponent].code() ) { + if ( !_ampN.isFinished() ) { + _ampN.feed(record); + if ( _ampN.status() == InProgress ) + setStatus(WaveformProcessor::InProgress, _ampN.statusValue()); + else if ( _ampN.isFinished() && _ampE.isFinished() ) { + if ( !isFinished() ) { + if ( _ampE.status() == Finished ) + setStatus(_ampN.status(), _ampN.statusValue()); + else + setStatus(_ampE.status(), _ampE.statusValue()); + } + } + } + } + else if ( record->channelCode() == _streamConfig[SecondHorizontalComponent].code() ) { + if ( !_ampE.isFinished() ) { + _ampE.feed(record); + if ( _ampE.status() == InProgress ) + setStatus(WaveformProcessor::InProgress, _ampE.statusValue()); + else if ( _ampE.isFinished() && _ampN.isFinished() ) { + if ( !isFinished() ) { + if ( _ampN.status() == Finished ) + setStatus(_ampE.status(), _ampE.statusValue()); + else + setStatus(_ampN.status(), _ampN.statusValue()); + } + } + } + } + + return true; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +void AmplitudeProcessor_mb_Lg2h::newAmplitude(const AmplitudeProcessor *proc, + const AmplitudeProcessor::Result &res) { + + if ( isFinished() ) return; + + int idx = 0; + + if ( proc == &_ampE ) { + idx = 0; + } + else if ( proc == &_ampN ) { + idx = 1; + } + + _results[idx] = ComponentResult(); + _results[idx]->value = res.amplitude; + _results[idx]->time = res.time; + _results[idx]->snr = res.snr; + _results[idx]->period = res.period; + + if ( _results[0] && _results[1] ) { + setStatus(Finished, 100.); + Result newRes; + newRes.record = res.record; + + switch ( _combiner ) { + case TakeAverage: + newRes.amplitude = average(_results[0]->value, _results[1]->value); + newRes.time = average(_results[0]->time, _results[1]->time); + newRes.snr = (_results[0]->snr + _results[1]->snr) * 0.5; + newRes.period = (_results[0]->period + _results[1]->period) * 0.5; + newRes.component = Horizontal; + + break; + case TakeGeometricMean: + newRes.amplitude = gmean(_results[0]->value, _results[1]->value); + newRes.time = average(_results[0]->time, _results[1]->time); + newRes.snr = (_results[0]->snr + _results[1]->snr) * 0.5; + newRes.period = (_results[0]->period + _results[1]->period) * 0.5; + newRes.component = Horizontal; + break; + case TakeMin: + if ( _results[0]->value.value <= _results[1]->value.value ) { + newRes.amplitude = _results[0]->value; + newRes.time = _results[0]->time; + newRes.snr = _results[0]->snr; + newRes.period = _results[0]->period; + newRes.component = _ampE.usedComponent(); + } + else { + newRes.amplitude = _results[1]->value; + newRes.time = _results[1]->time; + newRes.snr = _results[1]->snr; + newRes.period = _results[1]->period; + newRes.component = _ampN.usedComponent(); + } + break; + case TakeMax: + if ( _results[0]->value.value >= _results[1]->value.value ) { + newRes.amplitude = _results[0]->value; + newRes.time = _results[0]->time; + newRes.snr = _results[0]->snr; + newRes.period = _results[0]->period; + newRes.component = _ampE.usedComponent(); + } + else { + newRes.amplitude = _results[1]->value; + newRes.time = _results[1]->time; + newRes.snr = _results[1]->snr; + newRes.period = _results[1]->period; + newRes.component = _ampN.usedComponent(); + } + break; + }; + + // apply amplitude scaling + newRes.amplitude.value *= _amplitudeScale; + + if ( newRes.amplitude.upperUncertainty ) { + *newRes.amplitude.upperUncertainty *= _amplitudeScale; + } + + if ( newRes.amplitude.lowerUncertainty ) { + *newRes.amplitude.lowerUncertainty *= _amplitudeScale; + + } + + emitAmplitude(newRes); + } +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +const AmplitudeProcessor *AmplitudeProcessor_mb_Lg2h::componentProcessor(Component comp) const { + switch ( comp ) { + case FirstHorizontalComponent: + return &_ampN; + case SecondHorizontalComponent: + return &_ampE; + default: + break; + } + + return nullptr; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +const DoubleArray *AmplitudeProcessor_mb_Lg2h::processedData(Component comp) const { + switch ( comp ) { + case FirstHorizontalComponent: + return _ampN.processedData(comp); + case SecondHorizontalComponent: + return _ampE.processedData(comp); + default: + break; + } + + return nullptr; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +void AmplitudeProcessor_mb_Lg2h::reprocess(OPT(double) searchBegin, OPT(double) searchEnd) { + setStatus(WaitingForData, 0); + _ampN.setConfig(config()); + _ampE.setConfig(config()); + + _results[0] = _results[1] = Core::None; + + _ampN.reprocess(searchBegin, searchEnd); + _ampE.reprocess(searchBegin, searchEnd); + + if ( !isFinished() ) { + if ( _ampN.status() > Finished ) + setStatus(_ampN.status(), _ampN.statusValue()); + else + setStatus(_ampE.status(), _ampE.statusValue()); + } +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +} +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< diff --git a/libs/seiscomp/processing/amplitudes/mb_Lg.h b/libs/seiscomp/processing/amplitudes/mb_Lg.h new file mode 100644 index 00000000..a1ba64b4 --- /dev/null +++ b/libs/seiscomp/processing/amplitudes/mb_Lg.h @@ -0,0 +1,117 @@ +/* ################################################################################ +* # Copyright (C) 2024 by IGN Spain # +* # # +* # author: J. Barco, E. Suarez # +* # email: jbarco@transportes.gob.es , eadiaz@transportes.gob.es # +* # last modified: 2024-03-20 # +* # # +* # This program is free software; you can redistribute it and/or modify # +* # it under the terms of the GNU General Public License as published by # +* # the Free Software Foundation; either version 2 of the License, or # +* # (at your option) any later version. # +* # # +* # This program is distributed in the hope that it will be useful, # +* # but WITHOUT ANY WARRANTY; without even the implied warranty of # +* # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # +* # GNU General Public License for more details. # +* # # +* # You should have received a copy of the GNU General Public License # +* # along with this program; if not, write to the # +* # Free Software Foundation, Inc., # +* # 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. # +* ################################################################################ */ + + +#ifndef SEISCOMP_PROCESSING_AMPLITUDEPROCESSOR_mb_Lg_H +#define SEISCOMP_PROCESSING_AMPLITUDEPROCESSOR_mb_Lg_H + + +#include + + +namespace Seiscomp { +namespace Processing { + + +//! Wrapper class that allows access to protected methods for +//! the proxy (see below). +class SC_SYSTEM_CLIENT_API AmplitudeProcessor_mb_Lg : public AbstractAmplitudeProcessor_ML { + public: + AmplitudeProcessor_mb_Lg(); + + std::string _preFilter; + + friend class AmplitudeProcessor_MLc2h; +}; + + +class SC_SYSTEM_CLIENT_API AmplitudeProcessor_mb_Lg2h : public AmplitudeProcessor { + public: + AmplitudeProcessor_mb_Lg2h(); + + public: + int capabilities() const override; + IDList capabilityParameters(Capability cap) const override; + bool setParameter(Capability cap, const std::string &value) override; + std::string parameter(Capability cap) const override; + + void reset() override; + bool setup(const Settings &settings) override; + + void setTrigger(const Core::Time &trigger) override; + + void setEnvironment(const DataModel::Origin *hypocenter, + const DataModel::SensorLocation *receiver, + const DataModel::Pick *pick) override; + + void computeTimeWindow() override; + + void close() const override; + + bool feed(const Record *record) override; + + bool computeAmplitude(const DoubleArray &data, + size_t i1, size_t i2, + size_t si1, size_t si2, + double offset, + AmplitudeIndex *dt, AmplitudeValue *amplitude, + double *period, double *snr) override; + + const AmplitudeProcessor *componentProcessor(Component comp) const override; + const DoubleArray *processedData(Component comp) const override; + + void reprocess(OPT(double) searchBegin, OPT(double) searchEnd) override; + + + private: + void setDefaultConfiguration(); + + void newAmplitude(const AmplitudeProcessor *proc, + const AmplitudeProcessor::Result &res); + + struct ComponentResult { + AmplitudeValue value; + AmplitudeTime time; + double snr; + double period; + }; + + enum CombinerProc { + TakeMax, + TakeMin, + TakeAverage, + TakeGeometricMean + }; + + mutable AmplitudeProcessor_mb_Lg _ampE, _ampN; + CombinerProc _combiner{TakeMax}; + OPT(ComponentResult) _results[2]; + double _amplitudeScale{1.0}; +}; + + +} +} + + +#endif diff --git a/libs/seiscomp/processing/amplitudes/mb_VC.cpp b/libs/seiscomp/processing/amplitudes/mb_VC.cpp new file mode 100644 index 00000000..01d9d10e --- /dev/null +++ b/libs/seiscomp/processing/amplitudes/mb_VC.cpp @@ -0,0 +1,144 @@ +/* ################################################################################ +* # Copyright (C) 2024 by IGN Spain # +* # # +* # author: J. Barco, E. Suarez # +* # email: jbarco@transportes.gob.es , eadiaz@transportes.gob.es # +* # last modified: 2024-03-20 # +* # # +* # This program is free software; you can redistribute it and/or modify # +* # it under the terms of the GNU General Public License as published by # +* # the Free Software Foundation; either version 2 of the License, or # +* # (at your option) any later version. # +* # # +* # This program is distributed in the hope that it will be useful, # +* # but WITHOUT ANY WARRANTY; without even the implied warranty of # +* # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # +* # GNU General Public License for more details. # +* # # +* # You should have received a copy of the GNU General Public License # +* # along with this program; if not, write to the # +* # Free Software Foundation, Inc., # +* # 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. # +* ################################################################################ */ + + +#include +#include + +#include +#include + + +namespace Seiscomp { +namespace Processing { + + +REGISTER_AMPLITUDEPROCESSOR(AmplitudeProcessor_mB_VC, "mB_VC"); +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +AmplitudeProcessor_mB_VC::AmplitudeProcessor_mB_VC() +: AmplitudeProcessor("mB_VC") { + // 60 s should be OK for rupture durations up to ~100 s. + // This default value MUST NOT be increased because if the amplitudes are + // computed in the stream picker, which doesn't know of origins, it also + // doesn't know the S-P time. For distances of 5 degrees the S wave would + // leak into the time window thus contaminating the measurement. + setSignalEnd("min(D * 11.5, 60)"); + setMinDist(5); + setMaxDist(105); +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +void AmplitudeProcessor_mB_VC::finalizeAmplitude(DataModel::Amplitude *amplitude) const { + if ( !amplitude ) + return; + + try { + DataModel::TimeQuantity time(amplitude->timeWindow().reference()); + amplitude->setScalingTime(time); + } + catch ( ... ) { + } +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +bool AmplitudeProcessor_mB_VC::computeAmplitude(const DoubleArray &data, + size_t i1, size_t i2, + size_t si1, size_t si2, + double offset,AmplitudeIndex *dt, + AmplitudeValue *amplitude, + double *period, double *snr) { + /* + * Low-level signal amplitude computation. This is magnitude specific. + * + * Input: + * f double array of length n + * i1,i2 indices defining the measurement window, + * 0 <= i1 < i2 <= n + * offset this is subtracted from the samples in f before + * computation + * + * Output: + * dt Point at which the measurement was mad/completed. May + * be the peak time or end of integration. + * amplitude amplitude. This may be a peak amplitude as well as a + * sum or integral. + * period dominant period of the signal. Optional. If the period + * is not computed, set it to -1. + */ + int imax = find_absmax(data.size(), data.typedData(), si1, si2, offset); + double amax = fabs(data[imax] - offset); + double pmax = -1; + + // Bei Mwp bestimmt man amax zur Berechnung des SNR. Die eigentliche + // Amplitude ist aber was anderes! Daher ist SNR-Bestimmung auch + // magnitudenspezifisch! + + if ( *_noiseAmplitude == 0. ) + *snr = 1000000.0; + else + *snr = amax / *_noiseAmplitude; + + // SNR check + if (*snr < _config.snrMin) { + setStatus(LowSNR, *snr); + return false; + } + + dt->index = imax; + // Amplitudes are send in nanometers + *period = pmax; + + amplitude->value = amax; + + if ( _streamConfig[_usedComponent].gain != 0.0 ) + amplitude->value /= _streamConfig[_usedComponent].gain; + else { + setStatus(MissingGain, 0.0); + return false; + } + + // Convert m/s to nm/s + amplitude->value *= 1.E09; + + amplitude->value = std::abs(amplitude->value); + + return true; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + +} +} diff --git a/libs/seiscomp/processing/amplitudes/mb_VC.h b/libs/seiscomp/processing/amplitudes/mb_VC.h new file mode 100644 index 00000000..d008d57c --- /dev/null +++ b/libs/seiscomp/processing/amplitudes/mb_VC.h @@ -0,0 +1,60 @@ +/* ################################################################################ +* # Copyright (C) 2024 by IGN Spain # +* # # +* # author: J. Barco, E. Suarez # +* # email: jbarco@transportes.gob.es , eadiaz@transportes.gob.es # +* # last modified: 2024-03-20 # +* # # +* # This program is free software; you can redistribute it and/or modify # +* # it under the terms of the GNU General Public License as published by # +* # the Free Software Foundation; either version 2 of the License, or # +* # (at your option) any later version. # +* # # +* # This program is distributed in the hope that it will be useful, # +* # but WITHOUT ANY WARRANTY; without even the implied warranty of # +* # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # +* # GNU General Public License for more details. # +* # # +* # You should have received a copy of the GNU General Public License # +* # along with this program; if not, write to the # +* # Free Software Foundation, Inc., # +* # 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. # +* ################################################################################ */ + + + +#ifndef SEISCOMP_PROCESSING_AMPLITUDEPROCESSOR_mB_VC_H +#define SEISCOMP_PROCESSING_AMPLITUDEPROCESSOR_mB_VC_H + +#include + + +namespace Seiscomp { + +namespace Processing { + + +class SC_SYSTEM_CLIENT_API AmplitudeProcessor_mB_VC : public AmplitudeProcessor { + public: + AmplitudeProcessor_mB_VC(); + + public: + void finalizeAmplitude(DataModel::Amplitude *amplitude) const override; + + protected: + bool computeAmplitude(const DoubleArray &data, + size_t i1, size_t i2, + size_t si1, size_t si2, + double offset, + AmplitudeIndex *dt, + AmplitudeValue *amplitude, + double *period, double *snr) override; +}; + + +} + +} + + +#endif diff --git a/libs/seiscomp/processing/magnitudes/CMakeLists.txt b/libs/seiscomp/processing/magnitudes/CMakeLists.txt index b6073f87..aa30c918 100644 --- a/libs/seiscomp/processing/magnitudes/CMakeLists.txt +++ b/libs/seiscomp/processing/magnitudes/CMakeLists.txt @@ -9,6 +9,9 @@ SET( ML_idc.cpp MLc.cpp MLv.cpp + mb_VC.cpp + mb_Lg.cpp + MLcan.cpp Ms20.cpp msbb.cpp Mwp.cpp @@ -24,6 +27,9 @@ SET( ML.h MLc.h MLv.h + mb_VC.h + mb_Lg.h + MLcan.h Ms20.h msbb.h Mwp.h diff --git a/libs/seiscomp/processing/magnitudes/MLcan.cpp b/libs/seiscomp/processing/magnitudes/MLcan.cpp new file mode 100644 index 00000000..32245e77 --- /dev/null +++ b/libs/seiscomp/processing/magnitudes/MLcan.cpp @@ -0,0 +1,415 @@ +/* ################################################################################ +* # Copyright (C) 2024 by IGN Spain # +* # # +* # author: J. Barco, E. Suarez # +* # email: jbarco@transportes.gob.es , eadiaz@transportes.gob.es # +* # last modified: 2024-03-20 # +* # # +* # This program is free software; you can redistribute it and/or modify # +* # it under the terms of the GNU General Public License as published by # +* # the Free Software Foundation; either version 2 of the License, or # +* # (at your option) any later version. # +* # # +* # This program is distributed in the hope that it will be useful, # +* # but WITHOUT ANY WARRANTY; without even the implied warranty of # +* # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # +* # GNU General Public License for more details. # +* # # +* # You should have received a copy of the GNU General Public License # +* # along with this program; if not, write to the # +* # Free Software Foundation, Inc., # +* # 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. # +* ################################################################################ */ + +#define SEISCOMP_COMPONENT MLcan + +#include +#include +#include +#include +#include +#include + +using namespace std; + + +namespace Seiscomp { +namespace Processing { + + +namespace { + +std::string ExpectedAmplitudeUnit = "mm"; + +DEFINE_SMARTPOINTER(ExtraLocale); +class ExtraLocale : public Core::BaseObject { + public: + // general magnitude parameters + OPT(string) distanceMode; + OPT(string) calibrationType; + // parametric coefficients + OPT(double) c0_first; + OPT(double) c0_second; + OPT(double) c1; + OPT(double) c2; + OPT(double) c3; + OPT(double) c4; + OPT(double) c5; + // A0, non-parametric coefficients + OPT(LogA0) logA0; +}; + +} + + +REGISTER_MAGNITUDEPROCESSOR(MagnitudeProcessor_MLcan, "MLcan"); +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +MagnitudeProcessor_MLcan::MagnitudeProcessor_MLcan() +: Processing::MagnitudeProcessor("MLcan") {} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +std::string MagnitudeProcessor_MLcan::amplitudeType() const { + return MagnitudeProcessor::amplitudeType(); +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +bool MagnitudeProcessor_MLcan::setup(const Settings &settings) { + if ( !MagnitudeProcessor::setup(settings) ) { + return false; + } + + try {_maxDepth = settings.getDouble("magnitudes." + _type + ".maxDepth"); } + catch ( ... ) {} + + // distance constraints + try {_distanceMode = settings.getString("magnitudes." + _type + ".distMode"); } + catch ( ... ) {} + try {_minDistanceKm = Math::Geo::deg2km(settings.getDouble("magnitudes." + _type + ".minDist")); } + catch ( ... ) {} + try {_maxDistanceKm = Math::Geo::deg2km(settings.getDouble("magnitudes." + _type + ".maxDist")); } + catch ( ... ) {} + + // calibration function + try {_calibrationType = settings.getString("magnitudes." + _type + ".calibrationType"); } + catch ( ... ) {} + + if ( (_calibrationType != "A0") && (_calibrationType != "parametric") ) { + SEISCOMP_ERROR("%s: unrecognized calibration type: %s", _type.c_str(), + _calibrationType.c_str()); + return false; + } + + // parametric calibration function + try { _c0_first = settings.getDouble("magnitudes." + _type + ".parametric.c0_first"); } + catch ( ... ) {} + try { _c0_second = settings.getDouble("magnitudes." + _type + ".parametric.c0_second"); } + catch ( ... ) {} + try { _c1 = settings.getDouble("magnitudes." + _type + ".parametric.c1"); } + catch ( ... ) {} + try { _c2 = settings.getDouble("magnitudes." + _type + ".parametric.c2"); } + catch ( ... ) {} + try { _c3 = settings.getDouble("magnitudes." + _type + ".parametric.c3"); } + catch ( ... ) {} + try { _c4 = settings.getDouble("magnitudes." + _type + ".parametric.c4"); } + catch ( ... ) {} + try { _c5 = settings.getDouble("magnitudes." + _type + ".parametric.c5"); } + catch ( ... ) {} + + // A0, non-parametric calibration function + std::string defLogA0 = "0:-1.3,60:-2.8,100:-3.0,400:-4.5,1000:-5.85"; + bool logA0Default = true; + + try { + defLogA0 = settings.getString("magnitudes." + _type + ".A0.logA0"); + logA0Default = false; + } + catch ( ... ) {} + + // set distance-dependent intervals + if ( !_logA0.set(defLogA0) ) { + SEISCOMP_ERROR("%s: incorrect correction term log(A0): %s", _type.c_str(), + defLogA0.c_str()); + return false; + } + + if ( !logA0Default ) { + SEISCOMP_DEBUG("%s.%s: %s: logA0 from bindings = %s", + settings.networkCode, settings.stationCode, + type(), Core::toString(_logA0)); + } + + return true; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +bool MagnitudeProcessor_MLcan::initLocale(Locale *locale, + const Settings &settings, + const string &configPrefix) { + const Seiscomp::Config::Config *cfg = settings.localConfiguration; + ExtraLocalePtr extra = new ExtraLocale; + + // general + try { + extra->distanceMode = cfg->getString(configPrefix + "distMode"); + } + catch ( ... ) {} + try { + extra->calibrationType = cfg->getString(configPrefix + "calibrationType"); + } + catch ( ... ) {} + + // parametric + try { + extra->c0_first = cfg->getDouble(configPrefix + ".parametric.c0_first"); + } + catch ( ... ){} + try { + extra->c0_second = cfg->getDouble(configPrefix + ".parametric.c0_second"); + } + catch ( ... ) {}{} + try { + extra->c1 = cfg->getDouble(configPrefix + "parametric.c1"); + } + catch ( ... ) {} + try { + extra->c2 = cfg->getDouble(configPrefix + "parametric.c2"); + } + catch ( ... ) {} + try { + extra->c3 = cfg->getDouble(configPrefix + "parametric.c3"); + } + catch ( ... ) {} + try { + extra->c4 = cfg->getDouble(configPrefix + "parametric.c4"); + } + catch ( ... ) {} + try { + extra->c5 = cfg->getDouble(configPrefix + "parametric.c5"); + } + catch ( ... ) {} + + // A0, non-parametric + try { + auto logA0 = cfg->getStrings(configPrefix + "A0.logA0"); + if ( !logA0.empty() ) { + // If the first item contains a comma, then maybe it has been configured + // with double quotes. Raise an error. + if ( logA0[0].find(',') != string::npos ) { + SEISCOMP_ERROR("%sA0.logA0[0] contains a comma. Are the coefficients " + "enclosed with double quotes in the configuration?", + configPrefix); + return false; + } + + if ( logA0[0].find(';') != string::npos ) { + SEISCOMP_ERROR("%sA0.logA0 = %s contains semicolon. Supported format:" + " distance1:correction1,distance2:correction2, ...", + configPrefix, Core::toString(logA0).c_str()); + return false; + } + + extra->logA0 = LogA0(); + if ( !extra->logA0->set(logA0) ) { + SEISCOMP_ERROR("%s@%s: incorrect correction term log(A0)", + _type.c_str(), locale->name.c_str()); + return false; + } + } + } + catch ( ... ) {} + + locale->extra = extra; + + return true; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + +double MagnitudeProcessor_MLcan::getChanCorr(const double c0_first, const double c0_second, std::string ChanelCode){ + double chanCorr = 0; + string component = ""; + component = ChanelCode.back(); + if ((component == "1") || (component == "N")){ + chanCorr = c0_first; + } + if ((component == "2") || (component == "E")){ + chanCorr = c0_second; + } + else{ + chanCorr = (c0_first + c0_second)/2; + } + return chanCorr; + + }; + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +MagnitudeProcessor::Status MagnitudeProcessor_MLcan::computeMagnitude( + double amplitude, const std::string &unit, + double, double, double delta, double depth, + const DataModel::Origin *, + const DataModel::SensorLocation *receiver, + const DataModel::Amplitude *new_amplitude, + const Locale *locale, + double &value) { + + if ( amplitude <= 0 ) { + return AmplitudeOutOfRange; + } + + auto distanceMode = _distanceMode; + auto calibrationType = _calibrationType; + auto minimumDistanceKm = _minDistanceKm; + auto maximumDistanceKm = _maxDistanceKm; + auto maximumDepth = _maxDepth; + + ExtraLocale *extra = nullptr; + if ( locale ) { + extra = static_cast(locale->extra.get()); + if ( extra ) { + if ( extra->distanceMode ) { + distanceMode = *extra->distanceMode; + } + + if ( extra->calibrationType ) { + calibrationType = *extra->calibrationType; + } + } + + if ( locale->minimumDistance ) { + minimumDistanceKm = Math::Geo::deg2km(*locale->minimumDistance); + } + + if ( locale->maximumDistance ) { + maximumDistanceKm = Math::Geo::deg2km(*locale->maximumDistance); + } + } + else { + + SEISCOMP_DEBUG(" + maximum depth: %.3f km", maximumDepth); + if ( depth > maximumDepth ) { + return DepthOutOfRange; + } + SEISCOMP_DEBUG(" + minimum distance: %.3f km", minimumDistanceKm); + SEISCOMP_DEBUG(" + maximum distance: %.3f km", maximumDistanceKm); + } + + SEISCOMP_DEBUG(" + distance type: %s", distanceMode.c_str()); + + double hDistanceKm = Math::Geo::deg2km(delta); + double vDistanceKm = 0; + + if ( !receiver && distanceMode == "hypocentral" ) { + return MetaDataRequired; + } + + if ( receiver ) { + vDistanceKm = receiver->elevation() / 1000. + depth; + } + + double distanceKm; + if ( distanceMode == "hypocentral" ) { + distanceKm = sqrt(pow(hDistanceKm, 2) + pow(vDistanceKm, 2)); + } + else { + distanceKm = hDistanceKm; + } + + SEISCOMP_DEBUG(" + considered distance to station: %.3f km", distanceKm); + + if ( minimumDistanceKm >= 0 && distanceKm < minimumDistanceKm ) { + return DistanceOutOfRange; + } + + if ( maximumDistanceKm >= 0 && distanceKm > maximumDistanceKm ) { + return DistanceOutOfRange; + } + + if ( !convertAmplitude(amplitude, unit, ExpectedAmplitudeUnit) ) { + return InvalidAmplitudeUnit; + } + + SEISCOMP_DEBUG(" + %s distance: %.3f km (horiz. %.3f vert. %.3f)", + distanceMode.c_str(), distanceKm, hDistanceKm, vDistanceKm); + + double correction; + SEISCOMP_DEBUG(" + calibration type: %s", calibrationType.c_str()); + + if ( calibrationType == "parametric" ) { + // parametric calibration function + double c0 = 0; + // parametric calibration function + auto c0_first = (extra and extra->c0_first) ? *extra->c0_first : _c0_first; + auto c0_second = (extra and extra->c0_second) ? *extra->c0_second : _c0_second; + auto c1 = (extra and extra->c1) ? *extra->c1 : _c1; + auto c2 = (extra and extra->c2) ? *extra->c2 : _c2; + auto c3 = (extra and extra->c3) ? *extra->c3 : _c3; + auto c4 = (extra and extra->c4) ? *extra->c4 : _c4; + auto c5 = (extra and extra->c5) ? *extra->c5 : _c5; + + + + // extract the channel code and the location code from the amplitude + // as strings + std::string channelCode = new_amplitude->waveformID().channelCode(); + std::string locationCode = new_amplitude->waveformID().locationCode(); + + c0 = getChanCorr(c0_first, c0_second, channelCode); + SEISCOMP_DEBUG(" + Location: %s", locationCode); + SEISCOMP_DEBUG(" + Channel: %s", channelCode); + + + SEISCOMP_DEBUG(" + c0 - c5: %.5f %.5f %.5f %.5f %.5f %.5f", + c0, c1, c2, c3, c4, c5); + correction = c3 * log10(distanceKm / c5) + + c2 * (distanceKm + c4) + + c1 + + c0; + value = log10(amplitude) + correction; + } + else if ( calibrationType == "A0" ) { + // A0, non-parametric calibration function + try { + correction = -1.0 * (extra and extra->logA0 ? extra->logA0->at(distanceKm) : _logA0.at(distanceKm)); + + SEISCOMP_DEBUG(" + -log10(A0) at %.3f km: %.5f", distanceKm, correction); + value = log10(amplitude) + correction; + } + catch ( std::out_of_range & ) { + return DistanceOutOfRange; + } + } + + else { + return IncompleteConfiguration; + } + + SEISCOMP_DEBUG(" + amplitude: %.5f, magnitude: %.3f", amplitude, value); + + return OK; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +} +} diff --git a/libs/seiscomp/processing/magnitudes/MLcan.h b/libs/seiscomp/processing/magnitudes/MLcan.h new file mode 100644 index 00000000..e7786c83 --- /dev/null +++ b/libs/seiscomp/processing/magnitudes/MLcan.h @@ -0,0 +1,89 @@ +/* ################################################################################ +* # Copyright (C) 2024 by IGN Spain # +* # # +* # author: J. Barco, E. Suarez # +* # email: jbarco@transportes.gob.es , eadiaz@transportes.gob.es # +* # last modified: 2024-03-20 # +* # derived from Mezcua & Rueda (2021) # +* # https://doi.org/10.1007/s00445-022-01553-9 # +* # # +* # This program is free software; you can redistribute it and/or modify # +* # it under the terms of the GNU General Public License as published by # +* # the Free Software Foundation; either version 2 of the License, or # +* # (at your option) any later version. # +* # # +* # This program is distributed in the hope that it will be useful, # +* # but WITHOUT ANY WARRANTY; without even the implied warranty of # +* # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # +* # GNU General Public License for more details. # +* # # +* # You should have received a copy of the GNU General Public License # +* # along with this program; if not, write to the # +* # Free Software Foundation, Inc., # +* # 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. # +* ################################################################################ */ + + +#ifndef SEISCOMP_PROCESSING_MAGNITUDEPROCESSOR_MLcan_H +#define SEISCOMP_PROCESSING_MAGNITUDEPROCESSOR_MLcan_H + + +#include +#include +#include + + +namespace Seiscomp { +namespace Processing { + + +class SC_SYSTEM_CLIENT_API MagnitudeProcessor_MLcan : public MagnitudeProcessor { + public: + MagnitudeProcessor_MLcan(); + + + public: + bool setup(const Settings &settings) override; + + std::string amplitudeType() const override; + + + protected: + Status computeMagnitude(double amplitude, const std::string &unit, + double period, double snr, + double delta, double depth, + const DataModel::Origin *hypocenter, + const DataModel::SensorLocation *receiver, + const DataModel::Amplitude *, + const Locale *, + double &value) override; + + protected: + bool initLocale(Locale *locale, const Settings &settings, + const std::string &configPrefix) override; + + private: + double getChanCorr(const double c0_first, const double c0_second, std::string ChannelCode); + double _minDistanceKm{-1.0}; + double _maxDistanceKm{8.0 * KM_OF_DEGREE}; + double _maxDepth{80.0}; + std::string _distanceMode{"hypocentral"}; + std::string _calibrationType{"parametric"}; + // parameters for parametric magnitude calibration + double _c0_first{0.0}; + double _c0_second{0.0}; + double _c1{2.4449}; + double _c2{0.00142}; + double _c3{0.96657}; + double _c4{-40}; + double _c5{40}; + // parameters for non-parametric magnitude calibration + LogA0 _logA0; +}; + + +} +} + + +#endif diff --git a/libs/seiscomp/processing/magnitudes/descriptions/global_mb_VC.rst b/libs/seiscomp/processing/magnitudes/descriptions/global_mb_VC.rst new file mode 100644 index 00000000..ab3c7bf2 --- /dev/null +++ b/libs/seiscomp/processing/magnitudes/descriptions/global_mb_VC.rst @@ -0,0 +1,72 @@ +mb_VC is the standard body-wave magnitude. +Compare also with the :ref:`mb_VC magnitude `. + + +Amplitude +--------- + +mb_VC is defined on the amplitude of the first few cycles of the P-wave, +typically a time window of 20 s - 30 s. Only the first few cycles are used to +minimize the effects of radiation pattern and depth phases, which result in +complicate waveform signatures. +In |scname| mb_VC amplitudes are measured on vertical-component displacement seismograms +in a 30 s time window after simulation of a :term:`WWSSN_SP` short-period +seismometer. Amplitudes are used from stations with epicentral distances between +5° and 105° (configurable). The methods for measuring amplitudes are configurable +in the global bindings. + + +Station Magnitude +----------------- + +The general formula is + +.. math:: + + mb_VC = \log \left(\frac{A}{T}\right) + Q(h,\Delta) - 3.0 + +with A as the displacement amplitude in micrometers, T as the dominant period of +the signal in seconds, Q as a correction term for depth and distance. mb_VC is +usually determined at periods around 1s in adaptation to the use +of the World-Wide Standard Seismograph Network (WWSSN) short-period stations. +A scatter in the order of +/- 0.3 for the station magnitudes is usual. +Typically, mb_VC is determined for stations with distances larger than 5° to +have a distinct direct P-wave phase. A correction term for the distance has to +be determined empirically, which is quite complicate for distances smaller than 20°. +This reflects the complexity of the body waves that traverse only in the upper +mantle. mb_VC saturates at about magnitude 5.5 to 6.0 because the maximum amplitudes of larger +earthquakes occur at lower frequencies than the frequency range between 0.7 Hz - 2 Hz +used for the magnitude calculation. + +* Amplitude unit in |scname|: **nanometers** (nm) +* Time window: 30 s +* Default distance range: 5 - 105 deg, configurable: :confval:`magnitudes.mb.minDist`, + :confval:`magnitudes.mb.maxDist` +* Depth range: no limitation, for depth < 0 km, depth = 0 km is assumed + +.. note:: + + In 2013 the IASPEI commission (:cite:t:`iaspei-2013`) recommended a minimum distance of + 20 deg. However, :ref:`scautoloc` requires mb_VC amplitudes by default for + considering a pick. + For maintaining consistency, 5 deg is therefore kept as the default + for :confval:`magnitudes.mb.minDist`. + + +Network magnitude +----------------- + +By default, the trimmed mean is calculated from the station magnitudes to form +the :term:`network magnitude`. Outliers beyond the outer 12.5% percentiles are +removed before forming the mean. + + +Configuration +------------- + +Adjust the configurable parameters in global bindings in the mb_VC section or use +:file:`global.cfg` +as in :ref:`global_mlv`. Add mb_VC to the list of computed amplitudes and magnitudes +in the configuration of +:ref:`scamp` and :ref:`scmag` and in :ref:`scesv` or :ref:`scolv`/:ref:`scesv` +for visibility. diff --git a/libs/seiscomp/processing/magnitudes/descriptions/global_mb_VC.xml b/libs/seiscomp/processing/magnitudes/descriptions/global_mb_VC.xml new file mode 100644 index 00000000..dab31c04 --- /dev/null +++ b/libs/seiscomp/processing/magnitudes/descriptions/global_mb_VC.xml @@ -0,0 +1,35 @@ + + + + global + + Body wave magnitude at teleseismic distances + + + + + Body wave magnitude at teleseismic distances measured at 1 s period. + + + + + + Parameters for computing mb magnitudes from mb amplitudes. + + + + Minimum epicentral distance for computing mb. Note: According + to the IASPEI recommendations in 2013, the minimum distance + should be 20 deg. + + + + + Maximum epicentral distance for computing mb. + + + + + + + diff --git a/libs/seiscomp/processing/magnitudes/descriptions/global_mb_lg.rst b/libs/seiscomp/processing/magnitudes/descriptions/global_mb_lg.rst new file mode 100644 index 00000000..d26de044 --- /dev/null +++ b/libs/seiscomp/processing/magnitudes/descriptions/global_mb_lg.rst @@ -0,0 +1,169 @@ +mb_Lg is the custom magnitude for local events derived from the ML plugin. +The implementation is based on specifications by Spanish National Geographic Institute, +based on López, C. (2008). Nuevas fórmulas de magnitud para la Península Ibérica y su entorno. +Trabajo de investigación del Máster en Geofísica y Meteorología. +Departamento de Física de la Tierra, Astronomía y Astrofísica +Universidad Complutense de Madrid. Madrid. + + +The mb_Lg magnitude is very similar to the original :ref:`ML`, +except that by default + +* Amplitude pre-filtering is applied. +* A parametric :ref:`magnitude calibration ` function + applies. +* Hypocentral distance is used. + +Regionalization of magnitude computation is provided through global module +configuration. +Configuration of global bindings provides additional flexibility: + +* Amplitudes can be pre-filtered before simulating a :term:`Wood-Anderson seismometer` + (:confval:`amplitudes.mb_Lg.preFilter`), +* Wood-Anderson simulation is optional: + :confval:`amplitudes.mb_Lg.applyWoodAnderson`, +* Measured amplitudes can be scaled accounting for expected unit + (:confval:`amplitudes.mb_Lg.amplitudeScale`), +* A parametric or A0-based non-parametric + :ref:`magnitude calibration ` + function can be applied as controlled by + :confval:`magnitudes.mb_Lg.calibrationType`. +* Consider either hypocentral or epicentral distance for computing magnitudes + (:confval:`magnitudes.mb_Lg.distMode`). + +General (default) conditions apply, all values are configurable: + +* Amplitude pre-filtering: :ref:`BW(3,0.5,12) `. +* Amplitude unit in SeisComP: **millimeter** (mm) or as considered by the + configured calibration parameters. +* Optional amplitude scaling: 1. +* :term:`Wood-Anderson seismometer` simulation: yes. +* Time window: 150 s by :ref:`scautopick` or distance dependent + with :math:`endTime = min(R / 4, 150)`, configurable. +* Distance type: hypocentral. +* Distance range: 0 - 8 deg, measurements beyond 8 deg will be + strictly ignored. +* Depth range: <= 80 km. +* Magnitude calibration: parametric. + + +Amplitudes +---------- +mb_Lg amplitudes can be measured automatically by :ref:`scautopick` or :ref:`scamp` +or interactively by :ref:`scolv` very similarly to the original :ref:`ML`, +except that they can be pre-filtered and simulation of a :term:`Wood-Anderson seismometer` +is optional: :confval:`amplitudes.mb_Lg.preFilter`, +:confval:`amplitudes.mb_Lg.applyWoodAnderson`. +By default amplitudes are measured on both horizontal components where the +absolute maxima are taken. They are combined to a final measured amplitude by +taking the mean. The methods for measuring and combining the amplitudes are +configurable: +:confval:`amplitudes.mb_Lg.measureType`, :confval:`amplitudes.mb_Lg.combiner`. + +The Wood-Anderson simulation will convert input velocity data to ground +displacement in mm. The input data may be of a different unit after applying +:confval:`amplitudes.mb_Lg.preFilter`, e.g. when integration is applied, and / or +when Wood-Anderson simulation is disabled. Configure +:confval:`amplitudes.mb_Lg.amplitudeScale` for converting the unit of the +processed data to the unit expected by the +:ref:`station magnitude calibration ` for the measured +amplitude. + +.. note:: + + For comparing mb_Lg amplitudes with :ref:`ML amplitudes ` set the + global bindings parameters :: + + amplitudes.mb_Lg.preFilter = "" + amplitudes.mb_Lg.combiner = average + + +.. _mlc_station_magnitude: + +Station Magnitudes +------------------ + +Station magnitudes are computed from measured amplitudes automatically by +:ref:`scmag` +or interactively by :ref:`scolv`. By global bindings configuration mb_Lg considers + +* Hypocentral (default) or epicentral distance: :confval:`magnitudes.mb_Lg.distMode`. +* Distance range: :confval:`magnitudes.mb_Lg.minDist`, :confval:`magnitudes.mb_Lg.maxDist`. +* Events with depth up to :confval:`magnitudes.mb_Lg.maxDepth`. +* Parametric or non-parametric calibration functions + + * parametric when :confval:`magnitudes.mb_Lg.calibrationType` = "parametric"`: + + .. math:: + + mb_Lg = \log_{10}(A/(2^PI)) + c_0 * \log_{10}(r) + c_1 * (r) + c_2) + + where + + * *A*: displacement amplitude measured in unit of mm or as per configuration + * *r*: hypocentral (default) or epicentral distance + * *c0*, *c1*, *c2*: general calibration parameters + * *r*: Hypocentral (default) or epicentral distance as configured by + :confval:`magnitudes.mb_Lg.distMode`. + * *PI*: The value of Pi number (3.14159265359) + + * A0-based non-parametric when :confval:`magnitudes.mb_Lg.calibrationType` = "A0"`: + + .. math:: + + mb_Lg = \log_{10}(A) - \log_{10}(A_0) + + where + + * :math:`log_{10}(A_0)`: distance-dependent correction value. Read + :ref:`global_mlv` for the details. + +.. note:: + + The magnitude calibration function can regionalized by adjusting global module + configuration parameters in mb_Lg region profiles of + :confval:`magnitudes.mb_Lg.region.*` and in a *mb_Lg* Magnitude type profile e.g. + in :file:`global.cfg`. + + +Network Magnitude +----------------- + +The network magnitude is computed from station magnitudes automatically by +:ref:`scmag` or interactively by :ref:`scolv`. +Originally the median was computed from all station mb_Lg to form the +:term:`network magnitude` mb_Lg. Here, the trimmed mean is applied. Outliers +beyond the outer 12.5% percentiles are removed before forming the mean. The +method can be adjusted in :ref:`scmag` by :confval:`magnitudes.average`. + + + +Setup +===== + +#. **Set the configuration and calibration parameters** in the global bindings + similar + to :ref:`global_ml`. Instead of configuring lots of global bindings profiles + or station bindings one line per parameter can be added to the global module + configuration (:file:`global.cfg`) which takes the form + + .. code-block:: params + + module.trunk.NET.STA.amplitude.mb_Lg.preFilter = value + module.trunk.NET.STA.magnitude.mb_Lg.parametric.c0 = value + +#. Add mb_Lg to the list of default amplitudes and magnitudes if mb_Lg is to be + computed by automatic modules, e.g. of :ref:`scamp`, :ref:`scmag`. +#. Configure :ref:`scmag` (:confval:`magnitudes.average` in :file:`scmag.cfg`) + for choosing the method to form the + network magnitude from station magnitudes, e.g. + + .. code-block:: params + + magnitudes.average = mb_Lg:median + +#. Add mb_Lg to the list of magnitudes preferred by :ref:`scevent` + (:confval:`eventAssociation.magTypes` in :file:`scevent.cfg`) in order to let + mb_Lg become the preferred magnitude. +#. Set defaults/visibility of mb_Lg in :term:`GUI` modules, e.g. :ref:`scolv` + or :ref:`scesv`. diff --git a/libs/seiscomp/processing/magnitudes/descriptions/global_mb_lg.xml b/libs/seiscomp/processing/magnitudes/descriptions/global_mb_lg.xml new file mode 100644 index 00000000..740da4e8 --- /dev/null +++ b/libs/seiscomp/processing/magnitudes/descriptions/global_mb_lg.xml @@ -0,0 +1,256 @@ + + + + global + + Custom mbLg events measured on horizontal components with + correction for each component, from López, C. (2008). + + + + + + + Considered distance measure between source and receiver. + Possible values are + + hypocentral: hypocentral distance + + epicentral: epicentral distance + + + + + Type of magnitude calibration formula to be considered. + The calibration parameters are considered accordingly. + Currently supported are + + "parametric": consider parameters of parametric + configuration in parametric section + + "A0": consider parameters of non-parametric + configuration in A0 section. + + + + + Parameters for A0, non-parametric magnitude calibration. + + + + Overrides the calibration function log10(A0) + for computing mb_Lg per region. See logA0 + description in the bindings. + + + + + + + Parameters for parametric magnitude calibration: + mb_Lg = log10(A) + c0 * log10(r) + c1 * r + c2 + + + + + Overrides the calibration parameter c0 + + + + + + Overrides the calibration parameter c1 + for computing mb_Lg per region. See c1 + description in the bindings. + + + + + Overrides the calibration parameter c2 + for computing mb_Lg per region. See c2 + description in the bindings. + + + + + + + + + + + Custom magnitude for local events measured on horizontal components + + + + + + Parameters for measuring mb_Lg amplitudes. Add more parameters + by adding an amplitude type profile 'mb_Lg', + + + + The filter applied to raw records before applying + Wood-Anderson simulation. + + + + + Applying Wood-Anderson simulation. To achieve displacement + records without WA simulation, an integration filter can + be applied with the pre-filter. + + + + + Scaling value multiplied to the measured amplitudes to + match the amplitude units expected by the magnitude + calibration function. + + Expected amplitudes are + in units of mym but actual amplitudes provided from + Wood-Anderson-corrected seismograms are in units of mm: + amplitudeScale = 1000. + + If data are not corrected for WA, measured amplitudes + take the unit of gain-corrected data considering the + preFilter: + amplitudeScale converts between units of measured and + excpected amplitude. + + + + + Type for measuring amplitudes. Available: + + AbsMax: absolute maximum + + MinMax: half difference between absolute maximum and minimum + + PeakTrough: half difference between maximum and minimum + on a half cycle + + + + + Define how to combine the amplitudes measured on both + horizontals components: + + min: take the minimum + + max: take the maxium + + avgerage: form the average + + geometric_mean: form the geometric mean + + + + + + + + Parameters for computing mb_Lg magnitudes from mb_Lg amplitudes. + + + + Considered distance measure between source and receiver. + Possible values are + + hypocentral: hypocentral distance + + epicentral: epicentral distance + + + + + The minimum distance for computing magnitudes from amplitudes. + Negative values deactivate the check. + + + + + The maximum distance for computing magnitudes from amplitudes. + Negative values deactivate the check. + + + + + The maximum depth up to which magnitudes are computed. + + + + + Type of magnitude calibration formula to be considered. + The calibration parameters are considered accordingly. + Currently supported are + + "parametric": consider parameters of parametric + configuration in parametric section + + "A0": consider parameters of non-parametric + configuration in A0 section. + + + + + Parameters for A0, non-parametric magnitude calibration. + Considered if magnitude.mb_Lg.calibrationType = "A0". + + + + The non-parametric calibration function log10(A0). + + Format: any list of distance-value pairs separated by + comma. Values within pairs are separated by colon. + + Example: "0:-1.3,60:-2.8,100:-3.0,400:-4.5,1000:-5.85" + specifies 4 distance intervals from + 0...60, 60...100, 100...400 and 400...1000 km distance. + Within these intervals log10(A0) is interpolated linearly + between -1.3...-2.8, -2.8...-3.0, -3.0...-4.5 and -4.5...-5.8, + respectively. + + Note: The first and last distance samples limit the + maximum distance range for computing MLv. + + + + + + Parameters for parametric magnitude calibration: + mb_Lg = log10(A) + c0 * log10(r) + c1 * r + c2 + Considered if magnitude.mb_Lg.calibrationType = "parametric". + + + + Station correction. This is the calibration value 'c0' + applied in the magnitude calibration formula for the north component or + the closest component to the north. If amplitude is measured in the north + component, c0(station_horizontal_component) = c0_first + + mb_Lg = log10(A) + c0 * log10(r) + c1 * r + c2 + + + + + The calibration value 'c1' applied in the magnitude + calibration formula + + mb_Lg = log10(A) + c0 * log10(r) + c1 * r + c2 + + + + + The calibration value 'c2' applied in the + magnitude calibration formula + + mb_Lg = log10(A) + c0 * log10(r) + c1 * r + c2 + + + + + + + + diff --git a/libs/seiscomp/processing/magnitudes/descriptions/global_mlcan.rst b/libs/seiscomp/processing/magnitudes/descriptions/global_mlcan.rst new file mode 100644 index 00000000..db9808c5 --- /dev/null +++ b/libs/seiscomp/processing/magnitudes/descriptions/global_mlcan.rst @@ -0,0 +1,210 @@ +MLcan is the custom magnitude for local events provided by the mlcan plugin. +The implementation is based on specifications by Spanish National Geographic Institute, +based on the corrections obtained by Mezcua and Rueda (2022): +https://doi.org/10.1007/s00445-022-01553-9 +This magnitude changes the c0 station correction for the c0_first +and c0_second, which corresponds to channel correction for a given station. + +The MLcan magnitude is very similar to the original :ref:`ML`, +except that by default + +* Amplitude pre-filtering is applied. +* A parametric :ref:`magnitude calibration ` function + applies. +* Hypocentral distance is used. + +Regionalization of magnitude computation is provided through global module +configuration. +Configuration of global bindings provides additional flexibility: + +* Amplitudes can be pre-filtered before simulating a :term:`Wood-Anderson seismometer` + (:confval:`amplitudes.MLcan.preFilter`), +* Wood-Anderson simulation is optional: + :confval:`amplitudes.MLcan.applyWoodAnderson`, +* Measured amplitudes can be scaled accounting for expected unit + (:confval:`amplitudes.MLcan.amplitudeScale`), +* A parametric or A0-based non-parametric + :ref:`magnitude calibration ` + function can be applied as controlled by + :confval:`magnitudes.MLcan.calibrationType`. +* Consider either hypocentral or epicentral distance for computing magnitudes + (:confval:`magnitudes.MLcan.distMode`). + +General (default) conditions apply, all values are configurable: + +* Amplitude pre-filtering: :ref:`BW(3,0.5,12) `. +* Amplitude unit in SeisComP: **millimeter** (mm) or as considered by the + configured calibration parameters. +* Optional amplitude scaling: 1. +* :term:`Wood-Anderson seismometer` simulation: yes. +* Time window: 150 s by :ref:`scautopick` or distance dependent + with :math:`endTime = min(R / 4, 150)`, configurable. +* Distance type: hypocentral. +* Distance range: 0 - 8 deg, measurements beyond 8 deg will be + strictly ignored. +* Depth range: <= 80 km. +* Magnitude calibration: parametric. + + +Amplitudes +---------- + +MLcan amplitudes can be measured automatically by :ref:`scautopick` or :ref:`scamp` +or interactively by :ref:`scolv` very similarly to the original :ref:`ML`, +except that they can be pre-filtered and simulation of a :term:`Wood-Anderson seismometer` +is optional: :confval:`amplitudes.MLcan.preFilter`, +:confval:`amplitudes.MLcan.applyWoodAnderson`. +By default amplitudes are measured on both horizontal components where the +absolute maxima are taken. They are combined to a final measured amplitude by +taking the mean. The methods for measuring and combining the amplitudes are +configurable: +:confval:`amplitudes.MLcan.measureType`, :confval:`amplitudes.MLcan.combiner`. + +The Wood-Anderson simulation will convert input velocity data to ground +displacement in mm. The input data may be of a different unit after applying +:confval:`amplitudes.MLcan.preFilter`, e.g. when integration is applied, and / or +when Wood-Anderson simulation is disabled. Configure +:confval:`amplitudes.MLcan.amplitudeScale` for converting the unit of the +processed data to the unit expected by the +:ref:`station magnitude calibration ` for the measured +amplitude. + +.. note:: + + For comparing MLcan amplitudes with :ref:`ML amplitudes ` set the + global bindings parameters :: + + amplitudes.MLcan.preFilter = "" + amplitudes.MLcan.combiner = average + + +.. _mlc_station_magnitude: + +Station Magnitudes +------------------ + +Station magnitudes are computed from measured amplitudes automatically by +:ref:`scmag` +or interactively by :ref:`scolv`. By global bindings configuration MLcan considers + +* Hypocentral (default) or epicentral distance: :confval:`magnitudes.MLcan.distMode`. +* Distance range: :confval:`magnitudes.MLcan.minDist`, :confval:`magnitudes.MLcan.maxDist`. +* Events with depth up to :confval:`magnitudes.MLcan.maxDepth`. +* Parametric or non-parametric calibration functions + + * parametric when :confval:`magnitudes.MLcan.calibrationType` = "parametric"`: + + .. math:: + + MLcan = \log_{10}(A) + c_3 * \log_{10}(r/c_5) + c_2 * (r + c_4) + c_1 + c_0(station) + + where + + * *A*: displacement amplitude measured in unit of mm or as per configuration + * *r*: hypocentral (default) or epicentral distance + * *c1*, *c2*, *c3*, *c4*, *c5*: general calibration parameters + * *c0*: channel-station-specific correction + * *r*: Hypocentral (default) or epicentral distance as configured by + :confval:`magnitudes.MLcan.distMode`. + + c0, is the correction for the component used in amplitude measurement. In case of average or + geometric average measurement, the mean correction for c0_first and c0_second will be used. + c0_first is the north component or the closest to the nort. c0 second is the east component + or the closest to the east. + + * A0-based non-parametric when :confval:`magnitudes.MLcan.calibrationType` = "A0"`: + + .. math:: + + MLcan = \log_{10}(A) - \log_{10}(A_0) + + where + + * :math:`log_{10}(A_0)`: distance-dependent correction value. Read + :ref:`global_mlv` for the details. + +.. note:: + + The magnitude calibration function can regionalized by adjusting global module + configuration parameters in MLcan region profiles of + :confval:`magnitudes.MLcan.region.*` and in a *MLcan* Magnitude type profile e.g. + in :file:`global.cfg`. + + +Network Magnitude +----------------- + +The network magnitude is computed from station magnitudes automatically by +:ref:`scmag` or interactively by :ref:`scolv`. +Originally the median was computed from all station MLcan to form the +:term:`network magnitude` MLcan. Here, the trimmed mean is applied. Outliers +beyond the outer 12.5% percentiles are removed before forming the mean. The +method can be adjusted in :ref:`scmag` by :confval:`magnitudes.average`. + + +Examples +-------- + +The flexibility of the amplitude and magnitude processing allows to apply MLcan +in various use cases, e.g. + +* **Default:** Pre-filtered and gain-corrected amplitudes, Wood-Anderson + corrected and measured in mm for Southwestern Germany, :cite:t:`stange-2006`: + + .. math:: + + MLcan = \log_{10}(A) + 1.11 * \log_{10}(r) + 0.00095 * r + 0.69 + c_0 + +* Wood-Anderson-corrected displacement amplitudes measured in mm for + Southern California, :cite:t:`hutton-1987`: + + .. math:: + + MLcan = \log_{10}(A) + 1.110 * \log_{10}(r / 100) + 0.00189 * (r - 100) + 3.0 + +* Pre-filtered velocity amplitudes in units of mym/s (requiring to set + :confval:`amplitudes.MLcan.amplitudeScale`), no Wood-Anderson correction, + for West Bohemia, e.g. :cite:t:`hiemer-2012`: + + .. math:: + + MLcan = \log_{10}(A) - log_{10}(2\Pi) + 2.1 * \log_{10}(r) - 1.7 + c_0 + +.. figure:: media/magnitude-calibrations_MLcan_s_MLcan_hb.png + :align: center + :width: 18cm + + MLcan magnitudes for measured amplitude of 1 mm with default magnitude + calibration (*MLcan_s*, :cite:t:`stange-2006`) and calibration values for Southern + California (*MLcan_hb*, :cite:t:`hutton-1987`). + + +Setup +===== + +#. **Set the configuration and calibration parameters** in the global bindings + similar + to :ref:`global_ml`. Instead of configuring lots of global bindings profiles + or station bindings one line per parameter can be added to the global module + configuration (:file:`global.cfg`) which takes the form + + .. code-block:: params + + module.trunk.NET.STA.amplitude.MLcan.preFilter = value + module.trunk.NET.STA.magnitude.MLcan.parametric.c0 = value + +#. Add MLcan to the list of default amplitudes and magnitudes if MLcan is to be + computed by automatic modules, e.g. of :ref:`scamp`, :ref:`scmag`. +#. Configure :ref:`scmag` (:confval:`magnitudes.average` in :file:`scmag.cfg`) + for choosing the method to form the + network magnitude from station magnitudes, e.g. + + .. code-block:: params + + magnitudes.average = MLcan:median + +#. Add MLcan to the list of magnitudes preferred by :ref:`scevent` + (:confval:`eventAssociation.magTypes` in :file:`scevent.cfg`) in order to let + MLcan become the preferred magnitude. +#. Set defaults/visibility of MLcan in :term:`GUI` modules, e.g. :ref:`scolv` + or :ref:`scesv`. diff --git a/libs/seiscomp/processing/magnitudes/descriptions/global_mlcan.xml b/libs/seiscomp/processing/magnitudes/descriptions/global_mlcan.xml new file mode 100644 index 00000000..04100bae --- /dev/null +++ b/libs/seiscomp/processing/magnitudes/descriptions/global_mlcan.xml @@ -0,0 +1,322 @@ + + + + global + + Custom magnitude for local events measured on horizontal components with + correction for each component. If the amplitude is measured as "mean" or + "geometric_mean", the correction will be the mean of each correction. If + no correction is specified, a c0_first and second will be 0 + + + + + + + Considered distance measure between source and receiver. + Possible values are + + hypocentral: hypocentral distance + + epicentral: epicentral distance + + + + + Type of magnitude calibration formula to be considered. + The calibration parameters are considered accordingly. + Currently supported are + + "parametric": consider parameters of parametric + configuration in parametric section + + "A0": consider parameters of non-parametric + configuration in A0 section. + + + + + Parameters for A0, non-parametric magnitude calibration. + + + + Overrides the calibration function log10(A0) + for computing MLcan per region. See logA0 + description in the bindings. + + + + + + + Parameters for parametric magnitude calibration: + MLcan = log10(A) + c3 * log10(r/c5) + c2 * (r + c4) + c1 + c0(station_horizontal_component) + where c0(station_horizontal_component) is c0_first, c0_second or their mean, according to which + component has been used to measure the amplitude + + + + Overrides the calibration parameter c0_first + for computing MLcan per region. See c0_first + description in the bindings. + + + + + Overrides the calibration parameter c0_second + for computing MLcan per region. See c0_second + description in the bindings. + + + + + Overrides the calibration parameter c1 + for computing MLcan per region. See c1 + description in the bindings. + + + + + Overrides the calibration parameter c2 + for computing MLcan per region. See c2 + description in the bindings. + + + + + Overrides the calibration parameter c3 + for computing MLcan per region. See c3 + description in the bindings. + + + + + Overrides the calibration parameter c4 + for computing MLcan per region. See c4 + description in the bindings. + + + + + Overrides the calibration parameter c5 + for computing MLcan per region. See c5 + description in the bindings. + + + + + + + + + + Custom magnitude for local events measured on horizontal components + + + + + + Parameters for measuring MLcan amplitudes. Add more parameters + by adding an amplitude type profile 'MLcan', + + + + The filter applied to raw records before applying + Wood-Anderson simulation. + + + + + Applying Wood-Anderson simulation. To achieve displacement + records without WA simulation, an integration filter can + be applied with the pre-filter. + + + + + Scaling value multiplied to the measured amplitudes to + match the amplitude units expected by the magnitude + calibration function. + + Expected amplitudes are + in units of mym but actual amplitudes provided from + Wood-Anderson-corrected seismograms are in units of mm: + amplitudeScale = 1000. + + If data are not corrected for WA, measured amplitudes + take the unit of gain-corrected data considering the + preFilter: + amplitudeScale converts between units of measured and + excpected amplitude. + + + + + Type for measuring amplitudes. Available: + + AbsMax: absolute maximum + + MinMax: half difference between absolute maximum and minimum + + PeakTrough: half difference between maximum and minimum + on a half cycle + + + + + Define how to combine the amplitudes measured on both + horizontals components: + + min: take the minimum + + max: take the maxium + + avgerage: form the average + + geometric_mean: form the geometric mean + + + + + + + + Parameters for computing MLcan magnitudes from MLcan amplitudes. + + + + Considered distance measure between source and receiver. + Possible values are + + hypocentral: hypocentral distance + + epicentral: epicentral distance + + + + + The minimum distance for computing magnitudes from amplitudes. + Negative values deactivate the check. + + + + + The maximum distance for computing magnitudes from amplitudes. + Negative values deactivate the check. + + + + + The maximum depth up to which magnitudes are computed. + + + + + Type of magnitude calibration formula to be considered. + The calibration parameters are considered accordingly. + Currently supported are + + "parametric": consider parameters of parametric + configuration in parametric section + + "A0": consider parameters of non-parametric + configuration in A0 section. + + + + + Parameters for A0, non-parametric magnitude calibration. + Considered if magnitude.MLcan.calibrationType = "A0". + + + + The non-parametric calibration function log10(A0). + + Format: any list of distance-value pairs separated by + comma. Values within pairs are separated by colon. + + Example: "0:-1.3,60:-2.8,100:-3.0,400:-4.5,1000:-5.85" + specifies 4 distance intervals from + 0...60, 60...100, 100...400 and 400...1000 km distance. + Within these intervals log10(A0) is interpolated linearly + between -1.3...-2.8, -2.8...-3.0, -3.0...-4.5 and -4.5...-5.8, + respectively. + + Note: The first and last distance samples limit the + maximum distance range for computing MLv. + + + + + + Parameters for parametric magnitude calibration: + MLcan = log10(A) + c3 * log10(r/c5) + c2 * (r + c4) + c1 + c0(station_horizontal_component) + + Considered if magnitude.MLcan.calibrationType = "parametric". + + + + Station correction. This is the calibration value 'c0' + applied in the magnitude calibration formula for the north component or + the closest component to the north. If amplitude is measured in the north + component, c0(station_horizontal_component) = c0_first + + MLcan = c0(station_horizontal_component) + c1 + c2 * (r + c4) + c3 * log(r/c5) + log10(A) + + + + + Station correction. This is the calibration value 'c0' + applied in the magnitude calibration formula for the east component or + the closest component to the east. If amplitude is measured in the north + component, c0(station_horizontal_component) = c0_second + + MLcan = c0(station_horizontal_component) + c1 + c2 * (r + c4) + c3 * log(r/c5) + log10(A) + + + + + The calibration value 'c1' applied in the magnitude + calibration formula + + MLcan = c0(station_horizontal_component) + c1 + c2 * (r + c4) + c3 * log(r/c5) + log10(A) + + + + + The calibration value 'c2' applied in the + magnitude calibration formula + + MLcan = c0(station_horizontal_component) + c1 + c2 * (r + c4) + c3 * log(r/c5) + log10(A) + + + + + The calibration value 'c3' applied in the + magnitude calibration formula + + MLcan = c0(station_horizontal_component) + c1 + c2 * (r + c4) + c3 * log(r/c5) + log10(A) + + + + + The calibration value 'c4' applied in the + magnitude calibration formula + + MLcan = c0(station_horizontal_component) + c1 + c2 * (r + c4) + c3 * log(r/c5) + log10(A) + + + + + The calibration value 'c4' applied in the + magnitude calibration formula + + MLcan = c0(station_horizontal_component) + c1 + c2 * (r + c4) + c3 * log(r/c5) + log10(A) + + + + + + + + diff --git a/libs/seiscomp/processing/magnitudes/mb_Lg.cpp b/libs/seiscomp/processing/magnitudes/mb_Lg.cpp new file mode 100644 index 00000000..353d8842 --- /dev/null +++ b/libs/seiscomp/processing/magnitudes/mb_Lg.cpp @@ -0,0 +1,358 @@ +/* ################################################################################ +* # Copyright (C) 2024 by IGN Spain # +* # # +* # author: J. Barco, E. Suarez # +* # email: jbarco@transportes.gob.es, eadiaz@transportes.gob.es # +* # last modified: 2024-03-20 # +* # # +* # This program is free software; you can redistribute it and/or modify # +* # it under the terms of the GNU General Public License as published by # +* # the Free Software Foundation; either version 2 of the License, or # +* # (at your option) any later version. # +* # # +* # This program is distributed in the hope that it will be useful, # +* # but WITHOUT ANY WARRANTY; without even the implied warranty of # +* # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # +* # GNU General Public License for more details. # +* # # +* # You should have received a copy of the GNU General Public License # +* # along with this program; if not, write to the # +* # Free Software Foundation, Inc., # +* # 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. # +* ################################################################################ */ + + +#define SEISCOMP_COMPONENT mb_Lg + +#include +#include +#include +#include +#include + +using namespace std; + + +namespace Seiscomp { +namespace Processing { + + +namespace { + +std::string ExpectedAmplitudeUnit = "mm"; + +DEFINE_SMARTPOINTER(ExtraLocale); +class ExtraLocale : public Core::BaseObject { + public: + // general magnitude parameters + OPT(string) distanceMode; + OPT(string) calibrationType; + // parametric coefficients + OPT(double) c0; + OPT(double) c1; + OPT(double) c2; + // A0, non-parametric coefficients + OPT(LogA0) logA0; +}; + +} + + +REGISTER_MAGNITUDEPROCESSOR(MagnitudeProcessor_mb_Lg, "mb_Lg"); +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +MagnitudeProcessor_mb_Lg::MagnitudeProcessor_mb_Lg() +: Processing::MagnitudeProcessor("mb_Lg") {} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +std::string MagnitudeProcessor_mb_Lg::amplitudeType() const { + return MagnitudeProcessor::amplitudeType(); +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +bool MagnitudeProcessor_mb_Lg::setup(const Settings &settings) { + if ( !MagnitudeProcessor::setup(settings) ) { + return false; + } + + try {_maxDepth = settings.getDouble("magnitudes." + _type + ".maxDepth"); } + catch ( ... ) {} + + // distance constraints + try {_distanceMode = settings.getString("magnitudes." + _type + ".distMode"); } + catch ( ... ) {} + try {_minDistanceKm = Math::Geo::deg2km(settings.getDouble("magnitudes." + _type + ".minDist")); } + catch ( ... ) {} + try {_maxDistanceKm = Math::Geo::deg2km(settings.getDouble("magnitudes." + _type + ".maxDist")); } + catch ( ... ) {} + + // calibration function + try {_calibrationType = settings.getString("magnitudes." + _type + ".calibrationType"); } + catch ( ... ) {} + + if ( (_calibrationType != "A0") && (_calibrationType != "parametric") ) { + SEISCOMP_ERROR("%s: unrecognized calibration type: %s", _type.c_str(), + _calibrationType.c_str()); + return false; + } + + // parametric calibration function + try { _c0 = settings.getDouble("magnitudes." + _type + ".parametric.c0"); } + catch ( ... ) {} + try { _c1 = settings.getDouble("magnitudes." + _type + ".parametric.c1"); } + catch ( ... ) {} + try { _c2 = settings.getDouble("magnitudes." + _type + ".parametric.c2"); } + catch ( ... ) {} + + // A0, non-parametric calibration function + std::string defLogA0 = "0:-1.3,60:-2.8,100:-3.0,400:-4.5,1000:-5.85"; + bool logA0Default = true; + + try { + defLogA0 = settings.getString("magnitudes." + _type + ".A0.logA0"); + logA0Default = false; + } + catch ( ... ) {} + + // set distance-dependent intervals + if ( !_logA0.set(defLogA0) ) { + SEISCOMP_ERROR("%s: incorrect correction term log(A0): %s", _type.c_str(), + defLogA0.c_str()); + return false; + } + + if ( !logA0Default ) { + SEISCOMP_DEBUG("%s.%s: %s: logA0 from bindings = %s", + settings.networkCode, settings.stationCode, + type(), Core::toString(_logA0)); + } + + return true; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +bool MagnitudeProcessor_mb_Lg::initLocale(Locale *locale, + const Settings &settings, + const string &configPrefix) { + const Seiscomp::Config::Config *cfg = settings.localConfiguration; + ExtraLocalePtr extra = new ExtraLocale; + + // general + try { + extra->distanceMode = cfg->getString(configPrefix + "distMode"); + } + catch ( ... ) {} + try { + extra->calibrationType = cfg->getString(configPrefix + "calibrationType"); + } + catch ( ... ) {} + + // parametric + try { + // extra->c0 = cfg->getDouble(configPrefix + ".parametric.c0"); + extra->c0 = cfg->getDouble(configPrefix + "parametric.c0"); + } + catch ( ... ) {} + + try { + extra->c1 = cfg->getDouble(configPrefix + "parametric.c1"); + } + catch ( ... ) {} + + try { + extra->c2 = cfg->getDouble(configPrefix + "parametric.c2"); + } + catch ( ... ) {} + + + // A0, non-parametric + try { + auto logA0 = cfg->getStrings(configPrefix + "A0.logA0"); + if ( !logA0.empty() ) { + // If the first item contains a comma, then maybe it has been configured + // with double quotes. Raise an error. + if ( logA0[0].find(',') != string::npos ) { + SEISCOMP_ERROR("%sA0.logA0[0] contains a comma. Are the coefficients " + "enclosed with double quotes in the configuration?", + configPrefix); + return false; + } + + if ( logA0[0].find(';') != string::npos ) { + SEISCOMP_ERROR("%sA0.logA0 = %s contains semicolon. Supported format:" + " distance1:correction1,distance2:correction2, ...", + configPrefix, Core::toString(logA0).c_str()); + return false; + } + + extra->logA0 = LogA0(); + if ( !extra->logA0->set(logA0) ) { + SEISCOMP_ERROR("%s@%s: incorrect correction term log(A0)", + _type.c_str(), locale->name.c_str()); + return false; + } + } + } + catch ( ... ) {} + + locale->extra = extra; + + return true; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +MagnitudeProcessor::Status MagnitudeProcessor_mb_Lg::computeMagnitude( + double amplitude, const std::string &unit, + double, double, double delta, double depth, + const DataModel::Origin *, const DataModel::SensorLocation *receiver, + const DataModel::Amplitude *, + const Locale *locale, + double &value) { + + if ( amplitude <= 0 ) { + return AmplitudeOutOfRange; + } + + auto distanceMode = _distanceMode; + auto calibrationType = _calibrationType; + auto minimumDistanceKm = _minDistanceKm; + auto maximumDistanceKm = _maxDistanceKm; + auto maximumDepth = _maxDepth; + + ExtraLocale *extra = nullptr; + if ( locale ) { + extra = static_cast(locale->extra.get()); + if ( extra ) { + if ( extra->distanceMode ) { + distanceMode = *extra->distanceMode; + } + + if ( extra->calibrationType ) { + calibrationType = *extra->calibrationType; + } + } + + if ( locale->minimumDistance ) { + minimumDistanceKm = Math::Geo::deg2km(*locale->minimumDistance); + } + + if ( locale->maximumDistance ) { + maximumDistanceKm = Math::Geo::deg2km(*locale->maximumDistance); + } + } + else { + + SEISCOMP_DEBUG(" + maximum depth: %.3f km", maximumDepth); + if ( depth > maximumDepth ) { + return DepthOutOfRange; + } + SEISCOMP_DEBUG(" + minimum distance: %.3f km", minimumDistanceKm); + SEISCOMP_DEBUG(" + maximum distance: %.3f km", maximumDistanceKm); + } + + SEISCOMP_DEBUG(" + distance type: %s", distanceMode.c_str()); + + double hDistanceKm = Math::Geo::deg2km(delta); + double vDistanceKm = 0; + + if ( !receiver && distanceMode == "hypocentral" ) { + return MetaDataRequired; + } + + if ( receiver ) { + vDistanceKm = receiver->elevation() / 1000. + depth; + } + + double distanceKm; + if ( distanceMode == "hypocentral" ) { + distanceKm = sqrt(pow(hDistanceKm, 2) + pow(vDistanceKm, 2)); + } + else { + distanceKm = hDistanceKm; + } + + SEISCOMP_DEBUG(" + considered distance to station: %.3f km", distanceKm); + + if ( minimumDistanceKm >= 0 && distanceKm < minimumDistanceKm ) { + return DistanceOutOfRange; + } + + if ( maximumDistanceKm >= 0 && distanceKm > maximumDistanceKm ) { + return DistanceOutOfRange; + } + + if ( !convertAmplitude(amplitude, unit, ExpectedAmplitudeUnit) ) { + return InvalidAmplitudeUnit; + } + + SEISCOMP_DEBUG(" + %s distance: %.3f km (horiz. %.3f vert. %.3f)", + distanceMode.c_str(), distanceKm, hDistanceKm, vDistanceKm); + + double correction; + SEISCOMP_DEBUG(" + calibration type: %s", calibrationType.c_str()); + + if ( calibrationType == "parametric" ) { + // parametric calibration function + //double c0 = 0; + // parametric calibration function + auto c0 = (extra and extra->c0) ? *extra->c0 : _c0; + auto c1 = (extra and extra->c1) ? *extra->c1 : _c1; + auto c2 = (extra and extra->c2) ? *extra->c2 : _c2; + + SEISCOMP_DEBUG(" + c0: %.5f", c0); + SEISCOMP_DEBUG(" + c1: %.5f", c1); + SEISCOMP_DEBUG(" + c2: %.5f", c2); + + correction = c0*log10(distanceKm) + c1 * (distanceKm) + c2; + double PI = 3.14159265359; + value = log10(amplitude/(2*PI)) + correction; + } + else if ( calibrationType == "A0" ) { + // A0, non-parametric calibration function + try { + correction = -1.0 * (extra and extra->logA0 ? extra->logA0->at(distanceKm) : _logA0.at(distanceKm)); + + SEISCOMP_DEBUG(" + -log10(A0) at %.3f km: %.5f", distanceKm, correction); + value = log10(amplitude) + correction; + } + catch ( std::out_of_range & ) { + return DistanceOutOfRange; + } + } + + else { + return IncompleteConfiguration; + } + + SEISCOMP_DEBUG(" + amplitude: %.5f, magnitude: %.3f", amplitude, value); + + return OK; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +} +} diff --git a/libs/seiscomp/processing/magnitudes/mb_Lg.h b/libs/seiscomp/processing/magnitudes/mb_Lg.h new file mode 100644 index 00000000..4be97e94 --- /dev/null +++ b/libs/seiscomp/processing/magnitudes/mb_Lg.h @@ -0,0 +1,86 @@ +/* ################################################################################ +* # Copyright (C) 2024 by IGN Spain # +* # # +* # author: J. Barco, E. Suarez # +* # email: jbarco@transportes.gob.es , eadiaz@transportes.gob.es # +* # last modified: 2024-03-20 # +* # based on Lopez, C. (2008). Nuevas formulas de magnitud para la # +* # Peninsula Iberica y su entorno. Trabajo de investigación del # +* # Master en Geofisica y Meteorologia. Departamento de Física de la Tierra, # +* # Astronomia y Astrofisica. Universidad Complutense de Madrid. Madrid. # +* # # +* # This program is free software; you can redistribute it and/or modify # +* # it under the terms of the GNU General Public License as published by # +* # the Free Software Foundation; either version 2 of the License, or # +* # (at your option) any later version. # +* # # +* # This program is distributed in the hope that it will be useful, # +* # but WITHOUT ANY WARRANTY; without even the implied warranty of # +* # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # +* # GNU General Public License for more details. # +* # # +* # You should have received a copy of the GNU General Public License # +* # along with this program; if not, write to the # +* # Free Software Foundation, Inc., # +* # 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. # +* ################################################################################ */ + + +#ifndef SEISCOMP_PROCESSING_MAGNITUDEPROCESSOR_mb_Lg_H +#define SEISCOMP_PROCESSING_MAGNITUDEPROCESSOR_mb_Lg_H + + +#include +#include +#include + + +namespace Seiscomp { +namespace Processing { + + +class SC_SYSTEM_CLIENT_API MagnitudeProcessor_mb_Lg : public MagnitudeProcessor { + public: + MagnitudeProcessor_mb_Lg(); + + + public: + bool setup(const Settings &settings) override; + + std::string amplitudeType() const override; + + + protected: + Status computeMagnitude(double amplitude, const std::string &unit, + double period, double snr, + double delta, double depth, + const DataModel::Origin *hypocenter, + const DataModel::SensorLocation *receiver, + const DataModel::Amplitude *, + const Locale *, + double &value) override; + + protected: + bool initLocale(Locale *locale, const Settings &settings, + const std::string &configPrefix) override; + + private: + double _minDistanceKm{-1.0}; + double _maxDistanceKm{8.0 * KM_OF_DEGREE}; + double _maxDepth{80.0}; + std::string _distanceMode{"hypocentral"}; + std::string _calibrationType{"parametric"}; + // parameters for parametric magnitude calibration + double _c0{1.17}; + double _c1{0.0012}; + double _c2{0.67}; + // parameters for non-parametric magnitude calibration + LogA0 _logA0; +}; + + +} +} + + +#endif diff --git a/libs/seiscomp/processing/magnitudes/mb_VC.cpp b/libs/seiscomp/processing/magnitudes/mb_VC.cpp new file mode 100644 index 00000000..a6852ade --- /dev/null +++ b/libs/seiscomp/processing/magnitudes/mb_VC.cpp @@ -0,0 +1,141 @@ +/* ################################################################################ +* # Copyright (C) 2024 by IGN Spain # +* # # +* # author: J. Barco, E. Suarez # +* # email: jbarco@transportes.gob.es , eadiaz@transportes.gob.es # +* # last modified: 2024-03-20 # +* # # +* # This program is free software; you can redistribute it and/or modify # +* # it under the terms of the GNU General Public License as published by # +* # the Free Software Foundation; either version 2 of the License, or # +* # (at your option) any later version. # +* # # +* # This program is distributed in the hope that it will be useful, # +* # but WITHOUT ANY WARRANTY; without even the implied warranty of # +* # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # +* # GNU General Public License for more details. # +* # # +* # You should have received a copy of the GNU General Public License # +* # along with this program; if not, write to the # +* # Free Software Foundation, Inc., # +* # 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. # +* ################################################################################ */ + + + +#include +#include +#include + +namespace Seiscomp { +namespace Processing { + + +namespace { + +std::string ExpectedAmplitudeUnit = "nm/s"; + +} + + +IMPLEMENT_SC_CLASS_DERIVED(MagnitudeProcessor_mB_VC, MagnitudeProcessor, "MagnitudeProcessor_mB_VC"); +REGISTER_MAGNITUDEPROCESSOR(MagnitudeProcessor_mB_VC, "mB_VC"); +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +bool MagnitudeProcessor_mB_VC::MagnitudeProcessor_mB_VC::setup(const Settings &settings) { + if ( !MagnitudeProcessor::setup(settings) ) + return false; + + minDistanceDeg = 5.0; // default minimum distance + maxDistanceDeg = 105.0; // default maximum distance + + // distance range in degree + try { minDistanceDeg = settings.getDouble("magnitudes." + type() + ".minDist"); } + catch ( ... ) {} + + try { maxDistanceDeg = settings.getDouble("magnitudes." + type() + ".maxDist"); } + catch ( ... ) {} + + return true; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +MagnitudeProcessor_mB_VC::MagnitudeProcessor_mB_VC() + : MagnitudeProcessor("mB_VC") {} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +MagnitudeProcessor_mB_VC::MagnitudeProcessor_mB_VC(const std::string& type) + : MagnitudeProcessor(type) {} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +MagnitudeProcessor::Status MagnitudeProcessor_mB_VC::computeMagnitude( + double amplitude, const std::string &unit, + double, double, + double delta, double depth, + const DataModel::Origin *, const DataModel::SensorLocation *, + const DataModel::Amplitude *, const Locale *, double &value) { + // Clip depth to 0 + if ( depth < 0 ) depth = 0; + + if ( delta < minDistanceDeg || delta > maxDistanceDeg ) + return DistanceOutOfRange; + + if ( amplitude <= 0 ) + return AmplitudeOutOfRange; + + if ( !convertAmplitude(amplitude, unit, ExpectedAmplitudeUnit) ) + return InvalidAmplitudeUnit; + + bool status = Magnitudes::compute_mb_VC(amplitude*1.E-3, 2*M_PI, delta, depth+1, &value); + value -= 0.14; // HACK until we have an optimal calibration function + return status ? OK : Error; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +MagnitudeProcessor::Status MagnitudeProcessor_mB_VC::estimateMw( + const Config::Config *, + double magnitude, + double &Mw_estimate, + double &Mw_stdError) +{ + // Relationships Between Mw and Other Earthquake Size Parameters in the Spanish IGN Seismic + // Catalog + // Cabanas et al. (2015) + // DOI 10.1007/s00024-014-1025-2 + const double a=1.213, b=-1.528; + Mw_estimate = a * magnitude + b; + + Mw_stdError = 0.4; + + return OK; +} +// <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<< + + + + +// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> +} +} diff --git a/libs/seiscomp/processing/magnitudes/mb_VC.h b/libs/seiscomp/processing/magnitudes/mb_VC.h new file mode 100644 index 00000000..5fee433f --- /dev/null +++ b/libs/seiscomp/processing/magnitudes/mb_VC.h @@ -0,0 +1,69 @@ +/* ################################################################################ +* # Copyright (C) 2024 by IGN Spain # +* # # +* # author: J. Barco, E. Suarez # +* # email: jbarco@transportes.gob.es , eadiaz@transportes.gob.es # +* # last modified: 2024-03-20 # +* # # +* # This program is free software; you can redistribute it and/or modify # +* # it under the terms of the GNU General Public License as published by # +* # the Free Software Foundation; either version 2 of the License, or # +* # (at your option) any later version. # +* # # +* # This program is distributed in the hope that it will be useful, # +* # but WITHOUT ANY WARRANTY; without even the implied warranty of # +* # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # +* # GNU General Public License for more details. # +* # # +* # You should have received a copy of the GNU General Public License # +* # along with this program; if not, write to the # +* # Free Software Foundation, Inc., # +* # 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. # +* ################################################################################ */ + + +#ifndef SEISCOMP_PROCESSING_MAGNITUDEPROCESSOR_M_B_VC_H +#define SEISCOMP_PROCESSING_MAGNITUDEPROCESSOR_M_B_VC_H + + +#include + + +namespace Seiscomp { +namespace Processing { + + +class SC_SYSTEM_CLIENT_API MagnitudeProcessor_mB_VC : public MagnitudeProcessor { + DECLARE_SC_CLASS(MagnitudeProcessor_mB_VC) + + public: + MagnitudeProcessor_mB_VC(); + MagnitudeProcessor_mB_VC(const std::string& type); + + + bool setup(const Settings &settings) override; + + Status computeMagnitude(double amplitude, const std::string &unit, + double period, double snr, + double delta, double depth, + const DataModel::Origin *hypocenter, + const DataModel::SensorLocation *receiver, + const DataModel::Amplitude *, + const Locale *, + double &value) override; + + Status estimateMw(const Config::Config *config, + double magnitude, double &Mw_estimate, + double &Mw_stdError) override; + + private: + double minDistanceDeg; + double maxDistanceDeg; +}; + + +} +} + + +#endif diff --git a/libs/seiscomp/seismology/CMakeLists.txt b/libs/seiscomp/seismology/CMakeLists.txt index 7dfa7c4b..46d9b9b3 100644 --- a/libs/seiscomp/seismology/CMakeLists.txt +++ b/libs/seiscomp/seismology/CMakeLists.txt @@ -5,6 +5,7 @@ SET(SEISMOLOGY_SOURCES Mjma.cpp Mwp.cpp ttt.cpp + mb_VC.cpp ) SET(SEISMOLOGY_HEADERS diff --git a/libs/seiscomp/seismology/magnitudes.h b/libs/seiscomp/seismology/magnitudes.h index 2dd65638..35ac04d4 100644 --- a/libs/seiscomp/seismology/magnitudes.h +++ b/libs/seiscomp/seismology/magnitudes.h @@ -97,6 +97,23 @@ compute_mb_fromVelocity( double depth, // in kilometers double *mag); // resulting magnitude +SC_SYSTEM_CORE_API +bool +compute_mb_VC( + double amplitude, // in micrometers + double period, // in seconds + double delta, // in degrees + double depth, // in kilometers + double *mag); // resulting magnitude + + +SC_SYSTEM_CORE_API +bool +compute_mb_VC_fromVelocity( + double amplitude, // in micrometers/second + double delta, // in degrees + double depth, // in kilometers + double *mag); // resulting magnitude SC_SYSTEM_CORE_API bool diff --git a/libs/seiscomp/seismology/mb_VC.cpp b/libs/seiscomp/seismology/mb_VC.cpp new file mode 100644 index 00000000..43e548ab --- /dev/null +++ b/libs/seiscomp/seismology/mb_VC.cpp @@ -0,0 +1,210 @@ +/* ################################################################################ +* # Copyright (C) 2024 by IGN Spain # +* # # +* # author: J. Barco, E. Suarez # +* # email: jbarco@transportes.gob.es , eadiaz@transportes.gob.es # +* # last modified: 2024-03-20 # +* # # +* # This program is free software; you can redistribute it and/or modify # +* # it under the terms of the GNU General Public License as published by # +* # the Free Software Foundation; either version 2 of the License, or # +* # (at your option) any later version. # +* # # +* # This program is distributed in the hope that it will be useful, # +* # but WITHOUT ANY WARRANTY; without even the implied warranty of # +* # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the # +* # GNU General Public License for more details. # +* # # +* # You should have received a copy of the GNU General Public License # +* # along with this program; if not, write to the # +* # Free Software Foundation, Inc., # +* # 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. # +* ################################################################################ */ + + + +#include +#include + +/* See http://jclahr.com/science/software/magnitude/mb/index.html + * + * These are the original Gutenberg and Richter (1956) calibration functions + * as they are still (2007) used at the NEIC + * + * The Q values for a distance of 5 degrees are in the 4th column, i.e. + * at __qmb[x][3] + */ + +/* +static float __qmb[17][108] = { + { 5.6, 5.8, 6.1, 6.4, 6.5, 7.0, 7.0, 7.2, 7.3, 7.2, 7.1, 7.0, 6.6, 6.3, 5.9, 5.9, 5.9, 6.0, 6.1, 6.1, 6.2, 6.3, 6.4, 6.5, 6.5, 6.5, 6.6, 6.6, 6.6, 6.7, 6.7, 6.7, 6.7, 6.6, 6.6, 6.5, 6.5, 6.4, 6.4, 6.5, 6.5, 6.5, 6.6, 6.7, 6.8, 6.9, 6.9, 6.8, 6.7, 6.7, 6.7, 6.7, 6.8, 6.8, 6.8, 6.8, 6.8, 6.9, 6.9, 6.9, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 7.0, 6.9, 6.9, 6.9, 6.9, 6.8, 6.8, 6.9, 6.9, 6.9, 6.8, 6.7, 6.8, 6.9, 7.0, 7.0, 7.0, 6.9, 7.0, 7.1, 7.0, 7.0, 7.1, 7.1, 7.2, 7.1, 7.2, 7.3, 7.4, 7.5, 7.5, 7.3, 7.4, 7.4, 7.5, 7.6, 7.7, 7.8, 7.8, 7.9, 8.0 }, + { 0.0, 0.0, 0.0, 6.3, 6.5, 6.8, 7.0, 7.0, 7.1, 7.0, 7.0, 6.9, 6.5, 6.1, 5.9, 5.9, 5.9, 6.0, 6.1, 6.2, 6.2, 6.3, 6.3, 6.4, 6.4, 6.4, 6.5, 6.5, 6.6, 6.6, 6.7, 6.7, 6.7, 6.7, 6.7, 6.6, 6.6, 6.5, 6.5, 6.5, 6.5, 6.5, 6.6, 6.7, 6.7, 6.8, 6.8, 6.8, 6.8, 6.7, 6.7, 6.7, 6.8, 6.8, 6.8, 6.8, 6.8, 6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 6.8, 6.8, 6.8, 6.8, 6.8, 6.8, 6.8, 6.8, 6.8, 6.8, 6.8, 6.9, 7.0, 7.0, 7.0, 7.0, 7.1, 7.1, 7.0, 7.1, 7.2, 7.2, 7.2, 7.2, 7.2, 7.3, 7.3, 7.3, 7.3, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.8, 7.9, 8.0 }, + { 0.0, 0.0, 0.0, 6.0, 6.2, 6.5, 6.6, 6.8, 6.9, 6.9, 6.8, 6.7, 6.5, 6.0, 5.9, 5.9, 5.9, 6.0, 6.1, 6.1, 6.2, 6.2, 6.3, 6.3, 6.3, 6.4, 6.4, 6.4, 6.5, 6.5, 6.6, 6.6, 6.7, 6.7, 6.7, 6.7, 6.7, 6.6, 6.6, 6.5, 6.5, 6.5, 6.5, 6.6, 6.7, 6.7, 6.8, 6.8, 6.8, 6.8, 6.8, 6.8, 6.8, 6.8, 6.8, 6.8, 6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 6.8, 6.8, 6.8, 6.8, 6.8, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.8, 6.8, 6.7, 6.7, 6.7, 6.8, 6.8, 6.8, 6.9, 7.0, 7.0, 7.0, 7.1, 7.1, 7.2, 7.2, 7.2, 7.2, 7.2, 7.3, 7.3, 7.3, 7.3, 7.3, 7.4, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 7.9, 8.0 }, + { 0.0, 0.0, 0.0, 6.0, 6.1, 6.4, 6.5, 6.6, 6.8, 6.7, 6.7, 6.5, 6.1, 6.0, 5.9, 6.0, 6.0, 6.0, 6.1, 6.1, 6.2, 6.2, 6.2, 6.2, 6.3, 6.3, 6.4, 6.4, 6.5, 6.5, 6.6, 6.6, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.6, 6.6, 6.6, 6.6, 6.6, 6.7, 6.7, 6.7, 6.8, 6.8, 6.8, 6.8, 6.8, 6.8, 6.8, 6.8, 6.9, 6.9, 6.9, 6.9, 6.9, 6.9, 6.8, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.8, 6.8, 6.8, 6.9, 6.9, 7.0, 7.0, 7.1, 7.2, 7.2, 7.2, 7.2, 7.3, 7.3, 7.3, 7.3, 7.4, 7.4, 7.5, 7.5, 7.7, 7.8, 7.8, 7.9, 7.9, 8.0 }, + { 0.0, 0.0, 0.0, 6.0, 6.0, 6.3, 6.4, 6.6, 6.7, 6.6, 6.5, 6.4, 6.0, 6.0, 6.0, 6.0, 6.0, 6.1, 6.1, 6.1, 6.1, 6.1, 6.2, 6.2, 6.3, 6.3, 6.4, 6.4, 6.5, 6.5, 6.5, 6.6, 6.6, 6.7, 6.7, 6.7, 6.7, 6.6, 6.6, 6.6, 6.6, 6.6, 6.6, 6.6, 6.6, 6.6, 6.6, 6.7, 6.8, 6.8, 6.8, 6.8, 6.8, 6.8, 6.8, 6.8, 6.8, 6.9, 6.9, 6.8, 6.8, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.6, 6.6, 6.6, 6.6, 6.6, 6.6, 6.6, 6.6, 6.6, 6.6, 6.6, 6.6, 6.6, 6.6, 6.7, 6.7, 6.7, 6.8, 6.8, 6.8, 6.9, 6.9, 7.0, 7.1, 7.1, 7.2, 7.2, 7.3, 7.3, 7.3, 7.4, 7.4, 7.4, 7.5, 7.6, 7.7, 7.8, 7.8, 7.9, 7.9, 8.0 }, + { 0.0, 0.0, 0.0, 6.0, 6.0, 6.2, 6.4, 6.4, 6.4, 6.4, 6.3, 6.2, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.2, 6.2, 6.2, 6.2, 6.3, 6.3, 6.4, 6.4, 6.5, 6.5, 6.5, 6.5, 6.6, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.5, 6.5, 6.5, 6.5, 6.5, 6.6, 6.6, 6.6, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.6, 6.6, 6.6, 6.6, 6.6, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.6, 6.6, 6.6, 6.6, 6.7, 6.8, 6.8, 6.8, 6.9, 6.9, 7.0, 7.0, 7.1, 7.2, 7.2, 7.3, 7.3, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 7.9, 8.0, 8.0 }, + { 0.0, 0.0, 0.0, 5.9, 6.0, 6.0, 6.0, 6.1, 6.1, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.3, 6.3, 6.3, 6.3, 6.4, 6.4, 6.4, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.3, 6.3, 6.3, 6.3, 6.3, 6.2, 6.2, 6.2, 6.2, 6.2, 6.3, 6.3, 6.4, 6.4, 6.4, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.6, 6.6, 6.6, 6.6, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.6, 6.6, 6.6, 6.7, 6.7, 6.8, 6.8, 6.9, 6.9, 7.0, 7.1, 7.1, 7.2, 7.2, 7.2, 7.3, 7.4, 7.5, 7.7, 7.8, 7.9, 7.9, 8.0, 8.0 }, + { 0.0, 0.0, 0.0, 5.7, 5.7, 5.8, 5.8, 5.8, 5.9, 5.9, 6.0, 6.0, 6.1, 6.2, 6.2, 6.2, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.3, 6.3, 6.3, 6.3, 6.3, 6.2, 6.2, 6.2, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.2, 6.2, 6.2, 6.2, 6.2, 6.3, 6.3, 6.3, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.6, 6.6, 6.6, 6.6, 6.6, 6.7, 6.7, 6.7, 6.8, 6.8, 6.9, 7.0, 7.0, 7.1, 7.1, 7.2, 7.2, 7.3, 7.3, 7.4, 7.6, 7.8, 7.9, 8.0, 8.0, 8.0 }, + { 0.0, 0.0, 0.0, 5.6, 5.7, 5.7, 5.8, 5.8, 5.9, 6.0, 6.0, 6.1, 6.1, 6.1, 6.2, 6.2, 6.2, 6.1, 6.1, 6.1, 6.1, 6.2, 6.2, 6.2, 6.2, 6.2, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.2, 6.2, 6.1, 6.1, 6.0, 6.0, 6.0, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.3, 6.3, 6.3, 6.3, 6.4, 6.4, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.6, 6.6, 6.6, 6.6, 6.6, 6.6, 6.6, 6.6, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.8, 6.9, 6.9, 7.0, 7.0, 7.1, 7.2, 7.2, 7.3, 7.3, 7.4, 7.6, 7.8, 7.9, 8.0, 8.0, 8.0 }, + { 0.0, 0.0, 0.0, 5.7, 5.7, 5.8, 5.9, 6.0, 6.1, 6.2, 6.2, 6.2, 6.2, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.2, 6.2, 6.2, 6.2, 6.2, 6.3, 6.3, 6.3, 6.4, 6.4, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.6, 6.6, 6.6, 6.6, 6.6, 6.6, 6.6, 6.6, 6.6, 6.7, 6.7, 6.7, 6.7, 6.7, 6.7, 6.8, 6.8, 6.9, 6.9, 7.0, 7.0, 7.1, 7.1, 7.2, 7.3, 7.3, 7.4, 7.6, 7.8, 7.9, 8.0, 8.0, 8.0 }, + { 0.0, 0.0, 0.0, 5.7, 5.8, 5.9, 6.0, 6.1, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.1, 6.1, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.2, 6.2, 6.2, 6.2, 6.3, 6.3, 6.4, 6.4, 6.4, 6.5, 6.5, 6.5, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.6, 6.6, 6.6, 6.6, 6.6, 6.6, 6.7, 6.7, 6.8, 6.9, 6.9, 6.9, 7.0, 7.1, 7.1, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.9, 8.0, 8.0 }, + { 0.0, 0.0, 0.0, 5.8, 5.9, 5.9, 6.1, 6.1, 6.2, 6.2, 6.2, 6.2, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.2, 6.2, 6.2, 6.2, 6.3, 6.3, 6.3, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.5, 6.5, 6.5, 6.5, 6.5, 6.6, 6.6, 6.7, 6.7, 6.8, 6.8, 6.9, 6.9, 7.0, 7.0, 7.1, 7.2, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.8, 7.8, 7.8, 7.9 }, + { 0.0, 0.0, 0.0, 5.8, 5.8, 5.9, 6.0, 6.1, 6.2, 6.2, 6.2, 6.3, 6.3, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.2, 6.2, 6.2, 6.2, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.2, 6.2, 6.2, 6.3, 6.3, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.2, 6.2, 6.2, 6.2, 6.3, 6.3, 6.3, 6.4, 6.4, 6.5, 6.5, 6.6, 6.6, 6.7, 6.7, 6.8, 6.8, 6.9, 6.9, 7.0, 7.1, 7.1, 7.2, 7.2, 7.3, 7.4, 7.5, 7.6, 7.7, 7.7, 7.7, 7.8, 7.8 }, + { 0.0, 0.0, 0.0, 5.8, 5.8, 5.9, 5.9, 6.0, 6.1, 6.1, 6.2, 6.2, 6.2, 6.3, 6.3, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.2, 6.2, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.1, 6.2, 6.2, 6.2, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.3, 6.3, 6.3, 6.4, 6.4, 6.5, 6.5, 6.6, 6.7, 6.7, 6.8, 6.8, 6.9, 6.9, 7.0, 7.0, 7.1, 7.1, 7.2, 7.2, 7.3, 7.4, 7.5, 7.6, 7.6, 7.7, 7.7, 7.7, 7.7 }, + { 0.0, 0.0, 0.0, 5.7, 5.8, 5.8, 5.8, 5.9, 5.8, 5.9, 6.0, 6.0, 6.1, 6.1, 6.2, 6.2, 6.3, 6.3, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.2, 6.2, 6.2, 6.1, 6.1, 6.1, 6.0, 6.0, 6.0, 6.0, 6.1, 6.1, 6.2, 6.2, 6.2, 6.2, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.2, 6.2, 6.2, 6.2, 6.3, 6.3, 6.3, 6.4, 6.4, 6.5, 6.5, 6.6, 6.7, 6.7, 6.8, 6.8, 6.8, 7.0, 7.0, 7.0, 7.1, 7.1, 7.2, 7.2, 7.3, 7.3, 7.4, 7.5, 7.6, 7.6, 7.6, 7.7, 7.7 }, + { 0.0, 0.0, 0.0, 5.7, 5.7, 5.7, 5.8, 5.8, 5.8, 5.8, 5.8, 5.9, 5.9, 5.9, 6.0, 6.1, 6.1, 6.1, 6.2, 6.2, 6.3, 6.3, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.3, 6.3, 6.3, 6.3, 6.3, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.1, 6.1, 6.1, 6.1, 6.1, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.1, 6.1, 6.1, 6.1, 6.2, 6.2, 6.2, 6.2, 6.2, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.4, 6.4, 6.4, 6.5, 6.5, 6.5, 6.6, 6.7, 6.7, 6.8, 6.9, 6.9, 7.0, 7.0, 7.0, 7.1, 7.1, 7.2, 7.2, 7.3, 7.3, 7.4, 7.5, 7.5, 7.5, 7.6, 7.6 }, + { 0.0, 0.0, 0.0, 5.7, 5.7, 5.7, 5.7, 5.7, 5.7, 5.8, 5.8, 5.8, 5.8, 5.8, 5.9, 5.9, 6.0, 6.0, 6.0, 6.0, 6.1, 6.1, 6.1, 6.2, 6.2, 6.3, 6.3, 6.3, 6.3, 6.3, 6.4, 6.4, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.2, 6.2, 6.2, 6.2, 6.2, 6.2, 6.1, 6.1, 6.1, 6.1, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.0, 6.1, 6.1, 6.1, 6.1, 6.2, 6.2, 6.2, 6.2, 6.2, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.3, 6.4, 6.4, 6.4, 6.4, 6.5, 6.5, 6.6, 6.7, 6.7, 6.8, 6.9, 6.9, 7.0, 7.0, 7.0, 7.1, 7.1, 7.2, 7.2, 7.3, 7.4, 7.4, 7.4, 7.5, 7.5 } +}; + +/* Q-factor of P waves according to Veith-Clawson (1972) */ +static float depths[11] = {0., 15., 40., 100., 200., 300., 400., 500., 600., 700., 800.}; +static float qmb[11][181] = { + { // 0 km +0.301, 1.191, 2.371, 2.501, 2.851, 3.061, 3.201, 3.321, 3.401, 3.451, 3.491, 3.531, 3.551, 3.561, 3.561, 3.551, 3.511, 3.401, 3.281, 3.091, 3.071, 3.101, 3.151, 3.241, 3.341, 3.451, 3.551, 3.651, 3.721, 3.741, 3.721, 3.681, 3.661, 3.661, 3.651, 3.641, 3.641, 3.641, 3.631, 3.631, 3.621, 3.621, 3.621, 3.631, 3.631, 3.641, 3.641, 3.651, 3.661, 3.661, 3.671, 3.671, 3.681, 3.691, 3.691, 3.701, 3.701, 3.711, 3.721, 3.721, 3.731, 3.741, 3.741, 3.751, 3.751, 3.761, 3.761, 3.771, 3.781, 3.781, 3.791, 3.801, 3.801, 3.811, 3.811, 3.821, 3.831, 3.831, 3.841, 3.841, 3.851, 3.861, 3.871, 3.881, 3.891, 3.911, 3.941, 3.961, 3.981, 4.021, 4.061, 4.101, 4.151, 4.211, 4.281, 4.361, 4.441, 4.521, 4.601, 4.681, 4.761, 4.75, 4.74, 4.72, 4.71, 4.70, 4.70, 4.70, 4.70, 4.70, 4.70, 4.50, 4.50, 4.50, 4.30, 4.30, 4.30, 4.30, 4.30, 4.30, 4.25, 4.25, 4.25, 4.25, 4.20, 4.20, 4.20, 4.20, 4.20, 4.20, 4.20, 4.20, 4.20, 4.20, 4.20, 4.20, 4.20, 4.20, 4.20, 4.20, 4.20, 4.20, 3.50, 3.50, 3.50, 3.50, 3.50, 3.50, 3.50, 3.50, 3.50, 3.50, 3.60, 3.65, 3.80, 3.87, 3.93, 4.00, 4.02, 4.05, 4.07, 4.10, 4.13, 4.16, 4.19, 4.22, 4.25, 4.28, 4.31, 4.34, 4.37, 4.41, 4.45, 4.49, 4.52, 4.55, 4.58, 4.61, 4.64, 4.67, 4.70 + }, + { // 15 km +-1.519, 3.221, 2.321, 2.531, 2.831, 3.031, 3.161, 3.271, 3.341, 3.391, 3.431, 3.461, 3.481, 3.491, 3.491, 3.471, 3.411, 3.301, 3.151, 3.011, 3.011, 3.051, 3.121, 3.191, 3.291, 3.391, 3.491, 3.571, 3.631, 3.651, 3.631, 3.611, 3.591, 3.581, 3.571, 3.561, 3.551, 3.551, 3.541, 3.541, 3.541, 3.541, 3.541, 3.541, 3.541, 3.551, 3.551, 3.561, 3.561, 3.571, 3.581, 3.581, 3.591, 3.601, 3.601, 3.611, 3.621, 3.621, 3.631, 3.641, 3.641, 3.651, 3.651, 3.661, 3.661, 3.671, 3.681, 3.681, 3.691, 3.691, 3.701, 3.701, 3.711, 3.721, 3.721, 3.731, 3.731, 3.741, 3.741, 3.751, 3.751, 3.761, 3.761, 3.771, 3.791, 3.811, 3.831, 3.851, 3.881, 3.921, 3.961, 4.001, 4.051, 4.111, 4.191, 4.271, 4.351, 4.431, 4.511, 4.591, 4.671, 4.66, 4.65, 4.63, 4.62, 4.61, 4.61, 4.61, 4.61, 4.61, 4.61, 4.41, 4.41, 4.41, 4.21, 4.21, 4.21, 4.21, 4.21, 4.21, 4.16, 4.16, 4.16, 4.16, 4.11, 4.11, 4.11, 4.11, 4.11, 4.11, 4.11, 4.11, 4.11, 4.11, 4.11, 4.11, 4.11, 4.11, 4.11, 4.11, 4.11, 4.11, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.51, 3.56, 3.71, 3.78, 3.84, 3.91, 3.93, 3.96, 3.98, 4.01, 4.04, 4.07, 4.10, 4.13, 4.16, 4.19, 4.22, 4.25, 4.28, 4.32, 4.36, 4.40, 4.43, 4.46, 4.49, 4.52, 4.55, 4.58, 4.61 + }, + { // 40 km +-1.139, 2.461, 2.591, 2.831, 3.021, 3.171, 3.291, 3.381, 3.431, 3.471, 3.491, 3.501, 3.511, 3.511, 3.511, 3.471, 3.381, 3.241, 3.051, 2.931, 2.921, 2.951, 3.011, 3.091, 3.171, 3.271, 3.361, 3.431, 3.481, 3.491, 3.501, 3.501, 3.491, 3.481, 3.461, 3.451, 3.441, 3.431, 3.421, 3.421, 3.411, 3.411, 3.411, 3.411, 3.411, 3.411, 3.421, 3.421, 3.431, 3.441, 3.441, 3.451, 3.461, 3.461, 3.471, 3.481, 3.481, 3.491, 3.491, 3.501, 3.501, 3.511, 3.511, 3.521, 3.521, 3.531, 3.541, 3.541, 3.551, 3.551, 3.561, 3.561, 3.571, 3.571, 3.581, 3.581, 3.591, 3.591, 3.601, 3.601, 3.611, 3.611, 3.621, 3.631, 3.641, 3.661, 3.681, 3.701, 3.741, 3.781, 3.821, 3.861, 3.911, 3.971, 4.051, 4.131, 4.211, 4.291, 4.371, 4.451, 4.531, 4.52, 4.51, 4.49, 4.48, 4.47, 4.47, 4.47, 4.47, 4.47, 4.47, 4.27, 4.27, 4.27, 4.07, 4.07, 4.07, 4.07, 4.07, 4.07, 4.02, 4.02, 4.02, 4.02, 3.97, 3.97, 3.97, 3.97, 3.97, 3.97, 3.97, 3.97, 3.97, 3.97, 3.97, 3.97, 3.97, 3.97, 3.97, 3.97, 3.97, 3.97, 3.27, 3.27, 3.27, 3.27, 3.27, 3.27, 3.27, 3.27, 3.27, 3.27, 3.37, 3.42, 3.57, 3.64, 3.70, 3.77, 3.79, 3.82, 3.84, 3.87, 3.90, 3.93, 3.96, 3.99, 4.02, 4.05, 4.08, 4.11, 4.14, 4.18, 4.22, 4.26, 4.29, 4.32, 4.35, 4.38, 4.41, 4.44, 4.47 + }, + { // 100 km +-0.499, -0.009, 0.691, 1.241, 1.661, 1.981, 2.251, 2.481, 2.671, 2.821, 2.931, 3.021, 3.091, 3.141, 3.151, 3.091, 2.981, 2.841, 2.761, 2.751, 2.781, 2.841, 2.911, 2.991, 3.081, 3.181, 3.271, 3.341, 3.371, 3.381, 3.381, 3.371, 3.361, 3.351, 3.341, 3.331, 3.321, 3.311, 3.301, 3.301, 3.291, 3.291, 3.291, 3.291, 3.291, 3.301, 3.301, 3.311, 3.321, 3.321, 3.331, 3.341, 3.341, 3.351, 3.361, 3.361, 3.371, 3.381, 3.381, 3.391, 3.401, 3.401, 3.411, 3.421, 3.421, 3.431, 3.441, 3.441, 3.451, 3.461, 3.461, 3.461, 3.471, 3.471, 3.481, 3.481, 3.491, 3.491, 3.501, 3.501, 3.501, 3.501, 3.511, 3.531, 3.551, 3.571, 3.601, 3.621, 3.641, 3.681, 3.721, 3.761, 3.811, 3.871, 3.951, 4.031, 4.111, 4.191, 4.271, 4.351, 4.431, 4.42, 4.41, 4.39, 4.38, 4.37, 4.37, 4.37, 4.37, 4.37, 4.37, 4.17, 4.17, 4.17, 3.97, 3.97, 3.97, 3.97, 3.97, 3.97, 3.92, 3.92, 3.92, 3.92, 3.87, 3.87, 3.87, 3.87, 3.87, 3.87, 3.87, 3.87, 3.87, 3.87, 3.87, 3.87, 3.87, 3.87, 3.87, 3.87, 3.87, 3.87, 3.17, 3.17, 3.17, 3.17, 3.17, 3.17, 3.17, 3.17, 3.17, 3.17, 3.27, 3.32, 3.47, 3.54, 3.60, 3.67, 3.69, 3.72, 3.74, 3.77, 3.80, 3.83, 3.86, 3.89, 3.92, 3.95, 3.98, 4.01, 4.04, 4.08, 4.12, 4.16, 4.19, 4.22, 4.25, 4.28, 4.31, 4.34, 4.37 + }, + { // 200 km +-0.059, 0.111, 0.471, 0.851, 1.191, 1.491, 1.761, 1.991, 2.181, 2.331, 2.441, 2.511, 2.551, 2.561, 2.531, 2.481, 2.451, 2.461, 2.491, 2.541, 2.601, 2.681, 2.761, 2.861, 2.961, 3.041, 3.101, 3.131, 3.141, 3.151, 3.151, 3.141, 3.131, 3.121, 3.111, 3.101, 3.101, 3.091, 3.091, 3.091, 3.081, 3.081, 3.091, 3.091, 3.101, 3.111, 3.111, 3.121, 3.131, 3.131, 3.141, 3.151, 3.161, 3.161, 3.171, 3.181, 3.181, 3.191, 3.201, 3.211, 3.211, 3.221, 3.231, 3.231, 3.241, 3.251, 3.251, 3.261, 3.271, 3.281, 3.281, 3.291, 3.301, 3.301, 3.311, 3.311, 3.321, 3.321, 3.321, 3.331, 3.331, 3.331, 3.341, 3.351, 3.361, 3.381, 3.401, 3.441, 3.481, 3.521, 3.561, 3.611, 3.671, 3.741, 3.821, 3.901, 3.981, 4.061, 4.141, 4.221, 4.301, 4.29, 4.28, 4.26, 4.25, 4.24, 4.24, 4.24, 4.24, 4.24, 4.24, 4.04, 4.04, 4.04, 3.84, 3.84, 3.84, 3.84, 3.84, 3.84, 3.79, 3.79, 3.79, 3.79, 3.74, 3.74, 3.74, 3.74, 3.74, 3.74, 3.74, 3.74, 3.74, 3.74, 3.74, 3.74, 3.74, 3.74, 3.74, 3.74, 3.74, 3.74, 3.04, 3.04, 3.04, 3.04, 3.04, 3.04, 3.04, 3.04, 3.04, 3.04, 3.14, 3.19, 3.34, 3.41, 3.47, 3.54, 3.56, 3.59, 3.61, 3.64, 3.67, 3.70, 3.73, 3.76, 3.79, 3.82, 3.85, 3.88, 3.91, 3.95, 3.99, 4.03, 4.06, 4.09, 4.12, 4.15, 4.18, 4.21, 4.24 + }, + { // 300 km +0.201, 0.291, 0.491, 0.741, 0.991, 1.231, 1.441, 1.621, 1.771, 1.891, 1.991, 2.061, 2.091, 2.121, 2.161, 2.201, 2.251, 2.311, 2.371, 2.441, 2.521, 2.601, 2.691, 2.781, 2.861, 2.921, 2.961, 2.981, 2.991, 2.991, 2.981, 2.971, 2.961, 2.961, 2.951, 2.941, 2.941, 2.941, 2.931, 2.931, 2.931, 2.941, 2.941, 2.951, 2.951, 2.961, 2.971, 2.971, 2.981, 2.991, 3.001, 3.001, 3.011, 3.021, 3.031, 3.031, 3.041, 3.051, 3.061, 3.061, 3.071, 3.081, 3.091, 3.101, 3.101, 3.111, 3.121, 3.131, 3.141, 3.141, 3.151, 3.161, 3.161, 3.171, 3.181, 3.181, 3.191, 3.191, 3.201, 3.201, 3.201, 3.211, 3.221, 3.241, 3.261, 3.291, 3.311, 3.331, 3.371, 3.411, 3.451, 3.501, 3.561, 3.641, 3.721, 3.801, 3.881, 3.961, 4.041, 4.121, 4.201, 4.19, 4.18, 4.16, 4.15, 4.14, 4.14, 4.14, 4.14, 4.14, 4.14, 3.94, 3.94, 3.94, 3.74, 3.74, 3.74, 3.74, 3.74, 3.74, 3.69, 3.69, 3.69, 3.69, 3.64, 3.64, 3.64, 3.64, 3.64, 3.64, 3.64, 3.64, 3.64, 3.64, 3.64, 3.64, 3.64, 3.64, 3.64, 3.64, 3.64, 3.64, 2.94, 2.94, 2.94, 2.94, 2.94, 2.94, 2.94, 2.94, 2.94, 2.94, 3.04, 3.09, 3.24, 3.31, 3.37, 3.44, 3.46, 3.49, 3.51, 3.54, 3.57, 3.60, 3.63, 3.66, 3.69, 3.72, 3.75, 3.78, 3.81, 3.85, 3.89, 3.93, 3.96, 3.99, 4.02, 4.05, 4.08, 4.10, 4.14 + }, + { // 400 km +0.381, 0.421, 0.551, 0.711, 0.901, 1.081, 1.251, 1.401, 1.531, 1.651, 1.751, 1.841, 1.921, 1.981, 2.041, 2.101, 2.171, 2.251, 2.321, 2.401, 2.481, 2.571, 2.651, 2.731, 2.801, 2.841, 2.861, 2.871, 2.861, 2.841, 2.831, 2.821, 2.811, 2.811, 2.811, 2.811, 2.811, 2.811, 2.811, 2.811, 2.821, 2.821, 2.831, 2.841, 2.841, 2.851, 2.861, 2.871, 2.871, 2.881, 2.891, 2.901, 2.911, 2.921, 2.931, 2.941, 2.941, 2.951, 2.961, 2.971, 2.981, 2.981, 2.991, 3.001, 3.011, 3.021, 3.031, 3.041, 3.041, 3.051, 3.061, 3.071, 3.071, 3.081, 3.091, 3.091, 3.101, 3.101, 3.111, 3.111, 3.111, 3.131, 3.151, 3.171, 3.191, 3.211, 3.231, 3.251, 3.291, 3.331, 3.381, 3.431, 3.501, 3.581, 3.661, 3.741, 3.821, 3.901, 3.981, 4.061, 4.141, 4.13, 4.12, 4.10, 4.09, 4.08, 4.08, 4.08, 4.08, 4.08, 4.08, 3.88, 3.88, 3.88, 3.68, 3.68, 3.68, 3.68, 3.68, 3.68, 3.63, 3.63, 3.63, 3.63, 3.58, 3.58, 3.58, 3.58, 3.58, 3.58, 3.58, 3.58, 3.58, 3.58, 3.58, 3.58, 3.58, 3.58, 3.58, 3.58, 3.58, 3.58, 2.88, 2.88, 2.88, 2.88, 2.88, 2.88, 2.88, 2.88, 2.88, 2.88, 2.98, 3.03, 3.18, 3.25, 3.31, 3.38, 3.40, 3.43, 3.45, 3.48, 3.51, 3.54, 3.57, 3.60, 3.63, 3.66, 3.69, 3.72, 3.75, 3.79, 3.83, 3.87, 3.90, 3.93, 3.96, 3.99, 4.02, 4.05, 4.08 + }, + { // 500 km +0.521, 0.551, 0.631, 0.751, 0.881, 1.031, 1.171, 1.301, 1.421, 1.541, 1.641, 1.731, 1.821, 1.901, 1.981, 2.061, 2.141, 2.221, 2.301, 2.371, 2.451, 2.521, 2.581, 2.631, 2.671, 2.691, 2.691, 2.681, 2.681, 2.681, 2.671, 2.671, 2.671, 2.671, 2.671, 2.671, 2.671, 2.681, 2.681, 2.691, 2.701, 2.701, 2.711, 2.721, 2.731, 2.731, 2.741, 2.751, 2.761, 2.771, 2.781, 2.791, 2.801, 2.811, 2.821, 2.831, 2.841, 2.841, 2.851, 2.861, 2.871, 2.881, 2.891, 2.901, 2.911, 2.921, 2.921, 2.931, 2.941, 2.951, 2.961, 2.971, 2.981, 2.981, 2.991, 3.001, 3.001, 3.001, 3.011, 3.011, 3.021, 3.031, 3.041, 3.061, 3.081, 3.101, 3.131, 3.171, 3.211, 3.251, 3.301, 3.361, 3.431, 3.511, 3.591, 3.671, 3.751, 3.831, 3.911, 3.991, 4.071, 4.06, 4.05, 4.03, 4.02, 4.01, 4.01, 4.01, 4.01, 4.01, 4.01, 3.81, 3.81, 3.81, 3.61, 3.61, 3.61, 3.61, 3.61, 3.61, 3.56, 3.56, 3.56, 3.56, 3.51, 3.51, 3.51, 3.51, 3.51, 3.51, 3.51, 3.51, 3.51, 3.51, 3.51, 3.51, 3.51, 3.51, 3.51, 3.51, 3.51, 3.51, 2.81, 2.81, 2.81, 2.81, 2.81, 2.81, 2.81, 2.81, 2.81, 2.81, 2.91, 2.96, 3.11, 3.18, 3.24, 3.31, 3.33, 3.36, 3.38, 3.41, 3.44, 3.47, 3.50, 3.53, 3.56, 3.59, 3.62, 3.65, 3.68, 3.72, 3.76, 3.80, 3.83, 3.86, 3.89, 3.92, 3.95, 3.98, 4.01 + }, + { // 600 km +0.641, 0.661, 0.721, 0.801, 0.911, 1.031, 1.141, 1.251, 1.361, 1.471, 1.571, 1.671, 1.771, 1.861, 1.941, 2.021, 2.101, 2.171, 2.241, 2.311, 2.371, 2.431, 2.471, 2.491, 2.501, 2.501, 2.491, 2.491, 2.481, 2.481, 2.481, 2.491, 2.501, 2.511, 2.521, 2.531, 2.541, 2.541, 2.551, 2.561, 2.571, 2.581, 2.591, 2.601, 2.611, 2.621, 2.631, 2.641, 2.651, 2.661, 2.671, 2.681, 2.691, 2.701, 2.711, 2.721, 2.731, 2.741, 2.751, 2.761, 2.771, 2.781, 2.791, 2.801, 2.811, 2.821, 2.831, 2.841, 2.851, 2.861, 2.861, 2.871, 2.881, 2.891, 2.891, 2.901, 2.911, 2.911, 2.921, 2.921, 2.931, 2.941, 2.961, 2.981, 3.001, 3.031, 3.061, 3.101, 3.141, 3.181, 3.241, 3.301, 3.381, 3.461, 3.541, 3.621, 3.701, 3.781, 3.861, 3.941, 4.021, 4.01, 4.00, 3.98, 3.97, 3.96, 3.96, 3.96, 3.96, 3.96, 3.96, 3.76, 3.76, 3.76, 3.56, 3.56, 3.56, 3.56, 3.56, 3.56, 3.51, 3.51, 3.51, 3.51, 3.46, 3.46, 3.46, 3.46, 3.46, 3.46, 3.46, 3.46, 3.46, 3.46, 3.46, 3.46, 3.46, 3.46, 3.46, 3.46, 3.46, 3.46, 2.76, 2.76, 2.76, 2.76, 2.76, 2.76, 2.76, 2.76, 2.76, 2.76, 2.86, 2.91, 3.06, 3.13, 3.19, 3.26, 3.28, 3.31, 3.33, 3.36, 3.39, 3.42, 3.45, 3.48, 3.51, 3.54, 3.57, 3.60, 3.63, 3.67, 3.71, 3.75, 3.78, 3.81, 3.84, 3.87, 3.90, 3.93, 3.96 + }, + { // 700 km +0.741, 0.761, 0.801, 0.871, 0.951, 1.051, 1.141, 1.241, 1.341, 1.441, 1.531, 1.621, 1.711, 1.791, 1.871, 1.941, 2.011, 2.081, 2.141, 2.181, 2.211, 2.231, 2.241, 2.241, 2.251, 2.261, 2.271, 2.281, 2.291, 2.301, 2.311, 2.331, 2.341, 2.351, 2.361, 2.371, 2.391, 2.401, 2.411, 2.421, 2.431, 2.451, 2.461, 2.471, 2.481, 2.491, 2.511, 2.521, 2.531, 2.541, 2.561, 2.571, 2.581, 2.591, 2.611, 2.621, 2.631, 2.641, 2.651, 2.671, 2.681, 2.691, 2.701, 2.711, 2.721, 2.731, 2.741, 2.751, 2.761, 2.771, 2.781, 2.791, 2.791, 2.801, 2.801, 2.811, 2.821, 2.821, 2.821, 2.831, 2.841, 2.861, 2.881, 2.901, 2.931, 2.961, 2.991, 3.031, 3.071, 3.121, 3.181, 3.251, 3.331, 3.411, 3.491, 3.571, 3.651, 3.731, 3.811, 3.891, 3.971, 3.96, 3.95, 3.93, 3.92, 3.91, 3.91, 3.91, 3.91, 3.91, 3.91, 3.71, 3.71, 3.71, 3.51, 3.51, 3.51, 3.51, 3.51, 3.51, 3.46, 3.46, 3.46, 3.46, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 2.71, 2.71, 2.71, 2.71, 2.71, 2.71, 2.71, 2.71, 2.71, 2.71, 2.81, 2.86, 3.01, 3.08, 3.14, 3.21, 3.23, 3.26, 3.28, 3.31, 3.34, 3.37, 3.40, 3.43, 3.46, 3.49, 3.52, 3.55, 3.58, 3.62, 3.66, 3.70, 3.73, 3.76, 3.79, 3.82, 3.85, 3.88, 3.91 + }, + { // 800 km +0.821, 0.831, 0.861, 0.921, 0.981, 1.061, 1.141, 1.221, 1.301, 1.391, 1.471, 1.551, 1.621, 1.691, 1.751, 1.801, 1.851, 1.901, 1.941, 1.981, 2.011, 2.041, 2.061, 2.081, 2.101, 2.121, 2.141, 2.161, 2.181, 2.201, 2.221, 2.241, 2.261, 2.271, 2.291, 2.301, 2.321, 2.341, 2.351, 2.371, 2.381, 2.391, 2.411, 2.421, 2.431, 2.451, 2.461, 2.471, 2.481, 2.491, 2.511, 2.521, 2.531, 2.541, 2.551, 2.561, 2.571, 2.591, 2.601, 2.611, 2.621, 2.631, 2.641, 2.651, 2.661, 2.671, 2.681, 2.691, 2.701, 2.711, 2.721, 2.731, 2.741, 2.751, 2.761, 2.771, 2.771, 2.781, 2.781, 2.791, 2.801, 2.821, 2.851, 2.881, 2.911, 2.951, 2.991, 3.031, 3.081, 3.141, 3.211, 3.271, 3.331, 3.411, 3.491, 3.571, 3.651, 3.731, 3.811, 3.891, 3.971, 3.96, 3.95, 3.93, 3.92, 3.91, 3.91, 3.91, 3.91, 3.91, 3.91, 3.71, 3.71, 3.71, 3.51, 3.51, 3.51, 3.51, 3.51, 3.51, 3.46, 3.46, 3.46, 3.46, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 3.41, 2.71, 2.71, 2.71, 2.71, 2.71, 2.71, 2.71, 2.71, 2.71, 2.71, 2.81, 2.86, 3.01, 3.08, 3.14, 3.21, 3.23, 3.26, 3.28, 3.31, 3.34, 3.37, 3.40, 3.43, 3.46, 3.49, 3.52, 3.55, 3.58, 3.62, 3.66, 3.70, 3.73, 3.76, 3.79, 3.82, 3.85, 3.88, 3.91 + } +}; + + + + +#define max(x,y) ((x)>(y)?(x):(y)) +#define min(x,y) ((x)<(y)?(x):(y)) + +static double +bmagnz(double amplitude, double period, double delta, double depth, int *err) +{ + int j, k; + double Q, Q1, Q2, mb; + + *err = 0; + + if (depth < 0. || depth > 800. || delta < 0.1 || delta > 180) { + *err = -1; + return 0.; + } + + if (depth < 15) + k = 0; + else if (depth < 40) + k = 1; + else if (depth < 100) + k = 2; + else if (depth < 200) + k = 3; + else if (depth < 300) + k = 4; + else if (depth < 400) + k = 5; + else if (depth < 500) + k = 6; + else if (depth < 600) + k = 7; + else if (depth < 700) + k = 8; + else if (depth < 800) + k = 9; + else + k = 10; + + + j = min((int)(delta), 181); + // Get the attenuation factor. If the (depth, delta) value is not in the table + // we perform a linear interpolation + if (depths[k] == depth) { + if (j == delta) + Q = qmb[k][j]; + else { + Q = qmb[k][j+1] - qmb[k][j]; + Q *= (delta - j); + Q += qmb[k][j]; + } + } + else { + Q1 = qmb[k][j+1] - qmb[k][j]; + Q1 *= (delta - j); + Q1 += qmb[k][j]; + + Q2 = qmb[k+1][j+1] - qmb[k+1][j]; + Q2 *= (delta - j); + Q2 += qmb[k+1][j]; + + Q = Q2 - Q1; + Q /= (depths[k+1] - depths[k]); + Q *= (depth - depths[k]); + Q += Q1; + } + + amplitude /= 2; //p2p -> z2p + mb = log10(amplitude/period) + Q; + + return mb; +} + +namespace Seiscomp { +namespace Magnitudes { + +bool +compute_mb_VC( + double amplitude, // in micrometers + double period, // in seconds + double delta, // in degrees + double depth, // in kilometers + double *mag) +{ + int err; + double m = bmagnz(amplitude, period, delta, depth, &err); + if (err) return false; + + *mag = m; + + return true; +} + +bool +compute_mb_VC_fromVelocity( + double amplitude, // in micrometers/second + double delta, // in degrees + double depth, // in kilometers + double *mag) +{ + int err; + double m = bmagnz(amplitude, 2*M_PI, delta, depth, &err); + if (err) return false; + + *mag = m; + + return true; +} + +} // namespace Magnitudes +} // namespace Seiscomp