BLonD v1.5.0. See README.md for details.

Author: htimko

.gitignore

@@ -20,3 +20,9 @@ stats.html
 # _doc/buil folder #
 #######################
 /__doc/build/
+
+# Figure folder #
+#################
+/__EXAMPLE_MAIN_FILES/fig/
+/__EXAMPLE_MAIN_FILES/*.h5
+/__EXAMPLE_MAIN_FILES/null

README.md

@@ -9,11 +9,11 @@ longitudinal beam dynamics with collective effects.
 
 Developers:
 
-Theodoros Argyropoulos (Theodoros.Argyropoulos@cern.ch)
-Alexandre Lasheen (alexandre.lasheen@cern.ch)
-Juan Esteban Muller (juan.fem@cern.ch)
-Danilo Quartullo (danilo.quartullo@cern.ch)
-Helga Timko (Helga.Timko@cern.ch)
+Theodoros Argyropoulos (Theodoros.Argyropoulos (at) cern.ch)
+Alexandre Lasheen (alexandre.lasheen (at) cern.ch)
+Juan Esteban Muller (juan.fem (at) cern.ch)
+Danilo Quartullo (danilo.quartullo (at) cern.ch)
+Helga Timko (Helga.Timko (at) cern.ch)
 
 
 The structure is as follows:
@@ -37,6 +37,13 @@ The structure is as follows:
 VERSION CONTENTS
 ==========
 
+2014-10-13
+v1.5.0 - Phase loop implemented and tested
+	   - New feedback features in plotting and monitoring
+	   - Removed average values (beta etc.) from input_parameters
+	   - Plotting: formatting setup is separate, files made "uniform"
+	   - Calculated synchronous phase now has the same array length as momentum
+
 2014-10-13
 v1.4.0 - The name of the code has been changed from PyLongitudinal to BLond and 
          a new repository has been created especially for this code.

__EXAMPLE_MAIN_FILES/_Acceleration.py

@@ -16,7 +16,7 @@
 
 # Simulation parameters --------------------------------------------------------
 # Bunch parameters
-N_b = 1.e9           # Intensity
+N_b = 1.e9           # Intensity
 N_p = 50000         # Macro-particles
 tau_0 = 0.4          # Initial bunch length, 4 sigma [ns]
 
@@ -55,7 +55,7 @@
 # Define beam and distribution
 beam = Beam(general_params, N_p, N_b)
 longitudinal_bigaussian(general_params, rf_params, beam, tau_0/4, 
-                              xunit = 'ns', reinsertion = 'on')
+                        xunit = 'ns', reinsertion = 'on')
 
 
 print "Beam set and distribution generated..."
@@ -94,7 +94,7 @@
         print "   Gaussian bunch length %.4e rad" %beam.bl_gauss
         print ""
         # In plots, you can choose following units: rad, ns, m  
-        plot_long_phase_space(beam, general_params, rf_params, 0, 0.0001763, -400, 400, separatrix_plot = True)
+        plot_long_phase_space(beam, general_params, rf_params, 0, 0.0001763, -400, 400, separatrix_plot = True)
 
     # Track
     for m in map_:

__EXAMPLE_MAIN_FILES/_Main_long_ps_booster.py

@@ -166,8 +166,7 @@
 
         plot_bunch_length_evol(my_beam, 'beam', general_params, n_turns)
 
-        plot_position_evol(i+1, my_beam, 'beam', general_params, unit = None, dirname = 'fig')
-
+        plot_position_evol(my_beam, 'beam', general_params, i+1, unit = None, dirname = 'fig')
 
 print "Done!"
 

__EXAMPLE_MAIN_FILES/_Phase_loop.py

@@ -0,0 +1,145 @@
+# Example input for longitudinal simulation with phase loop (LHC)
+# No intensity effects
+
+
+import time 
+import numpy as np
+
+from llrf.RF_noise import *
+from input_parameters.general_parameters import *
+from input_parameters.rf_parameters import *
+from trackers.longitudinal_tracker import *
+from llrf.feedbacks import *
+from beams.beams import *
+from beams.longitudinal_distributions import *
+from beams.slices import *
+from monitors.monitors import *
+from longitudinal_plots.plot_settings import *
+from longitudinal_plots.plot_beams import *
+from longitudinal_plots.plot_llrf import *
+from longitudinal_plots.plot_slices import *
+
+
+
+# Simulation parameters --------------------------------------------------------
+# Bunch parameters
+N_b = 1.e9           # Intensity
+N_p = 10001          # Macro-particles
+tau_0 = 0.4          # Initial bunch length, 4 sigma [ns]
+
+# Machine and RF parameters
+C = 26658.883        # Machine circumference [m]
+p_s = 450.e9         # Synchronous momentum [eV]
+h = 35640            # Harmonic number
+V = 6.e6             # RF voltage [eV]
+dphi = 0             # Phase modulation/offset
+gamma_t = 55.759505  # Transition gamma
+alpha = 1./gamma_t/gamma_t        # First order mom. comp. factor
+
+# Tracking details
+N_t = 1001           # Number of turns to track
+dt_plt = 200         # Time steps between plots
+dt_mon = 1           # Time steps between monitoring
+
+
+# Simulation setup -------------------------------------------------------------
+print "Setting up the simulation..."
+print ""
+
+# Define general parameters
+general_params = GeneralParameters(N_t, C, alpha, p_s, 'proton')
+
+# Define RF station parameters, phase loop, and corresponding tracker
+RF_params = RFSectionParameters(general_params, 1, h, V, dphi)
+PL_gain = 1./(5.*general_params.t_rev[0])
+print "PL gain is %.4e 1/s, Trev = %.4e s" %(PL_gain, general_params.t_rev[0])
+PL = PhaseLoop(general_params, RF_params, PL_gain, sampling_frequency = 1, 
+               machine = 'LHC')
+long_tracker = RingAndRFSection(RF_params, PhaseLoop=PL)
+
+print "General and RF parameters set..."
+
+
+# Define beam and distribution
+beam = Beam(general_params, N_p, N_b)
+# Generate new distribution
+longitudinal_bigaussian(general_params, RF_params, beam, tau_0/4, seed=1234,
+                        xunit = 'ns', reinsertion = 'on')
+print "Beam set and distribution generated..."
+
+
+# Need slices for the Gaussian fit; slice for the first plot
+slice_beam = Slices(beam, 100, slicing_coord = 'theta')#, #fit_option = 'gaussian', 
+                    #slice_immediately = 'on')
+
+# Define what to save in file
+bunchmonitor = BunchMonitor('output_data', N_t+1, "Longitudinal", PhaseLoop=PL)#slice_beam, PL)
+
+print "Statistics set..."
+
+# Initialize plots
+PlotSettings().set_plot_format()
+
+# Accelerator map
+map_ = [slice_beam] + [bunchmonitor] +[long_tracker] 
+print "Map set"
+print ""
+
+
+# Initial injection kick/error
+beam.theta += 1.e-5
+#beam.dE *= -1.
+bunchmonitor.track(beam)
+slice_beam.track(beam)
+
+
+# Tracking ---------------------------------------------------------------------
+for i in range(N_t):
+    t0 = time.clock()
+    
+   
+    # Plot has to be done before tracking (at least for cases with separatrix)
+    if (i % dt_plt) == 0:
+        print "Outputting at time step %d..." %i
+        print "   Beam momentum %.6e eV" %beam.momentum
+        print "   Beam gamma %3.3f" %beam.gamma_r
+        print "   Beam beta %3.3f" %beam.beta_r
+        print "   Beam energy %.6e eV" %beam.energy
+        print "   Four-times r.m.s. bunch length %.4e rad" %(4.*beam.sigma_theta)
+        #print "   Gaussian bunch length %.4e rad" %beam.bl_gauss
+        print ""
+        # In plots, you can choose following units: rad, ns, m  
+        plot_long_phase_space(beam, general_params, RF_params, 0, 0.0001763, 
+                              -450, 450, separatrix_plot = True)
+        plot_beam_profile(i, general_params, slice_beam)
+        plot_COM_motion(beam, general_params, RF_params, 'output_data', 
+                        0.8e-4, 1.1e-4, -75, 75, separatrix_plot = False)
+
+
+    # Track
+    slice_beam.track(beam)
+    #slice_beam.mean_theta = beam.theta[0] # For single particle
+    long_tracker.track(beam)
+    bunchmonitor.track(beam)
+        
+    # These plots have to be done after the tracking
+    if (i % dt_plt) == 0 and i > dt_mon:
+        plot_bunch_length_evol(beam, 'output_data', general_params, i, 
+                               output_freq=dt_mon, unit='ns')
+        #plot_bunch_length_evol_gaussian(beam, 'output_data', general_params, 
+        #                                slice_beam, i, output_freq=dt_mon, unit='ns')
+        plot_position_evol(beam, 'output_data', general_params, i,
+                           output_freq=dt_mon, unit = None, style = '.') 
+        plot_PL_phase_corr(PL, 'output_data', i, output_freq=dt_mon)
+        plot_PL_freq_corr(PL, 'output_data', i, output_freq=dt_mon)
+
+    
+    # Define losses according to separatrix and/or longitudinal position
+    #beam.losses_separatrix(general_params, rf_params)
+    #beam.losses_longitudinal_cut(0.28e-4, 0.75e-4)
+
+
+bunchmonitor.h5file.close()
+print "Done!"
+print ""
+

__EXAMPLE_MAIN_FILES/_RFnoise.py

@@ -42,7 +42,8 @@
 # Pre-processing: RF phase noise -----------------------------------------------
 f = np.arange(0, 5.6227612455e+03, 1.12455000e-02)
 spectrum = np.concatenate((1.11100000e-07 * np.ones(4980), np.zeros(495021)))
-noise_t, noise_dphi = PhaseNoise(f, spectrum, seed1=1234, seed2=7564).spectrum_to_phase_noise()
+RFnoise = PhaseNoise(f, spectrum, seed1=1234, seed2=7564)
+RFnoise.spectrum_to_phase_noise()
 
 # Hermitian vs complex FFT (gives the same result)
 # plot_noise_spectrum(f, spectrum, sampling=100)
@@ -62,10 +63,10 @@
 # os.rename('temp/phase_noise.png', 'temp/phase_noise_c.png')
 
 plot_noise_spectrum(f, spectrum, sampling=100)
-plot_phase_noise(noise_t, noise_dphi, sampling=100)
-#plot_phase_noise(noise_t[0:10000], noise_dphi[0:10000], sampling=1)
-print "   Sigma of RF noise is %.4e" %np.std(noise_dphi)
-print "   Time step of RF noise is %.4e" %noise_t[1]
+plot_phase_noise(RFnoise.t, RFnoise.dphi, sampling=100)
+#plot_phase_noise(RFnoise.t[0:10000], RFnoise.dphi[0:10000], sampling=1)
+print "   Sigma of RF noise is %.4e" %np.std(RFnoise.dphi)
+print "   Time step of RF noise is %.4e" %RFnoise.t[1]
 print ""
 
 
@@ -74,11 +75,10 @@
 print ""
 
 # Define general parameters
-general_params = GeneralParameters(N_t, C, alpha, p_s, 
-                                   'proton')
+general_params = GeneralParameters(N_t, C, alpha, p_s,'proton')
 
 # Define RF station parameters and corresponding tracker
-rf_params = RFSectionParameters(general_params, 1, h, V, dphi)
+rf_params = RFSectionParameters(general_params, 1, h, V, RFnoise.dphi[0:N_t+1])
 long_tracker = RingAndRFSection(rf_params)
 
 print "General and RF parameters set..."
@@ -96,7 +96,8 @@
 
 
 # Need slices for the Gaussian fit; slice for the first plot
-slice_beam = Slices(beam, 100, slicing_coord = 'theta', fit_option = 'gaussian', slice_immediately = 'on')
+slice_beam = Slices(beam, 100, slicing_coord = 'theta', fit_option = 'gaussian', 
+                    slice_immediately = 'on')
 
 # Define what to save in file
 bunchmonitor = BunchMonitor('output_data', N_t+1, "Longitudinal", slice_beam)

beams/beams.py

@@ -8,7 +8,7 @@
 import numpy as np
 import warnings
 from scipy.constants import c, e, m_p
-
+import cython_functions.stats as cp
 from scipy.optimize import curve_fit
 from trackers.longitudinal_utilities import is_in_separatrix
 from scipy import ndimage

beams/longitudinal_distributions.py

@@ -677,6 +677,7 @@ def longitudinal_bigaussian(GeneralParameters, RFSectionParameters, beam,
         warnings.warn("longitudinal_bigaussian for multiple RF is not yet implemented")
 
     counter = RFSectionParameters.counter[0]
+    
     harmonic = RFSectionParameters.harmonic[0,counter]
     energy = RFSectionParameters.energy[counter]
     beta = RFSectionParameters.beta_r[counter]

beams/slices.py

@@ -73,7 +73,7 @@ def __init__(self, Beam, n_slices, n_sigma = None, cut_left = None,
             raise RuntimeError('The slicing_coord is not recognized')
 
         #: *Number of macroparticles per slice (~profile).*
-        self.n_macroparticles = np.zeros(n_slices, dtype='uint32')
+        self.n_macroparticles = np.zeros(n_slices)
 
         #: *Edges positions of the slicing*
         self.edges = np.zeros(n_slices + 1)
@@ -311,8 +311,7 @@ def beam_spectrum_generation(self, n_sampling_fft, filter_option = None, only_rf
         self.beam_spectrum_freq = rfftfreq(n_sampling_fft, time_step)
 
         if not only_rfft:
-            self.beam_spectrum = rfft(self.n_macroparticles, n_sampling_fft)
-             
+            self.beam_spectrum = rfft(self.n_macroparticles, n_sampling_fft)             
 
 
     def beam_profile_derivative(self, mode = 'gradient', coord = 'theta'):      

impedances/longitudinal_impedance.py

@@ -333,8 +333,7 @@ def track(self, Beam):
         self.induced_voltage_generation(Beam)
 
         induced_voltage_kick = np.interp(Beam.tau, self.slices.bins_centers, self.induced_voltage)
-        Beam.dE += induced_voltage_kick
-        
+        Beam.dE += induced_voltage_kick        
 
 
 class InductiveImpedance(object):