Increased slightly tolerances, working on radiation in multipole

Author: giadarol

tests/test_full_rings.py

@@ -23,8 +23,8 @@ def test_full_rings(element_by_element=False):
             ]
 
      tolerances_10_turns = [
-                    (1e-9, 1e-11),
-                    (1e-9, 1e-11),
+                    (1e-9, 3e-11),
+                    (1e-9, 3e-11),
                     (2e-8, 7e-9)]
 
      for icase, fname_line_particles in enumerate(test_fnames):

xtrack/beam_elements/elements.py

@@ -142,6 +142,7 @@ class Multipole(BeamElement):
         'length': xo.Float64,
         'hxl': xo.Float64,
         'hyl': xo.Float64,
+        'radiation_flag': xo.Int64,
         'bal': xo.Float64[:],
         }
 
@@ -209,7 +210,8 @@ def ksl(self):
         #return self.bal[idx + 1] * factorial(idx // 2, exact=True)
 
 Multipole.XoStruct.extra_sources = [
-    _pkg_root.joinpath('headers/rng.h'),
+    _pkg_root.joinpath('random_number_generator/rng_src/base_rng.h'),
+    _pkg_root.joinpath('random_number_generator/rng_src/local_particle_rng.h'),
     _pkg_root.joinpath('beam_elements/elements_src/multipole.h')]
 
 

xtrack/beam_elements/elements_src/multipole.h

@@ -4,6 +4,8 @@
 /*gpufun*/
 void Multipole_track_local_particle(MultipoleData el, LocalParticle* part0){
 
+   int64_t radiation_flag = MultipoleData_get_radiation_flag(el);
+
    //start_per_particle_block (part0->part)
     	int64_t order = MultipoleData_get_order(el);
     	int64_t index_x = 2 * order;
@@ -36,6 +38,7 @@ void Multipole_track_local_particle(MultipoleData el, LocalParticle* part0){
 
 	// compute the average energy loss by synchrotron radiation
 	double const length = MultipoleData_get_length(el); // m
+        //if (radiation_flag>0 && length!=0.0) {
 	if (length!=0.0) {
 	  double const h = hypot(dpx, dpy) / length; // 1/m, 1/rho, curvature
 	  double const p0c = LocalParticle_get_p0c(part); // eV
@@ -47,27 +50,14 @@ void Multipole_track_local_particle(MultipoleData el, LocalParticle* part0){
 	  double const q0 = LocalParticle_get_q0(part); // e
 	  // e^2 / 4 pi epsilon0 eV = (1 / 694461541.7756249) m
 	  double const classical_radius = q0*q0 / m0 / 694461541.7756249; // m, classical electromagnetic radius
-	  double const eloss = 0 * 2.0 / 3.0 * classical_radius*length * beta_gamma*beta_gamma*beta_gamma * h*h * energy; // eV
+	  double const eloss = 2.0 / 3.0 * classical_radius*length * beta_gamma*beta_gamma*beta_gamma * h*h * energy; // eV
 
 	  // apply the energy kick
-	  LocalParticle_add_to_energy(part, -eloss);
+	  //LocalParticle_add_to_energy(part, -eloss);
+	  LocalParticle_add_to_energy(part, 0.);
 
-	  // random number generator test
-	  uint32_t s1, s2, s3;
-	  rng_set(&s1, &s2, &s3, 1234);
+	  double r = LocalParticle_generate_random_double(part); 
 
-	  LocalParticle_set___rng_s1(part, s1);
-	  LocalParticle_set___rng_s2(part, s2);
-	  LocalParticle_set___rng_s3(part, s3);
-	  
-	  /* uint32_t */ s1 = LocalParticle_get___rng_s1(part);
-	  /* uint32_t */ s2 = LocalParticle_get___rng_s2(part);
-	  /* uint32_t */ s3 = LocalParticle_get___rng_s3(part);
-
-	  for (int i=0; i<3; i++) {
-	    double r = rng_get(&s1, &s2, &s3);
-	    printf("rnd = %g\n", r);
-	  }
 	}
 
     	if( ( hxl > 0) || ( hyl > 0) || ( hxl < 0 ) || ( hyl < 0 ) )