Minimal working example

src/Classic/AbsBeamline/Ring.cpp

@@ -120,7 +120,7 @@ bool Ring::apply(const size_t &id, const double &t, Vector_t &E,
 
 bool Ring::apply(const Vector_t &R, const Vector_t &P,
                  const double &t, Vector_t &E, Vector_t &B) {
-    std::cerr << "Ring::apply " << this << std::endl;
+    std::cerr << "Ring::apply " << this << " " << R << std::endl;
     B = Vector_t(0.0, 0.0, 0.0);
     E = Vector_t(0.0, 0.0, 0.0);
     std::vector<RingSection*> sections = getSectionsAt(R); // I think this doesn't actually use R -DW

src/Classic/AbsBeamline/VerticalFFAMagnet.cpp

@@ -29,6 +29,8 @@ VerticalFFAMagnet::VerticalFFAMagnet(const VerticalFFAMagnet &right)
         }
         std::cerr << "VerticalFFAMagnet::CopyCtor 0 " << this << std::endl;
         std::cerr << "   right k " << right.k_m << " this k " << k_m << std::endl;
+        std::cerr << "   right ef " << right.endField_m->function(0.1, 1) 
+                  << " this ef " << endField_m->function(0.1, 1) << std::endl;
         RefPartBunch_m = right.RefPartBunch_m;
         std::cerr << "VerticalFFAMagnet::CopyCtor 1" << std::endl;
 }
