From 7050b34f7158393f77f4e0d080ee4cb78bc37e88 Mon Sep 17 00:00:00 2001 From: Alan Salcedo Gomez Date: Fri, 6 Mar 2026 10:20:46 -0600 Subject: [PATCH] Reimplement xpol --- Detector.cc | 191 +++++++++++++--------- Detector.h | 31 +++- RayTrace.h | 4 +- Report.cc | 451 ++++++++++++++++++++++++++++++++++++++-------------- Report.h | 39 ++++- Settings.cc | 129 ++++++++------- Settings.h | 3 +- 7 files changed, 579 insertions(+), 269 deletions(-) diff --git a/Detector.cc b/Detector.cc index 110f44db..1f821b44 100644 --- a/Detector.cc +++ b/Detector.cc @@ -2254,10 +2254,15 @@ inline void Detector::ReadAllAntennaGains(Settings *settings1){ std::string VgainFile; std::string VgainTopFile; std::string HgainFile; - std::string TxgainFile; + std::string TxgainFile; + std::string TxgainFileCross; + std::string VgainCrossFile; + std::string VgainTopCrossFile; + std::string HgainCrossFile; //Adding step to read Tx gain. Will hardcode to PVA gain for now. - TxgainFile = string(getenv("ARA_SIM_DIR"))+"/data/antennas/realizedGain/PVA_RealizedGainAndPhase_Copol_Kansas2024.txt"; + TxgainFile = string(getenv("ARA_SIM_DIR"))+"/data/antennas/realizedGain/PVA_RealizedGainAndPhase_Copol_Kansas2024.txt"; + TxgainFileCross = string(getenv("ARA_SIM_DIR"))+"/data/antennas/realizedGain/PVA_RealizedGainAndPhase_Crosspol_Kansas2024.txt"; if (settings1->ANTENNA_MODE == 0){ // use the orignal Vpol/Hpol gains @@ -2307,7 +2312,12 @@ inline void Detector::ReadAllAntennaGains(Settings *settings1){ VgainTopFile = string(getenv("ARA_SIM_DIR"))+"/data/antennas/realizedGain/ARA_TVpol_RealizedGainAndPhase_Copol_Custom.txt"; HgainFile = string(getenv("ARA_SIM_DIR"))+"/data/antennas/realizedGain/ARA_Hpol_RealizedGainAndPhase_Copol_Custom.txt"; } - + + // Add cross-pol gain files + VgainCrossFile = string(getenv("ARA_SIM_DIR"))+"/data/antennas/realizedGain/ARA_BVpol_RealizedGainAndPhase_Crosspol_Kansas2024.txt"; + VgainTopCrossFile = string(getenv("ARA_SIM_DIR"))+"/data/antennas/realizedGain/ARA_BVpol_RealizedGainAndPhase_Crosspol_Kansas2024.txt"; //Should be TV but it's not working due to a negative SWR (first SWR entry) -ASG 12/09/24 + HgainCrossFile = string(getenv("ARA_SIM_DIR"))+"/data/antennas/realizedGain/ARA_Hpol_RealizedGainAndPhase_Crosspol_Kansas2024.txt"; + // Check for ALL_ANT_V_ON, then set all antennas to VPol if true if (settings1->ALL_ANT_V_ON == 1) { HgainFile = VgainFile; @@ -2320,12 +2330,18 @@ inline void Detector::ReadAllAntennaGains(Settings *settings1){ freq_width = -1; freq_init = -1; - //Read in antenna gain files. + //Read co-pol Rx gain files. ReadAntennaGain(VgainFile, settings1, eVPol); ReadAntennaGain(VgainTopFile, settings1, eVPolTop); ReadAntennaGain(HgainFile, settings1, eHPol); ReadAntennaGain(TxgainFile, settings1, eTx); + // Read cross-pol Rx gain files + ReadAntennaGain(VgainCrossFile, settings1, eVPolCross); + ReadAntennaGain(VgainTopCrossFile, settings1, eVPolTopCross); + ReadAntennaGain(HgainCrossFile, settings1, eHPolCross); + ReadAntennaGain(TxgainFileCross, settings1, eTxCross); + // update parameters to reflect what was read-in params.freq_step = freq_step; params.ang_step = ang_step; @@ -2378,34 +2394,57 @@ inline void Detector::ReadAntennaGain(string filename, Settings *settings1, EAnt // make sure dummy variables point to the right variables switch(type) { - case(eVPol) : + case(eVPol): freq = &Freq; gain = &Vgain; phase = &Vphase; transAnt_databin = &transV_databin; source_n = &antenna_source_medium_n; break; - case(eVPolTop) : + case(eVPolTop): freq = &Freq; gain = &VgainTop; phase = &VphaseTop; transAnt_databin = &transVTop_databin; source_n = &antenna_source_medium_n; break; - case(eHPol) : + case(eHPol): freq = &Freq; gain = &Hgain; phase = &Hphase; transAnt_databin = &transH_databin; source_n = &antenna_source_medium_n; break; - case(eTx) : + case(eVPolCross): // Cross-pol VPol + freq = &FreqCross; + gain = &VgainCross; + phase = &VphaseCross; + transAnt_databin = &transVCross_databin; + break; + case(eVPolTopCross): // Cross-pol VPol Top + freq = &FreqCross; + gain = &VgainTopCross; + phase = &VphaseTopCross; + transAnt_databin = &transVTopCross_databin; + break; + case(eHPolCross): // Cross-pol HPol + freq = &FreqCross; + gain = &HgainCross; + phase = &HphaseCross; + transAnt_databin = &transHCross_databin; + break; + case(eTx): freq = &TxFreq; gain = &Txgain; phase = &Txphase; source_n = &Txantenna_source_medium_n; break; - default : + case(eTxCross): + freq = &TxFreqCross; + gain = &TxgainCross; + phase = &TxphaseCross; + break; + default: throw runtime_error("Unknown antenna type!"); } @@ -2454,9 +2493,10 @@ inline void Detector::ReadAntennaGain(string filename, Settings *settings1, EAnt // clear vector in case there's any lingering data gain->clear(); - phase->clear(); - if(freq_step == -1 || type == eTx) { // only reset if it hasn't been read-in yet - freq->clear(); + phase->clear(); + + if(freq_step == -1 || type == eTx || type == eVPolCross || type == eVPolTopCross || type == eHPolCross || eTxCross) { // only reset if it hasn't been read-in yet + freq->clear(); } // check the file opened successfully @@ -2488,7 +2528,7 @@ inline void Detector::ReadAntennaGain(string filename, Settings *settings1, EAnt double buff_n = stof(words[2]); // always assign if this is the first read-in or its the transmitter - if(freq_step == -1 || type == eTx) { + if(freq_step == -1 || type == eTx || type == eVPolCross || type == eVPolTopCross || type == eHPolCross || eTxCross ) { *source_n = buff_n; } else { // if not the first read-in, check for consistency @@ -2616,30 +2656,34 @@ inline void Detector::ReadAntennaGain(string filename, Settings *settings1, EAnt return; } - // set parameter values if this is the first read-in - if(freq_step == -1) { - freq_step = (int)freq->size(); + int this_freq_step = (int)freq->size(); + // set parameter values for the smallest freq_step + if(freq_step == -1 || this_freq_step < freq_step) { + freq_step = this_freq_step; ang_step = (int)gain->back().size(); freq_width = freq->at(1)-freq->at(0); - freq_init = freq->at(0); + freq_init = freq->at(0); } - // check things look sensible - if(Transm.size() != freq_step) { - throw runtime_error("Transm has an unexpected length! "+filename); - } - if(gain->size() != freq_step) { - throw runtime_error("gain has an unexpected length! "+filename); - } - if(phase->size() != freq_step) { - throw runtime_error("phase has an unexpected length! "+filename); - } - for(int i = 0; i < freq_step; ++i) { - if(gain->at(i).size() != ang_step) { - throw runtime_error("gain vectors have inconsistent length! "+filename); + // Only enforce check for non-Tx and non-crosspol antennas + if(type != eTx && type != eTxCross && type != eVPolCross && type != eVPolTopCross && type != eHPolCross) { + // check things look sensible + if(Transm.size() != freq_step) { + throw runtime_error("Transm has an unexpected length! "+filename); + } + if(gain->size() != freq_step) { + throw runtime_error("gain has an unexpected length! "+filename); } - if(phase->at(i).size() != ang_step) { - throw runtime_error("gain vectors have inconsistent length! "+filename); + if(phase->size() != freq_step) { + throw runtime_error("phase has an unexpected length! "+filename); + } + for(int i = 0; i < freq_step; ++i) { + if(gain->at(i).size() != ang_step) { + throw runtime_error("gain vectors have inconsistent length! "+filename); + } + if(phase->at(i).size() != ang_step) { + throw runtime_error("gain vectors have inconsistent length! "+filename); + } } } @@ -3095,7 +3139,7 @@ double Detector::GetAntPhase( double freq, double theta, double phi, int ant_m, } -double Detector::GetGain_1D_OutZero( double freq, double theta, double phi, int ant_m, double antenna_target_medium_n, int string_number, int ant_number, bool useInTransmitterMode) { +double Detector::GetGain_1D_OutZero( double freq, double theta, double phi, int ant_m, double antenna_target_medium_n, int string_number, int ant_number, bool useInTransmitterMode, bool useCrossPol) { /* The purpose of this function is to interpolate the globally defined gain arrays (Vgain, VgainTop, Hgain, Txgain) @@ -3128,56 +3172,57 @@ double Detector::GetGain_1D_OutZero( double freq, double theta, double phi, int //VPol Rx if ( Detector_mode == 5 ){ // Phased Array mode if ( useInTransmitterMode ) { - tempGain = &Txgain; // Transmitter mode + tempGain = useCrossPol ? &TxgainCross : &Txgain; // Transmitter mode } else if ( ant_m == 1 ) { - tempGain = &Hgain; // PA Hpols + tempGain = useCrossPol ? &HgainCross : &Hgain; // HPol // PA Hpols } else { if ( string_number == 0 ) { - tempGain = &Vgain; // PA Vpols + tempGain = useCrossPol ? &VgainCross : &Vgain; // PA Vpols } else { if ( ant_number == 1 ) { - tempGain = &VgainTop; // A5 Top VPols + tempGain = useCrossPol ? &VgainTopCross : &VgainTop; // A5 Top VPols } else { - tempGain = &Vgain; // A5 Bottom Vpols + useCrossPol ? &VgainCross : &Vgain; // A5 Bottom Vpols } } } } else { // Traditional Station mode - //Tx - if (useInTransmitterMode) { - tempGain = &Txgain; - } - else if (ant_m == 0) { - if (ant_number == 0) { - tempGain = &Vgain; + if (!useInTransmitterMode) { + if (ant_m == 0) { // VPol + if (ant_number == 0) { + tempGain = useCrossPol ? &VgainCross : &Vgain; + } + else if (ant_number == 2) { + tempGain = useCrossPol ? &VgainTopCross : &VgainTop; + } } - else if (ant_number == 2) { - tempGain = &VgainTop; + else if (ant_m == 1) { // HPol + tempGain = useCrossPol ? &HgainCross : &Hgain; } } - //HPol Rx - else if (ant_m == 1) { - tempGain = &Hgain; - } - else { - throw runtime_error("In GetGain_1D_OutZero: No appropriate gain model for this simulation setup."); + else { // Tx mode + tempGain = useCrossPol ? &TxgainCross : &Txgain; } } + if (!tempGain) { + throw runtime_error("In GetGain_1D_OutZero: No appropriate gain model for this simulation setup."); + } + double thisFreq_init; double thisFreq_width; if(useInTransmitterMode) { - F = &TxFreq; + F = useCrossPol ? &TxFreqCross : &TxFreq; thisFreq_init = Tx_freq_init; thisFreq_width = Tx_freq_width; - } + } else { - F = &Freq; + F = useCrossPol ? &FreqCross : &Freq; thisFreq_init = freq_init; thisFreq_width = freq_width; } @@ -3281,7 +3326,7 @@ double Detector::GetImpedance( double freq, int ant_m, int ant_number, bool useI } -double Detector::GetAntPhase_1D( double freq, double theta, double phi, int ant_m, double antenna_target_medium_n, bool useInTransmitterMode ) { +double Detector::GetAntPhase_1D( double freq, double theta, double phi, int ant_m, double antenna_target_medium_n, bool useInTransmitterMode, bool useCrossPol ) { // check that target index of refraction is sensible if(antenna_target_medium_n < 1.0) { @@ -3302,34 +3347,34 @@ double Detector::GetAntPhase_1D( double freq, double theta, double phi, int ant_ //Creating tempPhase array to make this function more dynamic for Rx and Tx mode. vector > *tempPhase = nullptr; - vector * F; - - //Tx - if (useInTransmitterMode) { - tempPhase = &Txphase; + vector * F; + + // Assign tempPhase based on polarization type and cross-pol flag + if (!useInTransmitterMode) { + if (ant_m == 0) { // VPol + tempPhase = useCrossPol ? &VphaseCross : &Vphase; + } + else if (ant_m == 1) { // HPol + tempPhase = useCrossPol ? &HphaseCross : &Hphase; + } } - //VPol Rx - else if (ant_m == 0) { - tempPhase = &Vphase; + else { // Transmitter mode + tempPhase = useCrossPol ? &TxphaseCross : &Txphase; } - //HPol Rx - else if (ant_m == 1) { - tempPhase = &Hphase; + if (!tempPhase) { + throw runtime_error("In GetAntPhase_1D: No appropriate phase model for this simulation setup."); } - else { - throw runtime_error("In GetAntPhase_1D: No appropriate gain model for this simulation setup."); - } double thisFreq_init; double thisFreq_width; int thisFreq_step; if(useInTransmitterMode) { - F = &TxFreq; + F = useCrossPol ? &TxFreqCross : &TxFreq; thisFreq_init = Tx_freq_init; thisFreq_width = Tx_freq_width; } else { - F = &Freq; + F = useCrossPol ? &FreqCross : &Freq; thisFreq_init = freq_init; thisFreq_width = freq_width; } diff --git a/Detector.h b/Detector.h index 0b40a0a7..43d155e8 100644 --- a/Detector.h +++ b/Detector.h @@ -217,7 +217,11 @@ enum EAntennaType { eVPol, // (bottom) Vpol eVPolTop, // top Vpol eHPol, // Hpol - eTx // transmitter + eTx, // transmitter + eTxCross, // transmitter + eVPolCross, // Cross-pol VPol + eVPolTopCross, // Cross-pol VPol Top + eHPolCross // Cross-pol HPol }; class Detector { @@ -233,6 +237,15 @@ class Detector { vector > Hgain; vector > Hphase; vector Freq; + vector FreqCross; + + // Cross-pol vectors + vector > VgainCross; + vector > VphaseCross; + vector > VgainTopCross; + vector > VphaseTopCross; + vector > HgainCross; + vector > HphaseCross; double antenna_source_medium_n; //Define impedance and gain for receiving antenna @@ -247,11 +260,14 @@ class Detector { //Define impedance and gain for transmitting antenna vector RealImpedanceTx; vector ImagImpedanceTx; - int Tx_freq_init; - int Tx_freq_width; + double Tx_freq_init; + double Tx_freq_width; vector TxFreq; + vector TxFreqCross; vector > Txgain; vector > Txphase; + vector > TxgainCross; + vector > TxphaseCross; double Txantenna_source_medium_n; void ReadImpedance(string filename, vector *TempRealImpedance, vector *TempImagImpedance); void ReadAllAntennaImpedance(Settings *settings1); @@ -348,6 +364,11 @@ class Detector { vector transVTop_databin; vector transH_databin; + // Cross-pol data bins + vector transVCross_databin; + vector transVTopCross_databin; + vector transHCross_databin; + void ReadAmplifierNoiseFigure(Settings *settings1); vector< vector > amplifierNoiseFig_ch; @@ -384,7 +405,7 @@ class Detector { double GetGain(double freq, double theta, double phi, int ant_m, int ant_o, double antenna_target_medium_n); //read antenna gain at certain angle, certain type, and certain orientation double GetGain(double freq, double theta, double phi, int ant_m, double antenna_target_medium_n); //read antenna gain at certain angle, certain type. (orientation : default) - double GetGain_1D_OutZero(double freq, double theta, double phi, int ant_m, double antenna_target_medium_n, int string_number=0, int ant_number=0, bool useInTransmitterMode=false); //read antenna gain at certain angle, certain type. (orientation : default) and use 1-D interpolation to get gain, if freq bigger than freq range, return 0 gain + double GetGain_1D_OutZero(double freq, double theta, double phi, int ant_m, double antenna_target_medium_n, int string_number=0, int ant_number=0, bool useInTransmitterMode=false, bool useCrossPol=false); //read antenna gain at certain angle, certain type. (orientation : default) and use 1-D interpolation to get gain, if freq bigger than freq range, return 0 gain //Creating function to interpolate antenna impedance to frequency binning. double GetImpedance(double freq, int ant_m=0, int ant_number=0, bool useInTransmitterMode=false); @@ -394,7 +415,7 @@ class Detector { double GetAntPhase(double freq, double theta, double phi, int ant_m, double antenna_target_medium_n); // return antenna phase with 2-D interpolation - double GetAntPhase_1D(double freq, double theta, double phi, int ant_m, double antenna_target_medium_n, bool useInTransmitterMode=false); // return antenna phase with 1-D interpolation + double GetAntPhase_1D(double freq, double theta, double phi, int ant_m, double antenna_target_medium_n, bool useInTransmitterMode=false, bool useCrossPol=false); // return antenna phase with 1-D interpolation double GetFilterGain(int bin) { return FilterGain[bin]; } // same bin with Vgain, Hgain diff --git a/RayTrace.h b/RayTrace.h index 18f975d1..f49ece31 100644 --- a/RayTrace.h +++ b/RayTrace.h @@ -421,11 +421,11 @@ namespace RayTrace{ struct positionRecordingWrapper : public positionType{ rkStepRecord rec; - positionRecordingWrapper(){ + positionRecordingWrapper(){ rec.length=0.0; std::fill(&rec.z[0],&rec.z[6],0.0); } - positionRecordingWrapper(const positionType& pos):positionType(pos){ + positionRecordingWrapper(const positionType& pos):positionType(pos){ rec.length=0.0; std::fill(&rec.z[0],&rec.z[6],0.0); } diff --git a/Report.cc b/Report.cc index 5bd56548..88dd7f2e 100644 --- a/Report.cc +++ b/Report.cc @@ -165,6 +165,8 @@ void Antenna_r::Prepare_Outputs(int n_interactions) { Pol_vector.resize(n_interactions); vmmhz.resize(n_interactions); Heff.resize(n_interactions); + Heff_copol.resize(n_interactions); + Heff_crosspol.resize(n_interactions); Mag.resize(n_interactions); Fresnel.resize(n_interactions); Pol_factor.resize(n_interactions); @@ -211,6 +213,8 @@ void Antenna_r::clear() { Pol_vector.clear(); vmmhz.clear(); Heff.clear(); + Heff_copol.clear(); + Heff_crosspol.clear(); Mag.clear(); Fresnel.clear(); Pol_factor.clear(); @@ -280,6 +284,8 @@ void Antenna_r::clear_useless(Settings *settings1) { // Clear EM signal and waveform data for each interaction and ray vmmhz.clear(); Heff.clear(); + Heff_copol.clear(); + Heff_crosspol.clear(); Vfft.clear(); Vfft_noise.clear(); V.clear(); @@ -1187,7 +1193,7 @@ void Report::ModelRay( PropagateSignal( dT_forfft, outbin, Tarray_vector, Earray_vector, T_forint, - interaction_idx, ray_idx, ray_output, launch_vector, time_diff_birefringence, + interaction_idx, ray_idx, ray_output, launch_vector, receive_vector, fresnel, time_diff_birefringence, Pol_vector_src, Pol_vector, n_trg_slappy, n_trg_pokey, antenna_r, antenna_d, gain_ch_no, j, k, birefringence, detector, event, icemodel, settings); @@ -1236,7 +1242,7 @@ void Report::ModelRay( PropagateSignal( dT_forfft, waveform_bin, ArbitraryWaveform_T_vector, ArbitraryWaveform_V_vector, T_forint, - interaction_idx, ray_idx, ray_output, launch_vector, time_diff_birefringence, + interaction_idx, ray_idx, ray_output, launch_vector, receive_vector, fresnel, time_diff_birefringence, Pol_vector_src, Pol_vector, n_trg_slappy, n_trg_pokey, antenna_r, antenna_d, gain_ch_no, j, k, birefringence, detector, event, icemodel, settings); @@ -1286,7 +1292,7 @@ void Report::ModelRay( PropagateSignal( dT_forfft, waveform_bin, signal->PulserWaveform_T, signal->PulserWaveform_V, T_forint, - interaction_idx, ray_idx, ray_output, launch_vector, time_diff_birefringence, + interaction_idx, ray_idx, ray_output, launch_vector, receive_vector, fresnel, time_diff_birefringence, Pol_vector_src, Pol_vector, n_trg_slappy, n_trg_pokey, antenna_r, antenna_d, gain_ch_no, j, k, birefringence, detector, event, icemodel, settings); @@ -1319,21 +1325,14 @@ void Report::ModelRay( double newPol_vectorZ = cos(psi)*sin(theta); Vector Pol_vector = Vector(newPol_vectorX, newPol_vectorY, newPol_vectorZ); - //Apply Fresnel factors for magnification and 1/r dependence - icemodel->GetFresnel( - ray_output[1][ray_idx], // launch_angle - ray_output[2][ray_idx], // rec_angle - ray_output[3][ray_idx], // reflect_angle - event->Nu_Interaction[interaction_idx].posnu, - launch_vector, - receive_vector, - settings, - fresnel, - Pol_vector); // input src Pol and return Pol at trg + std::vector scaled_waveform; + scaled_waveform.reserve(signal->InputVoltage_V.size()); + for (const auto& v : signal->InputVoltage_V) + scaled_waveform.push_back(v * 1e-3); PropagateSignal( - dT_forfft, waveform_bin, signal->InputVoltage_T, signal->InputVoltage_V, T_forint, - interaction_idx, ray_idx, ray_output, launch_vector, time_diff_birefringence, + dT_forfft, waveform_bin, signal->InputVoltage_T, scaled_waveform, T_forint, + interaction_idx, ray_idx, ray_output, launch_vector, receive_vector, fresnel, time_diff_birefringence, Pol_vector_src, Pol_vector, n_trg_slappy, n_trg_pokey, antenna_r, antenna_d, gain_ch_no, j, k, birefringence, detector, event, icemodel, settings); @@ -1356,7 +1355,7 @@ void Report::ModelRay( PropagateSignal( dT_forfft, CP_bin, detector->CalPulserWF_ns, detector->CalPulserWF_V, T_forint, - interaction_idx, ray_idx, ray_output, launch_vector, time_diff_birefringence, + interaction_idx, ray_idx, ray_output, launch_vector, receive_vector, fresnel, time_diff_birefringence, Pol_vector_src, Pol_vector, n_trg_slappy, n_trg_pokey, antenna_r, antenna_d, gain_ch_no, j, k, birefringence, detector, event, icemodel, settings); @@ -1404,6 +1403,8 @@ void Report::GetRayParameters( antenna_r->reflect_ang[interaction_idx].push_back(ray_output[3][ray_idx]); antenna_r->vmmhz[interaction_idx].resize(ray_idx + 1); antenna_r->Heff[interaction_idx].resize(ray_idx + 1); + antenna_r->Heff_copol[interaction_idx].resize(ray_idx + 1); + antenna_r->Heff_crosspol[interaction_idx].resize(ray_idx + 1); antenna_r->Vm_zoom[interaction_idx].resize(ray_idx + 1); antenna_r->Vm_zoom_T[interaction_idx].resize(ray_idx + 1); antenna_r->Vfft[interaction_idx].resize(ray_idx + 1); @@ -1492,7 +1493,7 @@ void Report::DetermineWFBins( void Report::PropagateSignal( double dT_forfft, int efield_length, vector< double > efield_time, vector< double > efield, double *T_forint, - int interaction_idx, int ray_idx, vector > ray_output, Position launch_vector, double time_diff_birefringence, + int interaction_idx, int ray_idx, vector > ray_output, Position launch_vector, Position receive_vector, double fresnel, double time_diff_birefringence, Vector Pol_vector_src, Vector Pol_vector, Vector n_trg_slappy, Vector n_trg_pokey, Antenna_r *antenna_r, Antenna *antenna_d, int gain_ch_no, int j, int k, Birefringence *birefringence, Detector *detector, Event *event, IceModel *icemodel, Settings *settings @@ -1571,41 +1572,69 @@ void Report::PropagateSignal( // Get ant gain with 2-D interpolation double n_eff = GetEffectiveIndex(icemodel->GetN(*antenna_d)); - double heff_lastbin = GaintoHeight( - detector->GetGain_1D_OutZero( - freq_tmp *1.E-6, antenna_theta, antenna_phi, - antenna_d->type, n_eff, j, k), - freq_tmp, + + // Co-pol effective height for the last bin + double heff_copol_lastbin = GaintoHeight( + detector->GetGain_1D_OutZero(freq_tmp * 1.E-6, antenna_theta, antenna_phi, + antenna_d->type, n_eff, j, k, false, false), + freq_tmp, + icemodel->GetN(*antenna_d), + detector->GetImpedance(freq_tmp * 1.E-6, antenna_d->type, k)); + + // Cross-pol effective height for the last bin + double heff_crosspol_lastbin = GaintoHeight( + detector->GetGain_1D_OutZero(freq_tmp * 1.E-6, antenna_theta, antenna_phi, + antenna_d->type, n_eff, j, k, false, true), + freq_tmp, icemodel->GetN(*antenna_d), - detector->GetImpedance(freq_tmp*1.E-6, antenna_d->type, k)); + detector->GetImpedance(freq_tmp * 1.E-6, antenna_d->type, k)); // Tx effective height for last bin. Currently locked to standard ARA Vpol and HPol antennas. // Need to add selection mode. // Only used in EVENT_MODE == 12 double Tx_theta = 0; double Tx_phi = 0; - double heff_Tx_lastbin = 0; - if (settings->EVENT_TYPE == 11 || settings->EVENT_TYPE == 12){ // Pulser simulations + double heff_Tx_copol_lastbin=0; + double heff_Tx_crosspol_lastbin=0; + if (settings->EVENT_TYPE == 12){ // Pulser simulations - //Defining polarization at the source (using launch_vector (a unit vector)) + //Defining transmission angle at the source (using launch_vector (a unit vector)) double theta = acos(launch_vector[2]); double phi = atan2(launch_vector[1],launch_vector[0]); - double Tx_theta = theta*180/PI; - double Tx_phi = phi*180/PI; + Tx_theta = theta*180/PI; + Tx_phi = phi*180/PI; double n_eff = GetEffectiveIndex(icemodel->GetN(*antenna_d)); - double heff_Tx_lastbin = GaintoHeight( - detector->GetGain_1D_OutZero(freq_tmp *1.E-6, Tx_theta, Tx_phi, 0, n_eff, 0, 0, true), - freq_tmp, + + // Tx effective height for co-pol and cross-pol for the last bin + heff_Tx_copol_lastbin = GaintoHeight( + detector->GetGain_1D_OutZero(freq_tmp * 1.E-6, Tx_theta, Tx_phi, 0, n_eff, 0, 0, true, false), + freq_tmp, + icemodel->GetN(*antenna_d), + detector->GetImpedance(freq_tmp * 1.E-6, 0, 0, true)); + + + heff_Tx_crosspol_lastbin = GaintoHeight( + detector->GetGain_1D_OutZero(freq_tmp * 1.E-6, Tx_theta, Tx_phi, 0, n_eff, 0, 0, true, true), + freq_tmp, icemodel->GetN(*antenna_d), - detector->GetImpedance(freq_tmp*1.E-6, 0, 0, true)); + detector->GetImpedance(freq_tmp * 1.E-6, 0, 0, true)); + // End Tx effective height for last bin } + else if (event->IsCalpulser > 0){ // Calibration Pulser simulations Tx_theta = ray_output[1][ray_idx] *DEGRAD; // from 0 to 180 - double n_eff = GetEffectiveIndex(icemodel->GetN(event->Nu_Interaction[interaction_idx].posnu)); - heff_Tx_lastbin = GaintoHeight( - detector->GetGain_1D_OutZero(freq_tmp *1.E-6, Tx_theta, antenna_phi, antenna_d->type, n_eff, j, k), - freq_tmp, icemodel->GetN(event->Nu_Interaction[interaction_idx].posnu)); + double n_eff = GetEffectiveIndex(icemodel->GetN(event->Nu_Interaction[interaction_idx].posnu)); + + // Tx effective height for co-pol and cross-pol for the last bin + heff_Tx_copol_lastbin = GaintoHeight( + detector->GetGain_1D_OutZero(freq_tmp * 1.E-6, Tx_theta, antenna_phi, antenna_d->type, n_eff, j, k, true, false), + freq_tmp, icemodel->GetN(event->Nu_Interaction[interaction_idx].posnu)); + + heff_Tx_crosspol_lastbin = GaintoHeight( + detector->GetGain_1D_OutZero(freq_tmp * 1.E-6, Tx_theta, antenna_phi, antenna_d->type, n_eff, j, k, true, true), + freq_tmp, icemodel->GetN(event->Nu_Interaction[interaction_idx].posnu)); + // End Tx effective height for last bin if (event->IsCalpulser == 1) { Pol_vector = n_trg_slappy; @@ -1629,12 +1658,23 @@ void Report::PropagateSignal( freq_tmp = dF_Nnew *((double) n + 0.5); // in Hz 0.5 to place the middle of the bin and avoid zero freq double n_eff = GetEffectiveIndex(icemodel->GetN(*antenna_d)); - double heff = GaintoHeight( - detector->GetGain_1D_OutZero(freq_tmp *1.E-6, antenna_theta, antenna_phi, antenna_d->type, n_eff, j, k), - freq_tmp, + + // Co-pol effective height + double heff_copol = GaintoHeight( + detector->GetGain_1D_OutZero(freq_tmp * 1.E-6, antenna_theta, antenna_phi, antenna_d->type, n_eff, j, k, false, false), + freq_tmp, + icemodel->GetN(*antenna_d), + detector->GetImpedance(freq_tmp * 1.E-6, antenna_d->type, k)); + + // Cross-pol effective height + double heff_crosspol = GaintoHeight( + detector->GetGain_1D_OutZero(freq_tmp * 1.E-6, antenna_theta, antenna_phi, antenna_d->type, n_eff, j, k, false, true), + freq_tmp, icemodel->GetN(*antenna_d), - detector->GetImpedance(freq_tmp*1.E-6, antenna_d->type, k)); - antenna_r->Heff[interaction_idx][ray_idx].push_back(heff); + detector->GetImpedance(freq_tmp * 1.E-6, antenna_d->type, k)); + + antenna_r->Heff_copol[interaction_idx][ray_idx].push_back(heff_copol); + antenna_r->Heff_crosspol[interaction_idx][ray_idx].push_back(heff_crosspol); //apply freq dependent attenuation model if in neutrino mode if (settings->EVENT_TYPE == 0 && settings->USE_ARA_ICEATTENU == 2) { @@ -1666,56 +1706,106 @@ void Report::PropagateSignal( // Apply transmitter effective height if in PVA Pulser mode if (settings->EVENT_TYPE == 12){ double n_eff = GetEffectiveIndex(icemodel->GetN(*antenna_d)); - double heff_Tx = GaintoHeight( - detector->GetGain_1D_OutZero(freq_tmp *1.E-6, Tx_theta, Tx_phi, 0, n_eff, 0, 0, true), - freq_tmp, + + // Co-pol effective height Tx + double heff_Tx_copol = GaintoHeight( + detector->GetGain_1D_OutZero(freq_tmp *1.E-6, Tx_theta, Tx_phi, 0, n_eff, 0, 0, true, false), + freq_tmp, icemodel->GetN(*antenna_d), - detector->GetImpedance(freq_tmp*1.E-6, 0, 0, true)); - if (n > 0) { - ApplyAntFactors_Tdomain( - detector->GetAntPhase_1D(freq_tmp *1.e-6, Tx_theta, Tx_phi, 0, n_eff), - heff_Tx, Pol_vector, 0, Pol_factor, V_forfft[2 *n], V_forfft[2 *n + 1], - settings, Tx_theta, Tx_phi, freq_tmp, detector->GetImpedance(freq_tmp*1.E-6, 0, 0, true), true); + detector->GetImpedance(freq_tmp*1.E-6, 0, 0, true)); + + // Cross-pol effective height Tx + double heff_Tx_crosspol = GaintoHeight( + detector->GetGain_1D_OutZero(freq_tmp *1.E-6, Tx_theta, Tx_phi, 0, n_eff, 0, 0, true, true), + freq_tmp, + icemodel->GetN(*antenna_d), + detector->GetImpedance(freq_tmp*1.E-6, 0, 0, true)); + + // Retrieve co-pol and cross-pol phases + double phase_copol_tx = detector->GetAntPhase_1D(freq_tmp *1.e-6, Tx_theta, Tx_phi, 0, n_eff, false, false); //ASG change when fixing phases + double phase_crosspol_tx = detector->GetAntPhase_1D(freq_tmp *1.e-6, Tx_theta, Tx_phi, 0, n_eff, false, true); //ASG change when fixing phases + + if (n > 0) { + ApplyAntFactors_Tdomain(phase_copol_tx, phase_crosspol_tx, heff_Tx_copol, heff_Tx_crosspol, Pol_vector, 0, Pol_factor, + V_forfft[2 * n], V_forfft[2 * n + 1], settings, Tx_theta, Tx_phi, freq_tmp, true, false); } else { ApplyAntFactors_Tdomain_FirstTwo( - heff_Tx, heff_Tx_lastbin, Pol_vector, 0, Pol_factor, - V_forfft[2 *n], V_forfft[2 *n + 1], Tx_theta, Tx_phi, freq_tmp); + heff_Tx_copol, heff_Tx_copol_lastbin, heff_Tx_crosspol, heff_Tx_crosspol_lastbin, Pol_vector, 0, + Pol_factor, V_forfft[2 * n], V_forfft[2 * n + 1], settings, Tx_theta, Tx_phi, freq_tmp, true, false); } + + //Apply Fresnel factors for magnification and 1/r dependence + icemodel->GetFresnel( + ray_output[1][ray_idx], // launch_angle + ray_output[2][ray_idx], // rec_angle + ray_output[3][ray_idx], // reflect_angle + event->Nu_Interaction[interaction_idx].posnu, + launch_vector, + receive_vector, + settings, + fresnel, + Pol_vector); // input src Pol and return Pol at trg + } else if (event->IsCalpulser > 0){ // apply ant factors (transmitter ant) double n_eff = GetEffectiveIndex(icemodel->GetN(*antenna_d)); - heff = GaintoHeight( + + // Co-pol effective height + double heff_Tx_copol = GaintoHeight( + detector->GetGain_1D_OutZero( + freq_tmp * 1.E-6, Tx_theta, antenna_phi, antenna_d->type, n_eff, j, k, false, false), + freq_tmp, icemodel->GetN(event->Nu_Interaction[interaction_idx].posnu)); + + // Cross-pol effective height + double heff_Tx_crosspol = GaintoHeight( detector->GetGain_1D_OutZero( freq_tmp *1.E-6, // to MHz - Tx_theta, antenna_phi, antenna_d->type, n_eff, j, k), + Tx_theta, antenna_phi, antenna_d->type, n_eff, j, k, false, true), freq_tmp, icemodel->GetN(event->Nu_Interaction[interaction_idx].posnu)); + + // Retrieve co-pol and cross-pol phases + double phase_copol_tx = detector->GetAntPhase_1D(freq_tmp *1.e-6, Tx_theta, antenna_phi, antenna_d->type, n_eff, true, false); + double phase_crosspol_tx = detector->GetAntPhase_1D(freq_tmp *1.e-6, Tx_theta, antenna_phi, antenna_d->type, n_eff, true, true); if (n > 0) { - ApplyAntFactors_Tdomain( - detector->GetAntPhase_1D(freq_tmp *1.e-6, Tx_theta, antenna_phi, antenna_d->type, n_eff), - heff, Pol_vector, antenna_d->type, Pol_factor, - V_forfft[2 *n], V_forfft[2 *n + 1], settings, true, antenna_theta, antenna_phi, freq_tmp); + ApplyAntFactors_Tdomain(phase_copol_tx, phase_crosspol_tx, heff_Tx_copol, heff_Tx_crosspol, Pol_vector, antenna_d->type, Pol_factor, + V_forfft[2 * n], V_forfft[2 * n + 1], settings, Tx_theta, Tx_phi, freq_tmp, true, false); } else { ApplyAntFactors_Tdomain_FirstTwo( - heff, heff_Tx_lastbin, Pol_vector, antenna_d->type, Pol_factor, - V_forfft[2 *n], V_forfft[2 *n + 1], antenna_theta, antenna_phi, freq_tmp); + heff_Tx_copol, heff_Tx_copol_lastbin, heff_Tx_crosspol, heff_Tx_crosspol_lastbin, Pol_vector, antenna_d->type, Pol_factor, + V_forfft[2 * n], V_forfft[2 * n + 1], settings, Tx_theta, Tx_phi, freq_tmp, true, false); } + + //Apply Fresnel factors for magnification and 1/r dependence + icemodel->GetFresnel( + ray_output[1][ray_idx], // launch_angle + ray_output[2][ray_idx], // rec_angle + ray_output[3][ray_idx], // reflect_angle + event->Nu_Interaction[interaction_idx].posnu, + launch_vector, + receive_vector, + settings, + fresnel, + Pol_vector); // input src Pol and return Pol at trg + } + double phase_copol = detector->GetAntPhase_1D(freq_tmp * 1.e-6, antenna_theta, antenna_phi, + antenna_d->type, n_eff, false, false); + double phase_crosspol = detector->GetAntPhase_1D(freq_tmp * 1.e-6, antenna_theta, antenna_phi, + antenna_d->type, n_eff, false, true); + // apply ant factors if (n > 0) { - ApplyAntFactors_Tdomain( - detector->GetAntPhase_1D(freq_tmp *1.e-6, antenna_theta, antenna_phi, antenna_d->type, n_eff), - heff, Pol_vector, antenna_d->type, Pol_factor, - V_forfft[2 *n], V_forfft[2 *n + 1], settings, antenna_theta, antenna_phi, freq_tmp); + ApplyAntFactors_Tdomain(phase_copol, phase_crosspol, heff_copol, heff_crosspol, Pol_vector, antenna_d->type, Pol_factor, + V_forfft[2 * n], V_forfft[2 * n + 1], settings, antenna_theta, antenna_phi, freq_tmp); } else { - ApplyAntFactors_Tdomain_FirstTwo( - heff, heff_lastbin, Pol_vector, antenna_d->type, Pol_factor, - V_forfft[2 *n], V_forfft[2 *n + 1], antenna_theta, antenna_phi, freq_tmp); + ApplyAntFactors_Tdomain_FirstTwo(heff_copol, heff_copol_lastbin, heff_crosspol, heff_crosspol_lastbin, Pol_vector, antenna_d->type, Pol_factor, + V_forfft[2 * n], V_forfft[2 * n + 1], settings, antenna_theta, antenna_phi, freq_tmp); } // apply entire elect chain gain, phase @@ -3564,12 +3654,30 @@ void Report::ApplyAntFactors( } +void Report::ApplyAntFactors_Tdomain( + double AntPhase, double heff, Vector &Pol_vector, int ant_type, double &pol_factor, double &vm_real, double &vm_img, + Settings *settings1, double antenna_theta, double antenna_phi, double freq, bool useInTransmitterMode, bool applyInverse +) { + double heff_copol = heff; + double heff_crosspol = 0.0; + double phase_copol = AntPhase; + double phase_crosspol = 0.0; + + ApplyAntFactors_Tdomain( + phase_copol, phase_crosspol, + heff_copol, heff_crosspol, + Pol_vector, ant_type, + pol_factor, vm_real, vm_img, + settings1, antenna_theta, antenna_phi, + freq, useInTransmitterMode, applyInverse); +} + -void Report::ApplyAntFactors_Tdomain ( - double AntPhase, double heff, Vector &Pol_vector, - int ant_type, double &pol_factor, double &vm_real, double &vm_img, Settings *settings1, double antenna_theta, double antenna_phi, - double freq, bool useInTransmitterMode, bool applyInverse -) { +void Report::ApplyAntFactors_Tdomain(double phase_copol, double phase_crosspol, double heff_copol, double heff_crosspol, + Vector &Pol_vector, int ant_type, double &pol_factor, double &vm_real, double &vm_img, + Settings *settings1, double antenna_theta, double antenna_phi, double freq, + bool useInTransmitterMode, bool applyInverse +) { /* Report::ApplyAntFactors_Tdomain() Purpose: Multiply (or divide) voltage in Fourier space by the antenna gain and phase. @@ -3602,7 +3710,22 @@ void Report::ApplyAntFactors_Tdomain ( if(applyInverse==true){ phaseSign*=-1.; amplitudeSign*=-1;}; //Calculate the polarization factor, which is essentially the vector component of the Electric field that is projected onto the antenna. - pol_factor = calculatePolFactor(Pol_vector, ant_type, antenna_theta, antenna_phi); + //Do both co-pol and cross-pol polarization factors + double pol_factor_copol = calculatePolFactor(Pol_vector, ant_type, antenna_theta, antenna_phi); + double pol_factor_crosspol = calculatePolFactor(Pol_vector, 1 - ant_type, antenna_theta, antenna_phi); + + // Calculate amplitude amplifications from co-pol and cross-pol + double v_amplification_copol = heff_copol * pol_factor_copol; + double v_amplification_crosspol = 0.0; + + //turn on cross-pol for receiver + if (settings1->CROSSPOL_RX){ + v_amplification_crosspol = heff_crosspol * pol_factor_crosspol; + } + + // Add the contributions linearly for Rx + double v_amplification = v_amplification_copol + v_amplification_crosspol; + if ( settings1->PHASE_SKIP_MODE != 1 ) { double phase_current; if ( vm_real != 0. ) { @@ -3619,32 +3742,72 @@ void Report::ApplyAntFactors_Tdomain ( else phase_current = 0.; } //Calculate amplitude via the real and imaginary components. - double v_amp = sqrt(vm_real*vm_real + vm_img*vm_img); - //Apply effective height and polarization factors to ampltitude. - v_amp *= pow(heff*pol_factor, amplitudeSign); + double v_amp_in = sqrt(vm_real*vm_real + vm_img*vm_img); + double v_amp = v_amp_in; //save the incomming voltage + // Apply amplitude sign for inversion if necessary + v_amp *= pow(v_amplification, amplitudeSign); //If in transmitter mode, we must apply additional frequency and impedance terms to the amplitude. if (useInTransmitterMode==true){ phase_current += PI/2; - // The factors of two in this line are still under consideration, see GitHub Issue #175 - v_amp *= pow(freq/CLIGHT*(Z0/Zr)/4/sqrt(2.), amplitudeSign); + double psi = 0.0; + double delta_psi = TMath::ATan(heff_crosspol / heff_copol); // Tx cross-pol tilt + double theta = antenna_theta*PI/180.0; + double phi = antenna_phi*PI/180.0; + + //turn on cross-pol Tx + if (settings1->CROSSPOL_TX == 0){ + heff_crosspol = 0.0; //need to be turned off for the calculation of amplitude v_amp + delta_psi = 0.0; // turn off cross-pol tx tilt + } + + //Adjust polarization vector + double newPol_vectorX = -cos(psi+delta_psi)*cos(theta)*cos(phi) + sin(psi+delta_psi)*sin(phi); + double newPol_vectorY = -cos(psi+delta_psi)*cos(theta)*sin(phi) - sin(psi+delta_psi)*cos(phi); + double newPol_vectorZ = cos(psi+delta_psi)*sin(theta); + + Pol_vector = Vector(newPol_vectorX, newPol_vectorY, newPol_vectorZ); + + //copol and cross-pol add quadratically in E-field + v_amp *= pow(freq/CLIGHT*(Z0/Zr)/4/sqrt(2.)*(1/v_amplification)*(sqrt(heff_crosspol* heff_crosspol + heff_copol*heff_copol )), amplitudeSign); } - //Calculate the real and imaginary terms using the new ampltitude and phase. - vm_real = v_amp * cos( phase_current + (phaseSign * AntPhase*RADDEG) ); - vm_img = v_amp * sin( phase_current + (phaseSign * AntPhase*RADDEG) ); + // Calculate the combined real and imaginary terms from co-pol and cross-pol + double vm_real_copol = v_amp_in * (v_amplification_copol * cos(phase_current + (phaseSign * phase_copol * RADDEG))); + double vm_img_copol = v_amp_in * (v_amplification_copol * sin(phase_current + (phaseSign * phase_copol * RADDEG))); + + double vm_real_crosspol = v_amp_in * (v_amplification_crosspol * cos(phase_current + (phaseSign * phase_crosspol * RADDEG))); + double vm_img_crosspol = v_amp_in * (v_amplification_crosspol * sin(phase_current + (phaseSign * phase_crosspol * RADDEG))); + + // Combine co-pol and cross-pol real and imaginary parts + vm_real = vm_real_copol + vm_real_crosspol; + vm_img = vm_img_copol + vm_img_crosspol; } else { // only amplitude - vm_real *= pow(heff * pol_factor, amplitudeSign); - vm_img *= pow(heff * pol_factor, amplitudeSign); + vm_real *= pow(v_amplification, amplitudeSign); + vm_img *= pow(v_amplification, amplitudeSign); } } +void Report::ApplyAntFactors_Tdomain_FirstTwo( + double heff, double heff_lastbin, Vector &Pol_vector, int ant_type, double &pol_factor, double &vm_bin0, double &vm_bin1, + Settings *settings1, double antenna_theta, double antenna_phi, double freq, bool useInTransmitterMode, bool applyInverse +) { + // Default to no cross-pol + double heff_crosspol = 0.0; + double heff_crosspol_lastbin = 0.0; + + ApplyAntFactors_Tdomain_FirstTwo( + heff, heff_lastbin, heff_crosspol, heff_crosspol_lastbin, + Pol_vector, ant_type, pol_factor, vm_bin0, vm_bin1, + settings1, antenna_theta, antenna_phi, freq, useInTransmitterMode, applyInverse); +} -void Report::ApplyAntFactors_Tdomain_FirstTwo ( - double heff, double heff_lastbin, Vector &Pol_vector, - int ant_type, double &pol_factor, double &vm_bin0, double &vm_bin1, double antenna_theta, double antenna_phi, - double freq, bool useInTransmitterMode, bool applyInverse +void Report::ApplyAntFactors_Tdomain_FirstTwo(double heff_copol, double heff_copol_lastbin, double heff_crosspol, + double heff_crosspol_lastbin, Vector &Pol_vector, int ant_type, + double &pol_factor, double &vm_bin0, double &vm_bin1, Settings *settings1, + double antenna_theta, double antenna_phi, double freq, bool useInTransmitterMode, + bool applyInverse ) { /* Report::ApplyAntFactors_Tdomain_FirstTwo() Purpose: @@ -3660,45 +3823,103 @@ void Report::ApplyAntFactors_Tdomain_FirstTwo ( //double pol_factor; pol_factor = calculatePolFactor(Pol_vector, ant_type, antenna_theta, antenna_phi); - //First step in the voltage calculation. Both Tx and Rx mode use this. - vm_bin0 *= pow(heff * pol_factor, amplitudeSign); - vm_bin1 *= pow(heff_lastbin * pol_factor, amplitudeSign); - - if (useInTransmitterMode) { - // The factors of two in these two lines are currently up for debate, see GitHub Issue #175 - vm_bin0 *= pow(freq/CLIGHT*(Z0/(Zr))/4/sqrt(2.), amplitudeSign); - vm_bin1 *= pow(freq/CLIGHT*(Z0/(Zr))/4/sqrt(2.), amplitudeSign); + // Initialize cross-pol amplification factors + double pol_factor_crosspol = 0.0; + double v_amplification_copol_bin0 = heff_copol * pol_factor; + double v_amplification_copol_bin1 = heff_copol_lastbin * pol_factor; + double v_amplification_crosspol_bin0 = 0.0; + double v_amplification_crosspol_bin1 = 0.0; + + // Handle cross-pol factors for receiver mode + if (!useInTransmitterMode && settings1->CROSSPOL_RX) { + pol_factor_crosspol = calculatePolFactor(Pol_vector, 1 - ant_type, antenna_theta, antenna_phi); + v_amplification_crosspol_bin0 = heff_crosspol * pol_factor_crosspol; + v_amplification_crosspol_bin1 = heff_crosspol_lastbin * pol_factor_crosspol; } - else { - vm_bin0 *= pow(heff * pol_factor, amplitudeSign); - vm_bin1 *= pow(heff_lastbin * pol_factor, amplitudeSign); + + // Combine co-pol and cross-pol amplifications + double v_amplification_bin0 = v_amplification_copol_bin0 + v_amplification_crosspol_bin0; + double v_amplification_bin1 = v_amplification_copol_bin1 + v_amplification_crosspol_bin1; + + // Apply amplification for receiver mode + vm_bin0 *= pow(v_amplification_bin0, amplitudeSign); + vm_bin1 *= pow(v_amplification_bin1, amplitudeSign); + + if (useInTransmitterMode) { + // Handle cross-pol factors for transmitter mode + if (settings1->CROSSPOL_TX) { + pol_factor_crosspol = calculatePolFactor(Pol_vector, 1 - ant_type, antenna_theta, antenna_phi); + v_amplification_crosspol_bin0 = heff_crosspol * pol_factor_crosspol; + v_amplification_crosspol_bin1 = heff_crosspol_lastbin * pol_factor_crosspol; + } + + // Co-pol and cross-pol add quadratically in E-field + double tx_amplification_bin0 = freq / CLIGHT * (Z0 / Zr) / 4 / sqrt(2.0) * (1/v_amplification_bin0) * + (sqrt(heff_crosspol * heff_crosspol + heff_copol * heff_copol)); + double tx_amplification_bin1 = freq / CLIGHT * (Z0 / Zr) / 4 / sqrt(2.0) * (1/v_amplification_bin1) * + (sqrt(heff_crosspol_lastbin * heff_crosspol_lastbin + heff_copol_lastbin * heff_copol_lastbin)); + + vm_bin0 *= pow(tx_amplification_bin0, amplitudeSign); + vm_bin1 *= pow(tx_amplification_bin1, amplitudeSign); } } -void Report::InvertAntFactors_Tdomain ( - double AntPhase, double heff, Vector &Pol_vector, - int ant_type, double &pol_factor, double &vm_real, double &vm_img, Settings *settings1, double antenna_theta, double antenna_phi, - double freq, bool useInTransmitterMode -) { - /* Report::InvertAntFactors_Tdomain() - Purpose: Inverts the antenna factors in a convenient way by simply calling ApplyAntFactors with the boolean applyInverse enabled. +void Report::InvertAntFactors_Tdomain( + double AntPhase, double heff, Vector &Pol_vector, int ant_type, double &pol_factor, double &vm_real, double &vm_img, Settings *settings1, + double antenna_theta, double antenna_phi, double freq, bool useInTransmitterMode +) { + ApplyAntFactors_Tdomain( + AntPhase, heff, Pol_vector, ant_type, + pol_factor, vm_real, vm_img, + settings1, antenna_theta, antenna_phi, + freq, useInTransmitterMode, true); // applyInverse = true +} + +void Report::InvertAntFactors_Tdomain(double AntPhase_copol, double AntPhase_crosspol, + double heff_copol, double heff_crosspol, + Vector &Pol_vector, int ant_type, + double &pol_factor, double &vm_real, double &vm_img, + Settings *settings1, double antenna_theta, + double antenna_phi, double freq, + bool useInTransmitterMode) { + /* Report::InvertAntFactors_Tdomain_new() + Purpose: Inverts the antenna factors for both co-polarization and cross-polarization + by simply calling ApplyAntFactors_Tdomain_new with the boolean applyInverse enabled. */ - ApplyAntFactors_Tdomain (AntPhase, heff, Pol_vector, ant_type, pol_factor, vm_real, vm_img, settings1, antenna_theta, antenna_phi, freq, useInTransmitterMode, true); + ApplyAntFactors_Tdomain(AntPhase_copol, AntPhase_crosspol, heff_copol, heff_crosspol, + Pol_vector, ant_type, pol_factor, vm_real, vm_img, + settings1, antenna_theta, antenna_phi, freq, useInTransmitterMode, true); + +} + +void Report::InvertAntFactors_Tdomain_FirstTwo( + double heff, double heff_lastbin, Vector &Pol_vector, int ant_type, double &pol_factor, double &vm_bin0, + double &vm_bin1, Settings *settings1, double antenna_theta, double antenna_phi, double freq, bool useInTransmitterMode +) { + ApplyAntFactors_Tdomain_FirstTwo( + heff, heff_lastbin, + Pol_vector, ant_type, pol_factor, vm_bin0, vm_bin1, + settings1, antenna_theta, antenna_phi, freq, + useInTransmitterMode, true); } -void Report::InvertAntFactors_Tdomain_FirstTwo ( - double heff, double heff_lastbin, Vector &Pol_vector, - int ant_type, double &pol_factor, double &vm_bin0, double &vm_bin1, double antenna_theta, double antenna_phi, - double freq, bool useInTransmitterMode -) { - /* Report::InvertAntFactors_Tdomain_FirstTwo() - Purpose: Inverts the antenna factors in a convenient way by simply calling ApplyAntFactors with the boolean applyInverse enabled. +void Report::InvertAntFactors_Tdomain_FirstTwo(double heff_copol, double heff_copol_lastbin, + double heff_crosspol, double heff_crosspol_lastbin, + Vector &Pol_vector, int ant_type, + double &pol_factor, double &vm_bin0, double &vm_bin1, + Settings *settings1, double antenna_theta, double antenna_phi, double freq, + bool useInTransmitterMode) { + /* Report::InvertAntFactors_Tdomain_FirstTwo_new() + Purpose: Inverts the antenna factors for both co-polarization and cross-polarization + by simply calling ApplyAntFactors_Tdomain_FirstTwo_new with the boolean applyInverse enabled. */ - ApplyAntFactors_Tdomain_FirstTwo (heff, heff_lastbin, Pol_vector, ant_type, pol_factor, vm_bin0, vm_bin1, antenna_theta, antenna_phi, freq, useInTransmitterMode, true); + ApplyAntFactors_Tdomain_FirstTwo(heff_copol, heff_copol_lastbin, heff_crosspol, heff_crosspol_lastbin, + Pol_vector, ant_type, pol_factor, vm_bin0, vm_bin1, + settings1, antenna_theta, antenna_phi, freq, useInTransmitterMode, true); } diff --git a/Report.h b/Report.h index 40b55049..b50d822f 100644 --- a/Report.h +++ b/Report.h @@ -78,6 +78,8 @@ class Antenna_r { vector < vector < vector > > vmmhz; //! // signal V/m/MHz for each freq bin vector < vector < vector > > Heff; //! // effective height for each freq bin + vector < vector < vector > > Heff_copol; // effective height for each freq bin + vector < vector < vector > > Heff_crosspol; // effective height for each freq bin vector < vector > Mag; //! // magnification factor vector < vector > Fresnel; //! // Fresnel factor vector < vector > Pol_factor; //! // Polarization factor @@ -257,8 +259,8 @@ class Report { IceModel *icemodel, Settings *settings1 ); void PropagateSignal( double dT_forfft, int efield_length, vector< double > efield_time, vector< double > efield, double *T_forint, - int interaction_idx, int ray_idx, vector > ray_output, Position launch_vector, double time_diff_birefringence, - Vector Pol_vector_src, Vector Pol_vector, Vector n_trg_slappy, Vector n_trg_pokey, + int interaction_idx, int ray_idx, vector > ray_output, Position launch_vector, Position receive_vector, double fresnel, + double time_diff_birefringence, Vector Pol_vector_src, Vector Pol_vector, Vector n_trg_slappy, Vector n_trg_pokey, Antenna_r *antenna_r, Antenna *antenna_d, int gain_ch_no, int j, int k, Birefringence *birefringence, Detector *detector, Event *event, IceModel *icemodel, Settings *settings); @@ -336,20 +338,41 @@ class Report { int ant_type, double &pol_factor, double &vm_real, double &vm_img, Settings *settings1, double antenna_theta, double antenna_phi, double freq, bool useInTransmitterMode=false, bool applyInverse=false ); + void ApplyAntFactors_Tdomain( + double phase_copol, double phase_crosspol, double heff_copol, double heff_crosspol, + Vector &Pol_vector, int ant_type, double &pol_factor, double &vm_real, double &vm_img, + Settings *settings1, double antenna_theta, double antenna_phi, double freq, bool useInTransmitterMode=false, bool applyInverse=false + ); void ApplyAntFactors_Tdomain_FirstTwo ( double heff, double heff_lastbin, Vector &Pol_vector, - int ant_type, double &pol_factor, double &vm_bin0, double &vm_bin1, double antenna_theta, double antenna_phi, + int ant_type, double &pol_factor, double &vm_bin0, double &vm_bin1, Settings *settings1, double antenna_theta, double antenna_phi, double freq, bool useInTransmitterMode=false, bool applyInverse=false ); + void ApplyAntFactors_Tdomain_FirstTwo( + double heff_copol, double heff_copol_lastbin, double heff_crosspol, double heff_crosspol_lastbin, + Vector &Pol_vector, int ant_type, double &pol_factor, double &vm_bin0, double &vm_bin1, + Settings *settings1, double antenna_theta, double antenna_phi, double freq, bool useInTransmitterMode=false, + bool applyInverse=false + ); void InvertAntFactors_Tdomain( - double AntPhase, double heff, Vector &Pol_vector, + double AntPhase, double heff, Vector &Pol_vector, int ant_type, double &pol_factor, double &vm_real, double &vm_img, Settings *settings1, double antenna_theta, double antenna_phi, double freq, bool useInTransmitterMode=false ); - void InvertAntFactors_Tdomain_FirstTwo ( - double heff, double heff_lastbin, Vector &Pol_vector, - int ant_type, double &pol_factor, double &vm_bin0, double &vm_bin1, double antenna_theta, double antenna_phi, - double freq, bool useInTransmitterMode=false + void InvertAntFactors_Tdomain( + double AntPhase_copol, double AntPhase_crosspol, double heff_copol, double heff_crosspol, Vector &Pol_vector, + int ant_type, double &pol_factor, double &vm_real, double &vm_img, Settings *settings1, double antenna_theta, + double antenna_phi, double freq, bool useInTransmitterMode=false + ); + void InvertAntFactors_Tdomain_FirstTwo ( + double heff, double heff_lastbin, Vector &Pol_vector, int ant_type, double &pol_factor, + double &vm_bin0, double &vm_bin1, Settings *settings1, double antenna_theta, double antenna_phi, double freq, + bool useInTransmitterMode=false + ); + void InvertAntFactors_Tdomain_FirstTwo( + double heff_copol, double heff_copol_lastbin, double heff_crosspol, double heff_crosspol_lastbin, + Vector &Pol_vector, int ant_type, double &pol_factor, double &vm_bin0, double &vm_bin1, Settings *settings1, + double antenna_theta, double antenna_phi, double freq, bool useInTransmitterMode=false ); void ApplyElect_Tdomain(double freq, Detector *detector, double &vm_real, double &vm_img, int gain_ch_no, Settings *settings1, bool applyInverse=false); diff --git a/Settings.cc b/Settings.cc index cacf7df7..52e3b7b8 100644 --- a/Settings.cc +++ b/Settings.cc @@ -325,7 +325,8 @@ outputdir="outputs"; // directory where outputs go // This setting is only used if DETECTOR > 3 && NOISE == 1 && CUSTOM_ELECTRONICS == 0 CLOCK_ANGLE=0; //Default: 0 -- Angle of polarization "on the clock". Angle of zero is pure thetaPol, whereas 90ยบ is pure phiPol. - + CROSSPOL_RX=0; //Default: 0 -- don't use cross-pol responses on receiving antennas + CROSSPOL_TX=0; //Default: 0 -- don't use cross-pol responses on transmitting antennas @@ -366,12 +367,12 @@ void Settings::ReadFile(string setupfile) { else if (label == "NOFZ") { NOFZ = atof( line.substr(line.find_first_of("=") + 1).c_str() ); } - else if (label == "BIREFRINGENCE"){ - BIREFRINGENCE = atof( line.substr(line.find_first_of("=") + 1).c_str() ); - } - else if (label == "BIAXIAL"){ - BIAXIAL = atof( line.substr(line.find_first_of("=") + 1).c_str() ); - } + else if (label == "BIREFRINGENCE"){ + BIREFRINGENCE = atof( line.substr(line.find_first_of("=") + 1).c_str() ); + } + else if (label == "BIAXIAL"){ + BIAXIAL = atof( line.substr(line.find_first_of("=") + 1).c_str() ); + } else if (label == "CONSTANTCRUST") { CONSTANTCRUST = atof( line.substr(line.find_first_of("=") + 1).c_str() ); } @@ -398,11 +399,11 @@ void Settings::ReadFile(string setupfile) { } else if (label == "DETECTOR_STATION") { DETECTOR_STATION = atof( line.substr(line.find_first_of("=") + 1).c_str() ); - DETECTOR_STATION_ARAROOT = DETECTOR_STATION; - if (DETECTOR_STATION == 100) { - DETECTOR_STATION = 1; - DETECTOR_STATION_ARAROOT = 100; - } + DETECTOR_STATION_ARAROOT = DETECTOR_STATION; + if (DETECTOR_STATION == 100) { + DETECTOR_STATION = 1; + DETECTOR_STATION_ARAROOT = 100; + } } else if (label == "DETECTOR_STATION_LIVETIME_CONFIG") { DETECTOR_STATION_LIVETIME_CONFIG = atof( line.substr(line.find_first_of("=") + 1).c_str() ); @@ -457,7 +458,7 @@ void Settings::ReadFile(string setupfile) { } else if(label == "TRIG_SCAN_MODE"){ TRIG_SCAN_MODE = atoi( line.substr(line.find_first_of("=") + 1).c_str() ); - } + } else if (label == "POWERTHRESHOLD") { POWERTHRESHOLD = atof( line.substr(line.find_first_of("=") + 1).c_str() ); } @@ -590,9 +591,6 @@ void Settings::ReadFile(string setupfile) { else if (label == "USE_CH_GAINOFFSET") { USE_CH_GAINOFFSET = atoi( line.substr(line.find_first_of("=") + 1).c_str() ); } - //else if (label == "GETCHORD_MODE") { - //GETCHORD_MODE = atoi( line.substr(line.find_first_of("=") + 1).c_str() ); - //} else if (label == "taumodes") { taumodes = atoi( line.substr(line.find_first_of("=") + 1).c_str() ); } @@ -689,42 +687,39 @@ void Settings::ReadFile(string setupfile) { else if (label == "AVZ_NORM_FACTOR_MODE") { AVZ_NORM_FACTOR_MODE = atoi( line.substr(line.find_first_of("=") + 1).c_str() ); } - else if (label == "number_of_stations") { - number_of_stations = atoi( line.substr(line.find_first_of("=") + 1).c_str() ); - } - else if (label == "RAY_TRACE_ICE_MODEL_PARAMS") { - RAY_TRACE_ICE_MODEL_PARAMS = atoi( line.substr(line.find_first_of("=") + 1).c_str() ); - } - else if (label == "WAVEFORM_LENGTH") { - WAVEFORM_LENGTH = atoi( line.substr(line.find_first_of("=") + 1).c_str() ); - } - else if (label == "WAVEFORM_CENTER") { - WAVEFORM_CENTER = atoi( line.substr(line.find_first_of("=") + 1).c_str() ); - } - else if (label == "DEADTIME") { - DEADTIME = atof( line.substr(line.find_first_of("=") + 1).c_str() ); - } - else if (label == "POSNU_R") { - POSNU_R = atof( line.substr(line.find_first_of("=") + 1).c_str() ); - } - else if (label == "POSNU_THETA") { - POSNU_THETA = atof( line.substr(line.find_first_of("=") + 1).c_str() ); - } - else if (label == "POSNU_PHI") { - POSNU_PHI = atof( line.substr(line.find_first_of("=") + 1).c_str() ); - } - else if (label == "ARBITRARY_EVENT_ATTENUATION") { - ARBITRARY_EVENT_ATTENUATION = atof( line.substr(line.find_first_of("=") + 1).c_str() ); - } - else if (label == "PICK_ABOVE_HEIGHT") { - PICK_ABOVE_HEIGHT = atof( line.substr(line.find_first_of("=") + 1).c_str() ); - } + else if (label == "number_of_stations") { + number_of_stations = atoi( line.substr(line.find_first_of("=") + 1).c_str() ); + } + else if (label == "RAY_TRACE_ICE_MODEL_PARAMS") { + RAY_TRACE_ICE_MODEL_PARAMS = atoi( line.substr(line.find_first_of("=") + 1).c_str() ); + } + else if (label == "WAVEFORM_LENGTH") { + WAVEFORM_LENGTH = atoi( line.substr(line.find_first_of("=") + 1).c_str() ); + } + else if (label == "WAVEFORM_CENTER") { + WAVEFORM_CENTER = atoi( line.substr(line.find_first_of("=") + 1).c_str() ); + } + else if (label == "DEADTIME") { + DEADTIME = atof( line.substr(line.find_first_of("=") + 1).c_str() ); + } + else if (label == "POSNU_R") { + POSNU_R = atof( line.substr(line.find_first_of("=") + 1).c_str() ); + } + else if (label == "POSNU_THETA") { + POSNU_THETA = atof( line.substr(line.find_first_of("=") + 1).c_str() ); + } + else if (label == "POSNU_PHI") { + POSNU_PHI = atof( line.substr(line.find_first_of("=") + 1).c_str() ); + } + else if (label == "ARBITRARY_EVENT_ATTENUATION") { + ARBITRARY_EVENT_ATTENUATION = atof( line.substr(line.find_first_of("=") + 1).c_str() ); + } + else if (label == "PICK_ABOVE_HEIGHT") { + PICK_ABOVE_HEIGHT = atof( line.substr(line.find_first_of("=") + 1).c_str() ); + } else if (label == "EVENT_GENERATION_MODE"){ EVENT_GENERATION_MODE = atoi(line.substr(line.find_first_of("=") + 1).c_str()); } - // else if (label == "EVENT_NUM"){ - // EVENT_NUM = atoi(line.substr(line.find_first_of("=") + 1).c_str()); - // } else if (label == "ANTENNA_MODE"){ ANTENNA_MODE = atoi(line.substr(line.find_first_of("=") + 1).c_str()); } @@ -749,22 +744,26 @@ void Settings::ReadFile(string setupfile) { else if (label == "ELECTRONICS_ANTENNA_CONSISTENCY"){ ELECTRONICS_ANTENNA_CONSISTENCY = atoi(line.substr(line.find_first_of("=") + 1).c_str()); } - else if (label == "CLOCK_ANGLE"){ - CLOCK_ANGLE = atof(line.substr(line.find_first_of("=") + 1).c_str()); - } - //Adding source easting, northing, and depth for INTERACTION_MODE=5. - else if (label == "SOURCE_LATITUDE"){ - SOURCE_LATITUDE = atof(line.substr(line.find_first_of("=") + 1).c_str()); - } - else if (label == "SOURCE_LONGITUDE"){ - SOURCE_LONGITUDE = atof(line.substr(line.find_first_of("=") + 1).c_str()); - } - else if (label == "SOURCE_DEPTH"){ - SOURCE_DEPTH = atof(line.substr(line.find_first_of("=") + 1).c_str()); - } - - - + else if (label == "CLOCK_ANGLE"){ + CLOCK_ANGLE = atof(line.substr(line.find_first_of("=") + 1).c_str()); + } + //Adding source easting, northing, and depth for INTERACTION_MODE=5. + else if (label == "SOURCE_LATITUDE"){ + SOURCE_LATITUDE = atof(line.substr(line.find_first_of("=") + 1).c_str()); + } + else if (label == "SOURCE_LONGITUDE"){ + SOURCE_LONGITUDE = atof(line.substr(line.find_first_of("=") + 1).c_str()); + } + else if (label == "SOURCE_DEPTH"){ + SOURCE_DEPTH = atof(line.substr(line.find_first_of("=") + 1).c_str()); + } + else if (label == "CROSSPOL_RX"){ + CROSSPOL_RX = atof(line.substr(line.find_first_of("=") + 1).c_str()); + } + else if (label == "CROSSPOL_TX"){ + CROSSPOL_TX = atof(line.substr(line.find_first_of("=") + 1).c_str()); + } + } } diff --git a/Settings.h b/Settings.h index cf46fec5..5e2ba028 100644 --- a/Settings.h +++ b/Settings.h @@ -347,7 +347,8 @@ class Settings // This setting is only used if DETECTOR > 3 && NOISE == 1 && CUSTOM_ELECTRONICS == 0 double CLOCK_ANGLE; // default: 0; Angle of polarization "on the clock" for use in pulser events (EVENT_TYPE=11) - + int CROSSPOL_RX=0; //Default: 0 -- don't use cross-pol responses on receiving antennas + int CROSSPOL_TX=0; //Default: 0 -- don't use cross-pol responses on transmitting antennas //arrays for saving read in event features in EVENT_GENERATION_MODE=1