zhangyiss/femtic
1
0
Fork 0
mirror of https://github.com/yoshiya-usui/femtic.git synced 2025-05-13 04:12:41 +08:00
Code Issues Actions Packages Projects Releases Wiki Activity
v4.1
femtic/src/ObservedDataStationMT.cpp

1453 lines
76 KiB
C++
Raw Permalink Normal View History

add new file
2021-11-09 00:06:52 +08:00
//-------------------------------------------------------------------------------------------------------
// 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);
}
}
Reference in New Issue Copy Permalink