diff --git a/Detector.cc b/Detector.cc index 7e509f80..641d6839 100644 --- a/Detector.cc +++ b/Detector.cc @@ -2256,11 +2256,16 @@ 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.gz"; + TxgainFileCross = string(getenv("ARA_SIM_DIR"))+"/data/antennas/realizedGain/PVA_RealizedGainAndPhase_Crosspol_Kansas2024.txt.gz"; + if (settings1->ANTENNA_MODE == 0){ // use the orignal Vpol/Hpol gains VgainFile = string(getenv("ARA_SIM_DIR"))+"/data/antennas/realizedGain/ARA_bicone6in_output.txt"; @@ -2309,31 +2314,42 @@ 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.gz"; + VgainTopCrossFile = string(getenv("ARA_SIM_DIR"))+"/data/antennas/realizedGain/ARA_TVpol_RealizedGainAndPhase_Crosspol_Kansas2024.txt.gz"; + HgainCrossFile = string(getenv("ARA_SIM_DIR"))+"/data/antennas/realizedGain/ARA_Hpol_RealizedGainAndPhase_Crosspol_Kansas2024.txt.gz"; + // Check for ALL_ANT_V_ON, then set all antennas to VPol if true if (settings1->ALL_ANT_V_ON == 1) { HgainFile = VgainFile; VgainTopFile = VgainFile; } - + // set parameters to "unset" values freq_step = -1; ang_step = -1; 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; params.freq_width = freq_width; params.freq_init = freq_init; - + return; } //Defining function that reads in TX antenna impedances @@ -2384,34 +2400,61 @@ 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; + source_n = &antenna_source_medium_n; + break; + case(eVPolTopCross): // Cross-pol VPol Top + freq = &FreqCross; + gain = &VgainTopCross; + phase = &VphaseTopCross; + transAnt_databin = &transVTopCross_databin; + source_n = &antenna_source_medium_n; + break; + case(eHPolCross): // Cross-pol HPol + freq = &FreqCross; + gain = &HgainCross; + phase = &HphaseCross; + transAnt_databin = &transHCross_databin; + source_n = &antenna_source_medium_n; + 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; + source_n = &Txantenna_source_medium_n; + break; + default: throw runtime_error("Unknown antenna type!"); } @@ -2460,10 +2503,9 @@ 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 == eTxCross) { freq->clear(); } // check the file opened successfully if (!NecOut ) { @@ -2494,7 +2536,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 || type == eTxCross ) { *source_n = buff_n; } else { // if not the first read-in, check for consistency @@ -2539,7 +2581,7 @@ inline void Detector::ReadAntennaGain(string filename, Settings *settings1, EAnt throw runtime_error("Non-finite frequency value found! "+filename); } - if(freq_step == -1 || type == eTx) { // add frequency if this is the first read-in + if(freq_step == -1 || type == eTx || type == eVPolCross || type == eTxCross) { // add frequency if this is the first read-in freq->push_back(thisFreq); } else { // otherwise make sure it matches the values already stored @@ -2621,36 +2663,50 @@ 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(); - ang_step = (int)gain->back().size(); - freq_width = freq->at(1)-freq->at(0); - freq_init = freq->at(0); + // Cross-pol Tx has its own frequency vector and does not need + // the Rx transmittance interpolation below + if(type == eTxCross) + { + return; } - // 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); + int this_freq_step = (int)freq->size(); + + // Only let the co-pol Rx models define the nominal frequency parameters + if(type == eVPol || type == eVPolTop || type == eHPol) { + + 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); + } } - 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(phase->at(i).size() != ang_step) { - throw runtime_error("gain vectors have inconsistent 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); + } + if(phase->at(i).size() != ang_step) { + throw runtime_error("gain vectors have inconsistent length! "+filename); + } } } // skip this for the transmitter case - if(type == eTx) { + if(type == eTx || type == eTxCross) { return; } @@ -3185,7 +3241,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) @@ -3218,59 +3274,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 + tempGain = 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; } } - double thisFreq_init; - double thisFreq_width; - if(useInTransmitterMode) { - F = &TxFreq; - thisFreq_init = Tx_freq_init; - thisFreq_width = Tx_freq_width; + if (!tempGain) { + throw runtime_error("In GetGain_1D_OutZero: No appropriate gain model for this simulation setup."); } + + if(useInTransmitterMode) { + F = useCrossPol ? &TxFreqCross : &TxFreq; + } else { - F = &Freq; - thisFreq_init = freq_init; - thisFreq_width = freq_width; - } + F = useCrossPol ? &FreqCross : &Freq; + } + + double thisFreq_init = F->at(0); + double thisFreq_width = F->at(1) - F->at(0); // check if angles range actually theta 0-180, phi 0-360 int i = (int)( (theta+2.5)/5. ); @@ -3371,7 +3425,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) { @@ -3392,38 +3446,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; - thisFreq_init = Tx_freq_init; - thisFreq_width = Tx_freq_width; + F = useCrossPol ? &TxFreqCross : &TxFreq; } else { - F = &Freq; - thisFreq_init = freq_init; - thisFreq_width = freq_width; + F = useCrossPol ? &FreqCross : &Freq; } - thisFreq_step = (int)F->size(); + + double thisFreq_init = F->at(0); + double thisFreq_width = F->at(1) - F->at(0); + int thisFreq_step = (int)F->size(); // check if angles range actually theta 0-180, phi 0-360 int i = (int)( (theta+2.5)/5. ); diff --git a/Detector.h b/Detector.h index 8a873681..0d8d9132 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 { @@ -236,6 +240,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 @@ -250,11 +263,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); @@ -351,6 +367,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; @@ -390,7 +411,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); @@ -400,7 +421,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/Primaries.cc b/Primaries.cc index 5bfb91e8..05a37573 100644 --- a/Primaries.cc +++ b/Primaries.cc @@ -2426,11 +2426,30 @@ void Interaction::PickExactGlobal(IceModel *antarctica, Detector *detector, Sett //Convert source vector into array coordinates TVector3 sourceStationVector = AraGeomTool::Instance()->convertArrayToStationCoords(stationId, sourceArrayVector); - + + // Calculate average antenna depth + double avgDepth = 0.; + int count = 0; + for (int chan = 0; chan < 16; chan++) { + double antDepth = AraGeomTool::Instance()->getStationInfo(settings1->DETECTOR_STATION_ARAROOT,settings1->DETECTOR_YEAR)->fAntInfo[chan].antLocation[2]; + avgDepth += antDepth; + count++; + } + + avgDepth /= double(count); + + + // SOURCE_DEPTH is referenced to the surface, + // while PickExact() later adds avgZ. + // Therefore shift only the Z coordinate. + double sourceX = sourceStationVector[0]; + double sourceY = sourceStationVector[1]; + double sourceZ = sourceStationVector[2] - avgDepth; + //Calculate posnu - double R = sqrt(pow(sourceStationVector[0],2) + pow(sourceStationVector[1],2) + pow(sourceStationVector[2],2)); - double phi = (360+(atan2(sourceStationVector[1], sourceStationVector[0]))*180/PI)*PI/180; - double theta = acos((sourceStationVector[2])/R); + double R = sqrt(pow(sourceX,2) + pow(sourceY,2) + pow(sourceZ,2)); + double phi = (360+(atan2(sourceY, sourceX))*180/PI)*PI/180; + double theta = acos((sourceZ)/R); settings1->POSNU_THETA = theta; settings1->POSNU_PHI = phi; settings1->POSNU_R = R; diff --git a/Report.cc b/Report.cc index 93e78a84..6b65827d 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,71 @@ void Report::PropagateSignal( // Get ant gain with 2-D interpolation double n_eff = icemodel->GetEffectiveN(*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, n_eff, - detector->GetImpedance(freq_tmp*1.E-6, antenna_d->type, k)); + 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, + n_eff, + 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 = icemodel->GetEffectiveN(*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, - n_eff, - detector->GetImpedance(freq_tmp*1.E-6, 0, 0, true)); + + // 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, + n_eff, + 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, + n_eff, + 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 = icemodel->GetEffectiveN(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, n_eff); + + // 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, false, false), + freq_tmp, n_eff); + + 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, false, true), + freq_tmp, n_eff); + // End Tx effective height for last bin if (event->IsCalpulser == 1) { Pol_vector = n_trg_slappy; @@ -1629,12 +1660,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 = icemodel->GetEffectiveN(*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, + n_eff, + 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, n_eff, - 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 +1708,106 @@ void Report::PropagateSignal( // Apply transmitter effective height if in PVA Pulser mode if (settings->EVENT_TYPE == 12){ double n_eff = icemodel->GetEffectiveN(*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, n_eff, - 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, + n_eff, + 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, true, 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, true, 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 = icemodel->GetEffectiveN(*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, n_eff); + + // 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, n_eff); + + // 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, false, false); + double phase_crosspol_tx = detector->GetAntPhase_1D(freq_tmp *1.e-6, Tx_theta, antenna_phi, antenna_d->type, n_eff, false, 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 @@ -3554,12 +3646,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. @@ -3586,17 +3696,38 @@ void Report::ApplyAntFactors_Tdomain ( */ double phaseSign = 1.; double amplitudeSign = 1.; + //If using in transmitter mode, the phase gets a minus sign since the signal is out-going from the antenna rather than incoming. if(useInTransmitterMode==true){ phaseSign*=-1.;}; - //If using this function to invert the antenna response, we apply a minus sign the phase in order to undo the phase applied by the antenna. We also divide the amplitude by the factor rather than multiply. To do this, we simply apply a -1 to the power of the factor applied to the amplitude so that it is divided out. + + //If using this function to invert the antenna response, we apply a minus sign the phase in order to undo the phase applied by the antenna. + //We also divide the amplitude by the factor rather than multiply. 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. ) { phase_current = atan( vm_img / vm_real ); + // phase in +-PI range if (vm_real<0.) { if (vm_img>0.) phase_current += PI; @@ -3608,33 +3739,114 @@ void Report::ApplyAntFactors_Tdomain ( else if (vm_img<0.) phase_current = -PI; 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 incoming voltage + //If in transmitter mode, we must apply additional frequency and impedance terms to the amplitude. - if (useInTransmitterMode==true){ + 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 = atan2(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; // turn 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( (0.25) * freq/CLIGHT * (Z0/Zr) * sqrt(heff_crosspol*heff_crosspol + heff_copol*heff_copol), // The 0.25 factor needs revision + amplitudeSign); + } + + // ========================================================= + // Transmitter + // ========================================================= + + if (useInTransmitterMode == true) { + + // Assume Tx co-pol and cross-pol have the same phase. + // Use co-pol phase as the common Tx antenna phase. + vm_real = v_amp * cos( phase_current + (phaseSign * phase_copol * RADDEG) ); + vm_img = v_amp * sin( phase_current + (phaseSign * phase_copol * RADDEG) ); + + } + + // ========================================================= + // Receiver + // ========================================================= + + else { + + // 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; } - //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) ); } else { // only amplitude - vm_real *= pow(heff * pol_factor, amplitudeSign); - vm_img *= pow(heff * pol_factor, amplitudeSign); + + // CHANGED: Tx needs the Tx amplitude conversion, not Rx projection + if (useInTransmitterMode == true) { + + if (settings1->CROSSPOL_TX == 0){ + heff_crosspol = 0.0; + } + + double heff_tx = sqrt( heff_crosspol*heff_crosspol + heff_copol*heff_copol ); + double tx_amplification = (0.25) * freq/CLIGHT * (Z0/Zr) * heff_tx; // The 0.25 factor needs revision + + vm_real *= pow(tx_amplification, amplitudeSign); + vm_img *= pow(tx_amplification, amplitudeSign); + } + + else { + + 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: @@ -3650,41 +3862,118 @@ 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); - + // 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; + } + + // 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 + if (!useInTransmitterMode) { + vm_bin0 *= pow(v_amplification_bin0, amplitudeSign); + vm_bin1 *= pow(v_amplification_bin1, 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); + + // 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 = (0.25) * freq / CLIGHT * (Z0 / Zr) * + sqrt(heff_crosspol * heff_crosspol + heff_copol * heff_copol); // The 0.25 factor needs revision + + double tx_amplification_bin1 = (0.25) * freq / CLIGHT * (Z0 / Zr) * + sqrt(heff_crosspol_lastbin * heff_crosspol_lastbin + + heff_copol_lastbin * heff_copol_lastbin); // The 0.25 factor needs revision + + vm_bin0 *= pow(tx_amplification_bin0, amplitudeSign); + vm_bin1 *= pow(tx_amplification_bin1, amplitudeSign); + } + + else { + + // Original co-pol-only transmitter behavior + vm_bin0 *= pow(v_amplification_copol_bin0, amplitudeSign); + vm_bin1 *= pow(v_amplification_copol_bin1, amplitudeSign); + + vm_bin0 *= pow((0.25) * freq / CLIGHT * (Z0 / Zr), amplitudeSign); // The 0.25 factor needs revision + vm_bin1 *= pow((0.25) * freq / CLIGHT * (Z0 / Zr), amplitudeSign); // The 0.25 factor needs revision + } } } -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 19e5d65a..5e3109df 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); @@ -335,20 +337,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 a10debfa..55da12b1 100644 --- a/Settings.cc +++ b/Settings.cc @@ -328,6 +328,9 @@ outputdir="outputs"; // directory where outputs go 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 + SYSTEMATICS_IceAttenuation=0; // 0=central (default), 1=up, 2=low SYSTEMATICS_AskaryanPercent=0.0; //Define the percentage (0% by default) @@ -406,11 +409,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() ); @@ -539,7 +542,6 @@ void Settings::ReadFile(string setupfile) { else if (label == "NNU_D_PHI") { NNU_D_PHI = atof( line.substr(line.find_first_of("=") + 1).c_str() ); } - else if (label == "Z_THIS_TOLERANCE") { Z_THIS_TOLERANCE = atof( line.substr(line.find_first_of("=") + 1).c_str() ); } @@ -764,7 +766,13 @@ void Settings::ReadFile(string setupfile) { 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()); + 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()); } else if (label == "SYSTEMATICS_IceAttenuation") { SYSTEMATICS_IceAttenuation = atoi(line.substr(line.find_first_of("=") + 1).c_str()); @@ -948,7 +956,6 @@ int Settings::CheckCompatibilitiesDetector(Detector *detector) { double heff1 = gain1 / freq1 / freq1; if(heff0 > heff1) { cerr << "Vpol antenna gain may not be falling fast enough! Effective height may grow and introduce power at low frequencies." << endl; - num_err++; } } @@ -962,7 +969,6 @@ int Settings::CheckCompatibilitiesDetector(Detector *detector) { double heff1 = gain1 / freq1 / freq1; if(heff0 > heff1) { cerr << "TVpol antenna gain may not be falling fast enough! Effective height may grow and introduce power at low frequencies." << endl; - num_err++; } } @@ -976,7 +982,6 @@ int Settings::CheckCompatibilitiesDetector(Detector *detector) { double heff1 = gain1 / freq1 / freq1; if(heff0 > heff1) { cerr << "Hpol antenna gain may not be falling fast enough! Effective height may grow and introduce power at low frequencies." << endl; - num_err++; } } } diff --git a/Settings.h b/Settings.h index 499ac14d..dba0d440 100644 --- a/Settings.h +++ b/Settings.h @@ -351,6 +351,9 @@ 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 int SYSTEMATICS_IceAttenuation; // 0=central, 1=up, 2=low diff --git a/data/antennas/realizedGain/ARA_BVpol_RealizedGainAndPhase_Crosspol_Kansas2024.txt b/data/antennas/realizedGain/ARA_BVpol_RealizedGainAndPhase_Crosspol_Kansas2024.txt index 88584774..2a7f86f4 100644 --- a/data/antennas/realizedGain/ARA_BVpol_RealizedGainAndPhase_Crosspol_Kansas2024.txt +++ b/data/antennas/realizedGain/ARA_BVpol_RealizedGainAndPhase_Crosspol_Kansas2024.txt @@ -1,3 +1,4 @@ +n : 1.0 freq : 83.33 MHz SWR : 11.85 Theta Phi Gain(dB) Gain Phase(deg) diff --git a/data/antennas/realizedGain/ARA_Hpol_RealizedGainAndPhase_Crosspol_Kansas2024.txt b/data/antennas/realizedGain/ARA_Hpol_RealizedGainAndPhase_Crosspol_Kansas2024.txt index 43ddc83a..36de658f 100644 --- a/data/antennas/realizedGain/ARA_Hpol_RealizedGainAndPhase_Crosspol_Kansas2024.txt +++ b/data/antennas/realizedGain/ARA_Hpol_RealizedGainAndPhase_Crosspol_Kansas2024.txt @@ -1,3 +1,4 @@ +n : 1.0 freq : 83.33 MHz SWR : 15.98 Theta Phi Gain(dB) Gain Phase(deg) diff --git a/data/antennas/realizedGain/ARA_TVpol_RealizedGainAndPhase_Crosspol_Kansas2024.txt b/data/antennas/realizedGain/ARA_TVpol_RealizedGainAndPhase_Crosspol_Kansas2024.txt index 3a63a272..7b22e08b 100644 --- a/data/antennas/realizedGain/ARA_TVpol_RealizedGainAndPhase_Crosspol_Kansas2024.txt +++ b/data/antennas/realizedGain/ARA_TVpol_RealizedGainAndPhase_Crosspol_Kansas2024.txt @@ -1,3 +1,4 @@ +n : 1.0 freq : 83.33 MHz SWR : -56.87 Theta Phi Gain(dB) Gain Phase(deg) diff --git a/data/antennas/realizedGain/PVA_RealizedGainAndPhase_Copol_Kansas2024.txt b/data/antennas/realizedGain/PVA_RealizedGainAndPhase_Copol_Kansas2024.txt index 720cacb3..47514170 100644 --- a/data/antennas/realizedGain/PVA_RealizedGainAndPhase_Copol_Kansas2024.txt +++ b/data/antennas/realizedGain/PVA_RealizedGainAndPhase_Copol_Kansas2024.txt @@ -1,3 +1,4 @@ +n : 1.0 freq : 83.33 MHz SWR : 8.64 Theta Phi Gain(dB) Gain Phase(deg) diff --git a/data/antennas/realizedGain/PVA_RealizedGainAndPhase_Crosspol_Kansas2024.txt b/data/antennas/realizedGain/PVA_RealizedGainAndPhase_Crosspol_Kansas2024.txt index 66137e31..ccb46fdc 100644 --- a/data/antennas/realizedGain/PVA_RealizedGainAndPhase_Crosspol_Kansas2024.txt +++ b/data/antennas/realizedGain/PVA_RealizedGainAndPhase_Crosspol_Kansas2024.txt @@ -1,3 +1,4 @@ +n : 1.0 freq : 83.33 MHz SWR : 8.64 Theta Phi Gain(dB) Gain Phase(deg)