Merge pull request #89 from PyCOMPLETE/feature/limit_distribution_generation

Feature/limit distribution generation

Author: aoeftiger

PyHEADTAIL/particles/generators.py

@@ -165,21 +165,30 @@ def _RF_bucket_dist(n_particles):
 
 
 def cut_distribution(distribution, is_accepted):
-    """ Generate coordinates according to some distribution inside the region
-    specified by is_accepted. (Wrapper for distributions, based on RF_cut..)
+    """Generate coordinates according to some distribution inside the
+    region specified by where the function is_accepted returns 1.
+    (Wrapper for distributions, based on RF_cut..)
     Args:
-        is_accepted: function taking two parameters (z,dp) and returns an
-                     array [0,1] specifying whether the coordinate lies
-                     inside the desired phase space volume. Possible sources
-                     are the rfbucket.make_is_accepted(...)
-        sigma_z: std of the normally distributed z coordinates
-        sigma_dp: std of the normally distributed dp coordinates
+        distribution: a function which takes the n_particles as a
+                      parameter and returns a list-like object
+                      containing a 2D phase space. result[0] should
+                      stand for the spatial, result[1] for the momentum
+                      coordinate
+        is_accepted: function taking two parameters (z, dp)
+                     [vectorised as arrays] and returning a boolean
+                     specifying whether the coordinate lies
+                     inside the desired phase space volume. A possible
+                     source to provide such an is_accepted function
+                     is the RFBucket.make_is_accepted or
+                     generators.make_is_accepted_within_n_sigma .
     Returns:
         A matcher with the specified bucket properties (closure)
     """
     def _cut_distribution(n_particles):
-        '''Removes all particles for which is_accepted(x,x') is false
-        and redistributes them until all lie inside the bucket
+        '''Regenerates all particles which fall outside a previously
+        specified phase space region (via the function is_accepted
+        in generators.cut_distribution) until all generated particles
+        have valid coordinates and momenta.
         '''
         z = np.zeros(n_particles)
         dp = np.zeros(n_particles)
@@ -196,6 +205,32 @@ def _cut_distribution(n_particles):
         return [z, dp]
     return _cut_distribution
 
+def make_is_accepted_within_n_sigma(epsn_rms, limit_n_rms, twiss_beta=1):
+    '''Closure creating an is_accepted function (e.g. for
+    cut_distribution). The is_accepted function will return whether
+    the canonical coordinate and momentum pair lies within the phase
+    space region limited by the action value limit_n_rms * epsn_rms.
+
+    Coordinate u and momentum up are assumed to be connected to the
+    amplitude J via the twiss_beta value,
+    J = sqrt(u^2 + twiss_beta^2 up^2) .
+    The amplitude is required to be below the limit to be accepted,
+    J < limit_n_rms * epsn_rms.
+    The usual use case will be generating u and up in normalised Floquet
+    space (i.e. before the normalised phase space coordinates
+    get matched to the optics or longitudinal eta and Qs).
+    In this case, twiss_beta takes the default value 1 in normalised
+    Floquet space. Consequently, the 1 sigma RMS reference value
+    epsn_rms corresponds to the normalised 1 sigma RMS emittance
+    (i.e. amounting to beam.epsn_x() and beam.epsn_y() in the transverse
+    plane, and beam.epsn_z()/4 in the longitudinal plane).
+    '''
+    threshold_amplitude_squared = (limit_n_rms * epsn_rms)**2
+    def is_accepted(u, up):
+        Jsq = u**2 + (twiss_beta * up)**2
+        return Jsq < threshold_amplitude_squared
+    return is_accepted
+
 
 class ParticleGenerator(Printing):
     '''Factory to generate Particle instances according to distributions

PyHEADTAIL/particles/rfbucket_matching.py

@@ -110,7 +110,7 @@ def error_from_target_epsn(ec):
         emittance = self._compute_emittance(self.rfbucket, self.psi)
         self.prints('--> Emittance: ' + str(emittance))
         sigma = self._compute_sigma(self.rfbucket, self.psi)
-        self.prints('--> Bunch length:' + str(sigma))
+        self.prints('--> Bunch length: ' + str(sigma))
 
     def psi_for_bunchlength_newton_method(self, sigma):
         # Maximum bunch length