mirror of
https://github.com/yoshiya-usui/femtic.git
synced 2025-05-12 19:51:12 +08:00
1453 lines
76 KiB
C++
1453 lines
76 KiB
C++
//-------------------------------------------------------------------------------------------------------
|
|
// The MIT License (MIT)
|
|
//
|
|
// Copyright (c) 2021 Yoshiya Usui
|
|
//
|
|
// Permission is hereby granted, free of charge, to any person obtaining a copy
|
|
// of this software and associated documentation files (the "Software"), to deal
|
|
// in the Software without restriction, including without limitation the rights
|
|
// to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
|
|
// copies of the Software, and to permit persons to whom the Software is
|
|
// furnished to do so, subject to the following conditions:
|
|
//
|
|
// The above copyright notice and this permission notice shall be included in all
|
|
// copies or substantial portions of the Software.
|
|
//
|
|
// THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
|
|
// IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
|
// FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
|
|
// AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
|
// LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
|
|
// OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
|
|
// SOFTWARE.
|
|
//-------------------------------------------------------------------------------------------------------
|
|
#include <fstream>
|
|
#include <iostream>
|
|
#include <string>
|
|
#include <algorithm>
|
|
#include <assert.h>
|
|
#include <iomanip>
|
|
|
|
#include "ObservedDataStationMT.h"
|
|
#include "ObservedDataStationPoint.h"
|
|
#include "OutputFiles.h"
|
|
#include "CommonParameters.h"
|
|
#include "ResistivityBlock.h"
|
|
|
|
// Constructer
|
|
ObservedDataStationMT::ObservedDataStationMT():
|
|
ObservedDataStationPoint(),
|
|
m_ZxxObserved(NULL),
|
|
m_ZxyObserved(NULL),
|
|
m_ZyxObserved(NULL),
|
|
m_ZyyObserved(NULL),
|
|
m_ZxxSD(NULL),
|
|
m_ZxySD(NULL),
|
|
m_ZyxSD(NULL),
|
|
m_ZyySD(NULL),
|
|
m_ZxxCalculated(NULL),
|
|
m_ZxyCalculated(NULL),
|
|
m_ZyxCalculated(NULL),
|
|
m_ZyyCalculated(NULL),
|
|
m_ZxxResidual(NULL),
|
|
m_ZxyResidual(NULL),
|
|
m_ZyxResidual(NULL),
|
|
m_ZyyResidual(NULL),
|
|
m_rhsVectorIDOfEx(0),
|
|
m_rhsVectorIDOfEy(0),
|
|
m_dataIDOfZxx(NULL),
|
|
m_dataIDOfZxy(NULL),
|
|
m_dataIDOfZyx(NULL),
|
|
m_dataIDOfZyy(NULL),
|
|
m_fixDistortionMatrix(true),
|
|
m_typeOfElectricField(AnalysisControl::USE_HORIZONTAL_ELECTRIC_FIELD),
|
|
m_arrayDistortionMatrixDifferences(NULL),
|
|
m_arrayGainsAndRotations(NULL)
|
|
{
|
|
|
|
for( int i = 0; i < 2; ++i ){
|
|
m_ExCalculated[i] = std::complex<double>(0.0,0.0);
|
|
m_EyCalculated[i] = std::complex<double>(0.0,0.0);
|
|
}
|
|
|
|
if( ( AnalysisControl::getInstance() )->getTypeOfDistortion() == AnalysisControl::DISTORTION_TYPE_UNDEFINED ){
|
|
OutputFiles::m_logFile << "Error : Type of distortion must be defined before instantiation of ObservedDataStationMT !!" << std::endl;
|
|
exit(1);
|
|
}
|
|
else if( ( AnalysisControl::getInstance() )->getTypeOfDistortion() == AnalysisControl::ESTIMATE_DISTORTION_MATRIX_DIFFERENCE ){
|
|
m_arrayDistortionMatrixDifferences = new DistortionMatrixDifferences;
|
|
for( int i = 0; i < 4; ++i ){
|
|
m_arrayDistortionMatrixDifferences->distortionMatrixDifferencePre[i] = 0.0;
|
|
m_arrayDistortionMatrixDifferences->distortionMatrixDifference[i] = 0.0;
|
|
m_arrayDistortionMatrixDifferences->distortionMatrixDifferenceUpdatedFull[i] = 0.0;
|
|
m_arrayDistortionMatrixDifferences->IDsOfDistortionMatrixDifference[i] = -1;
|
|
}
|
|
}
|
|
else if( ( AnalysisControl::getInstance() )->getTypeOfDistortion() == AnalysisControl::ESTIMATE_GAINS_AND_ROTATIONS ||
|
|
( AnalysisControl::getInstance() )->getTypeOfDistortion() == AnalysisControl::ESTIMATE_GAINS_ONLY ){
|
|
m_arrayGainsAndRotations = new GainsAndRotations;
|
|
for( int i = 0; i < 4; ++i ){
|
|
m_arrayGainsAndRotations->gainsAndRotationsPre[i] = 0.0;
|
|
m_arrayGainsAndRotations->gainsAndRotations[i] = 0.0;
|
|
m_arrayGainsAndRotations->gainsAndRotationsUpdatedFull[i] = 0.0;
|
|
m_arrayGainsAndRotations->IDsOfGainsAndRotations[i] = -1;
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
// Destructer
|
|
ObservedDataStationMT::~ObservedDataStationMT(){
|
|
|
|
if( m_ZxxObserved != NULL){
|
|
delete[] m_ZxxObserved;
|
|
m_ZxxObserved = NULL;
|
|
}
|
|
|
|
if( m_ZxyObserved != NULL){
|
|
delete[] m_ZxyObserved;
|
|
m_ZxyObserved = NULL;
|
|
}
|
|
|
|
if( m_ZyxObserved != NULL){
|
|
delete[] m_ZyxObserved;
|
|
m_ZyxObserved = NULL;
|
|
}
|
|
|
|
if( m_ZyyObserved != NULL){
|
|
delete[] m_ZyyObserved;
|
|
m_ZyyObserved = NULL;
|
|
}
|
|
|
|
if( m_ZxxSD != NULL){
|
|
delete[] m_ZxxSD;
|
|
m_ZxxSD = NULL;
|
|
}
|
|
|
|
if( m_ZxySD != NULL){
|
|
delete[] m_ZxySD;
|
|
m_ZxySD = NULL;
|
|
}
|
|
|
|
if( m_ZyxSD != NULL){
|
|
delete[] m_ZyxSD;
|
|
m_ZyxSD = NULL;
|
|
}
|
|
|
|
if( m_ZyySD != NULL){
|
|
delete[] m_ZyySD;
|
|
m_ZyySD = NULL;
|
|
}
|
|
|
|
if( m_ZxxCalculated != NULL){
|
|
delete[] m_ZxxCalculated;
|
|
m_ZxxCalculated = NULL;
|
|
}
|
|
|
|
if( m_ZxyCalculated != NULL){
|
|
delete[] m_ZxyCalculated;
|
|
m_ZxyCalculated = NULL;
|
|
}
|
|
|
|
if( m_ZyxCalculated != NULL){
|
|
delete[] m_ZyxCalculated;
|
|
m_ZyxCalculated = NULL;
|
|
}
|
|
|
|
if( m_ZyyCalculated != NULL){
|
|
delete[] m_ZyyCalculated;
|
|
m_ZyyCalculated = NULL;
|
|
}
|
|
|
|
if( m_ZxxResidual != NULL){
|
|
delete[] m_ZxxResidual;
|
|
m_ZxxResidual = NULL;
|
|
}
|
|
|
|
if( m_ZxyResidual != NULL){
|
|
delete[] m_ZxyResidual;
|
|
m_ZxyResidual = NULL;
|
|
}
|
|
|
|
if( m_ZyxResidual != NULL){
|
|
delete[] m_ZyxResidual;
|
|
m_ZyxResidual = NULL;
|
|
}
|
|
|
|
if( m_ZyyResidual != NULL){
|
|
delete[] m_ZyyResidual;
|
|
m_ZyyResidual = NULL;
|
|
}
|
|
|
|
if( m_dataIDOfZxx != NULL){
|
|
delete[] m_dataIDOfZxx;
|
|
m_dataIDOfZxx = NULL;
|
|
}
|
|
|
|
if( m_dataIDOfZxy != NULL){
|
|
delete[] m_dataIDOfZxy;
|
|
m_dataIDOfZxy = NULL;
|
|
}
|
|
|
|
if( m_dataIDOfZyx != NULL){
|
|
delete[] m_dataIDOfZyx;
|
|
m_dataIDOfZyx = NULL;
|
|
}
|
|
|
|
if( m_dataIDOfZyy != NULL){
|
|
delete[] m_dataIDOfZyy;
|
|
m_dataIDOfZyy = NULL;
|
|
}
|
|
|
|
if( m_arrayDistortionMatrixDifferences != NULL){
|
|
delete m_arrayDistortionMatrixDifferences;
|
|
m_arrayDistortionMatrixDifferences = NULL;
|
|
}
|
|
|
|
if( m_arrayGainsAndRotations != NULL){
|
|
delete m_arrayGainsAndRotations;
|
|
m_arrayGainsAndRotations = NULL;
|
|
}
|
|
|
|
}
|
|
|
|
// Read data from input file
|
|
void ObservedDataStationMT::inputObservedData( std::ifstream& inFile ){
|
|
|
|
inFile >> m_stationID;
|
|
|
|
inFile >> m_IDOfMagneticFieldStation;
|
|
|
|
OutputFiles::m_logFile << "# " << std::setw(15) << std::left << m_stationID << std::setw(18) << std::left << m_IDOfMagneticFieldStation;
|
|
|
|
const AnalysisControl* const pAnalysisControl = AnalysisControl::getInstance();
|
|
// Type of owner element
|
|
int ownerType(-1);
|
|
std::string ownerElemType;
|
|
if( pAnalysisControl->isTypeOfOwnerElementSetIndivisually() ){
|
|
// Owner element type of each site is specified
|
|
inFile >> ownerType;
|
|
}
|
|
else{
|
|
ownerType = pAnalysisControl->getTypeOfOwnerElement();
|
|
}
|
|
switch(ownerType){
|
|
case AnalysisControl::USE_LOWER_ELEMENT:
|
|
m_useUpperElementForInterpolationOfEMField = false;
|
|
ownerElemType = "Lower";
|
|
break;
|
|
case AnalysisControl::USE_UPPER_ELEMENT:
|
|
m_useUpperElementForInterpolationOfEMField = true;
|
|
ownerElemType = "Upper";
|
|
break;
|
|
default:
|
|
OutputFiles::m_logFile << std::endl << "Error : Unknown type of owner element : " << ownerType << std::endl;
|
|
exit(1);
|
|
break;
|
|
}
|
|
OutputFiles::m_logFile << std::setw(15) << std::left << ownerElemType;
|
|
|
|
// Type of electric field
|
|
if( pAnalysisControl->isTypeOfElectricFieldSetIndivisually() ){
|
|
// Electric field type of each site is specified
|
|
inFile >> m_typeOfElectricField;
|
|
}else{
|
|
m_typeOfElectricField = pAnalysisControl->getTypeOfElectricField();
|
|
}
|
|
if( pAnalysisControl->getTypeOfMesh() == MeshData::HEXA ){
|
|
if( m_typeOfElectricField != AnalysisControl::USE_HORIZONTAL_ELECTRIC_FIELD ){
|
|
OutputFiles::m_logFile << std::endl << "Warning : Horizontal electric field must be used for hexahedral mesh." << std::endl;
|
|
}
|
|
m_typeOfElectricField = AnalysisControl::USE_HORIZONTAL_ELECTRIC_FIELD;
|
|
}
|
|
std::string elecType;
|
|
switch( m_typeOfElectricField ){
|
|
case AnalysisControl::USE_TANGENTIAL_ELECTRIC_FIELD:
|
|
elecType = "Tangential";
|
|
break;
|
|
case AnalysisControl::USE_HORIZONTAL_ELECTRIC_FIELD:
|
|
elecType = "Horizontal";
|
|
break;
|
|
default:
|
|
OutputFiles::m_logFile << std::endl << "Error : Unknown type of the electric field : " << m_typeOfElectricField << std::endl;
|
|
exit(1);
|
|
break;
|
|
}
|
|
OutputFiles::m_logFile << std::left << elecType << std::endl;
|
|
|
|
double dbuf(0.0);
|
|
inFile >> dbuf;
|
|
m_location.X = dbuf * CommonParameters::convKilometerToMeter;
|
|
inFile >> dbuf;
|
|
m_location.Y = dbuf * CommonParameters::convKilometerToMeter;
|
|
inFile >> m_numOfFrequency;
|
|
|
|
const int nFreq = m_numOfFrequency;
|
|
|
|
if( nFreq > 0 ){
|
|
m_freq = new double[nFreq];
|
|
m_ZxxObserved = new std::complex<double>[nFreq];
|
|
m_ZxyObserved = new std::complex<double>[nFreq];
|
|
m_ZyxObserved = new std::complex<double>[nFreq];
|
|
m_ZyyObserved = new std::complex<double>[nFreq];
|
|
m_ZxxSD = new CommonParameters::DoubleComplexValues[nFreq];
|
|
m_ZxySD = new CommonParameters::DoubleComplexValues[nFreq];
|
|
m_ZyxSD = new CommonParameters::DoubleComplexValues[nFreq];
|
|
m_ZyySD = new CommonParameters::DoubleComplexValues[nFreq];
|
|
//m_ZxxCalculated = new std::complex<double>[nFreq];
|
|
//m_ZxyCalculated = new std::complex<double>[nFreq];
|
|
//m_ZyxCalculated = new std::complex<double>[nFreq];
|
|
//m_ZyyCalculated = new std::complex<double>[nFreq];
|
|
|
|
for( int i = 0; i < nFreq; ++i ){
|
|
inFile >> m_freq[i];
|
|
double dbuf1(0.0);
|
|
double dbuf2(0.0);
|
|
inFile >> dbuf1 >> dbuf2;
|
|
m_ZxxObserved[i] = std::complex<double>(dbuf1,dbuf2);
|
|
inFile >> dbuf1 >> dbuf2;
|
|
m_ZxyObserved[i] = std::complex<double>(dbuf1,dbuf2);
|
|
inFile >> dbuf1 >> dbuf2;
|
|
m_ZyxObserved[i] = std::complex<double>(dbuf1,dbuf2);
|
|
inFile >> dbuf1 >> dbuf2;
|
|
m_ZyyObserved[i] = std::complex<double>(dbuf1,dbuf2);
|
|
inFile >> m_ZxxSD[i].realPart;
|
|
inFile >> m_ZxxSD[i].imagPart;
|
|
inFile >> m_ZxySD[i].realPart;
|
|
inFile >> m_ZxySD[i].imagPart;
|
|
inFile >> m_ZyxSD[i].realPart;
|
|
inFile >> m_ZyxSD[i].imagPart;
|
|
inFile >> m_ZyySD[i].realPart;
|
|
inFile >> m_ZyySD[i].imagPart;
|
|
}
|
|
}
|
|
|
|
#ifdef _DEBUG_WRITE
|
|
std::cout << " MT " << m_stationID << " " << m_IDOfMagneticFieldStation << std::endl;
|
|
std::cout << m_location.X << " " << m_location.Y << std::endl;
|
|
std::cout << m_numOfFrequency << std::endl;
|
|
for( int i = 0; i < m_numOfFrequency; ++i ){
|
|
std::cout << m_freq[i] << " "
|
|
<< m_ZxxObserved[i] << " "
|
|
<< m_ZxyObserved[i] << " "
|
|
<< m_ZyxObserved[i] << " "
|
|
<< m_ZyyObserved[i] << " "
|
|
<< m_ZxxSD[i].realPart << " "
|
|
<< m_ZxxSD[i].imagPart << " "
|
|
<< m_ZxySD[i].realPart << " "
|
|
<< m_ZxySD[i].imagPart << " "
|
|
<< m_ZyxSD[i].realPart << " "
|
|
<< m_ZyxSD[i].imagPart << " "
|
|
<< m_ZyySD[i].realPart << " "
|
|
<< m_ZyySD[i].imagPart << std::endl;
|
|
}
|
|
#endif
|
|
|
|
}
|
|
|
|
// Calulate electric field
|
|
void ObservedDataStationMT::calculateElectricField( const Forward3D* const ptrForward3D, const int rhsVectorIDOfEx, const int rhsVectorIDOfEy ){
|
|
|
|
const int iPol = ptrForward3D->getPolarizationCurrent();
|
|
|
|
const int meshType = ( AnalysisControl::getInstance() )->getTypeOfMesh();
|
|
if( meshType == MeshData::TETRA ){// Tetra mesh
|
|
switch( getTypeOfElectricField() ){
|
|
case AnalysisControl::USE_HORIZONTAL_ELECTRIC_FIELD:
|
|
m_ExCalculated[iPol] = ptrForward3D->calcValueElectricFieldXDirection( m_elementIncludingStation, m_volumeCoordinateValues.coord1, m_volumeCoordinateValues.coord2, m_volumeCoordinateValues.coord3 );
|
|
m_EyCalculated[iPol] = ptrForward3D->calcValueElectricFieldYDirection( m_elementIncludingStation, m_volumeCoordinateValues.coord1, m_volumeCoordinateValues.coord2, m_volumeCoordinateValues.coord3 );
|
|
break;
|
|
case AnalysisControl::USE_TANGENTIAL_ELECTRIC_FIELD:
|
|
m_ExCalculated[iPol] = ptrForward3D->calcValueElectricFieldTangentialX( m_elementIncludingStation, m_faceIncludingStation, m_areaCoordinateValues.coord1, m_areaCoordinateValues.coord2 );
|
|
m_EyCalculated[iPol] = ptrForward3D->calcValueElectricFieldTangentialY( m_elementIncludingStation, m_faceIncludingStation, m_areaCoordinateValues.coord1, m_areaCoordinateValues.coord2 );
|
|
break;
|
|
default:
|
|
OutputFiles::m_logFile << "Error : Unknown type of the electric field." << std::endl;
|
|
exit(1);
|
|
break;
|
|
}
|
|
}else if( meshType == MeshData::HEXA ){// Hexa mesh
|
|
m_ExCalculated[iPol] = ptrForward3D->calcValueElectricFieldXDirection( m_elementIncludingStation, m_localCoordinateValues.X, m_localCoordinateValues.Y, m_localCoordinateValues.Z );
|
|
m_EyCalculated[iPol] = ptrForward3D->calcValueElectricFieldYDirection( m_elementIncludingStation, m_localCoordinateValues.X, m_localCoordinateValues.Y, m_localCoordinateValues.Z );
|
|
}else if( meshType == MeshData::NONCONFORMING_HEXA ){
|
|
switch( getTypeOfElectricField() ){
|
|
case AnalysisControl::USE_HORIZONTAL_ELECTRIC_FIELD:
|
|
m_ExCalculated[iPol] = ptrForward3D->calcValueElectricFieldXDirection( m_elementIncludingStation, m_localCoordinateValues.X, m_localCoordinateValues.Y, m_localCoordinateValues.Z );
|
|
m_EyCalculated[iPol] = ptrForward3D->calcValueElectricFieldYDirection( m_elementIncludingStation, m_localCoordinateValues.X, m_localCoordinateValues.Y, m_localCoordinateValues.Z );
|
|
break;
|
|
case AnalysisControl::USE_TANGENTIAL_ELECTRIC_FIELD:
|
|
m_ExCalculated[iPol] = ptrForward3D->calcValueElectricFieldTangentialXFromAllEdges( m_elementIncludingStation, m_faceIncludingStation, m_localCoordinateValues.X, m_localCoordinateValues.Y, m_localCoordinateValues.Z );
|
|
m_EyCalculated[iPol] = ptrForward3D->calcValueElectricFieldTangentialYFromAllEdges( m_elementIncludingStation, m_faceIncludingStation, m_localCoordinateValues.X, m_localCoordinateValues.Y, m_localCoordinateValues.Z );
|
|
break;
|
|
default:
|
|
OutputFiles::m_logFile << "Error : Unknown type of the electric field." << std::endl;
|
|
exit(1);
|
|
break;
|
|
}
|
|
}else{
|
|
OutputFiles::m_logFile << "Error : Wrong type of mesh : " << meshType << std::endl;
|
|
exit(1);
|
|
}
|
|
|
|
// For inversion
|
|
m_rhsVectorIDOfEx = rhsVectorIDOfEx;
|
|
m_rhsVectorIDOfEy = rhsVectorIDOfEy;
|
|
|
|
}
|
|
|
|
// Calulate Impedance tensor
|
|
void ObservedDataStationMT::calculateImpedanceTensor( const double freq, const ObservedDataStationPoint* const ptrStationOfMagneticField, int& icount ){
|
|
|
|
const int freqIDThisPEInSta = getFreqIDsAmongThisPE( freq );
|
|
|
|
if( freqIDThisPEInSta < 0 ){// Specified frequency is not the ones calculated by this PE in this station
|
|
return;
|
|
}
|
|
|
|
const int freqIDGlobalInSta = m_freqIDsAmongThisStationCalculatedByThisPE[ freqIDThisPEInSta ];
|
|
|
|
const std::complex<double> HxCalculated[2] = { ptrStationOfMagneticField->getHxCalculated( 0 ), ptrStationOfMagneticField->getHxCalculated( 1 ) };
|
|
const std::complex<double> HyCalculated[2] = { ptrStationOfMagneticField->getHyCalculated( 0 ), ptrStationOfMagneticField->getHyCalculated( 1 ) };
|
|
|
|
const std::complex<double> det = HxCalculated[0]*HyCalculated[1] - HxCalculated[1]*HyCalculated[0];
|
|
|
|
const std::complex<double> ZxxUndist = ( m_ExCalculated[0]*HyCalculated[1] - m_ExCalculated[1]*HyCalculated[0] ) / det;
|
|
const std::complex<double> ZxyUndist = ( m_ExCalculated[1]*HxCalculated[0] - m_ExCalculated[0]*HxCalculated[1] ) / det;
|
|
const std::complex<double> ZyxUndist = ( m_EyCalculated[0]*HyCalculated[1] - m_EyCalculated[1]*HyCalculated[0] ) / det;
|
|
const std::complex<double> ZyyUndist = ( m_EyCalculated[1]*HxCalculated[0] - m_EyCalculated[0]*HxCalculated[1] ) / det;
|
|
|
|
std::complex<double> cxx( 1.0, 0.0 );
|
|
std::complex<double> cxy( 0.0, 0.0 );
|
|
std::complex<double> cyx( 0.0, 0.0 );
|
|
std::complex<double> cyy( 1.0, 0.0 );
|
|
const int typeOfDistortion = ( AnalysisControl::getInstance() )->getTypeOfDistortion();
|
|
if( typeOfDistortion == AnalysisControl::ESTIMATE_DISTORTION_MATRIX_DIFFERENCE ){
|
|
assert( m_arrayDistortionMatrixDifferences != NULL );
|
|
cxx = std::complex<double>( m_arrayDistortionMatrixDifferences->distortionMatrixDifference[ObservedDataStationMT::COMPONENT_ID_CXX] + 1.0 , 0.0 );
|
|
cxy = std::complex<double>( m_arrayDistortionMatrixDifferences->distortionMatrixDifference[ObservedDataStationMT::COMPONENT_ID_CXY] , 0.0 );
|
|
cyx = std::complex<double>( m_arrayDistortionMatrixDifferences->distortionMatrixDifference[ObservedDataStationMT::COMPONENT_ID_CYX] , 0.0 );
|
|
cyy = std::complex<double>( m_arrayDistortionMatrixDifferences->distortionMatrixDifference[ObservedDataStationMT::COMPONENT_ID_CYY] + 1.0 , 0.0 );
|
|
}
|
|
else if( typeOfDistortion == AnalysisControl::ESTIMATE_GAINS_AND_ROTATIONS ){
|
|
assert( m_arrayGainsAndRotations != NULL );
|
|
const double gX = pow( 10.0, m_arrayGainsAndRotations->gainsAndRotations[ObservedDataStationMT::EX_GAIN] );
|
|
const double gY = pow( 10.0, m_arrayGainsAndRotations->gainsAndRotations[ObservedDataStationMT::EY_GAIN] );
|
|
const double betaX = m_arrayGainsAndRotations->gainsAndRotations[ObservedDataStationMT::EX_ROTATION];
|
|
const double betaY = m_arrayGainsAndRotations->gainsAndRotations[ObservedDataStationMT::EY_ROTATION];
|
|
|
|
cxx = std::complex<double>( gX * cos( betaX ), 0.0 );
|
|
cxy = std::complex<double>( - gY * sin( betaY ), 0.0 );
|
|
cyx = std::complex<double>( gX * sin( betaX ), 0.0 );
|
|
cyy = std::complex<double>( gY * cos( betaY ), 0.0 );
|
|
}
|
|
else if( typeOfDistortion == AnalysisControl::ESTIMATE_GAINS_ONLY ){
|
|
assert( m_arrayGainsAndRotations != NULL );
|
|
const double gX = pow( 10.0, m_arrayGainsAndRotations->gainsAndRotations[ObservedDataStationMT::EX_GAIN] );
|
|
const double gY = pow( 10.0, m_arrayGainsAndRotations->gainsAndRotations[ObservedDataStationMT::EY_GAIN] );
|
|
|
|
cxx = std::complex<double>( gX, 0.0 );
|
|
cxy = std::complex<double>( 0.0, 0.0 );
|
|
cyx = std::complex<double>( 0.0, 0.0 );
|
|
cyy = std::complex<double>( gY, 0.0 );
|
|
}
|
|
|
|
m_ZxxCalculated[freqIDThisPEInSta] = cxx * ZxxUndist + cxy * ZyxUndist;
|
|
m_ZxyCalculated[freqIDThisPEInSta] = cxx * ZxyUndist + cxy * ZyyUndist;
|
|
m_ZyxCalculated[freqIDThisPEInSta] = cyx * ZxxUndist + cyy * ZyxUndist;
|
|
m_ZyyCalculated[freqIDThisPEInSta] = cyx * ZxyUndist + cyy * ZyyUndist;
|
|
|
|
m_ZxxResidual[freqIDThisPEInSta].realPart = ( m_ZxxObserved[freqIDGlobalInSta].real() - m_ZxxCalculated[freqIDThisPEInSta].real() ) / m_ZxxSD[freqIDGlobalInSta].realPart;
|
|
m_ZxxResidual[freqIDThisPEInSta].imagPart = ( m_ZxxObserved[freqIDGlobalInSta].imag() - m_ZxxCalculated[freqIDThisPEInSta].imag() ) / m_ZxxSD[freqIDGlobalInSta].imagPart;
|
|
m_ZxyResidual[freqIDThisPEInSta].realPart = ( m_ZxyObserved[freqIDGlobalInSta].real() - m_ZxyCalculated[freqIDThisPEInSta].real() ) / m_ZxySD[freqIDGlobalInSta].realPart;
|
|
m_ZxyResidual[freqIDThisPEInSta].imagPart = ( m_ZxyObserved[freqIDGlobalInSta].imag() - m_ZxyCalculated[freqIDThisPEInSta].imag() ) / m_ZxySD[freqIDGlobalInSta].imagPart;
|
|
m_ZyxResidual[freqIDThisPEInSta].realPart = ( m_ZyxObserved[freqIDGlobalInSta].real() - m_ZyxCalculated[freqIDThisPEInSta].real() ) / m_ZyxSD[freqIDGlobalInSta].realPart;
|
|
m_ZyxResidual[freqIDThisPEInSta].imagPart = ( m_ZyxObserved[freqIDGlobalInSta].imag() - m_ZyxCalculated[freqIDThisPEInSta].imag() ) / m_ZyxSD[freqIDGlobalInSta].imagPart;
|
|
m_ZyyResidual[freqIDThisPEInSta].realPart = ( m_ZyyObserved[freqIDGlobalInSta].real() - m_ZyyCalculated[freqIDThisPEInSta].real() ) / m_ZyySD[freqIDGlobalInSta].realPart;
|
|
m_ZyyResidual[freqIDThisPEInSta].imagPart = ( m_ZyyObserved[freqIDGlobalInSta].imag() - m_ZyyCalculated[freqIDThisPEInSta].imag() ) / m_ZyySD[freqIDGlobalInSta].imagPart;
|
|
|
|
#ifdef _DEBUG_WRITE
|
|
std::cout << "freqIDThisPEInSta Zxx Zxy Zyx Zyy : " << freqIDThisPEInSta << " " << m_ZxxCalculated[freqIDThisPEInSta] << " " << m_ZxyCalculated[freqIDThisPEInSta] << " " << m_ZyxCalculated[freqIDThisPEInSta] << " " << m_ZyyCalculated[freqIDThisPEInSta] << std::endl;
|
|
#endif
|
|
|
|
// For inversion
|
|
m_dataIDOfZxx[freqIDThisPEInSta].realPart = icount++;
|
|
m_dataIDOfZxx[freqIDThisPEInSta].imagPart = icount++;
|
|
|
|
m_dataIDOfZxy[freqIDThisPEInSta].realPart = icount++;
|
|
m_dataIDOfZxy[freqIDThisPEInSta].imagPart = icount++;
|
|
|
|
m_dataIDOfZyx[freqIDThisPEInSta].realPart = icount++;
|
|
m_dataIDOfZyx[freqIDThisPEInSta].imagPart = icount++;
|
|
|
|
m_dataIDOfZyy[freqIDThisPEInSta].realPart = icount++;
|
|
m_dataIDOfZyy[freqIDThisPEInSta].imagPart = icount++;
|
|
|
|
#ifdef _DEBUG_WRITE
|
|
std::cout << "m_dataIDOfZxx[freqIDThisPEInSta].realPart : " << m_dataIDOfZxx[freqIDThisPEInSta].realPart << std::endl;
|
|
std::cout << "m_dataIDOfZxx[freqIDThisPEInSta].imagPart : " << m_dataIDOfZxx[freqIDThisPEInSta].imagPart << std::endl;
|
|
std::cout << "m_dataIDOfZxy[freqIDThisPEInSta].realPart : " << m_dataIDOfZxy[freqIDThisPEInSta].realPart << std::endl;
|
|
std::cout << "m_dataIDOfZxy[freqIDThisPEInSta].imagPart : " << m_dataIDOfZxy[freqIDThisPEInSta].imagPart << std::endl;
|
|
std::cout << "m_dataIDOfZyx[freqIDThisPEInSta].realPart : " << m_dataIDOfZyx[freqIDThisPEInSta].realPart << std::endl;
|
|
std::cout << "m_dataIDOfZyx[freqIDThisPEInSta].imagPart : " << m_dataIDOfZyx[freqIDThisPEInSta].imagPart << std::endl;
|
|
std::cout << "m_dataIDOfZyy[freqIDThisPEInSta].realPart : " << m_dataIDOfZyy[freqIDThisPEInSta].realPart << std::endl;
|
|
std::cout << "m_dataIDOfZyy[freqIDThisPEInSta].imagPart : " << m_dataIDOfZyy[freqIDThisPEInSta].imagPart << std::endl;
|
|
#endif
|
|
|
|
}
|
|
|
|
// Initialize electric field
|
|
void ObservedDataStationMT::initializeElectricField( const int iPol ){
|
|
|
|
m_ExCalculated[iPol] = std::complex<double>(0.0,0.0);
|
|
m_EyCalculated[iPol] = std::complex<double>(0.0,0.0);
|
|
|
|
}
|
|
|
|
// Initialize Impedance tensor and errors
|
|
void ObservedDataStationMT::initializeImpedanceTensorsAndErrors(){
|
|
|
|
//for( int i = 0; i < m_numOfFrequency; ++i ){
|
|
for( int i = 0; i < m_numOfFreqCalculatedByThisStaAndPE; ++i ){
|
|
m_ZxxResidual[i].realPart = 0.0;
|
|
m_ZxxResidual[i].imagPart = 0.0;
|
|
m_ZxyResidual[i].realPart = 0.0;
|
|
m_ZxyResidual[i].imagPart = 0.0;
|
|
m_ZyxResidual[i].realPart = 0.0;
|
|
m_ZyxResidual[i].imagPart = 0.0;
|
|
m_ZyyResidual[i].realPart = 0.0;
|
|
m_ZyyResidual[i].imagPart = 0.0;
|
|
|
|
m_ZxxCalculated[i] = std::complex<double>(0.0,0.0);
|
|
m_ZxyCalculated[i] = std::complex<double>(0.0,0.0);
|
|
m_ZyxCalculated[i] = std::complex<double>(0.0,0.0);
|
|
m_ZyyCalculated[i] = std::complex<double>(0.0,0.0);
|
|
}
|
|
|
|
}
|
|
|
|
// Allocate memory for the calculated Impedance and errors
|
|
void ObservedDataStationMT::allocateMemoryForCalculatedValues(){
|
|
|
|
if( m_ZxxCalculated != NULL){
|
|
delete[] m_ZxxCalculated;
|
|
m_ZxxCalculated = NULL;
|
|
}
|
|
|
|
if( m_ZxyCalculated != NULL){
|
|
delete[] m_ZxyCalculated;
|
|
m_ZxyCalculated = NULL;
|
|
}
|
|
|
|
if( m_ZyxCalculated != NULL){
|
|
delete[] m_ZyxCalculated;
|
|
m_ZyxCalculated = NULL;
|
|
}
|
|
|
|
if( m_ZyyCalculated != NULL){
|
|
delete[] m_ZyyCalculated;
|
|
m_ZyyCalculated = NULL;
|
|
}
|
|
|
|
if( m_ZxxResidual != NULL){
|
|
delete[] m_ZxxResidual;
|
|
m_ZxxResidual = NULL;
|
|
}
|
|
|
|
if( m_ZxyResidual != NULL){
|
|
delete[] m_ZxyResidual;
|
|
m_ZxyResidual = NULL;
|
|
}
|
|
|
|
if( m_ZyxResidual != NULL){
|
|
delete[] m_ZyxResidual;
|
|
m_ZyxResidual = NULL;
|
|
}
|
|
|
|
if( m_ZyyResidual != NULL){
|
|
delete[] m_ZyyResidual;
|
|
m_ZyyResidual = NULL;
|
|
}
|
|
|
|
if( m_dataIDOfZxx != NULL){
|
|
delete[] m_dataIDOfZxx;
|
|
m_dataIDOfZxx = NULL;
|
|
}
|
|
|
|
if( m_dataIDOfZxy != NULL){
|
|
delete[] m_dataIDOfZxy;
|
|
m_dataIDOfZxy = NULL;
|
|
}
|
|
|
|
if( m_dataIDOfZyx != NULL){
|
|
delete[] m_dataIDOfZyx;
|
|
m_dataIDOfZyx = NULL;
|
|
}
|
|
|
|
if( m_dataIDOfZyy != NULL){
|
|
delete[] m_dataIDOfZyy;
|
|
m_dataIDOfZyy = NULL;
|
|
}
|
|
|
|
if( m_numOfFreqCalculatedByThisStaAndPE > 0 ){
|
|
|
|
// If total number of frequencies calculated by this PE is not more than zero, do not allocate memory
|
|
m_ZxxCalculated = new std::complex<double>[m_numOfFreqCalculatedByThisStaAndPE];
|
|
m_ZxyCalculated = new std::complex<double>[m_numOfFreqCalculatedByThisStaAndPE];
|
|
m_ZyxCalculated = new std::complex<double>[m_numOfFreqCalculatedByThisStaAndPE];
|
|
m_ZyyCalculated = new std::complex<double>[m_numOfFreqCalculatedByThisStaAndPE];
|
|
m_ZxxResidual = new CommonParameters::DoubleComplexValues[m_numOfFreqCalculatedByThisStaAndPE];
|
|
m_ZxyResidual = new CommonParameters::DoubleComplexValues[m_numOfFreqCalculatedByThisStaAndPE];
|
|
m_ZyxResidual = new CommonParameters::DoubleComplexValues[m_numOfFreqCalculatedByThisStaAndPE];
|
|
m_ZyyResidual = new CommonParameters::DoubleComplexValues[m_numOfFreqCalculatedByThisStaAndPE];
|
|
m_dataIDOfZxx = new CommonParameters::InitComplexValues[m_numOfFreqCalculatedByThisStaAndPE];
|
|
m_dataIDOfZxy = new CommonParameters::InitComplexValues[m_numOfFreqCalculatedByThisStaAndPE];
|
|
m_dataIDOfZyx = new CommonParameters::InitComplexValues[m_numOfFreqCalculatedByThisStaAndPE];
|
|
m_dataIDOfZyy = new CommonParameters::InitComplexValues[m_numOfFreqCalculatedByThisStaAndPE];
|
|
|
|
// Initialize
|
|
for( int i = 0; i < m_numOfFreqCalculatedByThisStaAndPE; ++i ){
|
|
m_ZxxCalculated[i] = std::complex<double>(0.0,0.0);
|
|
m_ZxyCalculated[i] = std::complex<double>(0.0,0.0);
|
|
m_ZyxCalculated[i] = std::complex<double>(0.0,0.0);
|
|
m_ZyyCalculated[i] = std::complex<double>(0.0,0.0);
|
|
m_ZxxResidual[i].realPart = 0.0;
|
|
m_ZxxResidual[i].imagPart = 0.0;
|
|
m_ZxyResidual[i].realPart = 0.0;
|
|
m_ZxyResidual[i].imagPart = 0.0;
|
|
m_ZyxResidual[i].realPart = 0.0;
|
|
m_ZyxResidual[i].imagPart = 0.0;
|
|
m_ZyyResidual[i].realPart = 0.0;
|
|
m_ZyyResidual[i].imagPart = 0.0;
|
|
m_dataIDOfZxx[i].realPart = -1;
|
|
m_dataIDOfZxx[i].imagPart = -1;
|
|
m_dataIDOfZxy[i].realPart = -1;
|
|
m_dataIDOfZxy[i].imagPart = -1;
|
|
m_dataIDOfZyx[i].realPart = -1;
|
|
m_dataIDOfZyx[i].imagPart = -1;
|
|
m_dataIDOfZyy[i].realPart = -1;
|
|
m_dataIDOfZyy[i].imagPart = -1;
|
|
}
|
|
|
|
}
|
|
|
|
}
|
|
|
|
// Output calculated Impedance tensors
|
|
void ObservedDataStationMT::outputCalculatedValues() const{
|
|
|
|
int icount(0);
|
|
for( std::vector<int>::const_iterator itr = m_freqIDsAmongThisStationCalculatedByThisPE.begin(); itr != m_freqIDsAmongThisStationCalculatedByThisPE.end(); ++itr ){
|
|
fprintf( OutputFiles::m_csvFile, "%10d, %15.6e, %15.6e, %15.6e, %15.6e, %15.6e, %15.6e, %15.6e, %15.6e, %15.6e,",
|
|
m_stationID, m_freq[*itr], m_ZxxCalculated[icount].real(), m_ZxxCalculated[icount].imag(), m_ZxyCalculated[icount].real(), m_ZxyCalculated[icount].imag(),
|
|
m_ZyxCalculated[icount].real(), m_ZyxCalculated[icount].imag(), m_ZyyCalculated[icount].real(), m_ZyyCalculated[icount].imag());
|
|
fprintf( OutputFiles::m_csvFile, " %15.6e, %15.6e, %15.6e, %15.6e, %15.6e, %15.6e, %15.6e, %15.6e,",
|
|
m_ZxxResidual[icount].realPart, m_ZxxResidual[icount].imagPart, m_ZxyResidual[icount].realPart, m_ZxyResidual[icount].imagPart,
|
|
m_ZyxResidual[icount].realPart, m_ZyxResidual[icount].imagPart, m_ZyyResidual[icount].realPart, m_ZyyResidual[icount].imagPart );
|
|
fprintf( OutputFiles::m_csvFile, " %15.6e, %15.6e, %15.6e, %15.6e, %15.6e, %15.6e, %15.6e, %15.6e,",
|
|
m_ZxxObserved[*itr].real(),m_ZxxObserved[*itr].imag(),m_ZxyObserved[*itr].real(),m_ZxyObserved[*itr].imag(),
|
|
m_ZyxObserved[*itr].real(),m_ZyxObserved[*itr].imag(),m_ZyyObserved[*itr].real(),m_ZyyObserved[*itr].imag() );
|
|
fprintf( OutputFiles::m_csvFile, " %15.6e, %15.6e, %15.6e, %15.6e, %15.6e, %15.6e, %15.6e, %15.6e,\n",
|
|
m_ZxxSD[*itr].realPart,m_ZxxSD[*itr].imagPart,m_ZxySD[*itr].realPart,m_ZxySD[*itr].imagPart,
|
|
m_ZyxSD[*itr].realPart,m_ZyxSD[*itr].imagPart,m_ZyySD[*itr].realPart,m_ZyySD[*itr].imagPart );
|
|
|
|
++icount;
|
|
}
|
|
|
|
}
|
|
|
|
// Calulate interpolator vector of electric field
|
|
void ObservedDataStationMT::calcInterpolatorVectorOfElectricField( Forward3D* const ptrForward3D ){
|
|
|
|
const int meshType = ( AnalysisControl::getInstance() )->getTypeOfMesh();
|
|
if( meshType == MeshData::TETRA ){// Tetra mesh
|
|
switch( getTypeOfElectricField() ){
|
|
case AnalysisControl::USE_HORIZONTAL_ELECTRIC_FIELD:
|
|
ptrForward3D->calcInterpolatorVectorOfElectricFieldXDirection( m_elementIncludingStation, m_volumeCoordinateValues.coord1, m_volumeCoordinateValues.coord2, m_volumeCoordinateValues.coord3, m_rhsVectorIDOfEx );
|
|
ptrForward3D->calcInterpolatorVectorOfElectricFieldYDirection( m_elementIncludingStation, m_volumeCoordinateValues.coord1, m_volumeCoordinateValues.coord2, m_volumeCoordinateValues.coord3, m_rhsVectorIDOfEy );
|
|
break;
|
|
case AnalysisControl::USE_TANGENTIAL_ELECTRIC_FIELD:
|
|
ptrForward3D->calcInterpolatorVectorOfElectricFieldTangentialX( m_elementIncludingStation, m_faceIncludingStation, m_areaCoordinateValues.coord1, m_areaCoordinateValues.coord2, m_rhsVectorIDOfEx );
|
|
ptrForward3D->calcInterpolatorVectorOfElectricFieldTangentialY( m_elementIncludingStation, m_faceIncludingStation, m_areaCoordinateValues.coord1, m_areaCoordinateValues.coord2, m_rhsVectorIDOfEy );
|
|
break;
|
|
default:
|
|
OutputFiles::m_logFile << "Error : Unknown type of the electric field." << std::endl;
|
|
exit(1);
|
|
break;
|
|
}
|
|
}else if( meshType == MeshData::HEXA){// Hexa mesh
|
|
ptrForward3D->calcInterpolatorVectorOfElectricFieldXDirection( m_elementIncludingStation, m_localCoordinateValues.X, m_localCoordinateValues.Y, m_localCoordinateValues.Z, m_rhsVectorIDOfEx );
|
|
ptrForward3D->calcInterpolatorVectorOfElectricFieldYDirection( m_elementIncludingStation, m_localCoordinateValues.X, m_localCoordinateValues.Y, m_localCoordinateValues.Z, m_rhsVectorIDOfEy );
|
|
}else if( meshType == MeshData::NONCONFORMING_HEXA ){
|
|
switch( getTypeOfElectricField() ){
|
|
case AnalysisControl::USE_HORIZONTAL_ELECTRIC_FIELD:
|
|
ptrForward3D->calcInterpolatorVectorOfElectricFieldXDirection( m_elementIncludingStation, m_localCoordinateValues.X, m_localCoordinateValues.Y, m_localCoordinateValues.Z, m_rhsVectorIDOfEx );
|
|
ptrForward3D->calcInterpolatorVectorOfElectricFieldYDirection( m_elementIncludingStation, m_localCoordinateValues.X, m_localCoordinateValues.Y, m_localCoordinateValues.Z, m_rhsVectorIDOfEy );
|
|
break;
|
|
case AnalysisControl::USE_TANGENTIAL_ELECTRIC_FIELD:
|
|
ptrForward3D->calcInterpolatorVectorOfElectricFieldTangentialXFromAllEdges( m_elementIncludingStation, m_faceIncludingStation, m_localCoordinateValues.X, m_localCoordinateValues.Y, m_localCoordinateValues.Z, m_rhsVectorIDOfEx );
|
|
ptrForward3D->calcInterpolatorVectorOfElectricFieldTangentialYFromAllEdges( m_elementIncludingStation, m_faceIncludingStation, m_localCoordinateValues.X, m_localCoordinateValues.Y, m_localCoordinateValues.Z, m_rhsVectorIDOfEy );
|
|
break;
|
|
default:
|
|
OutputFiles::m_logFile << "Error : Unknown type of the electric field." << std::endl;
|
|
exit(1);
|
|
break;
|
|
}
|
|
}else{
|
|
OutputFiles::m_logFile << "Error : Wrong type of mesh : " << meshType << std::endl;
|
|
exit(1);
|
|
}
|
|
|
|
}
|
|
|
|
// Calulate sensitivity matrix of Impedance tensors
|
|
void ObservedDataStationMT::calculateSensitivityMatrix( const double freq, const int nModel,
|
|
const ObservedDataStationPoint* const ptrStationOfMagneticField,
|
|
const std::complex<double>* const derivativesOfEMFieldExPol,
|
|
const std::complex<double>* const derivativesOfEMFieldEyPol,
|
|
double* const sensitivityMatrix, const bool forceSDToOne ) const{
|
|
|
|
//const int freqID = getFreqIDs( freq );
|
|
//const int freqIDAmongThisStationAndPE = getFreqIDsAmongThisPE( freq );
|
|
|
|
const int freqIDThisPEInSta = getFreqIDsAmongThisPE( freq );
|
|
|
|
if( freqIDThisPEInSta < 0 ){// Specified frequency is not the ones calculated by this PE in this station
|
|
return;
|
|
}
|
|
|
|
const int freqIDGlobalInSta = m_freqIDsAmongThisStationCalculatedByThisPE[ freqIDThisPEInSta ];
|
|
|
|
const std::complex<double> HxCalculated[2] = { ptrStationOfMagneticField->getHxCalculated( 0 ), ptrStationOfMagneticField->getHxCalculated( 1 ) };
|
|
const std::complex<double> HyCalculated[2] = { ptrStationOfMagneticField->getHyCalculated( 0 ), ptrStationOfMagneticField->getHyCalculated( 1 ) };
|
|
|
|
const std::complex<double> det = HxCalculated[0]*HyCalculated[1] - HxCalculated[1]*HyCalculated[0];
|
|
const std::complex<double> divDet = std::complex<double>(1.0,0.0) / det;
|
|
const std::complex<double> divDet2 = divDet * divDet;
|
|
|
|
const long long rhsVectorIDOfHx = static_cast<long long>( ptrStationOfMagneticField->getRhsVectorIDOfHx() );
|
|
const long long rhsVectorIDOfHy = static_cast<long long>( ptrStationOfMagneticField->getRhsVectorIDOfHy() );
|
|
const long long nBlkNotFixed = static_cast<long long>( ( ResistivityBlock::getInstance() )->getNumResistivityBlockNotFixed() );
|
|
|
|
// //----- debug >>>>>
|
|
//#ifdef _DEBUG_WRITE
|
|
// std::cout << "divDet " << divDet << std::endl;
|
|
// for( int imdl = 0; imdl < nBlkNotFixed; ++imdl ){
|
|
// std::cout << "Ex imdl derivatives[0] " << imdl << " " << derivativesOfEMFieldExPol[ nBlkNotFixed*m_rhsVectorIDOfEx + imdl ] << std::endl;
|
|
// }
|
|
// for( int imdl = 0; imdl < nBlkNotFixed; ++imdl ){
|
|
// std::cout << "Ey imdl derivatives[0] " << imdl << " " << derivativesOfEMFieldExPol[ nBlkNotFixed*m_rhsVectorIDOfEy + imdl ] << std::endl;
|
|
// }
|
|
// for( int imdl = 0; imdl < nBlkNotFixed; ++imdl ){
|
|
// std::cout << "Hx imdl derivatives[0] " << imdl << " " << derivativesOfEMFieldExPol[ nBlkNotFixed*rhsVectorIDOfHx + imdl ] << std::endl;
|
|
// }
|
|
// for( int imdl = 0; imdl < nBlkNotFixed; ++imdl ){
|
|
// std::cout << "Hy imdl derivatives[0] " << imdl << " " << derivativesOfEMFieldExPol[ nBlkNotFixed*rhsVectorIDOfHy + imdl ] << std::endl;
|
|
// }
|
|
// for( int imdl = 0; imdl < nBlkNotFixed; ++imdl ){
|
|
// std::cout << "Ex imdl derivatives[1] " << imdl << " " << derivativesOfEMFieldEyPol[ nBlkNotFixed*m_rhsVectorIDOfEx + imdl ] << std::endl;
|
|
// }
|
|
// for( int imdl = 0; imdl < nBlkNotFixed; ++imdl ){
|
|
// std::cout << "Ey imdl derivatives[1] " << imdl << " " << derivativesOfEMFieldEyPol[ nBlkNotFixed*m_rhsVectorIDOfEy + imdl ] << std::endl;
|
|
// }
|
|
// for( int imdl = 0; imdl < nBlkNotFixed; ++imdl ){
|
|
// std::cout << "Hx imdl derivatives[1] " << imdl << " " << derivativesOfEMFieldEyPol[ nBlkNotFixed*rhsVectorIDOfHx + imdl ] << std::endl;
|
|
// }
|
|
// for( int imdl = 0; imdl < nBlkNotFixed; ++imdl ){
|
|
// std::cout << "Hy imdl derivatives[1] " << imdl << " " << derivativesOfEMFieldEyPol[ nBlkNotFixed*rhsVectorIDOfHy + imdl ] << std::endl;
|
|
// }
|
|
//#endif
|
|
// //----- debug <<<<<
|
|
|
|
//const double baseOfStaticShift = ( ObservedData::getInstance() )->getBaseOfStaticShift();
|
|
|
|
//const double factor( pow( 10.0, m_staticShiftFactor ) );
|
|
|
|
double cxx = 1.0;
|
|
double cxy = 0.0;
|
|
double cyx = 0.0;
|
|
double cyy = 1.0;
|
|
if( !m_fixDistortionMatrix ){ // Distortion matrix is not fixed
|
|
const int typeOfDistortion = ( AnalysisControl::getInstance() )->getTypeOfDistortion();
|
|
if( typeOfDistortion == AnalysisControl::ESTIMATE_DISTORTION_MATRIX_DIFFERENCE ){
|
|
assert( m_arrayDistortionMatrixDifferences != NULL );
|
|
cxx = m_arrayDistortionMatrixDifferences->distortionMatrixDifference[ObservedDataStationMT::COMPONENT_ID_CXX] + 1.0;
|
|
cxy = m_arrayDistortionMatrixDifferences->distortionMatrixDifference[ObservedDataStationMT::COMPONENT_ID_CXY];
|
|
cyx = m_arrayDistortionMatrixDifferences->distortionMatrixDifference[ObservedDataStationMT::COMPONENT_ID_CYX];
|
|
cyy = m_arrayDistortionMatrixDifferences->distortionMatrixDifference[ObservedDataStationMT::COMPONENT_ID_CYY] + 1.0;
|
|
}
|
|
else if( typeOfDistortion == AnalysisControl::ESTIMATE_GAINS_AND_ROTATIONS ){
|
|
assert( m_arrayGainsAndRotations != NULL );
|
|
const double gX = pow( 10.0, m_arrayGainsAndRotations->gainsAndRotations[ObservedDataStationMT::EX_GAIN] );
|
|
const double gY = pow( 10.0, m_arrayGainsAndRotations->gainsAndRotations[ObservedDataStationMT::EY_GAIN] );
|
|
const double betaX = m_arrayGainsAndRotations->gainsAndRotations[ObservedDataStationMT::EX_ROTATION];
|
|
const double betaY = m_arrayGainsAndRotations->gainsAndRotations[ObservedDataStationMT::EY_ROTATION];
|
|
cxx = gX * cos( betaX );
|
|
cxy = - gY * sin( betaY );
|
|
cyx = gX * sin( betaX );
|
|
cyy = gY * cos( betaY );
|
|
}
|
|
else if( typeOfDistortion == AnalysisControl::ESTIMATE_GAINS_ONLY ){
|
|
assert( m_arrayGainsAndRotations != NULL );
|
|
const double gX = pow( 10.0, m_arrayGainsAndRotations->gainsAndRotations[ObservedDataStationMT::EX_GAIN] );
|
|
const double gY = pow( 10.0, m_arrayGainsAndRotations->gainsAndRotations[ObservedDataStationMT::EY_GAIN] );
|
|
cxx = gX;
|
|
cxy = 0.0;
|
|
cyx = 0.0;
|
|
cyy = gY;
|
|
}
|
|
}
|
|
|
|
CommonParameters::DoubleComplexValues ZxxSD = m_ZxxSD[freqIDGlobalInSta];
|
|
CommonParameters::DoubleComplexValues ZxySD = m_ZxySD[freqIDGlobalInSta];
|
|
CommonParameters::DoubleComplexValues ZyxSD = m_ZyxSD[freqIDGlobalInSta];
|
|
CommonParameters::DoubleComplexValues ZyySD = m_ZyySD[freqIDGlobalInSta];
|
|
if( forceSDToOne ){
|
|
// force erro to one
|
|
ZxxSD.realPart = 1.0;
|
|
ZxySD.realPart = 1.0;
|
|
ZyxSD.realPart = 1.0;
|
|
ZyySD.realPart = 1.0;
|
|
ZxxSD.imagPart = 1.0;
|
|
ZxySD.imagPart = 1.0;
|
|
ZyxSD.imagPart = 1.0;
|
|
ZyySD.imagPart = 1.0;
|
|
}
|
|
|
|
for( long long imdl = 0; imdl < nBlkNotFixed; ++imdl ){
|
|
|
|
const std::complex<double> work1 = derivativesOfEMFieldExPol[ nBlkNotFixed * rhsVectorIDOfHx + imdl ] * HyCalculated[1]
|
|
+ derivativesOfEMFieldEyPol[ nBlkNotFixed * rhsVectorIDOfHy + imdl ] * HxCalculated[0]
|
|
- derivativesOfEMFieldEyPol[ nBlkNotFixed * rhsVectorIDOfHx + imdl ] * HyCalculated[0]
|
|
- derivativesOfEMFieldExPol[ nBlkNotFixed * rhsVectorIDOfHy + imdl ] * HxCalculated[1];
|
|
|
|
// dZxx/dm
|
|
const std::complex<double> workXX1 = derivativesOfEMFieldExPol[ nBlkNotFixed * static_cast<long long>(m_rhsVectorIDOfEx) + imdl ] * HyCalculated[1]
|
|
+ derivativesOfEMFieldEyPol[ nBlkNotFixed * rhsVectorIDOfHy + imdl ] * m_ExCalculated[0]
|
|
- derivativesOfEMFieldEyPol[ nBlkNotFixed * static_cast<long long>(m_rhsVectorIDOfEx) + imdl ] * HyCalculated[0]
|
|
- derivativesOfEMFieldExPol[ nBlkNotFixed * rhsVectorIDOfHy + imdl ] * m_ExCalculated[1];
|
|
|
|
const std::complex<double> workXX2 = m_ExCalculated[0]*HyCalculated[1] - m_ExCalculated[1]*HyCalculated[0];
|
|
|
|
const double dZxxRealUndist = ( workXX1 * divDet - work1 * workXX2 * divDet2 ).real();
|
|
const double dZxxImagUndist = ( workXX1 * divDet - work1 * workXX2 * divDet2 ).imag();
|
|
|
|
// dZxy/dm
|
|
const std::complex<double> workXY1 = derivativesOfEMFieldEyPol[ nBlkNotFixed * static_cast<long long>(m_rhsVectorIDOfEx) + imdl ] * HxCalculated[0]
|
|
+ derivativesOfEMFieldExPol[ nBlkNotFixed * rhsVectorIDOfHx + imdl ] * m_ExCalculated[1]
|
|
- derivativesOfEMFieldExPol[ nBlkNotFixed * static_cast<long long>(m_rhsVectorIDOfEx) + imdl ] * HxCalculated[1]
|
|
- derivativesOfEMFieldEyPol[ nBlkNotFixed * rhsVectorIDOfHx + imdl ] * m_ExCalculated[0];
|
|
|
|
const std::complex<double> workXY2 = m_ExCalculated[1]*HxCalculated[0] - m_ExCalculated[0]*HxCalculated[1];
|
|
|
|
const double dZxyRealUndist = ( workXY1 * divDet - work1 * workXY2 * divDet2 ).real();
|
|
const double dZxyImagUndist = ( workXY1 * divDet - work1 * workXY2 * divDet2 ).imag();
|
|
|
|
// dZyx/dm
|
|
const std::complex<double> workYX1 = derivativesOfEMFieldExPol[ nBlkNotFixed * static_cast<long long>(m_rhsVectorIDOfEy) + imdl ] * HyCalculated[1]
|
|
+ derivativesOfEMFieldEyPol[ nBlkNotFixed * rhsVectorIDOfHy + imdl ] * m_EyCalculated[0]
|
|
- derivativesOfEMFieldEyPol[ nBlkNotFixed * static_cast<long long>(m_rhsVectorIDOfEy) + imdl ] * HyCalculated[0]
|
|
- derivativesOfEMFieldExPol[ nBlkNotFixed * rhsVectorIDOfHy + imdl ] * m_EyCalculated[1];
|
|
|
|
const std::complex<double> workYX2 = m_EyCalculated[0]*HyCalculated[1] - m_EyCalculated[1]*HyCalculated[0];
|
|
|
|
const double dZyxRealUndist = ( workYX1 * divDet - work1 * workYX2 * divDet2 ).real();
|
|
const double dZyxImagUndist = ( workYX1 * divDet - work1 * workYX2 * divDet2 ).imag();
|
|
|
|
// dZyy/dm
|
|
const std::complex<double> workYY1 = derivativesOfEMFieldEyPol[ nBlkNotFixed * static_cast<long long>(m_rhsVectorIDOfEy) + imdl ] * HxCalculated[0]
|
|
+ derivativesOfEMFieldExPol[ nBlkNotFixed * rhsVectorIDOfHx + imdl ] * m_EyCalculated[1]
|
|
- derivativesOfEMFieldExPol[ nBlkNotFixed * static_cast<long long>(m_rhsVectorIDOfEy) + imdl ] * HxCalculated[1]
|
|
- derivativesOfEMFieldEyPol[ nBlkNotFixed * rhsVectorIDOfHx + imdl ] * m_EyCalculated[0];
|
|
|
|
const std::complex<double> workYY2 = m_EyCalculated[1]*HxCalculated[0] - m_EyCalculated[0]*HxCalculated[1];
|
|
|
|
const double dZyyRealUndist = ( workYY1 * divDet - work1 * workYY2 * divDet2 ).real();
|
|
const double dZyyImagUndist = ( workYY1 * divDet - work1 * workYY2 * divDet2 ).imag();
|
|
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxx[freqIDThisPEInSta].realPart) + imdl ] = ( cxx * dZxxRealUndist + cxy * dZyxRealUndist ) / ZxxSD.realPart;
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxy[freqIDThisPEInSta].realPart) + imdl ] = ( cxx * dZxyRealUndist + cxy * dZyyRealUndist ) / ZxySD.realPart;
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyx[freqIDThisPEInSta].realPart) + imdl ] = ( cyx * dZxxRealUndist + cyy * dZyxRealUndist ) / ZyxSD.realPart;
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyy[freqIDThisPEInSta].realPart) + imdl ] = ( cyx * dZxyRealUndist + cyy * dZyyRealUndist ) / ZyySD.realPart;
|
|
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxx[freqIDThisPEInSta].imagPart) + imdl ] = ( cxx * dZxxImagUndist + cxy * dZyxImagUndist ) / ZxxSD.imagPart;
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxy[freqIDThisPEInSta].imagPart) + imdl ] = ( cxx * dZxyImagUndist + cxy * dZyyImagUndist ) / ZxySD.imagPart;
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyx[freqIDThisPEInSta].imagPart) + imdl ] = ( cyx * dZxxImagUndist + cyy * dZyxImagUndist ) / ZyxSD.imagPart;
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyy[freqIDThisPEInSta].imagPart) + imdl ] = ( cyx * dZxyImagUndist + cyy * dZyyImagUndist ) / ZyySD.imagPart;
|
|
|
|
if( !m_fixDistortionMatrix ){ // Distortion matrix is not fixed
|
|
|
|
const std::complex<double> ZxxUndist = ( m_ExCalculated[0]*HyCalculated[1] - m_ExCalculated[1]*HyCalculated[0] ) / det;
|
|
const std::complex<double> ZxyUndist = ( m_ExCalculated[1]*HxCalculated[0] - m_ExCalculated[0]*HxCalculated[1] ) / det;
|
|
const std::complex<double> ZyxUndist = ( m_EyCalculated[0]*HyCalculated[1] - m_EyCalculated[1]*HyCalculated[0] ) / det;
|
|
const std::complex<double> ZyyUndist = ( m_EyCalculated[1]*HxCalculated[0] - m_EyCalculated[0]*HxCalculated[1] ) / det;
|
|
|
|
if( ( AnalysisControl::getInstance() )->getTypeOfDistortion() == AnalysisControl::ESTIMATE_DISTORTION_MATRIX_DIFFERENCE ){
|
|
const long long ID_Cxx = static_cast<long long>(m_arrayDistortionMatrixDifferences->IDsOfDistortionMatrixDifference[COMPONENT_ID_CXX]);
|
|
const long long ID_Cxy = static_cast<long long>(m_arrayDistortionMatrixDifferences->IDsOfDistortionMatrixDifference[COMPONENT_ID_CXY]);
|
|
const long long ID_Cyx = static_cast<long long>(m_arrayDistortionMatrixDifferences->IDsOfDistortionMatrixDifference[COMPONENT_ID_CYX]);
|
|
const long long ID_Cyy = static_cast<long long>(m_arrayDistortionMatrixDifferences->IDsOfDistortionMatrixDifference[COMPONENT_ID_CYY]);
|
|
if( ID_Cxx >= 0 ){
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxx[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_Cxx ] = ZxxUndist.real() / ZxxSD.realPart;// For Re(Zxx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxx[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_Cxx ] = ZxxUndist.imag() / ZxxSD.imagPart;// For Im(Zxx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxy[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_Cxx ] = ZxyUndist.real() / ZxySD.realPart;// For Re(Zxy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxy[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_Cxx ] = ZxyUndist.imag() / ZxySD.imagPart;// For Im(Zxy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyx[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_Cxx ] = 0.0;// For Re(Zyx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyx[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_Cxx ] = 0.0;// For Im(Zyx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyy[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_Cxx ] = 0.0;// For Re(Zyy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyy[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_Cxx ] = 0.0;// For Im(Zyy)
|
|
}
|
|
if( ID_Cxy >= 0 ){
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxx[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_Cxy ] = ZyxUndist.real() / ZxxSD.realPart;// For Re(Zxx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxx[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_Cxy ] = ZyxUndist.imag() / ZxxSD.imagPart;// For Im(Zxx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxy[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_Cxy ] = ZyyUndist.real() / ZxySD.realPart;// For Re(Zxy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxy[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_Cxy ] = ZyyUndist.imag() / ZxySD.imagPart;// For Im(Zxy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyx[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_Cxy ] = 0.0;// For Re(Zyx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyx[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_Cxy ] = 0.0;// For Im(Zyx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyy[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_Cxy ] = 0.0;// For Re(Zyy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyy[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_Cxy ] = 0.0;// For Im(Zyy)
|
|
}
|
|
if( ID_Cyx >= 0 ){
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxx[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_Cyx ] = 0.0;// For Re(Zxx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxx[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_Cyx ] = 0.0;// For Im(Zxx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxy[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_Cyx ] = 0.0;// For Re(Zxy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxy[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_Cyx ] = 0.0;// For Im(Zxy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyx[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_Cyx ] = ZxxUndist.real() / ZyxSD.realPart;// For Re(Zyx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyx[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_Cyx ] = ZxxUndist.imag() / ZyxSD.imagPart;// For Im(Zyx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyy[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_Cyx ] = ZxyUndist.real() / ZyySD.realPart;// For Re(Zyy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyy[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_Cyx ] = ZxyUndist.imag() / ZyySD.imagPart;// For Im(Zyy)
|
|
}
|
|
if( ID_Cyy >= 0 ){
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxx[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_Cyy ] = 0.0;// For Re(Zxx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxx[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_Cyy ] = 0.0;// For Im(Zxx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxy[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_Cyy ] = 0.0;// For Re(Zxy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxy[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_Cyy ] = 0.0;// For Im(Zxy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyx[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_Cyy ] = ZyxUndist.real() / ZyxSD.realPart;// For Re(Zyx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyx[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_Cyy ] = ZyxUndist.imag() / ZyxSD.imagPart;// For Im(Zyx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyy[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_Cyy ] = ZyyUndist.real() / ZyySD.realPart;// For Re(Zyy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyy[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_Cyy ] = ZyyUndist.imag() / ZyySD.imagPart;// For Im(Zyy)
|
|
}
|
|
}
|
|
else if( ( AnalysisControl::getInstance() )->getTypeOfDistortion() == AnalysisControl::ESTIMATE_GAINS_AND_ROTATIONS ){
|
|
|
|
const double ln10 = log(10.0);
|
|
|
|
const double gX = pow( 10.0, m_arrayGainsAndRotations->gainsAndRotations[ObservedDataStationMT::EX_GAIN] );
|
|
const double gY = pow( 10.0, m_arrayGainsAndRotations->gainsAndRotations[ObservedDataStationMT::EY_GAIN] );
|
|
const double betaX = m_arrayGainsAndRotations->gainsAndRotations[ObservedDataStationMT::EX_ROTATION];
|
|
const double betaY = m_arrayGainsAndRotations->gainsAndRotations[ObservedDataStationMT::EY_ROTATION];
|
|
|
|
const long long ID_GainX = static_cast<long long>(m_arrayGainsAndRotations->IDsOfGainsAndRotations[EX_GAIN]);
|
|
const long long ID_GainY = static_cast<long long>(m_arrayGainsAndRotations->IDsOfGainsAndRotations[EY_GAIN]);
|
|
const long long ID_RotX = static_cast<long long>(m_arrayGainsAndRotations->IDsOfGainsAndRotations[EX_ROTATION]);
|
|
const long long ID_RotY = static_cast<long long>(m_arrayGainsAndRotations->IDsOfGainsAndRotations[EY_ROTATION]);
|
|
|
|
if( ID_GainX >= 0 ){
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxx[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_GainX ] = gX * ln10 * cos( betaX ) * ZxxUndist.real() / ZxxSD.realPart;// For Re(Zxx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxx[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_GainX ] = gX * ln10 * cos( betaX ) * ZxxUndist.imag() / ZxxSD.imagPart;// For Im(Zxx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxy[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_GainX ] = gX * ln10 * cos( betaX ) * ZxyUndist.real() / ZxySD.realPart;// For Re(Zxy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxy[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_GainX ] = gX * ln10 * cos( betaX ) * ZxyUndist.imag() / ZxySD.imagPart;// For Im(Zxy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyx[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_GainX ] = gX * ln10 * sin( betaX ) * ZxxUndist.real() / ZyxSD.realPart;// For Re(Zyx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyx[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_GainX ] = gX * ln10 * sin( betaX ) * ZxxUndist.imag() / ZyxSD.imagPart;// For Im(Zyx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyy[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_GainX ] = gX * ln10 * sin( betaX ) * ZxyUndist.real() / ZyySD.realPart;// For Re(Zyy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyy[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_GainX ] = gX * ln10 * sin( betaX ) * ZxyUndist.imag() / ZyySD.imagPart;// For Im(Zyy)
|
|
}
|
|
if( ID_GainY >= 0 ){
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxx[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_GainY ] = - gY * ln10 * sin( betaY ) * ZyxUndist.real() / ZxxSD.realPart;// For Re(Zxx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxx[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_GainY ] = - gY * ln10 * sin( betaY ) * ZyxUndist.imag() / ZxxSD.imagPart;// For Im(Zxx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxy[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_GainY ] = - gY * ln10 * sin( betaY ) * ZyyUndist.real() / ZxySD.realPart;// For Re(Zxy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxy[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_GainY ] = - gY * ln10 * sin( betaY ) * ZyyUndist.imag() / ZxySD.imagPart;// For Im(Zxy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyx[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_GainY ] = gY * ln10 * cos( betaY ) * ZyxUndist.real() / ZyxSD.realPart;// For Re(Zyx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyx[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_GainY ] = gY * ln10 * cos( betaY ) * ZyxUndist.imag() / ZyxSD.imagPart;// For Im(Zyx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyy[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_GainY ] = gY * ln10 * cos( betaY ) * ZyyUndist.real() / ZyySD.realPart;// For Re(Zyy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyy[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_GainY ] = gY * ln10 * cos( betaY ) * ZyyUndist.imag() / ZyySD.imagPart;// For Im(Zyy)
|
|
}
|
|
if( ID_RotX >= 0 ){
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxx[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_RotX ] = - gX * sin( betaX ) * ZxxUndist.real() / ZxxSD.realPart;// For Re(Zxx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxx[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_RotX ] = - gX * sin( betaX ) * ZxxUndist.imag() / ZxxSD.imagPart;// For Im(Zxx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxy[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_RotX ] = - gX * sin( betaX ) * ZxyUndist.real() / ZxySD.realPart;// For Re(Zxy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxy[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_RotX ] = - gX * sin( betaX ) * ZxyUndist.imag() / ZxySD.imagPart;// For Im(Zxy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyx[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_RotX ] = gX * cos( betaX ) * ZxxUndist.real() / ZyxSD.realPart;// For Re(Zyx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyx[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_RotX ] = gX * cos( betaX ) * ZxxUndist.imag() / ZyxSD.imagPart;// For Im(Zyx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyy[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_RotX ] = gX * cos( betaX ) * ZxyUndist.real() / ZyySD.realPart;// For Re(Zyy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyy[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_RotX ] = gX * cos( betaX ) * ZxyUndist.imag() / ZyySD.imagPart;// For Im(Zyy)
|
|
}
|
|
if( ID_RotY >= 0 ){
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxx[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_RotY ] = - gY * cos( betaY ) * ZyxUndist.real() / ZxxSD.realPart;// For Re(Zxx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxx[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_RotY ] = - gY * cos( betaY ) * ZyxUndist.imag() / ZxxSD.imagPart;// For Im(Zxx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxy[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_RotY ] = - gY * cos( betaY ) * ZyyUndist.real() / ZxySD.realPart;// For Re(Zxy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxy[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_RotY ] = - gY * cos( betaY ) * ZyyUndist.imag() / ZxySD.imagPart;// For Im(Zxy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyx[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_RotY ] = - gY * sin( betaY ) * ZyxUndist.real() / ZyxSD.realPart;// For Re(Zyx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyx[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_RotY ] = - gY * sin( betaY ) * ZyxUndist.imag() / ZyxSD.imagPart;// For Im(Zyx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyy[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_RotY ] = - gY * sin( betaY ) * ZyyUndist.real() / ZyySD.realPart;// For Re(Zyy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyy[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_RotY ] = - gY * sin( betaY ) * ZyyUndist.imag() / ZyySD.imagPart;// For Im(Zyy)
|
|
}
|
|
}
|
|
else if( ( AnalysisControl::getInstance() )->getTypeOfDistortion() == AnalysisControl::ESTIMATE_GAINS_ONLY ){
|
|
|
|
const double ln10 = log(10.0);
|
|
|
|
const double gX = pow( 10.0, m_arrayGainsAndRotations->gainsAndRotations[ObservedDataStationMT::EX_GAIN] );
|
|
const double gY = pow( 10.0, m_arrayGainsAndRotations->gainsAndRotations[ObservedDataStationMT::EY_GAIN] );
|
|
|
|
const long long ID_GainX = static_cast<long long>(m_arrayGainsAndRotations->IDsOfGainsAndRotations[EX_GAIN]);
|
|
const long long ID_GainY = static_cast<long long>(m_arrayGainsAndRotations->IDsOfGainsAndRotations[EY_GAIN]);
|
|
if( ID_GainX >= 0 ){
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxx[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_GainX ] = gX * ln10 * ZxxUndist.real() / ZxxSD.realPart;// For Re(Zxx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxx[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_GainX ] = gX * ln10 * ZxxUndist.imag() / ZxxSD.imagPart;// For Im(Zxx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxy[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_GainX ] = gX * ln10 * ZxyUndist.real() / ZxySD.realPart;// For Re(Zxy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxy[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_GainX ] = gX * ln10 * ZxyUndist.imag() / ZxySD.imagPart;// For Im(Zxy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyx[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_GainX ] = 0.0;// For Re(Zyx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyx[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_GainX ] = 0.0;// For Im(Zyx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyy[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_GainX ] = 0.0;// For Re(Zyy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyy[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_GainX ] = 0.0;// For Im(Zyy)
|
|
}
|
|
if( ID_GainY >= 0 ){
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxx[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_GainY ] = 0.0;// For Re(Zxx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxx[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_GainY ] = 0.0;// For Im(Zxx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxy[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_GainY ] = 0.0;// For Re(Zxy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZxy[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_GainY ] = 0.0;// For Im(Zxy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyx[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_GainY ] = gY * ln10 * ZyxUndist.real() / ZyxSD.realPart;// For Re(Zyx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyx[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_GainY ] = gY * ln10 * ZyxUndist.imag() / ZyxSD.imagPart;// For Im(Zyx)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyy[freqIDThisPEInSta].realPart) + nBlkNotFixed + ID_GainY ] = gY * ln10 * ZyyUndist.real() / ZyySD.realPart;// For Re(Zyy)
|
|
sensitivityMatrix[ static_cast<long long>(nModel) * static_cast<long long>(m_dataIDOfZyy[freqIDThisPEInSta].imagPart) + nBlkNotFixed + ID_GainY ] = gY * ln10 * ZyyUndist.imag() / ZyySD.imagPart;// For Im(Zyy)
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
// Calculate data vector of this PE
|
|
void ObservedDataStationMT::calculateResidualVectorOfDataThisPE( const double freq, const int offset, double* vector ) const{
|
|
|
|
const int freqIDThisPEInSta = getFreqIDsAmongThisPE( freq );
|
|
|
|
if( freqIDThisPEInSta < 0 ){// Specified frequency is not the ones calculated by this PE in this station
|
|
return;
|
|
}
|
|
|
|
vector[ offset + m_dataIDOfZxx[freqIDThisPEInSta].realPart ] = m_ZxxResidual[freqIDThisPEInSta].realPart;
|
|
vector[ offset + m_dataIDOfZxx[freqIDThisPEInSta].imagPart ] = m_ZxxResidual[freqIDThisPEInSta].imagPart;
|
|
vector[ offset + m_dataIDOfZxy[freqIDThisPEInSta].realPart ] = m_ZxyResidual[freqIDThisPEInSta].realPart;
|
|
vector[ offset + m_dataIDOfZxy[freqIDThisPEInSta].imagPart ] = m_ZxyResidual[freqIDThisPEInSta].imagPart;
|
|
vector[ offset + m_dataIDOfZyx[freqIDThisPEInSta].realPart ] = m_ZyxResidual[freqIDThisPEInSta].realPart;
|
|
vector[ offset + m_dataIDOfZyx[freqIDThisPEInSta].imagPart ] = m_ZyxResidual[freqIDThisPEInSta].imagPart;
|
|
vector[ offset + m_dataIDOfZyy[freqIDThisPEInSta].realPart ] = m_ZyyResidual[freqIDThisPEInSta].realPart;
|
|
vector[ offset + m_dataIDOfZyy[freqIDThisPEInSta].imagPart ] = m_ZyyResidual[freqIDThisPEInSta].imagPart;
|
|
}
|
|
|
|
// Calulate L2 norm of misfit
|
|
double ObservedDataStationMT::calculateErrorSumOfSquaresThisPE() const{
|
|
|
|
double misfit(0.0);
|
|
|
|
for( int ifreq = 0; ifreq < m_numOfFreqCalculatedByThisStaAndPE; ++ifreq ){
|
|
misfit += pow( m_ZxxResidual[ifreq].realPart , 2 );
|
|
misfit += pow( m_ZxxResidual[ifreq].imagPart , 2 );
|
|
misfit += pow( m_ZxyResidual[ifreq].realPart , 2 );
|
|
misfit += pow( m_ZxyResidual[ifreq].imagPart , 2 );
|
|
misfit += pow( m_ZyxResidual[ifreq].realPart , 2 );
|
|
misfit += pow( m_ZyxResidual[ifreq].imagPart , 2 );
|
|
misfit += pow( m_ZyyResidual[ifreq].realPart , 2 );
|
|
misfit += pow( m_ZyyResidual[ifreq].imagPart , 2 );
|
|
}
|
|
|
|
return misfit;
|
|
|
|
}
|
|
|
|
// Copy current distortion parameters to previous ones
|
|
void ObservedDataStationMT::copyDistortionParamsCurToPre( const int iComp ){
|
|
|
|
const int type = ( AnalysisControl::getInstance() )->getTypeOfDistortion();
|
|
if( type == AnalysisControl::ESTIMATE_DISTORTION_MATRIX_DIFFERENCE ){
|
|
assert( m_arrayDistortionMatrixDifferences != NULL );
|
|
assert( iComp >= COMPONENT_ID_CXX );
|
|
assert( iComp <= COMPONENT_ID_CYY );
|
|
m_arrayDistortionMatrixDifferences->distortionMatrixDifferencePre[iComp] = m_arrayDistortionMatrixDifferences->distortionMatrixDifference[iComp];
|
|
}
|
|
else if( type == AnalysisControl::ESTIMATE_GAINS_AND_ROTATIONS ){
|
|
assert( m_arrayGainsAndRotations != NULL );
|
|
assert( iComp >= EX_GAIN );
|
|
assert( iComp <= EY_ROTATION );
|
|
m_arrayGainsAndRotations->gainsAndRotationsPre[iComp] = m_arrayGainsAndRotations->gainsAndRotations[iComp];
|
|
}
|
|
else if( type == AnalysisControl::ESTIMATE_GAINS_ONLY ){
|
|
m_arrayGainsAndRotations->gainsAndRotationsPre[iComp] = m_arrayGainsAndRotations->gainsAndRotations[iComp];
|
|
}else{
|
|
OutputFiles::m_logFile << "Error : Type of galvanic distortion is wrong. : " << type << std::endl;
|
|
exit(1);
|
|
}
|
|
|
|
}
|
|
|
|
/// Get flag specifing whether distortion matrix are fixed or not
|
|
bool ObservedDataStationMT::doesFixDistortionMatrix() const{
|
|
|
|
return m_fixDistortionMatrix;
|
|
|
|
}
|
|
|
|
// Get type of the electric field used to calculate response functions
|
|
int ObservedDataStationMT::getTypeOfElectricField() const{
|
|
|
|
return m_typeOfElectricField;
|
|
|
|
}
|
|
|
|
// Set flag specifing whether distortion matrix are fixed or not
|
|
void ObservedDataStationMT::setFixDistortionMatrix( const bool doesFix ){
|
|
|
|
m_fixDistortionMatrix = doesFix;
|
|
|
|
}
|
|
|
|
// Set type of the electric field used to calculate response functions
|
|
void ObservedDataStationMT::setTypeOfElectricField( const int type ){
|
|
|
|
m_typeOfElectricField = type;
|
|
|
|
}
|
|
|
|
// Set distortion parameters of previous iteration
|
|
void ObservedDataStationMT::setDistortionParamsPre( const int iComp, const double val ){
|
|
|
|
const int type = ( AnalysisControl::getInstance() )->getTypeOfDistortion();
|
|
if( type == AnalysisControl::ESTIMATE_DISTORTION_MATRIX_DIFFERENCE ){
|
|
assert( m_arrayDistortionMatrixDifferences != NULL );
|
|
assert( iComp >= COMPONENT_ID_CXX );
|
|
assert( iComp <= COMPONENT_ID_CYY );
|
|
m_arrayDistortionMatrixDifferences->distortionMatrixDifferencePre[iComp] = val;
|
|
}
|
|
else if( type == AnalysisControl::ESTIMATE_GAINS_AND_ROTATIONS ){
|
|
assert( m_arrayGainsAndRotations != NULL );
|
|
assert( iComp >= EX_GAIN );
|
|
assert( iComp <= EY_ROTATION );
|
|
m_arrayGainsAndRotations->gainsAndRotationsPre[iComp] = val;
|
|
}
|
|
else if( type == AnalysisControl::ESTIMATE_GAINS_ONLY ){
|
|
assert( m_arrayGainsAndRotations != NULL );
|
|
assert( iComp == EX_GAIN || iComp == EY_GAIN );
|
|
m_arrayGainsAndRotations->gainsAndRotationsPre[iComp] = val;
|
|
}
|
|
else{
|
|
OutputFiles::m_logFile << "Error : Type of galvanic distortion is wrong. : " << type << std::endl;
|
|
exit(1);
|
|
}
|
|
|
|
|
|
}
|
|
|
|
// Set distortion parameters
|
|
void ObservedDataStationMT::setDistortionParams( const int iComp, const double val ){
|
|
|
|
const int type = ( AnalysisControl::getInstance() )->getTypeOfDistortion();
|
|
if( type == AnalysisControl::ESTIMATE_DISTORTION_MATRIX_DIFFERENCE ){
|
|
assert( m_arrayDistortionMatrixDifferences != NULL );
|
|
assert( iComp >= COMPONENT_ID_CXX );
|
|
assert( iComp <= COMPONENT_ID_CYY );
|
|
m_arrayDistortionMatrixDifferences->distortionMatrixDifference[iComp] = val;
|
|
}
|
|
else if( type == AnalysisControl::ESTIMATE_GAINS_AND_ROTATIONS ){
|
|
assert( m_arrayGainsAndRotations != NULL );
|
|
assert( iComp >= EX_GAIN );
|
|
assert( iComp <= EY_ROTATION );
|
|
m_arrayGainsAndRotations->gainsAndRotations[iComp] = val;
|
|
}
|
|
else if( type == AnalysisControl::ESTIMATE_GAINS_ONLY ){
|
|
assert( m_arrayGainsAndRotations != NULL );
|
|
assert( iComp == EX_GAIN || iComp == EY_GAIN );
|
|
m_arrayGainsAndRotations->gainsAndRotations[iComp] = val;
|
|
}
|
|
else{
|
|
OutputFiles::m_logFile << "Error : Type of galvanic distortion is wrong. : " << type << std::endl;
|
|
exit(1);
|
|
}
|
|
|
|
}
|
|
|
|
// Set ID of distortion parameters
|
|
void ObservedDataStationMT::setIDOfDistortionParams( const int iComp, const int ID ){
|
|
|
|
const int type = ( AnalysisControl::getInstance() )->getTypeOfDistortion();
|
|
if( type == AnalysisControl::ESTIMATE_DISTORTION_MATRIX_DIFFERENCE ){
|
|
assert( m_arrayDistortionMatrixDifferences != NULL );
|
|
assert( iComp >= COMPONENT_ID_CXX );
|
|
assert( iComp <= COMPONENT_ID_CYY );
|
|
if( doesFixDistortionMatrix() ){
|
|
m_arrayDistortionMatrixDifferences->IDsOfDistortionMatrixDifference[iComp] = -1;
|
|
OutputFiles::m_logFile << "Warning : ID of distortion matrix must be -1 if distortion matrix is fixed." << std::endl;
|
|
}else{// Distortion matrix is not fixed
|
|
m_arrayDistortionMatrixDifferences->IDsOfDistortionMatrixDifference[iComp] = ID;
|
|
}
|
|
|
|
}
|
|
else if( type == AnalysisControl::ESTIMATE_GAINS_AND_ROTATIONS ){
|
|
assert( m_arrayGainsAndRotations != NULL );
|
|
assert( iComp >= EX_GAIN );
|
|
assert( iComp <= EY_ROTATION );
|
|
if( doesFixDistortionMatrix() ){
|
|
m_arrayGainsAndRotations->IDsOfGainsAndRotations[iComp] = -1;
|
|
OutputFiles::m_logFile << "Warning : ID of distortion matrix must be -1 if distortion matrix is fixed." << std::endl;
|
|
}else{// Distortion matrix is not fixed
|
|
m_arrayGainsAndRotations->IDsOfGainsAndRotations[iComp] = ID;
|
|
}
|
|
}
|
|
else if( type == AnalysisControl::ESTIMATE_GAINS_ONLY ){
|
|
assert( m_arrayGainsAndRotations != NULL );
|
|
assert( iComp == EX_GAIN || iComp == EY_GAIN );
|
|
if( doesFixDistortionMatrix() ){
|
|
m_arrayGainsAndRotations->IDsOfGainsAndRotations[iComp] = -1;
|
|
OutputFiles::m_logFile << "Warning : ID of distortion matrix must be -1 if distortion matrix is fixed." << std::endl;
|
|
}else{// Distortion matrix is not fixed
|
|
m_arrayGainsAndRotations->IDsOfGainsAndRotations[iComp] = ID;
|
|
}
|
|
}
|
|
else{
|
|
OutputFiles::m_logFile << "Error : Type of galvanic distortion is wrong. : " << type << std::endl;
|
|
exit(1);
|
|
}
|
|
|
|
}
|
|
|
|
// Set full updated value of distortion parameters
|
|
void ObservedDataStationMT::setDistortionParamsUpdatedFull( const int iComp, const double val ){
|
|
|
|
const int type = ( AnalysisControl::getInstance() )->getTypeOfDistortion();
|
|
if( type == AnalysisControl::ESTIMATE_DISTORTION_MATRIX_DIFFERENCE ){
|
|
assert( m_arrayDistortionMatrixDifferences != NULL );
|
|
assert( iComp >= COMPONENT_ID_CXX );
|
|
assert( iComp <= COMPONENT_ID_CYY );
|
|
if( doesFixDistortionMatrix() ){// Not change distortion matrix
|
|
m_arrayDistortionMatrixDifferences->distortionMatrixDifferenceUpdatedFull[iComp] = m_arrayDistortionMatrixDifferences->distortionMatrixDifferencePre[iComp];
|
|
}else{
|
|
m_arrayDistortionMatrixDifferences->distortionMatrixDifferenceUpdatedFull[iComp] = val;
|
|
}
|
|
}
|
|
else if( type == AnalysisControl::ESTIMATE_GAINS_AND_ROTATIONS ){
|
|
assert( m_arrayGainsAndRotations != NULL );
|
|
assert( iComp >= EX_GAIN );
|
|
assert( iComp <= EY_ROTATION );
|
|
if( doesFixDistortionMatrix() ){// Not change distortion matrix
|
|
m_arrayGainsAndRotations->gainsAndRotationsUpdatedFull[iComp] = m_arrayGainsAndRotations->gainsAndRotationsPre[iComp];
|
|
}else{
|
|
if( iComp == ObservedDataStationMT::EX_ROTATION || iComp == ObservedDataStationMT::EY_ROTATION ){
|
|
// Rotation angle is bounded in from -pi/2 to pi/2 radians
|
|
if( val > 0.5 * CommonParameters::PI ){
|
|
m_arrayGainsAndRotations->gainsAndRotationsUpdatedFull[iComp] = 0.5 * CommonParameters::PI;
|
|
}
|
|
else if( val < - 0.5 * CommonParameters::PI ){
|
|
m_arrayGainsAndRotations->gainsAndRotationsUpdatedFull[iComp] = - 0.5 * CommonParameters::PI;
|
|
}
|
|
else{
|
|
m_arrayGainsAndRotations->gainsAndRotationsUpdatedFull[iComp] = val;
|
|
}
|
|
}
|
|
else{
|
|
m_arrayGainsAndRotations->gainsAndRotationsUpdatedFull[iComp] = val;
|
|
}
|
|
}
|
|
|
|
}
|
|
else if( type == AnalysisControl::ESTIMATE_GAINS_ONLY ){
|
|
assert( m_arrayGainsAndRotations != NULL );
|
|
assert( iComp == EX_GAIN || iComp == EY_GAIN );
|
|
if( doesFixDistortionMatrix() ){// Not change distortion matrix
|
|
m_arrayGainsAndRotations->gainsAndRotationsUpdatedFull[iComp] = m_arrayGainsAndRotations->gainsAndRotationsPre[iComp];
|
|
}else{
|
|
m_arrayGainsAndRotations->gainsAndRotationsUpdatedFull[iComp] = val;
|
|
}
|
|
}
|
|
else{
|
|
OutputFiles::m_logFile << "Error : Type of galvanic distortion is wrong. : " << type << std::endl;
|
|
exit(1);
|
|
}
|
|
|
|
}
|
|
|
|
// Update distortion parameters
|
|
void ObservedDataStationMT::updateDistortionParams( const double dampingFactor ){
|
|
|
|
if( doesFixDistortionMatrix() ){// distortion matrix if fixed
|
|
return;
|
|
}
|
|
|
|
const int type = ( AnalysisControl::getInstance() )->getTypeOfDistortion();
|
|
|
|
if( type == AnalysisControl::ESTIMATE_DISTORTION_MATRIX_DIFFERENCE ){
|
|
assert( m_arrayDistortionMatrixDifferences != NULL );
|
|
|
|
for( int iComp = 0; iComp < 4; ++iComp ){
|
|
if( m_arrayDistortionMatrixDifferences->IDsOfDistortionMatrixDifference[iComp] >= 0 ){
|
|
m_arrayDistortionMatrixDifferences->distortionMatrixDifference[iComp] =
|
|
dampingFactor * m_arrayDistortionMatrixDifferences->distortionMatrixDifferenceUpdatedFull[iComp] + ( 1.0 - dampingFactor ) * m_arrayDistortionMatrixDifferences->distortionMatrixDifferencePre[iComp];
|
|
}
|
|
}
|
|
|
|
}
|
|
else if( type == AnalysisControl::ESTIMATE_GAINS_AND_ROTATIONS ){
|
|
assert( m_arrayGainsAndRotations != NULL );
|
|
|
|
for( int iComp = 0; iComp < 4; ++iComp ){
|
|
if( m_arrayGainsAndRotations->IDsOfGainsAndRotations[iComp] >= 0 ){
|
|
if( iComp == ObservedDataStationMT::EX_ROTATION || iComp == ObservedDataStationMT::EY_ROTATION ){
|
|
// Rotation angle is bounded in from -pi/2 to pi/2 radians
|
|
// y = tan(beta)
|
|
// beta = Arctan(y)
|
|
// dy/dbeta = 1 / cos(beta)^2
|
|
const double yUpdatedFull = tan(m_arrayGainsAndRotations->gainsAndRotationsUpdatedFull[iComp]);
|
|
const double yPre = tan(m_arrayGainsAndRotations->gainsAndRotationsPre[iComp]);
|
|
const double yNew = dampingFactor * yUpdatedFull + ( 1.0 - dampingFactor ) * yPre;
|
|
const double betaNew = atan(yNew);
|
|
m_arrayGainsAndRotations->gainsAndRotations[iComp] = betaNew;
|
|
}else{
|
|
m_arrayGainsAndRotations->gainsAndRotations[iComp] =
|
|
dampingFactor * m_arrayGainsAndRotations->gainsAndRotationsUpdatedFull[iComp] + ( 1.0 - dampingFactor ) * m_arrayGainsAndRotations->gainsAndRotationsPre[iComp];
|
|
}
|
|
}
|
|
}
|
|
|
|
}
|
|
else if( type == AnalysisControl::ESTIMATE_GAINS_ONLY ){
|
|
assert( m_arrayGainsAndRotations != NULL );
|
|
if( m_arrayGainsAndRotations->IDsOfGainsAndRotations[EX_GAIN] >= 0 ){
|
|
m_arrayGainsAndRotations->gainsAndRotations[EX_GAIN] =
|
|
dampingFactor * m_arrayGainsAndRotations->gainsAndRotationsUpdatedFull[EX_GAIN] + ( 1.0 - dampingFactor ) * m_arrayGainsAndRotations->gainsAndRotationsPre[EX_GAIN];
|
|
}
|
|
if( m_arrayGainsAndRotations->IDsOfGainsAndRotations[EY_GAIN] >= 0 ){
|
|
m_arrayGainsAndRotations->gainsAndRotations[EY_GAIN] =
|
|
dampingFactor * m_arrayGainsAndRotations->gainsAndRotationsUpdatedFull[EY_GAIN] + ( 1.0 - dampingFactor ) * m_arrayGainsAndRotations->gainsAndRotationsPre[EY_GAIN];
|
|
}
|
|
}
|
|
else{
|
|
OutputFiles::m_logFile << "Error : Type of galvanic distortion is wrong. : " << type << std::endl;
|
|
exit(1);
|
|
}
|
|
|
|
}
|
|
|
|
// Get distortion parameters of previous iteration
|
|
double ObservedDataStationMT::getDistortionParamsPre( const int iComp ) const{
|
|
|
|
const int type = ( AnalysisControl::getInstance() )->getTypeOfDistortion();
|
|
if( type == AnalysisControl::ESTIMATE_DISTORTION_MATRIX_DIFFERENCE ){
|
|
assert( m_arrayDistortionMatrixDifferences != NULL );
|
|
assert( iComp >= COMPONENT_ID_CXX );
|
|
assert( iComp <= COMPONENT_ID_CYY );
|
|
|
|
return m_arrayDistortionMatrixDifferences->distortionMatrixDifferencePre[iComp];
|
|
}
|
|
else if( type == AnalysisControl::ESTIMATE_GAINS_AND_ROTATIONS ){
|
|
assert( m_arrayGainsAndRotations != NULL );
|
|
assert( iComp >= EX_GAIN );
|
|
assert( iComp <= EY_ROTATION );
|
|
|
|
return m_arrayGainsAndRotations->gainsAndRotationsPre[iComp];
|
|
}
|
|
else if( type == AnalysisControl::ESTIMATE_GAINS_ONLY ){
|
|
assert( m_arrayGainsAndRotations != NULL );
|
|
assert( iComp == EX_GAIN || iComp == EY_GAIN );
|
|
|
|
return m_arrayGainsAndRotations->gainsAndRotationsPre[iComp];
|
|
}
|
|
else{
|
|
OutputFiles::m_logFile << "Error : Type of galvanic distortion is wrong. : " << type << std::endl;
|
|
exit(1);
|
|
}
|
|
|
|
}
|
|
|
|
// Get distortion parameters
|
|
double ObservedDataStationMT::getDistortionParams( const int iComp ) const{
|
|
|
|
const int type = ( AnalysisControl::getInstance() )->getTypeOfDistortion();
|
|
if( type == AnalysisControl::ESTIMATE_DISTORTION_MATRIX_DIFFERENCE ){
|
|
assert( m_arrayDistortionMatrixDifferences != NULL );
|
|
assert( iComp >= COMPONENT_ID_CXX );
|
|
assert( iComp <= COMPONENT_ID_CYY );
|
|
|
|
return m_arrayDistortionMatrixDifferences->distortionMatrixDifference[iComp];
|
|
}
|
|
else if( type == AnalysisControl::ESTIMATE_GAINS_AND_ROTATIONS ){
|
|
assert( m_arrayGainsAndRotations != NULL );
|
|
assert( iComp >= EX_GAIN );
|
|
assert( iComp <= EY_ROTATION );
|
|
|
|
return m_arrayGainsAndRotations->gainsAndRotations[iComp];
|
|
}
|
|
else if( type == AnalysisControl::ESTIMATE_GAINS_ONLY ){
|
|
assert( m_arrayGainsAndRotations != NULL );
|
|
assert( iComp == EX_GAIN || iComp == EY_GAIN );
|
|
|
|
return m_arrayGainsAndRotations->gainsAndRotations[iComp];
|
|
}
|
|
else{
|
|
OutputFiles::m_logFile << "Error : Type of galvanic distortion is wrong. : " << type << std::endl;
|
|
exit(1);
|
|
}
|
|
|
|
}
|
|
|
|
// Get ID of distortion parameters
|
|
int ObservedDataStationMT::getIDOfDistortionParams( const int iComp ) const{
|
|
|
|
const int type = ( AnalysisControl::getInstance() )->getTypeOfDistortion();
|
|
if( type == AnalysisControl::ESTIMATE_DISTORTION_MATRIX_DIFFERENCE ){
|
|
assert( m_arrayDistortionMatrixDifferences != NULL );
|
|
assert( iComp >= COMPONENT_ID_CXX );
|
|
assert( iComp <= COMPONENT_ID_CYY );
|
|
|
|
return m_arrayDistortionMatrixDifferences->IDsOfDistortionMatrixDifference[iComp];
|
|
}
|
|
else if( type == AnalysisControl::ESTIMATE_GAINS_AND_ROTATIONS ){
|
|
assert( m_arrayGainsAndRotations != NULL );
|
|
assert( iComp >= EX_GAIN );
|
|
assert( iComp <= EY_ROTATION );
|
|
|
|
return m_arrayGainsAndRotations->IDsOfGainsAndRotations[iComp];
|
|
}
|
|
else if( type == AnalysisControl::ESTIMATE_GAINS_ONLY ){
|
|
assert( m_arrayGainsAndRotations != NULL );
|
|
|
|
return m_arrayGainsAndRotations->IDsOfGainsAndRotations[iComp];
|
|
}
|
|
else{
|
|
OutputFiles::m_logFile << "Error : Type of galvanic distortion is wrong. : " << type << std::endl;
|
|
exit(1);
|
|
}
|
|
|
|
}
|
|
|
|
// Get full updated value of distortion parameters
|
|
double ObservedDataStationMT::getDistortionParamsUpdatedFull( const int iComp ) const{
|
|
|
|
const int type = ( AnalysisControl::getInstance() )->getTypeOfDistortion();
|
|
if( type == AnalysisControl::ESTIMATE_DISTORTION_MATRIX_DIFFERENCE ){
|
|
assert( m_arrayDistortionMatrixDifferences != NULL );
|
|
assert( iComp >= COMPONENT_ID_CXX );
|
|
assert( iComp <= COMPONENT_ID_CYY );
|
|
assert( iComp >= COMPONENT_ID_CXX );
|
|
assert( iComp <= COMPONENT_ID_CYY );
|
|
|
|
return m_arrayDistortionMatrixDifferences->distortionMatrixDifferenceUpdatedFull[iComp];
|
|
}
|
|
else if( type == AnalysisControl::ESTIMATE_GAINS_AND_ROTATIONS ){
|
|
assert( m_arrayGainsAndRotations != NULL );
|
|
assert( iComp >= EX_GAIN );
|
|
assert( iComp <= EY_ROTATION );
|
|
|
|
return m_arrayGainsAndRotations->gainsAndRotationsUpdatedFull[iComp];
|
|
}
|
|
else if( type == AnalysisControl::ESTIMATE_GAINS_ONLY ){
|
|
assert( m_arrayGainsAndRotations != NULL );
|
|
assert( iComp == EX_GAIN || iComp == EY_GAIN );
|
|
|
|
return m_arrayGainsAndRotations->gainsAndRotationsUpdatedFull[iComp];
|
|
}
|
|
else{
|
|
OutputFiles::m_logFile << "Error : Type of galvanic distortion is wrong. : " << type << std::endl;
|
|
exit(1);
|
|
}
|
|
|
|
}
|