diff --git a/Config/Project8_JustKass.xml b/Config/Project8_JustKass.xml
index 85b4ba24..5f9eca5d 100644
--- a/Config/Project8_JustKass.xml
+++ b/Config/Project8_JustKass.xml
@@ -21,8 +21,8 @@
-
-
+
+
@@ -361,8 +361,8 @@
-
-
+
+
@@ -442,7 +442,7 @@
-
+
@@ -464,7 +464,7 @@
-
+
@@ -192,7 +192,7 @@
-
+
@@ -201,7 +201,7 @@
-
+
diff --git a/Config/Tutorial/Project8Phase2_WithRoot_Template.xml b/Config/Tutorial/Project8Phase2_WithRoot_Template.xml
index 450fe741..9694d91d 100644
--- a/Config/Tutorial/Project8Phase2_WithRoot_Template.xml
+++ b/Config/Tutorial/Project8Phase2_WithRoot_Template.xml
@@ -3,22 +3,21 @@
-
-
+
+
-
-
-
-
-
+
+
+
+
-
+
@@ -30,7 +29,7 @@
-
+
@@ -46,19 +45,22 @@
-
+
-
+
+
+
+
@@ -200,16 +202,16 @@
-
+
-
-
-
+
+
+
-
+
@@ -219,6 +221,25 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
@@ -411,18 +432,19 @@
-
+
+
-
+
@@ -492,17 +514,18 @@
-
+
-
+
+
-
+
@@ -539,7 +562,7 @@
name="project8_simulation"
run="1"
seed="[seed]"
- events="10000"
+ events="25"
magnetic_field="field_electromagnet"
magnetic_field="field_magnetic_main"
space="space_world"
@@ -554,7 +577,7 @@
-
diff --git a/Config/Tutorial/Project8Phase3Geometry.xml b/Config/Tutorial/Project8Phase3Geometry.xml
new file mode 100644
index 00000000..ec328b22
--- /dev/null
+++ b/Config/Tutorial/Project8Phase3Geometry.xml
@@ -0,0 +1,367 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/Config/Tutorial/Project8Phase3LargeGeometry.xml b/Config/Tutorial/Project8Phase3LargeGeometry.xml
new file mode 100644
index 00000000..42315123
--- /dev/null
+++ b/Config/Tutorial/Project8Phase3LargeGeometry.xml
@@ -0,0 +1,453 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/Config/Tutorial/Project8Phase3LongGeometry.xml b/Config/Tutorial/Project8Phase3LongGeometry.xml
new file mode 100644
index 00000000..955df8b5
--- /dev/null
+++ b/Config/Tutorial/Project8Phase3LongGeometry.xml
@@ -0,0 +1,453 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/Config/Tutorial/Project8Phase3_WithRoot_Template.xml b/Config/Tutorial/Project8Phase3_WithRoot_Template.xml
index f13cab94..221eec34 100644
--- a/Config/Tutorial/Project8Phase3_WithRoot_Template.xml
+++ b/Config/Tutorial/Project8Phase3_WithRoot_Template.xml
@@ -4,33 +4,33 @@
-
+
-
+
-
+
-
+
-
+
-
+
@@ -45,7 +45,7 @@
-
+
@@ -109,7 +109,7 @@
-
+
@@ -201,16 +201,16 @@
-
+
-
+
-
+
@@ -406,15 +406,11 @@
-
+
-
-
-
-
-
-
+
+
@@ -423,7 +419,7 @@
-
+
@@ -486,7 +482,7 @@
-
+
@@ -498,22 +494,8 @@
-
-
-
-
-
-
-
-
-
-
-
-
-
-
@@ -523,7 +505,7 @@
-->
-
+
@@ -551,7 +533,7 @@
-
+
diff --git a/Config/Tutorial/katydid_new.json b/Config/Tutorial/katydid_new.json
index b8fe722b..b39af865 100644
--- a/Config/Tutorial/katydid_new.json
+++ b/Config/Tutorial/katydid_new.json
@@ -35,8 +35,12 @@
"signal": "to-ps1:psd",
"slot": "writer:psd"
},
-
-
+
+ {
+ "signal": "fft1:fft",
+ "slot": "writer:fs-fftw-phase"
+ },
+
{
"signal": "to-ps1:psd",
"slot": "waterfall-writer:psd"
@@ -61,10 +65,11 @@
"egg1":
{
- "filename": "/cmn/p8/locust_mc/cbuild/bin/locust_mc_test.egg",
+ "filename": "/home/slocum/locust_mc/cbuild/bin/locust_mc.egg",
"egg-reader": "egg3",
"number-of-slices": 500,
"slice-size": 8192,
+ "stride": 2048,
"dac":
{
"time-series-type": "fftw"
@@ -82,7 +87,7 @@
"waterfall-writer":
{
- "output-file": "/cmn/p8/locust_mc/cbuild/bin/katydidwaterfall.root",
+ "output-file": "/home/penny/locust_mc/cbuild/bin/katydidwaterfall.root",
"file-flag": "recreate",
"min-time": 0.0,
"max-time": 0.025,
@@ -94,7 +99,7 @@
"writer":
{
- "output-file": "/cmn/p8/locust_mc/cbuild/bin/basic.root",
+ "output-file": "/home/penny/locust_mc/cbuild/bin/basic.root",
"file-flag": "recreate"
}
diff --git a/Source/CMakeLists.txt b/Source/CMakeLists.txt
index 0f55b0c5..4b390e69 100644
--- a/Source/CMakeLists.txt
+++ b/Source/CMakeLists.txt
@@ -76,6 +76,7 @@ if (locust_mc_BUILD_WITH_KASSIOPEIA)
set( LOCUST_MC_HEADER_FILES ${LOCUST_MC_HEADER_FILES}
Generators/LMCFreeFieldSignalGenerator.hh
+ Generators/LMCPatchSignalGenerator.hh
Generators/LMCKassSignalGenerator.hh
Core/LMCGlobalsDeclaration.hh
@@ -93,6 +94,7 @@ if (locust_mc_BUILD_WITH_KASSIOPEIA)
set( LOCUST_MC_SOURCE_FILES ${LOCUST_MC_SOURCE_FILES}
Generators/LMCFreeFieldSignalGenerator.cc
+ Generators/LMCPatchSignalGenerator.cc
Generators/LMCKassSignalGenerator.cc
Core/LMCHFSSReader.cc
diff --git a/Source/Core/LMCGlobalsDeclaration.hh b/Source/Core/LMCGlobalsDeclaration.hh
index d101ea18..69401cea 100644
--- a/Source/Core/LMCGlobalsDeclaration.hh
+++ b/Source/Core/LMCGlobalsDeclaration.hh
@@ -7,8 +7,12 @@
#ifndef GLOBALSDECLARATION_HH_
#define GLOBALSDECLARATION_HH_
-#define CENTER_TO_SHORT 0.0760 // m
-#define CENTER_TO_ANTENNA 0.07685 // m
+
+#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 PATCH_RING_OFFSET 0.0 // m. Offset of entire array in X direction
+
#include
@@ -19,6 +23,8 @@
extern int Project8Phase; // 1, 2, or 3
+extern double CENTER_TO_SHORT;
+extern double CENTER_TO_ANTENNA;
extern double t_old;
extern double fDigitizerTimeStep;
diff --git a/Source/Core/LMCGlobalsDefinition.hh b/Source/Core/LMCGlobalsDefinition.hh
index 75014d1c..48b34d1f 100644
--- a/Source/Core/LMCGlobalsDefinition.hh
+++ b/Source/Core/LMCGlobalsDefinition.hh
@@ -9,6 +9,8 @@
#define GLOBALSDEFINITION_HH_
int Project8Phase = 0; // 1, 2, or 3, defined with the step modifier instance in the xml file.
+double CENTER_TO_SHORT = 0.0;
+double CENTER_TO_ANTENNA = 0.0;
double t_old = -99.;
diff --git a/Source/Core/LMCRunLengthCalculator.cc b/Source/Core/LMCRunLengthCalculator.cc
index 446bff66..a9bdc528 100644
--- a/Source/Core/LMCRunLengthCalculator.cc
+++ b/Source/Core/LMCRunLengthCalculator.cc
@@ -90,6 +90,11 @@ namespace locust
// nothing to see here, move along, please
return;
}
+ void RunLengthCalculator::Visit( const PatchSignalGenerator* )
+ {
+ // nothing to see here, move along, please
+ return;
+ }
void RunLengthCalculator::Visit( const TrappedElectronGenerator* )
{
diff --git a/Source/Core/LMCRunLengthCalculator.hh b/Source/Core/LMCRunLengthCalculator.hh
index 1658fee3..a945e1d1 100644
--- a/Source/Core/LMCRunLengthCalculator.hh
+++ b/Source/Core/LMCRunLengthCalculator.hh
@@ -90,6 +90,7 @@ namespace locust
private:
void Visit( const KassSignalGenerator* );
void Visit( const FreeFieldSignalGenerator* );
+ void Visit( const PatchSignalGenerator* );
void Visit( const TrappedElectronGenerator* );
void Visit( const GaussianNoiseGenerator* );
void Visit( const TestSignalGenerator* );
diff --git a/Source/Core/LMCVisitor.hh b/Source/Core/LMCVisitor.hh
index cb358583..0c8e8d14 100644
--- a/Source/Core/LMCVisitor.hh
+++ b/Source/Core/LMCVisitor.hh
@@ -19,6 +19,7 @@ namespace locust
class TrappedElectronGenerator;
class KassSignalGenerator;
class FreeFieldSignalGenerator;
+ class PatchSignalGenerator;
class DecimateSignalGenerator;
class GeneratorVisitor
@@ -31,6 +32,7 @@ namespace locust
virtual void Visit( const KassSignalGenerator* ) = 0;
virtual void Visit( const FreeFieldSignalGenerator* ) = 0;
+ virtual void Visit( const PatchSignalGenerator* ) = 0;
virtual void Visit( const GaussianNoiseGenerator* ) = 0;
virtual void Visit( const TrappedElectronGenerator* ) = 0;
virtual void Visit( const TestSignalGenerator* ) = 0;
diff --git a/Source/Generators/LMCDigitizer.cc b/Source/Generators/LMCDigitizer.cc
index d85af641..f95fea70 100644
--- a/Source/Generators/LMCDigitizer.cc
+++ b/Source/Generators/LMCDigitizer.cc
@@ -29,7 +29,7 @@ namespace locust
// get_calib_params( 8, 1, -3.e-6, 6.e-6, false, &fParams ); // if Gaussian noise is included.
- get_calib_params( 8, 1, -1.e-7, 2.e-7, false, &fParams ); // if Gaussian noise is not included.
+ get_calib_params( 8, 1, -1.e-8, 2.e-8, false, &fParams ); // if Gaussian noise is not included.
}
Digitizer::~Digitizer()
diff --git a/Source/Generators/LMCFreeFieldSignalGenerator.cc b/Source/Generators/LMCFreeFieldSignalGenerator.cc
index ba602337..915ea87d 100644
--- a/Source/Generators/LMCFreeFieldSignalGenerator.cc
+++ b/Source/Generators/LMCFreeFieldSignalGenerator.cc
@@ -98,20 +98,17 @@ namespace locust
return;
}
-
static bool ReceivedKassReady()
{
- if( !fKassEventReady)
- {
- std::unique_lock< std::mutex >tLock( fKassReadyMutex );
- fKassReadyCondition.wait( tLock );
- printf("LMC Got the fKassReadyCondition signal\n");
- }
+ 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");
return true;
}
-
double m(double angle, double x0, double y0)
{
angle = LMCConst::Pi() * angle / 180.;
diff --git a/Source/Generators/LMCGaussianNoiseGenerator.cc b/Source/Generators/LMCGaussianNoiseGenerator.cc
index f6e5efd9..3299adb7 100644
--- a/Source/Generators/LMCGaussianNoiseGenerator.cc
+++ b/Source/Generators/LMCGaussianNoiseGenerator.cc
@@ -8,6 +8,8 @@
#include "LMCGaussianNoiseGenerator.hh"
#include "logger.hh"
+#include "LMCSimulationController.hh"
+
using std::string;
@@ -22,6 +24,7 @@ namespace locust
fDoGenerateFunc( &GaussianNoiseGenerator::DoGenerateFreq ),
fMean( 0. ),
fSigma( 1. ),
+ fUniDist( 0., 360. ),
fNormDist( fMean, fSigma )
{
fRequiredSignalState = Signal::kFreq;
@@ -103,6 +106,7 @@ namespace locust
void GaussianNoiseGenerator::SetMeanAndSigma( double aMean, double aSigma )
{
fNormDist = std::normal_distribution< double >( aMean, aSigma );
+ fUniDist = std::uniform_real_distribution< double >(0.,360.);
fMean = aMean;
fSigma = aSigma;
return;
@@ -140,10 +144,20 @@ namespace locust
bool GaussianNoiseGenerator::DoGenerateTime( Signal* aSignal )
{
+ SimulationController SimulationController1;
+ const unsigned nchannels = SimulationController1.GetNChannels();
+ double phi = 0.; // voltage phase
+ double mag = 0.; // voltage mag
+
+ for (int ch=0; chTimeSize(); ++index )
{
- aSignal->SignalTimeComplex()[index][0] += fNormDist( fRNG );
- aSignal->SignalTimeComplex()[index][1] += fNormDist( fRNG );
+ phi = fUniDist( fRNG );
+ mag = fNormDist( fRNG );
+ aSignal->SignalTimeComplex()[ch*aSignal->TimeSize() + index][0] += sqrt(50.)* mag * cos(phi*LMCConst::Pi()/180.);
+ aSignal->SignalTimeComplex()[ch*aSignal->TimeSize() + index][1] += sqrt(50.)* mag * sin(phi*LMCConst::Pi()/180.);
+ }
}
return true;
diff --git a/Source/Generators/LMCGaussianNoiseGenerator.hh b/Source/Generators/LMCGaussianNoiseGenerator.hh
index f4a31683..20b28cbe 100644
--- a/Source/Generators/LMCGaussianNoiseGenerator.hh
+++ b/Source/Generators/LMCGaussianNoiseGenerator.hh
@@ -10,6 +10,8 @@
#include "LMCGenerator.hh"
#include "LMCRunLengthCalculator.hh"
+#include "LMCConst.hh"
+
#include
@@ -72,6 +74,8 @@ namespace locust
double fSigma;
mutable std::normal_distribution< double > fNormDist;
+ mutable std::uniform_real_distribution fUniDist;
+
};
diff --git a/Source/Generators/LMCKassSignalGenerator.cc b/Source/Generators/LMCKassSignalGenerator.cc
index 7dee26d0..e040f3fc 100644
--- a/Source/Generators/LMCKassSignalGenerator.cc
+++ b/Source/Generators/LMCKassSignalGenerator.cc
@@ -151,7 +151,17 @@ namespace locust
{
// initialize phases.
phi_t1 = 2.*LMCConst::Pi()*(CENTER_TO_ANTENNA - tPositionZ) / (tGroupVelocity / tDopplerFrequencyAntenna);
- phi_t2 = 2.*LMCConst::Pi()*(tPositionZ + 2.*CENTER_TO_SHORT + CENTER_TO_ANTENNA) / (tGroupVelocity / tDopplerFrequencyShort);
+ // printf("center_to_antenna is %f and tPositionZ is %f\n", CENTER_TO_ANTENNA, tPositionZ);
+ if (Project8Phase==1)
+ {
+ phi_t2 = LMCConst::Pi()/2. + 2.*LMCConst::Pi()*(CENTER_TO_SHORT + CENTER_TO_ANTENNA) /
+ (tGroupVelocity / tDopplerFrequencyShort); // phase of reflected field at antenna.
+ }
+ if (Project8Phase==2)
+ {
+ phi_t2 = 2.*LMCConst::Pi()*(tPositionZ + 2.*CENTER_TO_SHORT + CENTER_TO_ANTENNA) /
+ (tGroupVelocity / tDopplerFrequencyShort); // terminator in
+ }
EventStartTime = (double)index/RunLengthCalculator1.GetAcquisitionRate()/1.e6/aSignal->DecimationFactor();
EventToFile = false;
}
@@ -159,8 +169,7 @@ namespace locust
if ((tPitchAngle>0.)&&(EventToFile==false))
{
- fprintf(fp, "%10.4g %g\n", EventStartTime, tPitchAngle);
-// printf("start time %g and pitch angle %g\n", EventStartTime, tPitchAngle);
+ fprintf(fp, "%10.4g %g\n", EventStartTime, tPitchAngle);
EventToFile = true;
}
@@ -174,13 +183,13 @@ namespace locust
RealVoltage2 = cos( phi_t2 - phiLO_t ); // + cos( phi_t2 + phiLO_t ));
ImagVoltage2 = sin( phi_t2 - phiLO_t ); // + cos( phi_t2 + phiLO_t - PI/2.));
- //RealVoltage2 = 0.; // take out short?
- //ImagVoltage2 = 0.; // take out short?
if (Project8Phase == 2)
{
- aSignal->LongSignalTimeComplex()[ index ][0] += TE11ModeExcitation() * sqrt(tLarmorPower) * (RealVoltage1 + RealVoltage2);
- aSignal->LongSignalTimeComplex()[ index ][1] += TE11ModeExcitation() * sqrt(tLarmorPower) * (ImagVoltage1 + ImagVoltage2);
+ RealVoltage2 *= 0.03; // replace short with terminator.
+ ImagVoltage2 *= 0.03; // replace short with terminator.
+ aSignal->LongSignalTimeComplex()[ index ][0] += sqrt(50.)*TE11ModeExcitation() * sqrt(tLarmorPower) * (RealVoltage1 + RealVoltage2);
+ aSignal->LongSignalTimeComplex()[ index ][1] += sqrt(50.)*TE11ModeExcitation() * sqrt(tLarmorPower) * (ImagVoltage1 + ImagVoltage2);
}
else if (Project8Phase == 1)
{ // assume 50 ohm impedance
@@ -190,9 +199,10 @@ namespace locust
}
- /*
+ /*
printf("driving antenna, ModeExcitation is %g\n\n", TE11ModeExcitation());
printf("Realvoltage1 is %g and Realvoltage2 is %g\n", RealVoltage1, RealVoltage2);
+ printf("IMagVoltage1 is %g and ImagVoltage2 is %g\n", ImagVoltage1, ImagVoltage2);
printf("Locust says: signal %d is %g and zposition is %g and zvelocity is %g and sqrtLarmorPower is %g and "
" fcyc is %.10g and tDopplerFrequency is %g and GammaZ is %.10g\n\n\n",
index, aSignal->LongSignalTimeComplex()[ index ][0], tPositionZ, tVelocityZ, pow(tLarmorPower,0.5), tCyclotronFrequency, tDopplerFrequencyAntenna, tGammaZ);
@@ -200,7 +210,7 @@ namespace locust
printf("fLO_Frequency is %g\n", fLO_Frequency); getchar();
-
+
*/
t_old += fDigitizerTimeStep; // advance time here instead of in step modifier. This preserves the freefield sampling.
@@ -293,8 +303,10 @@ namespace locust
// 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));
diff --git a/Source/Generators/LMCLowPassFilterFFTGenerator.cc b/Source/Generators/LMCLowPassFilterFFTGenerator.cc
index 067344fe..cab535be 100644
--- a/Source/Generators/LMCLowPassFilterFFTGenerator.cc
+++ b/Source/Generators/LMCLowPassFilterFFTGenerator.cc
@@ -101,7 +101,7 @@ namespace locust
for( unsigned index = 0; index < windowsize; ++index )
{
- // normalize and take the real part of the reverse transform, for digitization.
+ // normalize
aSignal->LongSignalTimeComplex()[ ch*aSignal->TimeSize()*aSignal->DecimationFactor() + nwin*windowsize + index ][0] = SignalComplex[index][0]/norm;
aSignal->LongSignalTimeComplex()[ ch*aSignal->TimeSize()*aSignal->DecimationFactor() + nwin*windowsize + index ][1] = SignalComplex[index][1]/norm;
//if (index>=20000) {printf("filtered signal is %g\n", aSignal->SignalTime()[index]); getchar();}
diff --git a/Source/Generators/LMCPatchSignalGenerator.cc b/Source/Generators/LMCPatchSignalGenerator.cc
new file mode 100644
index 00000000..198f3b63
--- /dev/null
+++ b/Source/Generators/LMCPatchSignalGenerator.cc
@@ -0,0 +1,518 @@
+/*
+ * LMCPatchSignalGenerator.cc
+ *
+ * Created on: Feb 19, 2018
+ * Author: pslocum
+ */
+
+#include "LMCPatchSignalGenerator.hh"
+#include "LMCEventHold.hh"
+#include "LMCRunKassiopeia.hh"
+
+#include "logger.hh"
+#include
+#include
+
+#include
+#include
+
+#include "LMCGlobalsDeclaration.hh"
+#include "LMCHFSSReader.hh"
+#include "LMCSimulationController.hh"
+#include
+
+
+namespace locust
+{
+ LOGGER( lmclog, "PatchSignalGenerator" );
+
+ MT_REGISTER_GENERATOR(PatchSignalGenerator, "patch-signal");
+
+ PatchSignalGenerator::PatchSignalGenerator( const std::string& aName ) :
+ Generator( aName ),
+ fLO_Frequency( 0.),
+ gxml_filename("blank.xml")
+ {
+ fRequiredSignalState = Signal::kTime;
+ }
+
+ PatchSignalGenerator::~PatchSignalGenerator()
+ {
+ }
+
+ bool PatchSignalGenerator::Configure( const scarab::param_node* aParam )
+ {
+ 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" );
+ }
+
+ return true;
+ }
+
+ void PatchSignalGenerator::Accept( GeneratorVisitor* aVisitor ) const
+ {
+ aVisitor->Visit( this );
+ return;
+ }
+
+
+ static void* KassiopeiaInit(const std::string &aFile)
+ {
+ RunKassiopeia *RunKassiopeia1 = new RunKassiopeia;
+ RunKassiopeia1->Run(aFile);
+ delete RunKassiopeia1;
+
+ 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");
+
+ return true;
+ }
+
+
+
+ double PatchSignalGenerator::GetSpaceTimeInterval(const double &aParticleTime, const double &aReceiverTime, const LMCThreeVector &aParticlePosition, const LMCThreeVector &aReceiverPosition )
+ {
+ //return pow(aReceiverTime - aParticleTime,2.) - (aReceiverPosition - aParticlePosition).MagnitudeSquared() / pow(LMCConst::C() , 2.);
+ return aReceiverTime - aParticleTime - (aReceiverPosition - aParticlePosition).Magnitude() / LMCConst::C();
+ }
+
+ double PatchGetStepRoot(const locust::Particle aParticle, double aReceiverTime, LMCThreeVector aReceiverPosition, double aSpaceTimeInterval, const int aStepOrder = 0)
+ {
+ double tRetardedTime = aParticle.GetTime(true); //interpolate!!!
+
+ double c=LMCConst::C();
+
+ if(aStepOrder==0)
+ {
+ double tCorrection = sqrt(fabs(aSpaceTimeInterval));
+ double tSign;
+ (aSpaceTimeInterval > 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 GetMismatchFactor(double f)
+ {
+ // placeholder = 1 - mag(S11)
+ f /= 2.*LMCConst::Pi();
+ // fit to HFSS output
+ double MismatchFactor = 1. - (-5.39e16 / ((f-25.9141e9)*(f-25.9141e9) + 7.23e16) + 0.88);
+
+// printf("dopplerfrequency is %f and mismatchfactor is %g\n", f, MismatchFactor); getchar();
+// double MismatchFactor = 0.85; // punt.
+ return MismatchFactor;
+ }
+
+ double GetAOIFactor(LMCThreeVector IncidentKVector, double PatchPhi)
+ {
+ LMCThreeVector PatchNormalVector;
+ PatchNormalVector.SetComponents(cos(PatchPhi), sin(PatchPhi), 0.0);
+ double AOIFactor = fabs(IncidentKVector.Unit().Dot(PatchNormalVector));
+// printf("cos aoi is %f\n", AOIFactor);
+ return AOIFactor;
+ }
+
+
+ double GetVoltageAmplitude(LMCThreeVector IncidentElectricField, LMCThreeVector IncidentKVector, double PatchPhi, double DopplerFrequency)
+ {
+ double AntennaFactor = 1./600.; // 1/m, placeholder
+ 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();
+ return VoltageAmplitude;
+ }
+
+
+ void AddOnePatchVoltageToStripSum(Signal* aSignal, double VoltageAmplitude, double VoltagePhase, double phi_LO, unsigned channelindex)
+ {
+ 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)
+ {
+
+ locust::Particle tCurrentParticle = fParticleHistory.back();
+ int CurrentIndex;
+ HFSSReader HFRead;
+
+
+ SimulationController SimulationController1;
+ const unsigned nchannels = SimulationController1.GetNChannels();
+ const double dx = 0.00375; // m
+ const double dy = 0.002916; // m
+
+ double PatchPhi = 0.; // azimuthal angle of each antenna element.
+ static double tVoltagePhase[10000] = {0.}; // this is not resetting at the beginning of each event. big problem.
+ static double tVoltagePhaseMark[10000] = {0.}; // this is not resetting at the beginning of each event. big problem.
+
+ static double phi_LO = 0.;
+
+ const int signalSize = aSignal->TimeSize();
+ unsigned patchindex = 0;
+ unsigned channelindex = 0;
+
+ //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.;
+ LMCThreeVector tIncidentElectricField;
+ LMCThreeVector tIncidentKVector;
+ double tAverageDopplerFrequency = 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();
+
+ 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 for this patch.
+ tAverageDopplerFrequency = 0.; // initialize for this patch.
+ tIncidentElectricField.SetComponents(0.,0.,0.); // init
+ tIncidentKVector.SetComponents(0.,0.,0.); // init
+
+ patchindex = ch*NPATCHES_PER_STRIP + z_index; // which patch element (any strip (a.k.a. any channel).
+ channelindex = ch*signalSize*aSignal->DecimationFactor() + index; // which channel and which sample.
+
+// printf("zposition of the patch is %f\n", ZPositionPatch(z_index));
+
+
+
+ for(unsigned i=0;i dtStepSize)
+ {
+ CurrentIndex=FindNode(tRetardedTime,dtStepSize,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 * dy * tECrossH.Dot(tDirection.Unit()) / rReceiver.size() ;// * (fabs(tCurrentParticle.GetPosition(true).Z())<0.01);
+ tIncidentElectricField += tCurrentParticle.CalculateElectricField(rReceiver[i]) / rReceiver.size();
+ tIncidentKVector += tECrossH / rReceiver.size();
+
+ 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);
+ tAverageDopplerFrequency += tDopplerFrequency / rReceiver.size();
+
+ if (tRetardedTime > fDigitizerTimeStep) // if the signal has been present for longer than fDigitizerTimeStep
+ {
+ tVoltagePhase[patchindex]+= tDopplerFrequency * fDigitizerTimeStep / rReceiver.size();
+ }
+ else // if this is the first digitizer sample, the voltage phase doesn't advance for the full dt.
+ {
+// tVoltagePhase[patchindex]+= tDopplerFrequency * fDigitizerTimeStep / rReceiver.size();
+
+ tVoltagePhase[patchindex]+= // compressing
+ tDopplerFrequency * tRetardedTime / 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);
+// printf("after retarding, tVoltagePhase[%d] is %g\n", patchindex, tVoltagePhase[patchindex]); getchar();
+ }
+
+/*
+ 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_index, ch, tDopplerFrequency, tCurrentParticle.GetCyclotronFrequency(), tVelZ, tCosTheta);
+ printf("tVoltagePhase[%d] is %g\n", patchindex, tVoltagePhase[patchindex]);
+ printf("digitizer index is %d\n", index); getchar();
+ }
+*/
+
+
+ } // i, rReceiver.size() loop.
+
+// printf("average Doppler freq is %g\n", tAverageDopplerFrequency);
+// tVoltagePhaseMark[patchindex] = tAverageDopplerFrequency*2.*LMCConst::Pi()*(tReceiverPosition - tCurrentParticle.GetPosition(true)).Magnitude()/LMCConst::C();
+
+ double tVoltageAmplitude = GetVoltageAmplitude(tIncidentElectricField, tIncidentKVector, PatchPhi, tAverageDopplerFrequency);
+ AddOnePatchVoltageToStripSum(aSignal, tVoltageAmplitude, tVoltagePhase[patchindex], phi_LO, channelindex);
+
+// printf("writing phases to file. tVoltagePhase[%d] is %g\n", patchindex, tVoltagePhase[patchindex]);
+// fprintf(fp, "%g %g\n", tVoltagePhaseMark[patchindex], tVoltagePhase[patchindex]);
+
+
+
+ } // z_index waveguide element stepping loop.
+
+// printf("signal %d is %g\n", index, aSignal->LongSignalTimeComplex()[channelindex][0]); getchar();
+
+
+ } // nchannels loop.
+
+ t_old += fDigitizerTimeStep;
+
+ return 0;
+ }
+
+
+ //Return iterator of fParticleHistory particle closest to the time we are evaluating
+ int PatchSignalGenerator::FindNode(double tNew, double dtStepSize, int kIndexOld) const
+ {
+ int tHistorySize = fParticleHistory.size();
+
+ //Make sure we are not out of bounds of array!!!
+ kIndexOld = std::min( std::max(kIndexOld,0) , tHistorySize - 1 );
+
+ double tOld = fParticleHistory[ kIndexOld ].GetTime();
+
+ int kIndexMid=round((tNew-tOld)/dtStepSize) + kIndexOld;
+ kIndexMid = std::min( std::max(kIndexMid,0) , tHistorySize - 1 );
+
+ int kIndexSearchWidth;
+ int kIndexRange[2];
+ std::deque::iterator it;
+
+ for(int i = 0 ; i < 15 ; ++i){
+
+ kIndexSearchWidth = pow( 2 , i );
+ kIndexRange[0] = kIndexMid - kIndexSearchWidth;
+ kIndexRange[1] = kIndexMid + kIndexSearchWidth;
+
+ kIndexRange[0] = std::max(kIndexRange[0], 0 );
+ kIndexRange[1] = std::min(kIndexRange[1], tHistorySize - 1);
+
+ if( tNew >= fParticleHistory[ kIndexRange[0] ].GetTime() && tNew <= fParticleHistory[ kIndexRange[1] ].GetTime())
+ {
+ //Get iterator pointing to particle step closest to tNew
+ it = std::upper_bound( fParticleHistory.begin() , fParticleHistory.end() , tNew, [] (const double &a , const locust::Particle &b) { return a < b.GetTime();} );
+ break;
+ }
+ }
+
+ int tNodeIndex = it - fParticleHistory.begin();
+ //if(tNodeIndex < 0 || tNodeIndex > (tHistorySize - 1))
+ //{
+ // LERROR(lmclog, "OUT OF INDEX SEARCH");
+ //}
+
+ return tNodeIndex;
+ }
+
+ bool PatchSignalGenerator::DoGenerate( Signal* aSignal )
+ {
+ FILE *fp = fopen("phases.txt", "w");
+
+ //n samples for event spacing.
+ int PreEventCounter = 0;
+ const int NPreEventSamples = 150000;
+
+ std::thread Kassiopeia (KassiopeiaInit, gxml_filename); // spawn new thread
+ fRunInProgress = true;
+
+ for( unsigned index = 0; index < aSignal->DecimationFactor()*aSignal->TimeSize(); ++index )
+ {
+ if ((!fEventInProgress) && (fRunInProgress) && (!fPreEventInProgress))
+ {
+ if (ReceivedKassReady()) fPreEventInProgress = true;
+ }
+
+ 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.
+ }
+ }
+
+
+ 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
+ }
+ tLock.unlock();
+ }
+
+ } // for loop
+
+ // trigger any remaining events in Kassiopeia so that its thread can finish.
+ while (fRunInProgress)
+ {
+ if (fRunInProgress)
+ {
+ std::this_thread::sleep_for(std::chrono::milliseconds(100));
+ if (!fEventInProgress)
+ if (ReceivedKassReady())
+ WakeBeforeEvent();
+ }
+ }
+
+ Kassiopeia.join(); // wait for Kassiopeia to finish.
+
+ fclose(fp);
+
+
+ return true;
+ }
+
+} /* namespace locust */
diff --git a/Source/Generators/LMCPatchSignalGenerator.hh b/Source/Generators/LMCPatchSignalGenerator.hh
new file mode 100644
index 00000000..1d9f5998
--- /dev/null
+++ b/Source/Generators/LMCPatchSignalGenerator.hh
@@ -0,0 +1,69 @@
+/*
+ * LMCPatchSignalGenerator.hh
+ *
+ * Created on: Feb 19, 2018
+ * Author: pslocum
+ */
+
+#ifndef LMCPATCHSIGNALGENERATOR_HH_
+#define LMCPATCHSIGNALGENERATOR_HH_
+
+#include "LMCThreeVector.hh"
+#include "LMCGenerator.hh"
+
+namespace locust
+{
+
+ /*!
+ @class PatchSignalGenerator
+ @author P. L. Slocum
+
+ @brief Generate signal in free space(without wave guide) for phase III and detect with patch array.
+
+ @details
+ Operates in time space
+
+ Configuration name: "patch-signal"
+
+ Available configuration options:
+ - "param-name": type -- Description
+ - "lo-frequency" : double -- the special value tuned down by the local oscillator, e.g., the 24.something giga hertz.
+ - "xml-filename" : std::string -- the name of the xml locust config file.
+
+
+ */
+ class PatchSignalGenerator : public Generator
+ {
+ public:
+
+ PatchSignalGenerator( const std::string& aName = "patch-signal" );
+ virtual ~PatchSignalGenerator();
+
+ bool Configure( const scarab::param_node* aNode );
+
+ void Accept( GeneratorVisitor* aVisitor ) const;
+
+
+ private:
+ std::vector rReceiver; //Vector that contains 3D position of all points at which the fields are evaluated (ie. along receiver surface)
+ 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;
+
+
+ bool DoGenerate( Signal* aSignal );
+ void* DriveAntenna(int PreEventCounter, unsigned index, Signal* aSignal, FILE *fp);
+
+ int FindNode(double tNew, double dtStepSize, int IndexOld) const;
+ double GetSpaceTimeInterval(const double &aParticleTime, const double &aReceiverTime, const LMCThreeVector &aParticlePosition, const LMCThreeVector &aReceiverPosition );
+
+
+
+
+
+ };
+
+} /* namespace locust */
+
+#endif /* LMCPATCHSIGNALGENERATOR_HH_ */
diff --git a/Source/Kassiopeia/LMCCyclotronRadiationExtractor.cc b/Source/Kassiopeia/LMCCyclotronRadiationExtractor.cc
index ba792cad..7743a081 100644
--- a/Source/Kassiopeia/LMCCyclotronRadiationExtractor.cc
+++ b/Source/Kassiopeia/LMCCyclotronRadiationExtractor.cc
@@ -38,6 +38,16 @@ namespace locust
{
fP8Phase = P8Phase;
Project8Phase = P8Phase;
+ if (P8Phase==1)
+ {
+ CENTER_TO_SHORT = 0.047; // m
+ CENTER_TO_ANTENNA = 0.045; // m
+ }
+ if (P8Phase==2)
+ {
+ CENTER_TO_SHORT = 0.075; // m
+ CENTER_TO_ANTENNA = 0.075; // m
+ }
}
bool CyclotronRadiationExtractor::ExecutePreStepModification( KSParticle& anInitialParticle, KSParticleQueue& aQueue )
@@ -65,7 +75,7 @@ namespace locust
double CyclotronRadiationExtractor::GetGroupVelocityTE01(KSParticle& aFinalParticle) // Phase 1
{
- double SpeedOfLight = 2.99792458e8; // m/s
+ double SpeedOfLight = LMCConst::C(); // m/s
double CutOffFrequency = SpeedOfLight * LMCConst::Pi() / 10.668e-3; // a in m
double fcyc = aFinalParticle.GetCyclotronFrequency();
double GroupVelocity = SpeedOfLight * pow( 1. - pow(CutOffFrequency/(2.*LMCConst::Pi()*fcyc), 2.) , 0.5);
@@ -84,7 +94,7 @@ namespace locust
double r = sqrt( x * x + y * y);
double coupling = 119116./168.2 * 2./LMCConst::Pi() * 4./(2.*LMCConst::Pi()) / kc/2. * ( (j0(kc*r) - jn(2,kc*r)) +
(j0(kc*r) + jn(2, kc*r)) );
- return coupling;
+ return coupling*coupling;
}
double CyclotronRadiationExtractor::GetCouplingFactorTM01(KSParticle& aFinalParticle)
@@ -94,7 +104,7 @@ namespace locust
double y = aFinalParticle.GetPosition().GetY();
double r = sqrt(x*x + y*y);
double coupling = 146876.5/168.2 * 2./LMCConst::Pi() * 4./(2.*LMCConst::Pi()) / kc * j1(kc*r);
- return coupling;
+ return coupling*coupling;
}
double CyclotronRadiationExtractor::GetCouplingFactorTE01(KSParticle& aFinalParticle) // Phase 1
@@ -102,7 +112,7 @@ namespace locust
double dim1_wr42 = 10.668e-3; // a in m
double x = aFinalParticle.GetPosition().GetX() + dim1_wr42/2.;
double coupling = 0.63*sin(LMCConst::Pi()*x/dim1_wr42); // avg over cyclotron orbit.
- return coupling;
+ return coupling*coupling;
}
@@ -115,9 +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_short = 2.*LMCConst::Pi()*2.*(zPosition+CENTER_TO_SHORT)/(GroupVelocity/fprime_short);
-// double FieldFromShort = cos(phi_short); // no resonant enhancement.
- double FieldFromShort = cos(0.) + cos(phi_short); // yes resonant enhancement.
+ double phi_shortTE01 = LMCConst::Pi()/2. + 2.*LMCConst::Pi()*(zPosition+CENTER_TO_SHORT)/(GroupVelocity/fprime_short); // phase of reflected field at position of electron.
+// double FieldFromShort = cos(phi_shortTM01); // no resonant enhancement.
+ double FieldFromShort = cos(0.) + cos(phi_shortTE01); // yes resonant enhancement.
return FieldFromShort; // Phase 1
@@ -136,9 +146,9 @@ namespace locust
double TE11FieldAfterOneBounce = 0.;
double phi_shortTE11 = 0.;
- phi_shortTE11 = 2.*LMCConst::Pi()*2.*(zPosition+CENTER_TO_SHORT)/(GroupVelocity/fprime_short);
+ // Cu boundary condition gives PI/2 phase advancement to short.
+ phi_shortTE11 = LMCConst::Pi()/2. + 2.*LMCConst::Pi()*(zPosition+CENTER_TO_SHORT)/(GroupVelocity/fprime_short);
TE11FieldAfterOneBounce = cos(0.) + cos(phi_shortTE11);
- //printf("TE11FieldAfterOneBounce is %f\n", TE11FieldAfterOneBounce);
return TE11FieldAfterOneBounce;
}
@@ -162,61 +172,99 @@ namespace locust
double FieldFromShort=0.; // first doppler shift
double FieldFromPolarizer=0.; // other doppler shift
double TM01FieldAfterBounces = 0.;
- int nbounces = 20;
- double time_decay = 1.;
- double reflection_coefficient = 1.0;
+ int nbounces = 10; // even number please.
double phi_shortTM01 = 0.;
double phi_polarizerTM01 = 0.;
-
- //printf("TM01 l1 is %g and l2 is %g\n", GroupVelocity/fprime_short, GroupVelocity/fprime_polarizer); getchar();
-
-// if ((phi_shortTM01[0] == 0.)||(0==0)) // if the event has just started, or always.
-// {
- phi_shortTM01 = 2.* LMCConst::Pi() *2.*(tPositionZ+CENTER_TO_SHORT)/lambda_short; // starting phi after 0th bounce.
- phi_polarizerTM01 = 2.*LMCConst::Pi()*2.*(CENTER_TO_ANTENNA - tPositionZ)/lambda_polarizer + LMCConst::Pi(); // starting phi after 0th bounce.
+ double dphi = 0.;
// printf("phi_shortTM01[0] is %.10g and LMCConst::Pi() is %.10g and z is %.10g and lambda is %.10g\n", phi_shortTM01[0], LMCConst::Pi(), tPositionZ, lambda_short);
- FieldFromShort = cos(0.) + 1./1.4*reflection_coefficient*cos(phi_shortTM01); // starting field, after 0th bounce.
- FieldFromPolarizer = 1./1.4*reflection_coefficient*cos(phi_polarizerTM01); // starting field, after 0th bounce.
+ FieldFromShort = cos(0.); // starting field, after no reflections. divide by 2 later.
+ FieldFromPolarizer = cos(0.); // starting field, after no reflections. divide by 2 later.
- for (int i=0; i 0 or < 0.
- double DampingFactorTM01 = CouplingFactorTM01*(1. - TM01FieldAfterBounces*TM01FieldAfterBounces); // can be > 0 or < 0.
- double DampingFactor = DampingFactorTM01 + DampingFactorTE11;
+ // double DampingFactorTE11 = CouplingFactorTE11*(1. - TE11FieldFromShort*TE11FieldFromShort); // can be > 0 or < 0.
+ double DampingFactorTM01 = CouplingFactorTM01*(1. - TM01FieldWithTerminator*TM01FieldWithTerminator); // can be > 0 or < 0.
+ double DampingFactor = DampingFactorTM01;
return DampingFactor;
}
@@ -295,18 +343,11 @@ namespace locust
bool CyclotronRadiationExtractor::ExecutePostStepModification( KSParticle& anInitialParticle, KSParticle& aFinalParticle, KSParticleQueue& aQueue )
{
-// printf("fcyc before coupling is %.9g and Bz is %.10g\n\n", aFinalParticle.GetCyclotronFrequency(), aFinalParticle.GetMagneticField().GetZ());
-
- //printf("pre step kinetic energy - 4.84338e-15 is %g\n", anInitialParticle.GetKineticEnergy()- 4.84338e-15); //getchar();
- //printf("post step kinetic energy - 4.84338e-15 is %g\n", aFinalParticle.GetKineticEnergy()- 4.84338e-15); //getchar();
-
double DeltaE=0.;
-
if(fP8Phase==1)
{
// adjust power with reflections.
DeltaE = GetDampingFactorPhase1(anInitialParticle, aFinalParticle)*(aFinalParticle.GetKineticEnergy() - anInitialParticle.GetKineticEnergy());
- //printf("poststep says DeltaE is %g\n", DeltaE);
aFinalParticle.SetKineticEnergy((aFinalParticle.GetKineticEnergy() - DeltaE));
}
if(fP8Phase==2)
diff --git a/Source/Kassiopeia/LMCCyclotronRadiationExtractor.hh b/Source/Kassiopeia/LMCCyclotronRadiationExtractor.hh
index db8865ea..0095c946 100644
--- a/Source/Kassiopeia/LMCCyclotronRadiationExtractor.hh
+++ b/Source/Kassiopeia/LMCCyclotronRadiationExtractor.hh
@@ -61,6 +61,7 @@ namespace locust
double GetCouplingFactorTM01(Kassiopeia::KSParticle& aFinalParticle);
double GetCouplingFactorTE01(Kassiopeia::KSParticle& aFinalParticle);
double GetTM01FieldAfterBounces(Kassiopeia::KSParticle& anInitialParticle, Kassiopeia::KSParticle& aFinalParticle);
+ double GetTM01FieldWithTerminator(Kassiopeia::KSParticle& anInitialParticle, Kassiopeia::KSParticle& aFinalParticle);
double GetTE11FieldAfterOneBounce(Kassiopeia::KSParticle& anInitialParticle, Kassiopeia::KSParticle& aFinalParticle);
double GetTE01FieldAfterOneBounce(Kassiopeia::KSParticle& anInitialParticle, Kassiopeia::KSParticle& aFinalParticle);
locust::Particle ExtractKassiopeiaParticle( Kassiopeia::KSParticle &anInitialParticle, Kassiopeia::KSParticle &aFinalParticle);