/* * To change this license header, choose License Headers in Project Properties. * To change this template file, choose Tools | Templates * and open the template in the editor. */ /* * File: Shue.hh * Author: hacene * * Created on December 31, 2020, 10:14 AM */ #ifndef SHUE_HH #define SHUE_HH #include "Parameter.hh" #include "ParamData.hh" #include "DataTypeMath.hh" #include "Operation.hh" #include "shue_compute.hh" #include namespace AMDA { namespace Parameters { namespace Shue { #define AVERAGE_TIME 1200 // (seconds) #define MAX_GAP_SIZE 3600 // (seconds) #define DEFAULT_IMF_GSM_X 0. #define DEFAULT_IMF_GSM_Y 2. #define DEFAULT_IMF_GSM_Z -3. #define DEFAULT_PSW 3. class ShueBase : public Operation { public: ShueBase(Process& pProcess, ParamData& paramImfInput, ParamData& paramPswInput) : Operation(pProcess), _paramImfInput(dynamic_cast >&> (paramImfInput)), _paramPswInput(dynamic_cast&> (paramPswInput)) { } virtual ~ShueBase() { } void pushImfData(ParamDataIndexInfo &pParamDataIndexInfo) { for (unsigned int _index = pParamDataIndexInfo._startIndex; _index < pParamDataIndexInfo._startIndex + pParamDataIndexInfo._nbDataToProcess; ++_index) { double time = _paramImfInput.getTime(_index); std::vector inputElt = _paramImfInput.get(_index); _b_x_gse.push_back(std::pair(time, inputElt[0])); _b_y_gse.push_back(std::pair(time, inputElt[1])); _b_z_gse.push_back(std::pair(time, inputElt[2])); } } void pushPswData(ParamDataIndexInfo &pParamDataIndexInfo) { for (unsigned int _index = pParamDataIndexInfo._startIndex; _index < pParamDataIndexInfo._startIndex + pParamDataIndexInfo._nbDataToProcess; ++_index) { double time = _paramPswInput.getTime(_index); float inputElt = _paramPswInput.get(_index); _psw.push_back(std::pair(time, inputElt)); } } float getValue(std::vector >& input, double time) { double min_t = time - AVERAGE_TIME / 2.; double max_t = time + AVERAGE_TIME / 2.; std::vector > values_for_mean; std::pair prev_value(NAN, NAN); std::pair next_value(NAN, NAN); for (std::vector >::iterator it = input.begin(); it != input.end(); ++it) { if (isNAN(it->second)) continue; else if (it->first > max_t) { next_value = *it; break; } else if (it->first < min_t) { prev_value = *it; } else { values_for_mean.push_back(*it); } } float value = NAN; if (!values_for_mean.empty()) { //Compute mean float sum = 0; for (std::vector >::iterator it = values_for_mean.begin(); it != values_for_mean.end(); ++it) { sum += it->second; } value = sum / (float) values_for_mean.size(); } else { if (!isNAN(prev_value.first) && !isNAN(next_value.first) && (next_value.first - prev_value.first <= MAX_GAP_SIZE)) { //Compute interpolated value value = prev_value.second + (time - prev_value.first) / (next_value.first - prev_value.first) * (next_value.second - prev_value.second); } } return value; } void getImfData(double time, float& b_x, float& b_y, float& b_z) { b_x = getValue(_b_x_gse, time); b_y = getValue(_b_y_gse, time); b_z = getValue(_b_z_gse, time); } void getPswData(double time, float& p_sw) { p_sw = getValue(_psw, time); } private: ParamDataSpec >& _paramImfInput; ParamDataSpec& _paramPswInput; std::vector > _b_x_gse; std::vector > _b_y_gse; std::vector > _b_z_gse; std::vector > _psw; }; template class Shue : public ShueBase { public: Shue(Process& pProcess, ParamDataSpec >& paramInput, ParamData& paramImfInput, ParamData& paramPswInput, bool inGSE, bool model98) : ShueBase(pProcess, paramImfInput, paramPswInput), _paramInput(paramInput), _paramOutput(new TOutputParamData), _inGSE(inGSE), _model98(model98) { _paramDataOutput = _paramOutput; } void write(ParamDataIndexInfo &pParamDataIndexInfo) { for (unsigned int _index = pParamDataIndexInfo._startIndex; _index < pParamDataIndexInfo._startIndex + pParamDataIndexInfo._nbDataToProcess; ++_index) { double crtTime = _paramInput.getTime(_index); float b_x_gse, b_y_gse, b_z_gse; getImfData(crtTime, b_x_gse, b_y_gse, b_z_gse); float p_sw; getPswData(crtTime, p_sw); if (isNAN(p_sw)) { p_sw = DEFAULT_PSW; } std::vector inputElt = _paramInput.get(_index); time_t timestamp = crtTime; struct tm *tmp; tmp = gmtime(×tamp); std::vector ouputElt; ouputElt.resize(3); ouputElt << NotANumber(); //Init geopack with GSM frame ShueCompute::initInGSM(1900 + tmp->tm_year, 1 + tmp->tm_yday, tmp->tm_hour, tmp->tm_min, tmp->tm_sec); //Compute shue 98 // inout intialisation int transform_flag = -1; float sat_pos_X_GSE = inputElt[0]; float sat_pos_Y_GSE = inputElt[1]; float sat_pos_Z_GSE = inputElt[2]; float sat_pos_X_GSM = 0.; float sat_pos_Y_GSM = 0.; float sat_pos_Z_GSM = 0.; gswgse_08_(&sat_pos_X_GSM, &sat_pos_Y_GSM, &sat_pos_Z_GSM, &sat_pos_X_GSE, &sat_pos_Y_GSE, &sat_pos_Z_GSE, &transform_flag); //Compute Imf B field in GSM frame float b_x_gsm, b_y_gsm, b_z_gsm; if (!isNAN(b_x_gse) && !isNAN(b_y_gse) && !isNAN(b_z_gse)) { gswgse_08_(&b_x_gsm, &b_y_gsm, &b_z_gsm, &b_x_gse, &b_y_gse, &b_z_gse, &transform_flag); } else { b_x_gsm = DEFAULT_IMF_GSM_X; b_y_gsm = DEFAULT_IMF_GSM_Y; b_z_gsm = DEFAULT_IMF_GSM_Z; } float X_GSW, Y_GSW, Z_GSW; int flag_transToGSE = 1; float X_GSE = 0; float Y_GSE = 0; float Z_GSE = 0; float dst; int pos_id; bool computed = false; /** * for testing to reproduce JOURNAL OF GEOPHYSICAL RESEARCH, VOL. 103, NO. A8, PAGES 17,691-17,700, AUGUST 1, 1998 * fig 4 p_sw=14.0; b_z_gsm = 9.3; * */ if(_model98){ computed = ShueCompute::shue98(p_sw, b_z_gsm, sat_pos_X_GSM, sat_pos_Y_GSM, sat_pos_Z_GSM, X_GSW, Y_GSW, Z_GSW, dst, pos_id); }else{ computed = ShueCompute::shue97(p_sw, b_z_gsm, sat_pos_X_GSM, sat_pos_Y_GSM, sat_pos_Z_GSM, X_GSW, Y_GSW, Z_GSW, dst, pos_id); } if (computed) { dst = (float) dst*(-pos_id); if(_inGSE){ // transform results from GSM to GSE gswgse_08_(&X_GSW, &Y_GSW, &Z_GSW, &X_GSE, &Y_GSE, &Z_GSE, &flag_transToGSE); ouputElt[0] = X_GSE; ouputElt[1] = Y_GSE; ouputElt[2] = Z_GSE; }else{ ouputElt[0] = X_GSW; ouputElt[1] = Y_GSW; ouputElt[2] = Z_GSW; } } _paramOutput->pushTime(crtTime); pushData(ouputElt, dst, pos_id); } }; virtual void pushData(std::vector ouputElt, ElemType dst, int pos_id) = 0; protected: ParamDataSpec >& _paramInput; TOutputParamData* _paramOutput; bool _inGSE; bool _model98; }; template class ShuePos : public Shue>> { public: ShuePos(Process& pProcess, ParamDataSpec >& paramInput, ParamData& paramImfInput, ParamData& paramPswInput, bool inGSE, bool model98): Shue>>(pProcess, paramInput,paramImfInput,paramPswInput, inGSE, model98){ } virtual ~ShuePos(){ } void pushData(std::vector ouputElt, ElemType /*dst*/, int /*pos_id*/) { Shue>>::_paramOutput->getDataList().push_back(ouputElt); } }; template class ShueDst : public Shue> { public: ShueDst(Process& pProcess, ParamDataSpec >& paramInput, ParamData& paramImfInput, ParamData& paramPswInput, bool inGSE, bool model98): Shue>(pProcess, paramInput,paramImfInput,paramPswInput, inGSE, model98){ } virtual ~ShueDst(){ } void pushData(std::vector /*ouputElt*/, ElemType dst, int /*pos_id*/) { Shue>::_paramOutput->getDataList().push_back(dst); } }; template class ShuePos_flag : public Shue> { public: ShuePos_flag(Process& pProcess, ParamDataSpec >& paramInput, ParamData& paramImfInput, ParamData& paramPswInput, bool inGSE, bool model98): Shue>(pProcess, paramInput,paramImfInput,paramPswInput, inGSE, model98){ } virtual ~ShuePos_flag(){ } void pushData(std::vector /*ouputElt*/, ElemType /*dst*/, int pos_id) { Shue>::_paramOutput->getDataList().push_back(pos_id); } }; } // end Shue } // end Parameters } // end AMDA #endif /* SHUE_HH */