diff --git a/.gitmodules b/.gitmodules index 8afb5d0d..a7b48c82 100644 --- a/.gitmodules +++ b/.gitmodules @@ -8,4 +8,4 @@ [submodule "kassiopeia"] path = kassiopeia url = https://github.com/project8/kassiopeia.git - branch = feature/eventmod + branch = feature/eventmod diff --git a/Config/Tutorial/LocustFakeTrack.json b/Config/Tutorial/LocustFakeTrack.json new file mode 100644 index 00000000..5d9ae412 --- /dev/null +++ b/Config/Tutorial/LocustFakeTrack.json @@ -0,0 +1,43 @@ +{ + "generators": + [ + "fake-track", + "lpf-fft", + "decimate-signal", + "gaussian-noise", + "digitizer" + ], + + "fake-track": + { + "signal-power": 1.0e-15, + "start-frequency": 20.05e9, + "start-vphase": 0.0, + "slope": 0.6, + "start-time": 0.001, + "end-time": 0.005, + "lo-frequency": 20.0e9 + }, + + "simulation": + { + "egg-filename": "/home/penny/locust_mc/cbuild/bin/locust_mc.egg", + "n-records": 1, + "n-channels": 1, + "record-size": 1228800 + }, + + "gaussian-noise": + { + "noise-floor": 4.0e-21, + "domain": "time" + }, + + "digitizer": + { + "v-range": 1.0e-5, + "v-offset": -0.5e-5 + } + +} + diff --git a/Config/Tutorial/LocustPhase1Template.json b/Config/Tutorial/LocustPhase1Template.json index aa4ffff3..f16172a3 100644 --- a/Config/Tutorial/LocustPhase1Template.json +++ b/Config/Tutorial/LocustPhase1Template.json @@ -36,14 +36,14 @@ "gaussian-noise": { - "noise-floor": 4.0e-24, + "noise-floor": 2.0e-22, "domain": "time" }, "digitizer": { - "v-range": 1.0e-6, - "v-offset": -0.5e-6 + "v-range": 4.0e-5, + "v-offset": -2.0e-5 } } diff --git a/Config/Tutorial/LocustPhase3Template.json b/Config/Tutorial/LocustPhase3Template.json index 88156b71..351d821a 100644 --- a/Config/Tutorial/LocustPhase3Template.json +++ b/Config/Tutorial/LocustPhase3Template.json @@ -16,7 +16,8 @@ "patch-signal": { - "lo-frequency": 25.8781e9, + "lo-frequency": 26.0681e9, + "feed": "corporate", "xml-filename": "/home/hep/baker/ps48/locust_mc/Config/Tutorial/Project8Phase3_WithRoot_Template.xml" }, @@ -30,13 +31,13 @@ { "egg-filename": "/home/hep/baker/ps48/locust_mc/cbuild/bin/locust_mc.egg", "n-records": 1, - "record-size": 4194304, - "n-channels": 3 + "record-size": 81920, + "n-channels": 30 }, "gaussian-noise": { - "noise-floor": 4.0e-24, + "noise-floor": 2.0e-22, "domain": "time" }, diff --git a/Config/Tutorial/Project8Phase1_WithRoot_Template.xml b/Config/Tutorial/Project8Phase1_WithRoot_Template.xml index d6bd69f7..98e4f222 100644 --- a/Config/Tutorial/Project8Phase1_WithRoot_Template.xml +++ b/Config/Tutorial/Project8Phase1_WithRoot_Template.xml @@ -3,7 +3,7 @@ - diff --git a/Config/Tutorial/Project8Phase2_WithRoot_Template.xml b/Config/Tutorial/Project8Phase2_WithRoot_Template.xml index d15d8f3a..07d945b9 100644 --- a/Config/Tutorial/Project8Phase2_WithRoot_Template.xml +++ b/Config/Tutorial/Project8Phase2_WithRoot_Template.xml @@ -433,7 +433,7 @@ - + diff --git a/Config/Tutorial/Project8Phase3_WithRoot_Template.xml b/Config/Tutorial/Project8Phase3_WithRoot_Template.xml index 8ff62791..3e66ffeb 100644 --- a/Config/Tutorial/Project8Phase3_WithRoot_Template.xml +++ b/Config/Tutorial/Project8Phase3_WithRoot_Template.xml @@ -4,7 +4,7 @@ - + @@ -24,14 +24,18 @@ - - + + - + + + + + @@ -201,16 +205,16 @@ - + - + - + diff --git a/Source/CMakeLists.txt b/Source/CMakeLists.txt index 4b390e69..f7dddba7 100644 --- a/Source/CMakeLists.txt +++ b/Source/CMakeLists.txt @@ -38,9 +38,11 @@ set( LOCUST_MC_HEADER_FILES Generators/LMCDigitizer.hh Generators/LMCGaussianNoiseGenerator.hh + Generators/LMCFakeTrackGenerator.hh Generators/LMCTestSignalGenerator.hh Generators/LMCDecimateSignalGenerator.hh Generators/LMCLowPassFilterFFTGenerator.hh + Generators/LMCLocalOscillatorGenerator.hh Generators/LMCButterworthLPFGenerator.hh # Generators/LMCBasebandTrackGenerator.hh @@ -63,9 +65,11 @@ set( LOCUST_MC_SOURCE_FILES Generators/LMCDigitizer.cc Generators/LMCGaussianNoiseGenerator.cc + Generators/LMCFakeTrackGenerator.cc Generators/LMCTestSignalGenerator.cc Generators/LMCDecimateSignalGenerator.cc Generators/LMCLowPassFilterFFTGenerator.cc + Generators/LMCLocalOscillatorGenerator.cc Generators/LMCButterworthLPFGenerator.cc # Generators/LMCBasebandTrackGenerator.cc @@ -79,9 +83,12 @@ if (locust_mc_BUILD_WITH_KASSIOPEIA) Generators/LMCPatchSignalGenerator.hh Generators/LMCKassSignalGenerator.hh + Core/LMCChannel.hh Core/LMCGlobalsDeclaration.hh Core/LMCGlobalsDefinition.hh + Core/LMCPatchAntenna.hh Core/LMCParticle.hh + Core/LMCReceiver.hh Core/LMCHFSSReader.hh Kassiopeia/LMCCyclotronRadiationExtractor.hh @@ -97,8 +104,11 @@ if (locust_mc_BUILD_WITH_KASSIOPEIA) Generators/LMCPatchSignalGenerator.cc Generators/LMCKassSignalGenerator.cc + Core/LMCChannel.cc Core/LMCHFSSReader.cc Core/LMCParticle.cc + Core/LMCPatchAntenna.cc + Core/LMCReceiver.cc Kassiopeia/LMCCyclotronRadiationExtractor.cc Kassiopeia/LMCCyclotronRadiationExtractorBuilder.cc diff --git a/Source/Core/LMCChannel.cc b/Source/Core/LMCChannel.cc new file mode 100644 index 00000000..43f62e14 --- /dev/null +++ b/Source/Core/LMCChannel.cc @@ -0,0 +1,25 @@ +/* + * LMCChannel.cc + * + * Created on: Mar 5, 2018 + * Author: nbuzinsky + */ + +#include "LMCReceiver.hh" +#include "LMCChannel.hh" + +#include + +namespace locust +{ + + //Use lambda to add all of GetVoltages() from all receiver elements + template< class T> + double Channel::PhasedSum() + { + return std::accumulate(receiverElements.begin(), receiverElements.end(), 0, + [](double s ,T r) {return s + r->GetVoltage();}); + } + +} /* namespace locust */ + diff --git a/Source/Core/LMCChannel.hh b/Source/Core/LMCChannel.hh new file mode 100644 index 00000000..52c2b862 --- /dev/null +++ b/Source/Core/LMCChannel.hh @@ -0,0 +1,41 @@ +/* + * LMCChannel.hh + * + * Created on: Mar 1, 2018 + * Author: nbuzinsky + */ + +#ifndef LMCCHANNEL_HH_ +#define LMCCHANNEL_HH_ +#include + +namespace locust +{ + /*! + @class Channel + @author N. Buzinsky + @brief Used to conveniently group multiple receivers into a channel/ and do phased sum + @details + Available configuration options: + No input parameters + */ + template + class Channel + { + public: + + double PhasedSum(); + void AddReceiver(const T receiver){receiverElements.push_back(receiver);} + + T& operator[] (const int &index) {return receiverElements[index];} + std::size_t size() {return receiverElements.size();} + + private: + std::vector receiverElements; //Best way to add receivers w/o object slicing? + +}; + + +} /* namespace locust */ + +#endif /* LMCCHANNEL_HH_ */ diff --git a/Source/Core/LMCGlobalsDeclaration.hh b/Source/Core/LMCGlobalsDeclaration.hh index 69401cea..ef4a81db 100644 --- a/Source/Core/LMCGlobalsDeclaration.hh +++ b/Source/Core/LMCGlobalsDeclaration.hh @@ -8,9 +8,9 @@ #ifndef GLOBALSDECLARATION_HH_ #define GLOBALSDECLARATION_HH_ -#define NPATCHES_PER_STRIP 11 -#define PATCH_SPACING 0.0108 // m. Spacing along Z. -#define PATCH_RADIUS 0.0516 // m. Radius of patch ring. +#define NPATCHES_PER_STRIP 21 +#define PATCH_SPACING 0.0054 // m. Spacing along Z. +#define PATCH_RADIUS 0.045 // m. Radius of patch ring. #define PATCH_RING_OFFSET 0.0 // m. Offset of entire array in X direction diff --git a/Source/Core/LMCPatchAntenna.cc b/Source/Core/LMCPatchAntenna.cc new file mode 100644 index 00000000..e73d42e8 --- /dev/null +++ b/Source/Core/LMCPatchAntenna.cc @@ -0,0 +1,118 @@ +/* + * LMCPatchAntenna.cc + * + * Created on: Mar 1, 2018 + * Author: nbuzinsky + */ + +#include "LMCPatchAntenna.hh" + + +namespace locust +{ + + PatchAntenna::PatchAntenna(): + antennaFactorSpline(antennaFactor.begin(), antennaFactor.end(), lowerBoundFrequency, frequencySpacingSpline ), + gainSpline(gain.begin(), gain.end(), lowerBoundAngle, angularSpacingSpline), + copolarizationDirection(0,0,0), + normalDirection(0,0,0), + centerPosition(0,0,0), + incidentElectricField(0,0,0), + incidentMagneticField(0,0,0), + instantaneousFrequency(0.), + previousRetardedIndex(-99), + previousRetardedTime(-99), + timeDelay(0.) + { + } + + PatchAntenna::~PatchAntenna() + { + } + LMCThreeVector PatchAntenna::GetPosition() + { + return centerPosition; + } + + void PatchAntenna::SetIncidentElectricField(const LMCThreeVector &incomingElectricField) + { + incidentElectricField = incomingElectricField; + } + + void PatchAntenna::SetIncidentMagneticField(const LMCThreeVector &incomingMagneticField) + { + incidentMagneticField = incomingMagneticField; + } + void PatchAntenna::SetInstantaneousFrequency(const double &dopplerFrequency) + { + instantaneousFrequency = dopplerFrequency; + } + + int PatchAntenna::GetPreviousRetardedIndex() + { + return previousRetardedIndex; + } + + double PatchAntenna::GetPreviousRetardedTime() + { + return previousRetardedTime; + } + + double PatchAntenna::GetAnalogTimeDelay() + { + return timeDelay; + } + + void PatchAntenna::SetPreviousRetardedIndex(const int& index) + { + previousRetardedIndex = index; + } + void PatchAntenna::SetPreviousRetardedTime(const double &time) + { + previousRetardedTime = time; + } + + + double PatchAntenna::GetVoltage() + { + return GetCopolarizationFactor() * GetGainFactor() / GetAntennaFactor(); + } + + double PatchAntenna::GetAntennaFactor() + { + return antennaFactorSpline(instantaneousFrequency); + } + + double PatchAntenna::GetGainFactor() + { + LMCThreeVector waveVector = incidentElectricField.Cross(incidentMagneticField); + waveVector = waveVector.Unit(); //Normalize + double incidentAngle = acos(waveVector.Dot(normalDirection)); + if(incidentAngle > LMCConst::Pi() / 2.) + incidentAngle=LMCConst::Pi() - incidentAngle; + + return gainSpline(0.) / gainSpline(incidentAngle); + } + + double PatchAntenna::GetCopolarizationFactor() + { + return incidentElectricField.Dot(copolarizationDirection); + } + + void PatchAntenna::SetCenterPosition(const LMCThreeVector &newPosition) + { + centerPosition = newPosition; + } + void PatchAntenna::SetPolarizationDirection(const LMCThreeVector &copolDirection) + { + copolarizationDirection = copolDirection; + } + + void PatchAntenna::SetNormalDirection(const LMCThreeVector &normDirection) + { + normalDirection = normDirection; + } + + + +} /* namespace locust */ diff --git a/Source/Core/LMCPatchAntenna.hh b/Source/Core/LMCPatchAntenna.hh new file mode 100644 index 00000000..ed33fb58 --- /dev/null +++ b/Source/Core/LMCPatchAntenna.hh @@ -0,0 +1,94 @@ +/* + * LMCPatchAntenna.hh + * + * Created on: Mar 1, 2018 + * Author: nbuzinsky + */ + +#ifndef LMCPATCHANTENNA_HH_ +#define LMCPATCHANTENNA_HH_ + +#include +#include "LMCThreeVector.hh" +#include "LMCConst.hh" +#include "LMCReceiver.hh" + +namespace locust +{ + /*! + @class PatchAntenna + @author N. Buzinsky + @brief Derived class describing the behaviours of the patch antenna voltage, as a function of the E field + @details + Available configuration options: + No input parameters + */ + class PatchAntenna: public Receiver + { + + public: + PatchAntenna(); + PatchAntenna(const LMCThreeVector &patchPosition); + PatchAntenna(const LMCThreeVector &patchPosition, const double &timeDelay); + virtual ~PatchAntenna(); + + virtual double GetVoltage(); + virtual double GetAnalogTimeDelay(); + + + void SetIncidentElectricField(const LMCThreeVector &incomingElectricField); + void SetIncidentMagneticField(const LMCThreeVector &incomingMagneticField); + void SetInstantaneousFrequency(const double &dopplerFrequency); + + LMCThreeVector GetPosition(); + + int GetPreviousRetardedIndex(); + double GetPreviousRetardedTime(); + + void SetPreviousRetardedIndex(const int& index); + void SetPreviousRetardedTime(const double &time); + + void SetCenterPosition(const LMCThreeVector &newPosition); + void SetPolarizationDirection(const LMCThreeVector &copolDirection); + void SetNormalDirection(const LMCThreeVector &normDirection); + + + private: + double GetAntennaFactor(); + double GetGainFactor(); + double GetCopolarizationFactor(); + + LMCThreeVector copolarizationDirection; + LMCThreeVector normalDirection; + + LMCThreeVector centerPosition; + + LMCThreeVector incidentElectricField; + LMCThreeVector incidentMagneticField; + + int previousRetardedIndex; + double previousRetardedTime; + + double instantaneousFrequency; + + double timeDelay; + + const std::vector antennaFactor = {577.051, 571.497, 565.992, 560.537, 555.136, 549.789, 544.499, 539.269, 534.100, 528.995, 523.957, 518.988, 514.091, 509.268, 504.523, 499.857, 495.275, 490.779, 486.372, 482.058, 477.839, 473.719, 469.702, 465.790, 461.987, 458.296, 454.722, 451.268, 447.936, 444.732, 441.658, 438.718, 435.916, 433.255, 430.739, 428.371, 426.155, 424.093, 422.189, 420.446, 418.868, 417.456, 416.214, 415.144, 414.248, 413.529, 412.988, 412.627, 412.448, 412.451, 412.638, 413.010, 413.568, 414.310, 415.238, 416.352, 417.650, 419.133, 420.798, 422.646, 424.675, 426.884, 429.270, 431.831, 434.567, 437.474, 440.550, 443.792, 447.199, 450.768, 454.495, 458.378, 462.413, 466.599, 470.932, 475.409, 480.027, 484.784, 489.675, 494.698, 499.850, 505.129, 510.530, 516.052, 521.691, 527.445, 533.311, 539.286, 545.368, 551.553, 557.841, 564.227, 570.710, 577.287, 583.956, 590.715, 597.561, 604.494, 611.509, 618.607}; + + const double lowerBoundFrequency = 25.1e9; + const double frequencySpacingSpline = 1.e9 /99.; + + const std::vector gain= { 413.2552, 413.4029, 413.7279, 414.2307, 414.9122, 415.7736, 416.8165, 418.0426, 419.4540, 421.0532, 422.8427, 424.8257, 427.0053, 429.3850, 431.9687, 434.7602, 437.7639, 440.9843, 444.4262, 448.0945, 451.9944, 456.1315, 460.5115, 465.1403, 470.0241, 475.1694, 480.5829, 486.2715, 492.2426, 498.5036, 505.0624, 511.9270, 519.1060, 526.6080, 534.4421, 542.6178, 551.1447, 560.0330, 569.2932, 578.9360, 588.9725, 599.4145, 610.2737, 621.5624, 633.2932, 645.4791, 658.1333, 671.2694, 684.9012, 699.0430, 713.7091, 728.9140, 744.6726, 760.9996, 777.9100, 795.4189, 813.5410, 832.2914, 851.6845, 871.7349, 892.4566, 913.8632, 935.9679, 958.7832, 982.3209, 1006.591, 1031.605, 1057.372, 1083.897, 1111.189, 1139.251, 1168.086, 1197.697, 1228.081, 1259.237, 1291.159, 1323.840, 1357.270, 1391.437, 1426.325, 1461.917, 1498.192, 1535.126, 1572.693, 1610.862, 1649.602, 1688.878, 1728.650, 1768.877, 1809.517, 1850.522}; + + const double lowerBoundAngle = 0; + const double angularSpacingSpline = 1. * LMCConst::Pi() / 180.; //Angular spacing in radians + + boost::math::cubic_b_spline antennaFactorSpline; + boost::math::cubic_b_spline gainSpline; + + }; + + +} /* namespace locust */ + +#endif /* LMCPATCHANTENNA_HH_ */ diff --git a/Source/Core/LMCReceiver.cc b/Source/Core/LMCReceiver.cc new file mode 100644 index 00000000..9c54e9e0 --- /dev/null +++ b/Source/Core/LMCReceiver.cc @@ -0,0 +1,16 @@ +/* + * LMCReceiver.cc + * + * Created on: Mar 5, 2018 + * Author: nbuzinsky + */ + +#include "LMCReceiver.hh" + +namespace locust +{ + Receiver::Receiver() {} + Receiver::~Receiver() {} + +} /* namespace locust */ + diff --git a/Source/Core/LMCReceiver.hh b/Source/Core/LMCReceiver.hh new file mode 100644 index 00000000..636e5ea1 --- /dev/null +++ b/Source/Core/LMCReceiver.hh @@ -0,0 +1,37 @@ +/* + * LMCReceiver.hh + * + * Created on: Mar 1, 2018 + * Author: nbuzinsky + */ + +#ifndef LMCRECEIVER_HH_ +#define LMCRECEIVER_HH_ + +#include "LMCThreeVector.hh" + +namespace locust +{ + /*! + @class Receiver + @author N. Buzinsky + @brief Base class to characterize receiver elements (patches waveguides, etc.) + @details + Available configuration options: + No input parameters + */ + class Receiver + { + + public: + Receiver(); + virtual ~Receiver(); + + virtual double GetVoltage() = 0; + virtual double GetAnalogTimeDelay() = 0; +}; + + +} /* namespace locust */ + +#endif /* LMCCHANNEL_HH_ */ diff --git a/Source/Core/LMCRunLengthCalculator.cc b/Source/Core/LMCRunLengthCalculator.cc index 1a2e572c..8bdd403b 100644 --- a/Source/Core/LMCRunLengthCalculator.cc +++ b/Source/Core/LMCRunLengthCalculator.cc @@ -114,6 +114,12 @@ namespace locust return; } + void RunLengthCalculator::Visit( const LocalOscillatorGenerator* ) + { + // nothing to see here, move along, please + return; + } + void RunLengthCalculator::Visit( const GaussianNoiseGenerator* ) { // nothing to see here, move along, please @@ -132,6 +138,11 @@ namespace locust } + void RunLengthCalculator::Visit( const FakeTrackGenerator* ) + { + // nothing to see here, move along, please + return; + } void RunLengthCalculator::Visit( const TestSignalGenerator* ) { diff --git a/Source/Core/LMCRunLengthCalculator.hh b/Source/Core/LMCRunLengthCalculator.hh index 131d2de2..80655f4b 100644 --- a/Source/Core/LMCRunLengthCalculator.hh +++ b/Source/Core/LMCRunLengthCalculator.hh @@ -96,9 +96,11 @@ namespace locust void Visit( const PatchSignalGenerator* ); void Visit( const TrappedElectronGenerator* ); void Visit( const GaussianNoiseGenerator* ); + void Visit( const FakeTrackGenerator* ); void Visit( const TestSignalGenerator* ); void Visit( const BasebandTrackGenerator* ); void Visit( const LowPassFilterFFTGenerator* ); + void Visit( const LocalOscillatorGenerator* ); void Visit( const DecimateSignalGenerator* ); void Visit( const Digitizer* ); diff --git a/Source/Core/LMCSimulationController.cc b/Source/Core/LMCSimulationController.cc index f60bcb5b..a7919c20 100644 --- a/Source/Core/LMCSimulationController.cc +++ b/Source/Core/LMCSimulationController.cc @@ -19,6 +19,7 @@ namespace locust SimulationController::SimulationController() : fFirstGenerator( NULL ), fRunLengthCalc(), + fEggWriter() { SetRNGSeed(); diff --git a/Source/Core/LMCVisitor.hh b/Source/Core/LMCVisitor.hh index 0c8e8d14..326c0752 100644 --- a/Source/Core/LMCVisitor.hh +++ b/Source/Core/LMCVisitor.hh @@ -12,6 +12,7 @@ namespace locust { class Generator; class GaussianNoiseGenerator; + class FakeTrackGenerator; class TestSignalGenerator; class BasebandTrackGenerator; class LowPassFilterFFTGenerator; @@ -20,6 +21,7 @@ namespace locust class KassSignalGenerator; class FreeFieldSignalGenerator; class PatchSignalGenerator; + class LocalOscillatorGenerator; class DecimateSignalGenerator; class GeneratorVisitor @@ -35,9 +37,11 @@ namespace locust virtual void Visit( const PatchSignalGenerator* ) = 0; virtual void Visit( const GaussianNoiseGenerator* ) = 0; virtual void Visit( const TrappedElectronGenerator* ) = 0; + virtual void Visit( const FakeTrackGenerator* ) = 0; virtual void Visit( const TestSignalGenerator* ) = 0; virtual void Visit( const BasebandTrackGenerator* ) = 0; virtual void Visit( const LowPassFilterFFTGenerator* ) = 0; + virtual void Visit( const LocalOscillatorGenerator* ) = 0; virtual void Visit( const DecimateSignalGenerator* ) = 0; virtual void Visit( const Digitizer* ) = 0; }; diff --git a/Source/Generators/LMCDigitizer.cc b/Source/Generators/LMCDigitizer.cc index 328690a7..e6eaa6ab 100644 --- a/Source/Generators/LMCDigitizer.cc +++ b/Source/Generators/LMCDigitizer.cc @@ -74,47 +74,22 @@ namespace locust bool IQStream = true; // get rid of this. unsigned nchannels = fNChannels; - printf("digi says nchannels is %d\n", nchannels); unsigned signalSize = aSignal->TimeSize(); unsigned signalSizeComplex = 2*aSignal->TimeSize()*nchannels; double* analogData = aSignal->SignalTime(); // uint64_t* digitizedData = new uint64_t[ signalSize ]; -/* - std::complex* analogDataComplex; - analogDataComplex = (std::complex *) malloc(nchannels*signalSize * sizeof(std::complex)); - memcpy( analogDataComplex, aSignal->SignalTimeComplex(), nchannels*signalSize * sizeof( std::complex ) ); -*/ - if( fADCValuesSigned ) { - if (!IQStream) - { - int8_t* digitizedData = new int8_t[ signalSize ]; - for( unsigned index = 0; index < signalSize; ++index ) - { - digitizedData[ index ] = a2d< double, int8_t >( analogData[ index ], &fParams ); - if( index < 10 ) - { - LWARN( lmclog, "digitizing: " << index << ": " << analogData[ index ] << " --> " << (int) digitizedData[ index ] ); // pls added (int) - } - } - aSignal->ToDigital( digitizedData, signalSize ); - } // !IQStream - else - { int8_t* digitizedData = new int8_t[ signalSizeComplex ]; -// printf("signalSizeComplex is %d\n", signalSizeComplex); getchar(); for (unsigned ch = 0; ch < nchannels; ++ch) { for( unsigned index = 0; index < signalSize; ++index ) { -// printf("2*ch*signalSize+index*2 is %d\n", 2*ch*signalSize+index*2); -// digitizedData[2*ch*signalSize + index*2 ] = a2d< double, int8_t >( analogDataComplex[ch*signalSize + index ].real(), &fParams ); -// digitizedData[2*ch*signalSize + index*2+1 ] = a2d< double, int8_t >( analogDataComplex[ch*signalSize + index ].imag(), &fParams ); + digitizedData[2*ch*signalSize + index*2 ] = a2d< double, int8_t >( aSignal->SignalTimeComplex()[ch*signalSize + index ][0], &fParams ); digitizedData[2*ch*signalSize + index*2+1 ] = a2d< double, int8_t >( aSignal->SignalTimeComplex()[ch*signalSize + index ][1], &fParams ); @@ -124,60 +99,40 @@ namespace locust LWARN( lmclog, "digitizing channel " << ch << ": " << index << " I: " << aSignal->SignalTimeComplex()[ch*signalSize + index ][0] << " --> " << (int) digitizedData[2*ch*signalSize + index*2 ] ); // pls added (int) LWARN( lmclog, "digitizing channel " << ch << ": " << index << " Q: " << aSignal->SignalTimeComplex()[ch*signalSize + index ][1] << " --> " << (int) digitizedData[2*ch*signalSize + index*2+1 ] ); // pls added (int) } - } - } + } // signalsize + } // channels aSignal->ToDigital( digitizedData, signalSizeComplex ); - } // IQStream - } - else + } // fADCValuesSigned + else // unsigned { - if (!IQStream) - { - uint8_t* digitizedData = new uint8_t[ signalSize ]; - - for( unsigned index = 0; index < signalSize; ++index ) - { - digitizedData[ index ] = a2d< double, uint8_t >( analogData[ index ], &fParams ); - if( index < 10 ) - { - LWARN( lmclog, "digitizing: " << index << ": " << analogData[ index ] << " --> " << (unsigned) digitizedData[ index ] ); // pls added (int) - // printf("digitized data is %x\n", digitizedData[index]); - // getchar(); - } - } - aSignal->ToDigital( digitizedData, signalSize ); - } // !IQStream - else - { + uint8_t* digitizedData = new uint8_t[ signalSizeComplex ]; -// printf("signalSizeComplex is %d\n", signalSizeComplex); getchar(); + // FILE *fp = fopen("/home/hep/baker/ps48/data/Simulation/Phase3/timedata.txt", "w"); // write raw time series for testing. + for (unsigned ch = 0; ch < nchannels; ++ch) { + // fprintf(fp, "[%d]\n", ch); for( unsigned index = 0; index < signalSize; ++index ) { -// printf("ch is %d, index is %d, 2*ch*signalSize+index*2 is %d\n", ch, index, 2*ch*signalSize+index*2); -// digitizedData[2*ch*signalSize + index*2 ] = a2d< double, uint8_t >( analogDataComplex[ch*signalSize + index ].real(), &fParams ); -// digitizedData[2*ch*signalSize + index*2+1 ] = a2d< double, uint8_t >( analogDataComplex[ch*signalSize + index ].imag(), &fParams ); + + // fprintf(fp, "%g %g\n", aSignal->SignalTimeComplex()[ch*signalSize + index ][0], aSignal->SignalTimeComplex()[ch*signalSize + index ][1]); + digitizedData[2*ch*signalSize + index*2 ] = a2d< double, uint8_t >( aSignal->SignalTimeComplex()[ch*signalSize + index ][0], &fParams ); digitizedData[2*ch*signalSize + index*2+1 ] = a2d< double, uint8_t >( aSignal->SignalTimeComplex()[ch*signalSize + index ][1], &fParams ); - // fake data for debugging. -// digitizedData[2*ch*signalSize + index*2 ] = a2d< double, uint8_t >( 5.e-8, &fParams ); -// digitizedData[2*ch*signalSize + index*2+1 ] = a2d< double, uint8_t >( 5.e-8, &fParams ); - if( index < 20 ) { LWARN( lmclog, "digitizing channel " << ch << ": " << index << " I: " << aSignal->SignalTimeComplex()[ch*signalSize + index ][0] << " --> " << (int) digitizedData[2*ch*signalSize + index*2 ] ); // pls added (int) LWARN( lmclog, "digitizing channel " << ch << ": " << index << " Q: " << aSignal->SignalTimeComplex()[ch*signalSize + index ][1] << " --> " << (int) digitizedData[2*ch*signalSize + index*2+1 ] ); // pls added (int) } - } - } + } // signalsize + } // channels + // fclose(fp); aSignal->ToDigital( digitizedData, signalSizeComplex ); - } // IQStream - } + } // unsigned return true; } diff --git a/Source/Generators/LMCFakeTrackGenerator.cc b/Source/Generators/LMCFakeTrackGenerator.cc new file mode 100644 index 00000000..5e11393d --- /dev/null +++ b/Source/Generators/LMCFakeTrackGenerator.cc @@ -0,0 +1,254 @@ +/* + * LMCFakeTrackGenerator.cc + * + * Created on: Aug 8 2018 + * Author: plslocum + */ + +#include +#include "LMCFakeTrackGenerator.hh" +#include "LMCDigitizer.hh" +#include "logger.hh" +#include "LMCConst.hh" + + +using std::string; + +namespace locust +{ + LOGGER( lmclog, "FakeTrackGenerator" ); + + MT_REGISTER_GENERATOR(FakeTrackGenerator, "fake-track"); + + FakeTrackGenerator::FakeTrackGenerator( const std::string& aName ) : + Generator( aName ), + fDoGenerateFunc( &FakeTrackGenerator::DoGenerateTime ), + fSignalPower( 0. ), + fStartFrequency( 0. ), + fStartVPhase( 0. ), + fSlope( 0. ), + fStartTime( 0. ), + fEndTime( 0. ), + fLO_frequency( 0. ) + { + fRequiredSignalState = Signal::kTime; + } + + FakeTrackGenerator::~FakeTrackGenerator() + { + } + + + bool FakeTrackGenerator::Configure( const scarab::param_node* aParam ) + { + if( aParam == NULL) return true; + + if( aParam->has( "signal-power" ) ) + SetSignalPower( aParam->get_value< double >( "signal-power", fSignalPower ) ); + + if( aParam->has( "start-frequency" ) ) + SetStartFrequency( aParam->get_value< double >( "start-frequency", fStartFrequency ) ); + + if( aParam->has( "start-vphase" ) ) + SetStartVPhase( aParam->get_value< double >( "start-vphase", fStartVPhase ) ); + + if( aParam->has( "slope" ) ) + SetSlope( aParam->get_value< double >( "slope", fSlope ) ); + + if( aParam->has( "start-time" ) ) + SetStartTime( aParam->get_value< double >( "start-time", fStartTime ) ); + + if( aParam->has( "end-time" ) ) + SetEndTime( aParam->get_value< double >( "end-time", fEndTime ) ); + + if( aParam->has( "lo-frequency" ) ) + SetFrequency( aParam->get_value< double >( "lo-frequency", fLO_frequency ) ); + + + + if( aParam->has( "domain" ) ) + { + string domain = aParam->get_value( "domain" ); + if( domain == "time" ) + { + SetDomain( Signal::kTime ); + LDEBUG( lmclog, "Domain is equal to time."); + } + else if( domain == "freq" ) + { + SetDomain( Signal::kFreq ); + } + else + { + LERROR( lmclog, "Unable to use domain requested: <" << domain << ">" ); + return false; + } + } + + + return true; + + + + } + + + void FakeTrackGenerator::Accept( GeneratorVisitor* aVisitor ) const + { + aVisitor->Visit( this ); + return; + } + + double FakeTrackGenerator::GetSignalPower() const + { + return fSignalPower; + } + + void FakeTrackGenerator::SetSignalPower( double aPower ) + { + fSignalPower = aPower; + return; + } + + double FakeTrackGenerator::GetStartFrequency() const + { + return fStartFrequency; + } + + void FakeTrackGenerator::SetStartFrequency( double aFrequency ) + { + fStartFrequency = aFrequency; + return; + } + + double FakeTrackGenerator::GetStartVPhase() const + { + return fStartVPhase; + } + + void FakeTrackGenerator::SetStartVPhase( double aPhase ) + { + fStartVPhase = aPhase; + return; + } + + double FakeTrackGenerator::GetSlope() const + { + return fSlope; + } + + void FakeTrackGenerator::SetSlope( double aSlope ) + { + fSlope = aSlope; + return; + } + + double FakeTrackGenerator::GetStartTime() const + { + return fStartTime; + } + + void FakeTrackGenerator::SetStartTime( double aTime ) + { + fStartTime = aTime; + return; + } + + double FakeTrackGenerator::GetEndTime() const + { + return fEndTime; + } + + void FakeTrackGenerator::SetEndTime( double aTime ) + { + fEndTime = aTime; + return; + } + + double FakeTrackGenerator::GetFrequency() const + { + return fLO_frequency; + } + + void FakeTrackGenerator::SetFrequency( double aFrequency ) + { + fLO_frequency = aFrequency; + return; + } + + + + + Signal::State FakeTrackGenerator::GetDomain() const + { + return fRequiredSignalState; + } + + void FakeTrackGenerator::SetDomain( Signal::State aDomain ) + { + if( aDomain == fRequiredSignalState ) return; + fRequiredSignalState = aDomain; // pls changed == to =. + if( fRequiredSignalState == Signal::kTime ) + { + fDoGenerateFunc = &FakeTrackGenerator::DoGenerateTime; + } + else if( fRequiredSignalState == Signal::kFreq ) + { + fDoGenerateFunc = &FakeTrackGenerator::DoGenerateFreq; + } + else + { + LWARN( lmclog, "Unknown domain requested: " << aDomain ); + } + return; + } + + + bool FakeTrackGenerator::DoGenerate( Signal* aSignal ) + { + return (this->*fDoGenerateFunc)( aSignal ); + } + + bool FakeTrackGenerator::DoGenerateTime( Signal* aSignal ) + { + + RunLengthCalculator *RunLengthCalculator1 = new RunLengthCalculator; + + const unsigned nchannels = fNChannels; + double LO_phase = 0.; + double dt = 1./aSignal->DecimationFactor()/(RunLengthCalculator1->GetAcquisitionRate()*1.e6); + + for (unsigned ch = 0; ch < nchannels; ++ch) + { + double voltage_phase = fStartVPhase; + double track_frequency = fStartFrequency; + + for( unsigned index = 0; index < aSignal->TimeSize()*aSignal->DecimationFactor(); ++index ) + + { + double time = (double)index/aSignal->DecimationFactor()/(RunLengthCalculator1->GetAcquisitionRate()*1.e6); + + LO_phase += 2.*LMCConst::Pi()*fLO_frequency*dt; + + if ((time > fStartTime) && (time < fEndTime)) + { + track_frequency += fSlope*1.e6/1.e-3*dt; + voltage_phase += 2.*LMCConst::Pi()*track_frequency*(dt); + + + aSignal->LongSignalTimeComplex()[ch*aSignal->TimeSize()*aSignal->DecimationFactor() + index][0] += sqrt(50.)*sqrt(fSignalPower)*cos(voltage_phase-LO_phase); + aSignal->LongSignalTimeComplex()[ch*aSignal->TimeSize()*aSignal->DecimationFactor() + index][1] += sqrt(50.)*sqrt(fSignalPower)*cos(-LMCConst::Pi()/2. + voltage_phase-LO_phase); + } + + } + } + delete RunLengthCalculator1; + return true; + } + + bool FakeTrackGenerator::DoGenerateFreq( Signal* aSignal ) + { + return true; + } + +} /* namespace locust */ diff --git a/Source/Generators/LMCFakeTrackGenerator.hh b/Source/Generators/LMCFakeTrackGenerator.hh new file mode 100644 index 00000000..aed546d7 --- /dev/null +++ b/Source/Generators/LMCFakeTrackGenerator.hh @@ -0,0 +1,105 @@ +/* + * LMCFakeTrackGenerator.hh + * + * Created on: Aug 8 2018 + * Author: plslocum + */ + +#ifndef LMCFAKETRACKGENERATOR_HH_ +#define LMCFAKETRACKGENERATOR_HH_ + +#include "LMCGenerator.hh" +#include "LMCRunLengthCalculator.hh" +#include + + +namespace scarab +{ + class param_node; +} + +namespace locust +{ + class Digitizer; + + /*! + @class FakeTrackGenerator + @author P. L. Slocum + + @brief Add Sine Wave to the signal. + + @details + Operates in time. + + Configuration name: "fake-track" + + Available configuration options: + - "frequency": double -- Frequency of the sine wave. + - "amplitude": double -- Amplitude of the sine wave. + - "domain": string -- Determines whether the sinusoidal test signal is generated in the time + or frequency domain + + + */ + class FakeTrackGenerator : public Generator + { + public: + FakeTrackGenerator( const std::string& aName = "fake-track" ); + virtual ~FakeTrackGenerator(); + + bool Configure( const scarab::param_node* aNode ); + bool Configure2( const Digitizer* aDig ); + + + void Accept( GeneratorVisitor* aVisitor ) const; + + double GetFrequency() const; + void SetFrequency( double aFrequency ); + + double GetSignalPower() const; + void SetSignalPower( double aPower ); + + double GetStartFrequency() const; + void SetStartFrequency( double aFrequency ); + + double GetStartVPhase() const; + void SetStartVPhase( double aPhase ); + + double GetSlope() const; + void SetSlope( double aSlope ); + + double GetStartTime() const; + void SetStartTime( double aTime ); + + double GetEndTime() const; + void SetEndTime( double aTime ); + + + Signal::State GetDomain() const; + void SetDomain( Signal::State aDomain ); + + + private: + bool DoGenerate( Signal* aSignal ); + + bool DoGenerateTime( Signal* aSignal ); + bool DoGenerateFreq( Signal* aSignal ); + + bool (FakeTrackGenerator::*fDoGenerateFunc)( Signal* aSignal ); + + double fSignalPower; + double fStartFrequency; + double fStartVPhase; + double fSlope; + double fStartTime; + double fEndTime; + double fLO_frequency; + + + + }; + +} /* namespace locust */ + +#endif /* LMCFakeTrackGENERATOR_HH_ */ + diff --git a/Source/Generators/LMCFreeFieldSignalGenerator.cc b/Source/Generators/LMCFreeFieldSignalGenerator.cc index d1e2a7a3..246376c2 100644 --- a/Source/Generators/LMCFreeFieldSignalGenerator.cc +++ b/Source/Generators/LMCFreeFieldSignalGenerator.cc @@ -1,8 +1,8 @@ /* * LMCFreeFieldSignalGenerator.cc * - * Created on: Mar 12, 2014 - * Author: nsoblath + * Created on: Mar 12, 2017 + * Author: buzinsky */ #include "LMCFreeFieldSignalGenerator.hh" @@ -27,7 +27,7 @@ namespace locust FreeFieldSignalGenerator::FreeFieldSignalGenerator( const std::string& aName ) : Generator( aName ), - fWriteNFD(0.), + //fWriteNFD(0.), fLO_Frequency( 0.), gxml_filename("blank.xml") { @@ -42,34 +42,11 @@ namespace locust { if( aParam == NULL) return true; - if( aParam->has( "lo-frequency" ) ) - { - fLO_Frequency = aParam->get_value< double >( "lo-frequency" ); - } if( aParam->has( "xml-filename" ) ) { gxml_filename = aParam->get_value< std::string >( "xml-filename" ); } - if( aParam->has( "and-filename" ) ) - { - fWriteNFD = true; - fAND_filename = aParam->get_value< std::string >( "and-filename" ); - HFSSReader HFRead; - HFRead.ParseANDFile(fAND_filename); - rReceiver=HFRead.GetSurfacePoints(); - NFDFrequencies=HFRead.GetFrequencies(); - fNFD_filename=HFRead.GetNFDFilename(); - } - else - { - //If not, use existing code to generate plane receiver - HFSSReader HFRead; - rReceiver = HFRead.GeneratePlane({0.01,0.01},7);//Argumemts: Size, resolution - rReceiver = HFRead.RotateShift(rReceiver,{1.,0.,0.},{0.05,0.,0.});//Arguments Normal vector, Position (m) - } - PreviousTimes = std::vector >(rReceiver.size(),{-99.,-99.}); - return true; } @@ -89,18 +66,17 @@ namespace locust return 0; } - - static void WakeBeforeEvent() { fPreEventCondition.notify_one(); return; } - + static bool ReceivedKassReady() { printf("LMC about to wait ..\n"); + std::unique_lock< std::mutex >tLock( fKassReadyMutex); fKassReadyCondition.wait( tLock, [](){return fKassEventReady;} ); printf("LMC Got the fKassReadyCondition signal\n"); @@ -108,131 +84,15 @@ namespace locust return true; } - double m(double angle, double x0, double y0) - { - angle = LMCConst::Pi() * angle / 180.; - double xprime = cos(angle)*x0 - sin(angle)*y0; - double yprime = sin(angle)*x0 + cos(angle)*y0; - - double m = yprime/xprime; - - return m; - } - - double b(double m, double x0, double y0) - { - double b = y0 - m*x0; - return b; - } - - double directivity(double m1, double m2, double x0, double y0, double x, double y) - { - // cone shaped directivity with gain=1 and half angle OpeningAngle for antenna at x0,y0 - - double directivity = 0.; - - if (((y < m1*x + b(m1,x0,y0)) && (y > m2*x + b(m2,x0,y0))) | - ((y > m1*x + b(m1,x0,y0)) && (y < m2*x + b(m2,x0,y0)))) - { - if (fabs(atan(m1)-atan(m2)) < LMCConst::Pi()/2 ) - directivity = 1.; - } - else - { - if (fabs(atan(m1)-atan(m2)) > LMCConst::Pi()/2 ) - directivity = 1.; - } - return directivity; - - } - - void FreeFieldSignalGenerator::NFDWrite() - { - std::ofstream fNFDOutput; - fNFDOutput.open(fNFD_filename,std::ios::out | std::ios::trunc); - - fNFDOutput << "#Index, X, Y, Z, Ex(real, imag), Ey(real, imag), Ez(real, imag), Hx(real, imag), Hy(real, imag), Hz(real, imag)" < 0) ? tSign = 1. : tSign = -1.; - - return tRetardedTime + aSpaceTimeInterval; - //return tRetardedTime + tSign * aSpaceTimeInterval; - } - - LMCThreeVector tNewPosition = aParticle.GetPosition(true); - LMCThreeVector tNewVelocity = aParticle.GetVelocity(true); - - LMCThreeVector tReceiverVector = aReceiverPosition - tNewPosition; - double tReceiverDistance = tReceiverVector.Magnitude(); - - //Newtons Method X_{n+1} = X_{n} - f(X_{n}) / f'(X_{n}) - double fZero=pow((aReceiverTime - tRetardedTime),2.)-pow(tReceiverDistance,2.)/(c*c); - double fZeroPrime=2.*((tRetardedTime-aReceiverTime)-tNewVelocity.Dot(tNewPosition)/(c*c)+tNewVelocity.Dot(aReceiverPosition)/(c*c)); - double tNewtonRatio = fZero / fZeroPrime; - - if(aStepOrder==1) - { - return tRetardedTime-tNewtonRatio; - } - - //Householders Method - LMCThreeVector tNewAcceleration = aParticle.GetAcceleration(true); - double fZeroDoublePrime = 2. * (1. - tNewVelocity.Dot(tNewVelocity)/(c*c)-tNewAcceleration.Dot(tNewPosition-aReceiverPosition)/(c*c)); - - if(aStepOrder==2) - { - return tRetardedTime-tNewtonRatio* ( 1. + ( tNewtonRatio * fZeroDoublePrime) / ( 2. * fZeroPrime)); - } - - LERROR( lmclog, "Need to put root finding method with order 0-2!" ); - return 0; + double tRetardedTime = aParticle.GetTime(true); + return tRetardedTime + aSpaceTimeInterval; } @@ -241,208 +101,99 @@ namespace locust locust::Particle tCurrentParticle = fParticleHistory.back(); int CurrentIndex; - HFSSReader HFRead; - - - const unsigned nchannels = 1; - const double dx=0.01; - - const int nelementsZ = 9; - - double theta = 0.; // azimuthal angle of each amplifier channel. - double radius = 0.05; // radius of each amplifier channel in xy plane. - double OpeningAngle = 15.; //degrees, half angle fake directivity angle. - double DirectivityFactor = 0.; - - double m1 = 0.; - double m2 = 0.; - static double tVoltagePhase[10000] = {0.}; - static double phi_LO = 0.; - - static int fNFDIndex =- 1; - const int signalSize = aSignal->TimeSize(); + const double kassiopeiaTimeStep = fabs(fParticleHistory[0].GetTime() - fParticleHistory[1].GetTime()); + const int historySize = fParticleHistory.size(); + //Receiver Properties double tReceiverTime = t_old; - LMCThreeVector tReceiverPosition; double tRetardedTime = 0.; //Retarded time of particle corresponding to when emission occurs, reaching receiver at tReceiverTime - double tTotalPower=0.; double tSpaceTimeInterval=99.; double dtRetarded=0; double tTolerance=1e-23; - const double dtStepSize = fabs(fParticleHistory[0].GetTime() - fParticleHistory[1].GetTime()); - const int HistorySize = fParticleHistory.size(); + PatchAntenna *currentPatch; - phi_LO+= 2. * LMCConst::Pi() * fLO_Frequency * fDigitizerTimeStep; // this has to happen outside the signal generating loop. - - //int tAverageIterations=0; //Performance tracker. Count number of iterations to converge.... - - for (int ch=0; ch >(rReceiver.size(),{-99.,-99.}); // initialize - tTotalPower = 0.; // initialize - m1=m(-OpeningAngle*2.,radius*cos(theta), radius*sin(theta)); - m2=m(OpeningAngle*2.,radius*cos(theta), radius*sin(theta)); - - - for(unsigned i=0;iGetPosition() ) < 0 ) + { + //printf("Skipping! out of Bounds!: tReceiverTime=%e\n",tReceiverTime); + continue; + } + } - tCurrentParticle.Interpolate(tRetardedTime); + if(currentPatch->GetPreviousRetardedIndex() == -99.) + { + CurrentIndex=FindNode(tReceiverTime, kassiopeiaTimeStep, historySize-1); + tCurrentParticle = fParticleHistory[CurrentIndex]; + tRetardedTime = tReceiverTime - (tCurrentParticle.GetPosition() - currentPatch->GetPosition() ).Magnitude() / LMCConst::C(); + } + else + { + CurrentIndex = currentPatch->GetPreviousRetardedIndex(); + tRetardedTime = currentPatch->GetPreviousRetardedTime() + fDigitizerTimeStep; + } + + CurrentIndex = FindNode(tRetardedTime,kassiopeiaTimeStep,CurrentIndex); + CurrentIndex = std::min(std::max(CurrentIndex,0) , historySize - 1); - //Change the kassiopeia step we expand around if the interpolation time displacement is too large - if(fabs(tCurrentParticle.GetTime(true) - tCurrentParticle.GetTime(false)) > dtStepSize) - { - CurrentIndex=FindNode(tRetardedTime,dtStepSize,CurrentIndex); - tCurrentParticle=fParticleHistory[CurrentIndex]; + tCurrentParticle = fParticleHistory[CurrentIndex]; tCurrentParticle.Interpolate(tRetardedTime); - } - //printf("%e %e\n",tCurrentParticle.GetTimeDisplacement(),dtStepSize*0.75); - - tSpaceTimeInterval = GetSpaceTimeInterval(tCurrentParticle.GetTime(true), tReceiverTime, tCurrentParticle.GetPosition(true), tReceiverPosition); - - tOldSpaceTimeInterval = tSpaceTimeInterval; - } - - - PreviousTimes[i].first = CurrentIndex; - PreviousTimes[i].second = tRetardedTime; - - - LMCThreeVector tECrossH = tCurrentParticle.CalculateElectricField(rReceiver[i]).Cross(tCurrentParticle.CalculateMagneticField(rReceiver[i])); - LMCThreeVector tDirection = tReceiverPosition - tCurrentParticle.GetPosition(true); - - tTotalPower += dx * dx * tECrossH.Dot(tDirection.Unit()) / rReceiver.size() ;// * (fabs(tCurrentParticle.GetPosition(true).Z())<0.01); - - double tVelZ = tCurrentParticle.GetVelocity(true).Z(); - double tCosTheta = tVelZ * tDirection.Z() / tDirection.Magnitude() / fabs(tVelZ); - double tDopplerFrequency = tCurrentParticle.GetCyclotronFrequency() / ( 1. - fabs(tVelZ) / LMCConst::C() * tCosTheta); - double tTransitTime = (tReceiverPosition - tCurrentParticle.GetPosition(true)).Magnitude() / LMCConst::C(); - - if (tRetardedTime+tTransitTime > fDigitizerTimeStep) // if the signal has been present for longer than fDigitizerTimeStep - { - tVoltagePhase[ch*nelementsZ + z_position]+= tDopplerFrequency * fDigitizerTimeStep / rReceiver.size(); - } - else // if this is the first digitizer sample - { - tVoltagePhase[ch*nelementsZ + z_position]+= - tDopplerFrequency * (tReceiverTime - (tRetardedTime+tTransitTime)) / rReceiver.size(); -// printf("Retarding voltage phase: fDigitizerTimeStep is %g and tReceiverTime is %g and tRetardedTime is %g\n and t_old is %g\n", -// fDigitizerTimeStep, tReceiverTime, tRetardedTime, t_old); getchar(); - } + tSpaceTimeInterval = GetSpaceTimeInterval(tCurrentParticle.GetTime(true), tReceiverTime, tCurrentParticle.GetPosition(true), currentPatch->GetPosition()); - /* - if (i==0) // check receiver point 0 for each channel. It should be the same each time. - { - printf("rx point 0: zposition is %d and channel is %d, fcyc is %g and tVelZ is %g, dopplerfreq is %g, costheta is %f\n", - z_position, ch, tDopplerFrequency, tCurrentParticle.GetCyclotronFrequency(), tVelZ, tCosTheta); - printf("current index is %d\n", CurrentIndex); - } - */ + double tOldSpaceTimeInterval=99.; + //Converge to root + for(int j=0;j<25;++j) + { + //++tAverageIterations; - double tMinHFSS=1e-8; - const int nHFSSBins=2048; + tRetardedTime = GetStepRoot(tCurrentParticle, tReceiverTime, currentPatch->GetPosition(), tSpaceTimeInterval); + tCurrentParticle.Interpolate(tRetardedTime); + //Change the kassiopeia step we expand around if the interpolation time displacement is too large + if(fabs(tCurrentParticle.GetTime(true) - tCurrentParticle.GetTime(false)) > kassiopeiaTimeStep) + { + CurrentIndex=FindNode(tRetardedTime,kassiopeiaTimeStep,CurrentIndex); + tCurrentParticle=fParticleHistory[CurrentIndex]; + tCurrentParticle.Interpolate(tRetardedTime); + } - if( fWriteNFD && (fNFDIndex < nHFSSBins) && (tReceiverTime >= tMinHFSS) ) - { - LMCThreeVector tmpElectricField, tmpMagneticField; - tmpElectricField = tCurrentParticle.CalculateElectricField(rReceiver[i]); - tmpMagneticField = tCurrentParticle.CalculateMagneticField(rReceiver[i]); - - //If doing the first receiver - if(!i) ++fNFDIndex; - - for(int j=0;jGetPosition()); + tOldSpaceTimeInterval = tSpaceTimeInterval; + } - NFDElectricFieldFreq[j][i][k][0] += tmpElectricField[k] * ( DownConvert[0] * DFTFactor[0] - DownConvert[1] * DFTFactor[1] ) / nHFSSBins; - NFDElectricFieldFreq[j][i][k][1] += tmpElectricField[k] * ( DownConvert[0] * DFTFactor[1] + DownConvert[1] * DFTFactor[0] ) / nHFSSBins; + currentPatch->SetPreviousRetardedIndex(CurrentIndex); + currentPatch->SetPreviousRetardedTime(tRetardedTime); - NFDMagneticFieldFreq[j][i][k][0] += tmpMagneticField[k] * ( DownConvert[0] * DFTFactor[0] - DownConvert[1] * DFTFactor[1] ) / nHFSSBins; - NFDMagneticFieldFreq[j][i][k][1] += tmpMagneticField[k] * ( DownConvert[0] * DFTFactor[1] + DownConvert[1] * DFTFactor[0] ) / nHFSSBins; - } - } - } - else if(fWriteNFD && fNFDIndex >=nHFSSBins) - { - printf("Done! \n"); - } + LMCThreeVector tDirection = currentPatch->GetPosition() - tCurrentParticle.GetPosition(true); + double tVelZ = tCurrentParticle.GetVelocity(true).Z(); + double tCosTheta = tVelZ * tDirection.Z() / tDirection.Magnitude() / fabs(tVelZ); + double tDopplerFrequency = tCurrentParticle.GetCyclotronFrequency() / ( 1. - fabs(tVelZ) / LMCConst::C() * tCosTheta); - } // i, rReceiver.size() loop. + currentPatch->SetInstantaneousFrequency( tDopplerFrequency); + currentPatch->SetIncidentElectricField( tCurrentParticle.CalculateElectricField(currentPatch->GetPosition() )); + currentPatch->SetIncidentMagneticField( tCurrentParticle.CalculateMagneticField(currentPatch->GetPosition() )); + aSignal->LongSignalTimeComplex()[channelIndex*signalSize*aSignal->DecimationFactor() + index][0] += currentPatch->GetVoltage(); + aSignal->LongSignalTimeComplex()[channelIndex*signalSize*aSignal->DecimationFactor() + index][1] += currentPatch->GetVoltage(); - DirectivityFactor = directivity(m1, m2, radius*cos(theta), radius*sin(theta), - tCurrentParticle.GetPosition(true).GetX(), tCurrentParticle.GetPosition(true).GetY()); + } // z_position waveguide element stepping loop. - //if (fabs(tCurrentParticle.GetPosition(true).GetZ()-(double)(z_position-4)*0.01) < 0.005 ) - //{ - aSignal->LongSignalTimeComplex()[ch*signalSize*aSignal->DecimationFactor() + index][0] += DirectivityFactor * sqrt(50.)*sqrt(tTotalPower) * cos(tVoltagePhase[ch*nelementsZ + z_position] - phi_LO); - aSignal->LongSignalTimeComplex()[ch*signalSize*aSignal->DecimationFactor() + index][1] += DirectivityFactor * sqrt(50.)*sqrt(tTotalPower) * sin(tVoltagePhase[ch*nelementsZ + z_position] - phi_LO); - //} + //aSignal->LongSignalTimeComplex()[channelIndex*signalSize*aSignal->DecimationFactor() + index][0] += allChannels[channelIndex].PhasedSum(); - } // z_position waveguide element stepping loop. - } // nchannels loop. + } // nChannels loop. t_old += fDigitizerTimeStep; @@ -451,7 +202,8 @@ namespace locust //Return iterator of fParticleHistory particle closest to the time we are evaluating - int FreeFieldSignalGenerator::FindNode(double tNew, double dtStepSize, int kIndexOld) const + //Make simpler!!! + int FreeFieldSignalGenerator::FindNode(double tNew, double kassiopeiaTimeStep, int kIndexOld) const { int tHistorySize = fParticleHistory.size(); @@ -460,7 +212,7 @@ namespace locust double tOld = fParticleHistory[ kIndexOld ].GetTime(); - int kIndexMid=round((tNew-tOld)/dtStepSize) + kIndexOld; + int kIndexMid=round((tNew-tOld)/kassiopeiaTimeStep) + kIndexOld; kIndexMid = std::min( std::max(kIndexMid,0) , tHistorySize - 1 ); int kIndexSearchWidth; @@ -485,34 +237,50 @@ namespace locust } int tNodeIndex = it - fParticleHistory.begin(); - //if(tNodeIndex < 0 || tNodeIndex > (tHistorySize - 1)) - //{ - // LERROR(lmclog, "OUT OF INDEX SEARCH"); - //} return tNodeIndex; } - bool FreeFieldSignalGenerator::DoGenerate( Signal* aSignal ) + void FreeFieldSignalGenerator::InitializePatchArray() { - if(fWriteNFD) + + const unsigned nChannels = fNChannels; + const int nReceivers = 13; //Number of receivers per channel + + const double patchSpacingZ = 0.108; + const double patchRadius = 0.05; + double zPosition; + double theta; + const double dThetaArray = 2. * LMCConst::Pi() / nChannels; //Divide the circle into nChannels + + PatchAntenna modelPatch; + + allChannels.resize(nChannels); + + for(int channelIndex = 0; channelIndex < nChannels; ++channelIndex) { - //Initialize to correct sizes/ all zeros - NFDElectricFieldFreq.resize(NFDFrequencies.size()); - NFDMagneticFieldFreq.resize(NFDFrequencies.size()); - for(int i=0;i, 3> >(rReceiver.size(),{0.}); - NFDMagneticFieldFreq[i]=std::vector, 3> >(rReceiver.size(),{0.}); + zPosition = (receiverIndex - (nReceivers - 1.) /2.) * patchSpacingZ; + + modelPatch.SetCenterPosition({patchRadius * cos(theta) , patchRadius * sin(theta) , zPosition }); + modelPatch.SetPolarizationDirection({sin(theta), -cos(theta), 0.}); + modelPatch.SetNormalDirection({-cos(theta), -sin(theta), 0.}); //Say normals point inwards + allChannels[channelIndex].AddReceiver(modelPatch); } } + } + + bool FreeFieldSignalGenerator::DoGenerate( Signal* aSignal ) + { + InitializePatchArray(); //n samples for event spacing. int PreEventCounter = 0; const int NPreEventSamples = 150000; - //printf("fwritenfd is %d\n", fWriteNFD); getchar(); - std::thread Kassiopeia (KassiopeiaInit, gxml_filename); // spawn new thread fRunInProgress = true; @@ -526,12 +294,10 @@ namespace locust if (fPreEventInProgress) { PreEventCounter += 1; -// printf("preeventcounter is %d\n", PreEventCounter); if (PreEventCounter > NPreEventSamples) // finished noise samples. Start event. { fPreEventInProgress = false; // reset. fEventInProgress = true; - //printf("LMC about to wakebeforeevent\n"); WakeBeforeEvent(); // trigger Kass event. } } @@ -556,7 +322,7 @@ namespace locust } // for loop - if(fWriteNFD) NFDWrite(); + //if(fWriteNFD) NFDWrite(); //delete [] ImaginarySignal; diff --git a/Source/Generators/LMCFreeFieldSignalGenerator.hh b/Source/Generators/LMCFreeFieldSignalGenerator.hh index 17cb9c9d..71ea328b 100644 --- a/Source/Generators/LMCFreeFieldSignalGenerator.hh +++ b/Source/Generators/LMCFreeFieldSignalGenerator.hh @@ -11,6 +11,9 @@ #include "LMCThreeVector.hh" #include "LMCGenerator.hh" +#include "LMCChannel.hh" +#include "LMCPatchAntenna.hh" + namespace locust { @@ -46,24 +49,16 @@ namespace locust private: - std::vector rReceiver; //Vector that contains 3D position of all points at which the fields are evaluated (ie. along receiver surface) + std::vector< Channel > allChannels; //Vector that contains pointer to all channels + std::vector > PreviousTimes; //Cache the results from previous iteration. [0] is previous index, [1] is corresponding retarded time of previous solution double fLO_Frequency; // typically defined by a parameter in json file. std::string gxml_filename; - std::vector, 3 > > > NFDElectricFieldFreq; //Should use the LMCThreeVectors too..... - std::vector, 3 > > > NFDMagneticFieldFreq; - - bool fWriteNFD; - std::vector NFDFrequencies; - std::string fAND_filename; - std::string fNFD_filename; - bool DoGenerate( Signal* aSignal ); void* DriveAntenna(int PreEventCounter, unsigned index, Signal* aSignal); - - void NFDWrite(); + void InitializePatchArray(); int FindNode(double tNew, double dtStepSize, int IndexOld) const; diff --git a/Source/Generators/LMCKassSignalGenerator.cc b/Source/Generators/LMCKassSignalGenerator.cc index be69938a..c54fb31a 100644 --- a/Source/Generators/LMCKassSignalGenerator.cc +++ b/Source/Generators/LMCKassSignalGenerator.cc @@ -189,15 +189,15 @@ namespace locust { RealVoltage2 *= 0.03; // replace short with terminator. ImagVoltage2 *= 0.03; // replace short with terminator. - aSignal->LongSignalTimeComplex()[ index ][0] += gain*sqrt(50.)*TE11ModeExcitation() * sqrt(tLarmorPower) * (RealVoltage1 + RealVoltage2); - aSignal->LongSignalTimeComplex()[ index ][1] += gain*sqrt(50.)*TE11ModeExcitation() * sqrt(tLarmorPower) * (ImagVoltage1 + ImagVoltage2); + aSignal->LongSignalTimeComplex()[ index ][0] += gain*sqrt(50.)*TE11ModeExcitation() * sqrt(tLarmorPower/2.) * (RealVoltage1 + RealVoltage2); + aSignal->LongSignalTimeComplex()[ index ][1] += gain*sqrt(50.)*TE11ModeExcitation() * sqrt(tLarmorPower/2.) * (ImagVoltage1 + ImagVoltage2); } else if (Project8Phase == 1) { // assume 50 ohm impedance // RealVoltage2 *= 0.25; // some loss at short. - aSignal->LongSignalTimeComplex()[ index ][0] += gain*sqrt(50.) * TE01ModeExcitation() * sqrt(tLarmorPower) * (RealVoltage1 + RealVoltage2); - aSignal->LongSignalTimeComplex()[ index ][1] += gain*sqrt(50.) * TE01ModeExcitation() * sqrt(tLarmorPower) * (ImagVoltage1 + ImagVoltage2); + aSignal->LongSignalTimeComplex()[ index ][0] += gain*sqrt(50.) * TE01ModeExcitation() * sqrt(tLarmorPower/2.) * (RealVoltage1 + RealVoltage2); + aSignal->LongSignalTimeComplex()[ index ][1] += gain*sqrt(50.) * TE01ModeExcitation() * sqrt(tLarmorPower/2.) * (ImagVoltage1 + ImagVoltage2); } @@ -263,25 +263,33 @@ namespace locust std::thread Kassiopeia (KassiopeiaInit,gxml_filename); // spawn new thread fRunInProgress = true; fKassEventReady = false; + int StartEventTimer = 0; for( unsigned index = 0; index < aSignal->DecimationFactor()*aSignal->TimeSize(); ++index ) { if ((!fEventInProgress) && (fRunInProgress) && (!fPreEventInProgress)) { if (ReceivedKassReady()) fPreEventInProgress = true; - printf("LMC says it ReceivedKassReady()\n"); + printf("LMC says it ReceivedKassReady()\n"); + if (Project8Phase == 2) + { + if ((index-StartEventTimer) > 1e6) + { + break; + } + } } if (fPreEventInProgress) { - PreEventCounter += 1; -// printf("preeventcounter is %d\n", PreEventCounter); + PreEventCounter += 1; if (PreEventCounter > NPreEventSamples) // finished noise samples. Start event. { fPreEventInProgress = false; // reset. fEventInProgress = true; printf("LMC about to WakeBeforeEvent()\n"); + StartEventTimer = index; WakeBeforeEvent(); // trigger Kass event. } } @@ -289,13 +297,11 @@ namespace locust if (fEventInProgress) // fEventInProgress if (fEventInProgress) // check again. { - //printf("waiting for digitizer trigger ... index is %d\n", index); std::unique_lock< std::mutex >tLock( fMutexDigitizer, std::defer_lock ); tLock.lock(); fDigitizerCondition.wait( tLock ); if (fEventInProgress) { - //printf("about to drive antenna, PEV is %d\n", PreEventCounter); DriveAntenna(PreEventCounter, index, aSignal, fp); PreEventCounter = 0; // reset } @@ -303,24 +309,19 @@ namespace locust } } // for loop - fclose(fp); - - // trigger any remaining events in Kassiopeia so that its thread can finish. - fDoneWithSignalGeneration = true; printf("finished signal loop.\n"); - while (fRunInProgress) - { - std::this_thread::sleep_for(std::chrono::milliseconds(2000)); - if (fRunInProgress) - { - std::this_thread::sleep_for(std::chrono::milliseconds(100)); - if (!fEventInProgress) - if (ReceivedKassReady()) - WakeBeforeEvent(); - } - } + + fclose(fp); + fRunInProgress = false; // tell Kassiopeia to finish. + fDoneWithSignalGeneration = true; // tell LMCCyclotronRadExtractor + // if (fEventInProgress) + // if (ReceivedKassReady()) + WakeBeforeEvent(); + + Kassiopeia.join(); // wait for Kassiopeia to finish. + return true; } diff --git a/Source/Generators/LMCLocalOscillatorGenerator.cc b/Source/Generators/LMCLocalOscillatorGenerator.cc new file mode 100644 index 00000000..48cd0410 --- /dev/null +++ b/Source/Generators/LMCLocalOscillatorGenerator.cc @@ -0,0 +1,103 @@ +/* + * LMCLocalOscillatorGenerator.cc + * + * Created on: Apr 20, 2018 + * Author: buzinsky + */ + +#include "LMCLocalOscillatorGenerator.hh" + +#include "LMCRunLengthCalculator.hh" +#include "LMCSimulationController.hh" +#include "LMCConst.hh" + +#include "logger.hh" + +using std::string; + +namespace locust +{ + LOGGER( lmclog, "LocalOscillatorGenerator" ); + + MT_REGISTER_GENERATOR(LocalOscillatorGenerator, "local-oscillator"); + + LocalOscillatorGenerator::LocalOscillatorGenerator( const std::string& aName ) : + Generator( aName ), + fDoGenerateFunc( &LocalOscillatorGenerator::DoGenerateTime ) + { + fRequiredSignalState = Signal::kTime; + } + + LocalOscillatorGenerator::~LocalOscillatorGenerator() + { + } + + bool LocalOscillatorGenerator::Configure( const scarab::param_node* aParam ) + { + if( aParam == NULL) return true; + + if( aParam->has( "LO-Frequency" ) ) + { + fLOFrequency = aParam->get_value< double >( "LO-Frequency" ); + } + + return true; + } + + void LocalOscillatorGenerator::Accept( GeneratorVisitor* aVisitor ) const + { + aVisitor->Visit( this ); + return; + } + + + bool LocalOscillatorGenerator::DoGenerate( Signal* aSignal ) + { + return (this->*fDoGenerateFunc)( aSignal ); + } + + bool LocalOscillatorGenerator::DoGenerateTime( Signal* aSignal ) + { + + const unsigned nChannels = fNChannels; + + RunLengthCalculator *simpleCalculator = new RunLengthCalculator(); + const double acquisitionTimeStep = 1. / ( simpleCalculator->GetAcquisitionRate() * 1.e6); //Rate is in MHz + + const unsigned signalSize = aSignal->TimeSize(); + const unsigned decimationFactor = aSignal->DecimationFactor(); + double *downmixedSignal; + + + for(unsigned channelIndex = 0; channelIndex < nChannels; ++channelIndex) + { + double localOscillatorPhase = 0.; + + for( unsigned index = 0; index < decimationFactor*signalSize; ++index ) + { + localOscillatorPhase += 2. * LMCConst::Pi() * fLOFrequency * acquisitionTimeStep; + fftw_complex LOMixingPhase = {cos(localOscillatorPhase), -sin(localOscillatorPhase) }; + ComplexMultiplication( aSignal->LongSignalTimeComplex()[channelIndex*signalSize*decimationFactor + index] , LOMixingPhase); + + //Assign to aSignal + aSignal->LongSignalTimeComplex()[channelIndex*signalSize*decimationFactor + index][0] = downmixedSignal[0]; + aSignal->LongSignalTimeComplex()[channelIndex*signalSize*decimationFactor + index][1] = downmixedSignal[1]; + + } + } + + //Cleanup + delete []downmixedSignal; + + return true; + } + + double* LocalOscillatorGenerator::ComplexMultiplication(fftw_complex a, fftw_complex b) + { + double *product = new double[2]; + product[0] = a[0] * b[0] - a[1] * b[1]; + product[1] = a[1] * b[0] + a[0] * b[1]; + return product; + } + +} /* namespace locust */ diff --git a/Source/Generators/LMCLocalOscillatorGenerator.hh b/Source/Generators/LMCLocalOscillatorGenerator.hh new file mode 100644 index 00000000..f507e5eb --- /dev/null +++ b/Source/Generators/LMCLocalOscillatorGenerator.hh @@ -0,0 +1,52 @@ +/* + * LMCLocalOscillatorGenerator.hh + * + * Created on: Apr 20 2018 + * Author: buzinsky + */ + +#ifndef LMCLOCALOSCILLATORGENERATOR_HH_ +#define LMCLOCALOSCILLATORGENERATOR_HH_ + +#include "LMCGenerator.hh" + +namespace locust +{ + + /*! + @class LocalOscillatorGenerator + @author N. Buzinsky + + @brief Downmix the signal (multiplicative) by f_LO + + @details + Time domain only + + Configuration name: "local-oscillator" + + */ + +class LocalOscillatorGenerator : public Generator + { + public: + LocalOscillatorGenerator( const std::string& aName = "local-oscillator" ); + virtual ~LocalOscillatorGenerator(); + + bool Configure( const scarab::param_node* aNode ); + + void Accept( GeneratorVisitor* aVisitor ) const; + + private: + double fLOFrequency; + double* ComplexMultiplication(fftw_complex a, fftw_complex b); + + bool DoGenerate( Signal* aSignal ); + + bool DoGenerateTime( Signal* aSignal ); + + bool (LocalOscillatorGenerator::*fDoGenerateFunc)( Signal* aSignal ); + }; + +} /* namespace locust */ + +#endif /* LMCLOCALOSCILLATORGENERATOR_HH_ */ diff --git a/Source/Generators/LMCPatchSignalGenerator.cc b/Source/Generators/LMCPatchSignalGenerator.cc index 18712a07..0fb92874 100644 --- a/Source/Generators/LMCPatchSignalGenerator.cc +++ b/Source/Generators/LMCPatchSignalGenerator.cc @@ -31,6 +31,7 @@ namespace locust PatchSignalGenerator::PatchSignalGenerator( const std::string& aName ) : Generator( aName ), fLO_Frequency( 0.), + fCorporateFeed( 0 ), gxml_filename("blank.xml") { fRequiredSignalState = Signal::kTime; @@ -52,9 +53,16 @@ namespace locust { gxml_filename = aParam->get_value< std::string >( "xml-filename" ); } + if( aParam->has( "feed" ) ) + { + if (aParam->get_value< std::string >( "feed" ) == "corporate") + fCorporateFeed = 1; // default + else if (aParam->get_value< std::string >( "feed" ) == "series") + fCorporateFeed = 0; + else + fCorporateFeed = 1; // default - - + } return true; } @@ -168,43 +176,87 @@ namespace locust return AOIFactor; } + double GetVoltageAmpFromPlaneWave() + { + double AntennaFactor = 1./420.; + double epsilon0 = 8.854e-12; // C^2/N/m^2 + double c = 3.e8; // m/s + + // S = epsilon0 c E0^2 / 2. // power/area + // 0.6e-21 W/Hz * 24.e3 Hz / (0.00375*0.002916) = S = 1.3e-12 W/m^2 + // We should detect 0.6e-21 W/Hz * 24.e3 Hz in Katydid. + // E0 = sqrt(2.*S/epsilon0/c) + // effective patch area 0.00004583662 m^2 + + // double S = 0.6e-21*24.e3/(0.00004583662); // W/m^2, effective aperture. + // double E0 = pow(2.*S/epsilon0/c, 0.5); + double E0 = 1.0; // V/m, test case + double amplitude = E0*AntennaFactor; // volts + return amplitude; + + + } + double GetVoltageAmplitude(LMCThreeVector IncidentElectricField, LMCThreeVector IncidentKVector, double PatchPhi, double DopplerFrequency) { - double AntennaFactor = 1./400.; // 1/m, placeholder + // double AntennaFactor = 1.32*1.2*0.683*0.23*1./400.; // 1/m, placeholder + double AntennaFactor = 1./1607.5; double MismatchFactor = GetMismatchFactor(DopplerFrequency); double AOIFactor = GetAOIFactor(IncidentKVector, PatchPhi); // k dot patchnormal LMCThreeVector PatchPolarizationVector; PatchPolarizationVector.SetComponents(-sin(PatchPhi), cos(PatchPhi), 0.0); - double VoltageAmplitude = fabs( AntennaFactor * IncidentElectricField.Dot(PatchPolarizationVector) * MismatchFactor * AOIFactor); -// printf("IncidentElectricField.Dot(PatchPolarizationVector) is %g and VoltageAmplitude is %g\n", IncidentElectricField.Dot(PatchPolarizationVector), VoltageAmplitude); getchar(); + double VoltageAmplitude = fabs( AntennaFactor * IncidentElectricField.Dot(PatchPolarizationVector) * MismatchFactor * AOIFactor); + + // if (VoltageAmplitude>0.) {printf("IncidentElectricField.Dot(PatchPolarizationVector) is %g and VoltageAmplitude is %g\n", IncidentElectricField.Dot(PatchPolarizationVector), VoltageAmplitude); getchar();} return VoltageAmplitude; } + double ZPositionPatch(unsigned z_index) + { + double PatchOffset = 0.; + double ZPosition = 0.; + if (NPATCHES_PER_STRIP%2==1) + { + PatchOffset = -(NPATCHES_PER_STRIP-1.)/2. * PATCH_SPACING; // far left patch. + ZPosition = PatchOffset + (double)(z_index)*PATCH_SPACING; + } + else + { + PatchOffset = -(PATCH_SPACING)/2. - ((NPATCHES_PER_STRIP/2.) - 1.) * PATCH_SPACING; + ZPosition = PatchOffset + (double)(z_index)*PATCH_SPACING; + } + return ZPosition; + } + + + - void AddOnePatchVoltageToStripSum(Signal* aSignal, double VoltageAmplitude, double VoltagePhase, double phi_LO, unsigned channelindex) + double GetLinePhaseCorr(unsigned z_index, double DopplerFrequency) + { + double D = 0.007711*0.95; // m. 18.0 keV 90 degree electron, lambda in kapton. + double c_n = LMCConst::C()/1.5; // speed of light in Kapton. + double lambda = c_n/(DopplerFrequency/(2.*LMCConst::Pi())); + double dphi = 2.*LMCConst::Pi() * D * z_index / lambda; + return dphi; + } + + + + void PatchSignalGenerator::AddOnePatchVoltageToStripSum(Signal* aSignal, double VoltageAmplitude, double VoltagePhase, double phi_LO, unsigned channelindex, unsigned z_index, double DopplerFrequency) { + double LinePhaseCorr = 0.; + if (fCorporateFeed == false) + { + LinePhaseCorr = GetLinePhaseCorr(z_index, DopplerFrequency); + VoltagePhase += LinePhaseCorr; + } + + // if (VoltageAmplitude>0.) {printf("voltageamplitude is %g\n", VoltageAmplitude); getchar();} aSignal->LongSignalTimeComplex()[channelindex][0] += VoltageAmplitude * cos(VoltagePhase - phi_LO); aSignal->LongSignalTimeComplex()[channelindex][1] += VoltageAmplitude * sin(VoltagePhase - phi_LO); } - double ZPositionPatch(unsigned z_index) - { - double PatchOffset = 0.; - double ZPosition = 0.; - if (NPATCHES_PER_STRIP%2==1) - { - PatchOffset = -(NPATCHES_PER_STRIP-1.)/2. * PATCH_SPACING; // far left patch. - ZPosition = PatchOffset + (double)(z_index)*PATCH_SPACING; - } - else - { - PatchOffset = -(PATCH_SPACING)/2. - ((NPATCHES_PER_STRIP/2.) - 1.) * PATCH_SPACING; - ZPosition = PatchOffset + (double)(z_index)*PATCH_SPACING; - } - return ZPosition; - } - void* PatchSignalGenerator::DriveAntenna(int PreEventCounter, unsigned index, Signal* aSignal, FILE *fp) @@ -343,7 +395,8 @@ namespace locust LMCThreeVector tDirection = tReceiverPosition - tCurrentParticle.GetPosition(true); tTotalPower = dx * dy * tECrossH.Dot(tDirection.Unit()) / rReceiver.size() ; - // printf("total power hitting patch is %g\n", tTotalPower); getchar(); + // printf("total power hitting patch is %g\n", tTotalPower*rReceiver.size()); getchar(); + tIncidentElectricField += tCurrentParticle.CalculateElectricField(rReceiver[i]) / rReceiver.size(); tIncidentKVector += tECrossH / rReceiver.size(); @@ -359,8 +412,8 @@ namespace locust } else // if this is the first digitizer sample, the voltage phase doesn't advance for the full dt. { - tVoltagePhase[patchindex]+= // compressing - tDopplerFrequency * tRetardedTime / rReceiver.size(); + tVoltagePhase[patchindex]+= // compressing + tDopplerFrequency * tRetardedTime / rReceiver.size(); } /* @@ -376,13 +429,16 @@ namespace locust } // i, rReceiver.size() loop. - double tVoltageAmplitude = GetVoltageAmplitude(tIncidentElectricField, tIncidentKVector, PatchPhi, tAverageDopplerFrequency); - AddOnePatchVoltageToStripSum(aSignal, tVoltageAmplitude, tVoltagePhase[patchindex], phi_LO, channelindex); + double tVoltageAmplitude = GetVoltageAmplitude(tIncidentElectricField, tIncidentKVector, PatchPhi, tAverageDopplerFrequency); + // if (tAverageDopplerFrequency>0.) {printf("tAverageDopplerFrequency is %g\n", tAverageDopplerFrequency); getchar();} + // double tVoltageAmplitude = GetVoltageAmpFromPlaneWave(); + //printf("tvoltageamp is %g\n", tVoltageAmplitude); getchar(); + AddOnePatchVoltageToStripSum(aSignal, tVoltageAmplitude, tVoltagePhase[patchindex], phi_LO, channelindex, z_index, tAverageDopplerFrequency); } // z_index waveguide element stepping loop. - // printf("signal %d with frequency %g is %g\n", channelindex, tAverageDopplerFrequency, aSignal->LongSignalTimeComplex()[channelindex][0]); getchar(); + // printf("signal %d with frequency %g is %g\n", channelindex, tAverageDopplerFrequency, aSignal->LongSignalTimeComplex()[channelindex][0]); getchar(); } // nchannels loop. diff --git a/Source/Generators/LMCPatchSignalGenerator.hh b/Source/Generators/LMCPatchSignalGenerator.hh index 1d9f5998..bdbeaa78 100644 --- a/Source/Generators/LMCPatchSignalGenerator.hh +++ b/Source/Generators/LMCPatchSignalGenerator.hh @@ -42,6 +42,9 @@ namespace locust bool Configure( const scarab::param_node* aNode ); void Accept( GeneratorVisitor* aVisitor ) const; + + void AddOnePatchVoltageToStripSum(Signal* aSignal, double VoltageAmplitude, double VoltagePhase, double phi_LO, unsigned channelindex, unsigned z_index, double DopplerFrequency); + private: @@ -50,6 +53,7 @@ namespace locust double fLO_Frequency; // typically defined by a parameter in json file. std::string gxml_filename; + bool fCorporateFeed; bool DoGenerate( Signal* aSignal ); diff --git a/Source/Kassiopeia/LMCCyclotronRadiationExtractor.cc b/Source/Kassiopeia/LMCCyclotronRadiationExtractor.cc index f7d2e037..78fea82d 100644 --- a/Source/Kassiopeia/LMCCyclotronRadiationExtractor.cc +++ b/Source/Kassiopeia/LMCCyclotronRadiationExtractor.cc @@ -40,7 +40,7 @@ namespace locust Project8Phase = P8Phase; if (P8Phase==1) { - CENTER_TO_SHORT = 0.0485; // m, 0.047 is tuned. + CENTER_TO_SHORT = 0.0488; // m, 0.047 is tuned. CENTER_TO_ANTENNA = 0.045; // m } if (P8Phase==2) @@ -125,8 +125,9 @@ namespace locust double GammaZ = 1.0/pow(1.0-pow(zvelocity/GetGroupVelocityTE01(aFinalParticle),2.),0.5); double fprime_short = fcyc*GammaZ*(1.+zvelocity/GroupVelocity); - double phi_shortTE01 = LMCConst::Pi()/2. + 2.*LMCConst::Pi()*(fabs(zPosition)+CENTER_TO_SHORT)/(GroupVelocity/fprime_short); // phase of reflected field at position of electron. - double FieldFromShort = cos(0.) + cos(phi_shortTE01); // yes resonant enhancement. with reflection coefficient. + double phi_shortTE01 = LMCConst::Pi()/2. + 2.*LMCConst::Pi()*(fabs(zPosition) + CENTER_TO_SHORT)/(GroupVelocity/fprime_short); // phase of reflected field at position of electron. + + double FieldFromShort = cos(0.) + cos(phi_shortTE01); // printf("field sum at z=%g is %f with zvelocity %g\n", zPosition, FieldFromShort, zvelocity); // getchar(); @@ -243,23 +244,15 @@ namespace locust double CyclotronRadiationExtractor::GetTM01FieldWithTerminator(KSParticle& anInitialParticle, KSParticle& aFinalParticle) { - double dt = aFinalParticle.GetTime() - anInitialParticle.GetTime(); double tCyclotronFrequency = aFinalParticle.GetCyclotronFrequency(); double GroupVelocity = GetGroupVelocityTM01(aFinalParticle); double tVelocityZ = aFinalParticle.GetVelocity().GetZ(); double tPositionZ = aFinalParticle.GetPosition().GetZ(); double GammaZ = 1.0/sqrt(1.0-pow(tVelocityZ / GetGroupVelocityTM01(aFinalParticle),2.) ); - // printf("tcyc is %g and GammaZ is %g and tVelocityZ is %g\n", tCyclotronFrequency, GammaZ, tVelocityZ); double fprime_polarizer = tCyclotronFrequency*GammaZ*(1.-tVelocityZ/GroupVelocity); - // printf("GroupVelocity is %.10g and fprime_short is %.10g\n", GroupVelocity, fprime_short); - double lambda_polarizer = GroupVelocity/fprime_polarizer; - double FieldFromPolarizer=0.; // other doppler shift - double TM01FieldAfterBounces = 0.; - double phi_polarizerTM01 = 0.; - - phi_polarizerTM01 = 2.*LMCConst::Pi()*(2.*(CENTER_TO_ANTENNA-tPositionZ)) - /(GroupVelocity/fprime_polarizer); - TM01FieldAfterBounces = cos(0.) + cos(phi_polarizerTM01); + + double phi_polarizerTM01 = 2.*LMCConst::Pi()*(2.*(CENTER_TO_ANTENNA-fabs(tPositionZ)))/(GroupVelocity/fprime_polarizer); + double TM01FieldAfterBounces = cos(0.) + cos(phi_polarizerTM01); //printf("TE11FieldAfterOneBounce is %f\n", TE11FieldAfterOneBounce); return TM01FieldAfterBounces; @@ -283,7 +276,7 @@ namespace locust { double TM01FieldWithTerminator = GetTM01FieldWithTerminator(anInitialParticle, aFinalParticle); double A01squ = GetCouplingFactorTM01(aFinalParticle); - double DampingFactorTM01 = 1. - A01squ + A01squ*TM01FieldWithTerminator*TM01FieldWithTerminator; // can be > 0 or < 0. + double DampingFactorTM01 = 1. - A01squ + A01squ*TM01FieldWithTerminator*TM01FieldWithTerminator; // = P'/P double DampingFactor = DampingFactorTM01; return DampingFactor; @@ -342,17 +335,16 @@ namespace locust bool CyclotronRadiationExtractor::ExecutePostStepModification( KSParticle& anInitialParticle, KSParticle& aFinalParticle, KSParticleQueue& aQueue ) { double DeltaE=0.; - //printf("fcyc is %g\n", anInitialParticle.GetCyclotronFrequency()); //getchar(); + + // printf("dE/dt is %g\n", (aFinalParticle.GetKineticEnergy() - anInitialParticle.GetKineticEnergy())/(aFinalParticle.GetTime() - anInitialParticle.GetTime())); //getchar(); + if(fP8Phase==1) { - // adjust power with reflections. DeltaE = GetDampingFactorPhase1(anInitialParticle, aFinalParticle)*(aFinalParticle.GetKineticEnergy() - anInitialParticle.GetKineticEnergy()); - //aFinalParticle.SetKineticEnergy((aFinalParticle.GetKineticEnergy() - DeltaE)); aFinalParticle.SetKineticEnergy((anInitialParticle.GetKineticEnergy() + DeltaE)); } if(fP8Phase==2) { - // adjust power with reflections. DeltaE = GetDampingFactorPhase2(anInitialParticle, aFinalParticle)*(aFinalParticle.GetKineticEnergy() - anInitialParticle.GetKineticEnergy()); //printf("poststep says DeltaE is %g\n", DeltaE); aFinalParticle.SetKineticEnergy((anInitialParticle.GetKineticEnergy() + DeltaE)); @@ -361,7 +353,10 @@ namespace locust if (!fDoneWithSignalGeneration) // if Locust is still acquiring voltages. { - if (t_old == 0.) fPitchAngle = -99.; // new electron needs central pitch angle reset. + if (t_old == 0.) + { + fPitchAngle = -99.; // new electron needs central pitch angle reset. + } double t_poststep = aFinalParticle.GetTime(); fNewParticleHistory.push_back(ExtractKassiopeiaParticle(anInitialParticle, aFinalParticle)); @@ -375,7 +370,6 @@ namespace locust //Dont want to check .back() of history if it is empty! -> Segfault if(fParticleHistory.size() && (fNewParticleHistory.back().GetTime() < fParticleHistory.back().GetTime())) { -// printf("New Particle!, t_old is %g\n", t_old); getchar(); t_poststep = 0.; fParticleHistory.clear(); } @@ -389,7 +383,6 @@ namespace locust fParticleHistory.push_back(ExtractKassiopeiaParticle(anInitialParticle, tParticleCopy)); tHistoryMaxSize = 5; - } else { @@ -408,17 +401,16 @@ namespace locust fNewParticleHistory.clear(); //Purge fParticleHistory of overly old entries - while(t_poststep-fParticleHistory.front().GetTime()>1e-7 || fParticleHistory.size() > tHistoryMaxSize) - { - fParticleHistory.pop_front(); - } - // printf("done purging\n"); + while(t_poststep-fParticleHistory.front().GetTime()>1e-7 || fParticleHistory.size() > tHistoryMaxSize) + { + fParticleHistory.pop_front(); + } tLock.unlock(); fDigitizerCondition.notify_one(); // notify Locust after writing. } - } // fDoneWithSignalGeneration + } // DoneWithSignalGeneration return true; } diff --git a/kassiopeia b/kassiopeia index c0f761a3..fb2e5857 160000 --- a/kassiopeia +++ b/kassiopeia @@ -1 +1 @@ -Subproject commit c0f761a3400817b6289a717d7b5fa46765a14574 +Subproject commit fb2e585728fb48a3cfe174f9977a2e8f1b091bfe