@@ -39,10 +41,10 @@ VerticalFFAMagnet::~VerticalFFAMagnet() {
 
 ElementBase* VerticalFFAMagnet::clone() const {
     std::cerr << "VFFA Clone 0" << std::endl;
-    std::cerr << "VFFA " << this << std::endl;
-    std::cerr << "VFFA " << endField_m.get() << " " << boundingBox_m.get() << std::endl;
+    std::cerr << "VFFA this: " << this << std::endl;
+    std::cerr << "VFFA endfield: " << endField_m.get() << " bbox: " << boundingBox_m.get() << std::endl;
     VerticalFFAMagnet* magnet = new VerticalFFAMagnet(*this);
-    std::cerr << "VFFA Clone 0a " << std::endl;
+    std::cerr << "VFFA Clone 0a new: " << magnet << std::endl;
     magnet->initialise();
     std::cerr << "VFFA Clone 1" << std::endl;
     return magnet;
@@ -84,7 +86,6 @@ void VerticalFFAMagnet::accept(BeamlineVisitor& visitor) const {
 
 
 bool VerticalFFAMagnet::getFieldValue(const Vector_t &R, Vector_t &B) const {
-    std::cerr << "VFFA field 1 " << this << " " << R << std::endl;
     if (boundingBox_m.get() != nullptr) {
         if (boundingBox_m->isOutside(R)) {
             return true;
@@ -99,6 +100,7 @@ bool VerticalFFAMagnet::getFieldValue(const Vector_t &R, Vector_t &B) const {
     for (size_t i = 0; i < fringeDerivatives.size(); ++i) {
         fringeDerivatives[i] = endField_m->function(zRel, i); // d^i_phi f
     }
+    std::cerr << std::endl;
 
     std::vector<double> x_n(maxOrder_m+1); // x^n
     x_n[0] = 1.; // x^0
@@ -106,7 +108,6 @@ bool VerticalFFAMagnet::getFieldValue(const Vector_t &R, Vector_t &B) const {
         x_n[i] = x_n[i-1]*R[0];
     }
 
-    std::cerr << "VFFA field 2" << std::endl;
     std::vector<double> f_n(maxOrder_m+2, 0.);
     std::vector<double> dz_f_n(maxOrder_m+1, 0.);
     for (size_t n = 0; n < dfCoefficients_m.size(); ++n) {
@@ -114,10 +115,8 @@ bool VerticalFFAMagnet::getFieldValue(const Vector_t &R, Vector_t &B) const {
         for (size_t i = 0; i < coefficients.size(); ++i) {
             f_n[n] += coefficients[i]*fringeDerivatives[i];
             dz_f_n[n] += coefficients[i]*fringeDerivatives[i+1];
-            std::cerr << "    " << i << " " << n << " " << f_n[n] << " " << dz_f_n[n] << " " << coefficients[i] << " " << fringeDerivatives[i+1] << std::endl;
         }
     }
-    std::cerr << "VFFA field 3 " << std::endl;
     f_n[0] = fringeDerivatives[0];
     double bref = Bz_m*exp(k_m*R[1]);
     B[0] = 0.;
@@ -127,9 +126,7 @@ bool VerticalFFAMagnet::getFieldValue(const Vector_t &R, Vector_t &B) const {
         B[0] += bref*f_n[n+1]*(n+1)/k_m*x_n[n];
         B[1] += bref*f_n[n]*x_n[n];
         B[2] += bref*dz_f_n[n]/k_m*x_n[n];
-        std::cerr << "    " << B << " " << k_m << std::endl;
     }
-    std::cerr << "VFFA field 4" << std::endl;
     return false;
 }
 

src/Steppers/RK4.hpp

@@ -14,7 +14,7 @@ bool RK4<FieldFunction, Arguments ...>::doAdvance_m(PartBunchBase<double, 3>* bu
     double x[6];
     
     this->copyTo(bunch->R[i], bunch->P[i], &x[0]);
-    
+    std::cerr << "Entering RK4 " << bunch->R[i] << std::endl;   
     double  deriv1[6];
     double  deriv2[6];
     double  deriv3[6];
@@ -22,7 +22,7 @@ bool RK4<FieldFunction, Arguments ...>::doAdvance_m(PartBunchBase<double, 3>* bu
     double  xtemp[6];
 
     // Evaluate f1 = f(x,t).
-
+    std::cerr << "Calc deriv1" << std::endl;
     bool outOfBound = derivate_m(bunch, x, t, deriv1 , i, args ...);
     if (outOfBound)
         return false;
@@ -35,6 +35,7 @@ bool RK4<FieldFunction, Arguments ...>::doAdvance_m(PartBunchBase<double, 3>* bu
         xtemp[j] = x[j] + half_dt * deriv1[j];
     }
 
+    std::cerr << "Calc deriv2 " << half_dt << " " << xtemp[0] << " " << xtemp[1] << " " << xtemp[2] << " " << xtemp[3] << " " << xtemp[4] << " " << xtemp[5] <<std::endl;
     outOfBound = derivate_m(bunch, xtemp, t_half, deriv2 , i, args ...);
     if (outOfBound)
         return false;
@@ -59,8 +60,12 @@ bool RK4<FieldFunction, Arguments ...>::doAdvance_m(PartBunchBase<double, 3>* bu
     // Return x(t+dt) computed from fourth-order R-K.
     for(int j = 0; j < 6; ++j)
         x[j] += dt / 6.*(deriv1[j] + deriv4[j] + 2.*(deriv2[j] + deriv3[j]));
+
+    for(int j = 0; j < 6; ++j)
+        std::cerr << "Deriv " << j << " " << deriv1[j] << " " << deriv2[j] << " " << deriv3[j] << std::endl;
     
     this->copyFrom(bunch->R[i], bunch->P[i], &x[0]);
+    std::cerr << "Leaving RK4 " << bunch->R[i] << std::endl;   
 
     return true;
 }
@@ -86,11 +91,12 @@ bool RK4<FieldFunction, Arguments ...>::derivate_m(PartBunchBase<double, 3>* bun
         tempR(j) = y[j];
 
     bunch->R[i] = tempR;
+    std::cerr << "RK4 tempR " << tempR << std::endl;
 
     bool outOfBound = this->fieldfunc_m(t, i, externalE, externalB, args ...);
 
     double qtom = bunch->Q[i] / (bunch->M[i] * mass_coeff);   // m^2/s^2/GV
-
+    std::cerr << "RK4 qtom " << qtom << " " <<  bunch->M[i] << " " << mass_coeff << std::endl;
     double tempgamma = sqrt(1 + (y[3] * y[3] + y[4] * y[4] + y[5] * y[5]));
 
     yp[0] = c_mtns / tempgamma * y[3]; // [m/ns]

tests/opal_src/PyOpal/test_track_run.py

@@ -75,7 +75,10 @@ def make_line():
 
 print("Running track_run test")
 beam = PyOpal.beam.Beam()
+beam.mass = 0.938272
+beam.charge = 1.0
 beam.momentum = 0.1
+beam.beam_frequency = 1.0
 beam.number_of_slices = 10
 beam.register()