Skip to content
Snippets Groups Projects
Commit 38939b87 authored by Alexis GAMELIN's avatar Alexis GAMELIN
Browse files

Update plotting

parent c31f223a
No related branches found
No related tags found
1 merge request!1Update plotting
Show whitespace changes
Inline Side-by-side
Showing
with 224 additions and 155 deletions
growth_damp_utils.py
+ 224
− 155
View file @ 38939b87
import h5py as hp
import numpy as np
import os
from scipy.constants import pi, c, m_p, e
from matplotlib import pyplot as plt
from scipy.signal import hilbert
......@@ -33,7 +34,7 @@ class FBTPostmortemViewer():
# self.h_rf = h_rf
# print(filename, extension)
if extension == ".nxs":
with hp.File(filename, "r") as file:
with hp.File(self.filename, "r") as file:
print(file['post-mortem/rcm/ans/rf/bbfdataviewer.1/'].keys())
data = np.array(file['post-mortem/rcm/ans/rf/bbfdataviewer.1/clean_data'][:])
print('Shape of nexus raw data', data.shape)
......@@ -53,9 +54,9 @@ class FBTPostmortemViewer():
elif extension == ".hdf5":
with hp.File(self.filename) as f:
# print(f['V'].shape)
data_bunch_x = np.array(f['H'], dtype=np.int16)[:, :]
data_bunch_x = f["raw"][:]
# data_bunch_x -= np.mean(data_bunch_x, axis=1, dtype=np.int64, keepdims=True)
data_bunch_y = np.array(f['V'], dtype=np.int16)[:, :]
data_bunch_y = f["raw"][:]
# data_bunch_y -= np.mean(data_bunch_y, axis=1, dtype=np.int64, keepdims=True)
data_bunch_x = data_bunch_x.T
data_bunch_y = data_bunch_y.T
......@@ -91,10 +92,9 @@ class FBTPostmortemViewer():
# self.tunes_x = [nafflib.get_tune(self.data_bunch_x[:, i]) for i in range(self.h_rf)]
# self.tunes_y = [nafflib.get_tune(self.data_bunch_y[:, i]) for i in range(self.h_rf)]
self.fft_x = np.asarray([np.fft.rfft(self.data_bunch_x[:, i]) for i in range(self.h_rf)]).T
self.fft_y = np.asarray([np.fft.rfft(self.data_bunch_y[:, i]) for i in range(self.h_rf)]).T
self.fft_x = np.fft.rfft(self.data_bunch_x, axis=0).T
self.fft_y = np.fft.rfft(self.data_bunch_x, axis=1).T
self.fft_y = np.fft.rfft(self.data_bunch_y, axis=0).T
self.fftfreq = np.fft.rfftfreq(self.data_bunch_x.shape[0])
......@@ -124,36 +124,56 @@ class FBTPostmortemViewer():
def plot_partial_spectrum(self, ax=None, plane='y'):
plane = list(plane)
if ax == None:
_, ax = plt.subplots(1, 1)
_, ax = plt.subplots(4, 1, figsize=(10,8))
if 'x' in plane:
ax.plot(self.spectrum_x1[0], self.spectrum_x1[1], label='horizontal, ADC1')
ax.plot(self.spectrum_x2[0], self.spectrum_x2[1], label='horizontal, ADC2')
ax.plot(self.spectrum_x3[0], self.spectrum_x3[1], label='horizontal, ADC3')
ax.plot(self.spectrum_x4[0], self.spectrum_x4[1], label='horizontal, ADC4')
ax[0].plot(self.spectrum_x1[0], self.spectrum_x1[1], label='horizontal, ADC1')
ax[0].legend()
ax[1].plot(self.spectrum_x2[0], self.spectrum_x2[1], label='horizontal, ADC2')
ax[1].legend()
ax[2].plot(self.spectrum_x3[0], self.spectrum_x3[1], label='horizontal, ADC3')
ax[2].legend()
ax[3].plot(self.spectrum_x4[0], self.spectrum_x4[1], label='horizontal, ADC4')
ax[3].legend()
if 'y' in plane:
ax.plot(self.spectrum_y1[0], self.spectrum_y1[1], label='vertical, ADC1')
ax.plot(self.spectrum_y2[0], self.spectrum_y2[1], label='vertical, ADC2')
ax.plot(self.spectrum_y3[0], self.spectrum_y3[1], label='vertical, ADC3')
ax.plot(self.spectrum_y4[0], self.spectrum_y4[1], label='vertical, ADC4')
ax.legend()
ax[0].plot(self.spectrum_y1[0], self.spectrum_y1[1], label='vertical, ADC1')
ax[0].legend()
ax[1].plot(self.spectrum_y2[0], self.spectrum_y2[1], label='vertical, ADC2')
ax[1].legend()
ax[2].plot(self.spectrum_y3[0], self.spectrum_y3[1], label='vertical, ADC3')
ax[2].legend()
ax[3].plot(self.spectrum_y4[0], self.spectrum_y4[1], label='vertical, ADC4')
ax[3].legend()
def plot_offset(self, ax=None, plane='y'):
plane = list(plane)
if ax == None:
_, ax = plt.subplots(1, 1)
if 'y' in plane:
ax.plot(self.data_bunch_y.flatten())
ax.plot(self.data_bunch_y.mean(axis=1))
def plot_offset_hist(self, ax=None):
if ax == None:
_, ax = plt.subplots(1, 1)
ax.hist(self.data_bunch_y.flatten(), bins="auto")
ax.set_xlabel('Bunch c.m. offset (arb. units)')
ax.set_ylabel('Counts')
ax.set_title('Bunch c.m. offset histogram')
def plot_bunch_offset(self, bunch_numbers, axs=None, plane='y'):
plane = list(plane)
if axs == None:
fig, axs = plt.subplots(len(bunch_numbers), 1)
for i in range(len(bunch_numbers)):
if 'y' in plane:
if plane== "y":
axs[i].plot(self.data_bunch_y[:, bunch_numbers[i]], label=f'vertical, bunch {bunch_numbers[i]}')
if 'x' in plane:
axs[i].legend()
elif plane== "x":
axs[i].plot(self.data_bunch_x[:, bunch_numbers[i]], label=f'horizontal, bunch {bunch_numbers[i]}')
fig.legend()
axs[i].legend()
#This function does not do anything
def plot_fitexp(self, ax=None, plane='y'):
if ax == None:
_, ax = plt.subplots(1, 1)
......@@ -165,8 +185,8 @@ class FBTPostmortemViewer():
_, ax = plt.subplots(1, 1)
ax.plot(np.abs(self.filling_pattern) / self.filling_pattern_calibration)
ax.set_xlabel('Bunch number')
ax.set_ylabel('Bunch current (?mA?)')
ax.set_ylabel('Bunch current (mA)')
# this gives errors
def plot_tunes(self, ax=None):
if ax == None:
_, ax = plt.subplots(1, 1)
......@@ -191,184 +211,233 @@ class FBTPostmortemViewer():
pass
class GrowthDampDataAnalyzer():
def __init__(self, folder):
self.folder = folder
self.filelist = []
def read_bunch_data(self, filename, bunchnumber, N = 0):
with hp.File(self.folder + filename + '.hdf5') as f:
BPM_data = f[f'raw_{N:}'][:]
emitYC02before = np.array(f['emitZ_C02_before'])
emitYC02after = np.array(f['emitZ_C02_after'])
emitXC02before = np.array(f['emitX_C02_before'])
emitXC02after = np.array(f['emitX_C02_after'])
emitYC16before = np.array(f['emitZ_C16_before'])
emitYC16after = np.array(f['emitZ_C16_after'])
t0 = 1e-9*c/CIRCUMFERENCE*np.array(f['window_delay'])
t1 = 1e-6*c/CIRCUMFERENCE*np.array(f['window_width'])
Qxi = np.array(f['tuneX_before'])
Qxf = np.array(f['tuneX_after'])
Qyi = np.array(f['tuneZ_before'])
Qyf = np.array(f['tuneZ_after'])
filling_pattern_before = np.array(f['filling_pattern_before'])
filling_pattern_after = np.array(f['filling_pattern_after'])
return BPM_data[bunchnumber,:], t0, t1, emitYC02before, emitYC16before, emitYC02after, emitYC16after, emitXC02before, emitXC02after
def plot_bunch_offset_data(self, ax, bunchnumber, filename, N=0):
mean_y_bunch, t0, t1, emitYC02before, emitYC16before, emitYC02after, emitYC16after, emitXC02before, emitXC02after = read_bunch_data(bunchnumber, filename, N)
ax.plot(mean_y_bunch-np.mean(mean_y_bunch))
ax.axvline(t0, linestyle='dashed')
ax.axvline(t1, linestyle='dashed')
"""
A class to analyze growth and damping data .It processes data from an HDF5 file, containing raw signal data, emittance , tune etc parameters for a beam in the X and Z (vertical) planes.
"""
def __init__(self, filename):
"""
Parameters:
- filename (str): The path to the HDF5 file containing the beam data.
"""
self.filename= filename
with hp.File(filename, "r") as f:
self.mean_y_data = f[f'raw'][:] # Extract raw mean data
self.y_means= np.mean(self.mean_y_data, axis= 1) # Compute the mean of the data along axis=1
self.n_turns= self.mean_y_data.shape[1] # Determine the number of turns (columns) and the number of bunches (rows)
self.n_bunches= self.mean_y_data.shape[0]
self.mean_y_data-= self.y_means[:, np.newaxis] # Center the data by subtracting the mean offset for each bunch
self.emitYC02before = np.array(f['emitZ_C02_before'])
self.emitYC02after = np.array(f['emitZ_C02_after'])
self.emitXC02before = np.array(f['emitX_C02_before'])
self.emitXC02after = np.array(f['emitX_C02_after'])
self.emitYC16before = np.array(f['emitZ_C16_before'])
self.emitYC16after = np.array(f['emitZ_C16_after'])
self.t0 = 1e-9*c/CIRCUMFERENCE*np.array(f['window_delay'])
self.t1 = 1e-6*c/CIRCUMFERENCE*np.array(f['window_width'])
self.Qxi = np.array(f['tuneX_before'])
self.Qxf = np.array(f['tuneX_after'])
self.Qyi = np.array(f['tuneZ_before'])
self.Qyf = np.array(f['tuneZ_after'])
self.filling_pattern_before = np.array(f['filling_pattern_before'])
self.filling_pattern_after = np.array(f['filling_pattern_after'])
def plot_bunch_offset_data(self, bunchnumber, ax=None):
if ax == None:
_, ax = plt.subplots()
ax.plot(self.mean_y_data[bunchnumber])
ax.axvline(self.t0, linestyle='dashed')
ax.axvline(self.t1, linestyle='dashed')
ax.set_xlim(0,)
ax.set_xlabel('Time (turns)')
ax.set_ylabel('Bunch c.\,m. offset (arb. units)')
ax.title.set_text('Bunch \# {:} c.\,m. offset'.format(int(bunchnumber)))
return None
def plot_bunch_spectrum_data(self, ax, bunchnumber, filename, N=0):
mean_y_bunch, t0, t1, emitYC02before, emitYC16before, emitYC02after, emitYC16after, emitXC02before, emitXC02after = read_bunch_data(bunchnumber, filename, N)
mean_y_bunch-=np.mean(mean_y_bunch)
fftfreq = rfftfreq(mean_y_bunch.shape[0])
absfft = np.abs(rfft(mean_y_bunch))
ax.set_ylabel('Bunch c.,m. offset (arb. units)')
ax.title.set_text('Bunch # {:} c.,m. offset'.format(int(bunchnumber)))
def plot_bunch_spectrum_data(self, bunchnumber, ax= None):
if ax == None:
_, ax = plt.subplots()
mean_y= self.mean_y_data[bunchnumber]
fftfreq = rfftfreq(len(mean_y))
absfft = np.abs(rfft(mean_y))
ax.plot(fftfreq, absfft/np.max(absfft))
ax.set_xlim(0, .5)
ax.set_ylim(0, 1.05)
ax.set_xlabel('Coherent tune, $Q_\mathrm{coh}$')
ax.set_ylabel('Spectrum power (arb. units)')
ax.title.set_text('Bunch \# {:} spectrum'.format(int(bunchnumber)))
return None
def plot_bunch_risetime_fit(ax, bunchnumber, filename, folder, smoothing_window_size = 10, show_fitted_signal=False, N=0):
mean_y_bunch, t0, t1, emitYC02before, emitYC16before, emitYC02after, emitYC16after, emitXC02before, emitXC02after = read_bunch_data(bunchnumber, filename, folder, N)
mean_y_bunch -= np.mean(mean_y_bunch)
signal = np.abs(hilbert(mean_y_bunch))
ax.title.set_text('Bunch # {:} spectrum'.format(int(bunchnumber)))
def plot_bunch_risetime_fit(self,
bunchnumber,
avg_window=30,
min_level_select=50,
min_level_factor=1.8,
smoothing_window_size= 10,
min_points= 3,
min_n_risetimes= 0.5,
show_fitted_signal=False,
ax=None):
if ax == None:
_, ax = plt.subplots()
mean_y= self.mean_y_data[bunchnumber]
signal = np.abs(hilbert(mean_y))
if show_fitted_signal:
ax.plot(signal)
min_level = np.min((10.*np.mean(signal[:200]), 10.*np.mean(signal[-200:])))
ax.plot(signal, alpha= 0.5)
signal= np.convolve(signal, np.ones(avg_window)/avg_window, mode='valid')
min_level = np.min((np.mean(signal[:min_level_select]), np.mean(signal[-min_level_select:])))
risetime = fit_risetime(signal,
min_level=min_level,
min_level=min_level*min_level_factor,
smoothing_window_size=smoothing_window_size,
until= None,
start_from_0= False,
min_points=min_points,
min_n_risetimes=min_n_risetimes,
matplotlib_axis=ax)
ax.set_xlim(0,)
ax.set_xlim(0,len(signal))
ax.set_xlabel('Time (turns)')
ax.set_ylabel('Bunch c.\,m. offset (arb. units)')
ax.title.set_text('Bunch \# {:} c.\,m. offset'.format(int(bunchnumber)))
ax.title.set_text(r'Bunch $\#{}$: Rise time = {:.2f} turns'.format(bunchnumber, risetime))
return risetime
def plot_bunch_damping_time_fit(ax, bunchnumber, filename, folder, smoothing_window_size = 20, show_fitted_signal=False, N=0):
mean_y_bunch, t0, t1, emitYC02before, emitYC16before, emitYC02after, emitYC16after, emitXC02before, emitXC02after = read_bunch_data(bunchnumber, filename, folder, N)
mean_y_bunch -= np.mean(mean_y_bunch)
mean_y_bunch = np.flip(mean_y_bunch)
ax.plot(mean_y_bunch)
signal = np.abs(hilbert(mean_y_bunch))
def plot_bunch_damping_time_fit(self,
bunchnumber,
avg_window=30,
min_level_select=50,
min_level_factor=1.8,
smoothing_window_size= 10,
min_points= 3,
min_n_risetimes= 0.5,
show_fitted_signal=False,
ax=None):
if ax == None:
_, ax = plt.subplots()
mean_y= self.mean_y_data[bunchnumber]
mean_y= np.flip(mean_y)
signal = np.abs(hilbert(mean_y))
if show_fitted_signal:
ax.plot(signal)
min_level = np.min((15*np.mean(signal[:200]), 15*np.mean(signal[-200:])))
ax.plot(signal, alpha= 0.5)
signal= np.convolve(signal, np.ones(avg_window)/avg_window, mode='valid')
min_level = np.min((np.mean(signal[:min_level_select]), np.mean(signal[-min_level_select:])))
damping_time = fit_risetime(signal,
min_level=min_level,
min_level=min_level*min_level_factor,
smoothing_window_size=smoothing_window_size,
until= None,
start_from_0= False,
min_points=min_points,
min_n_risetimes=min_n_risetimes,
matplotlib_axis=ax)
ax.set_xlim(0,)
ax.set_xlabel('Reverse time (turns)')
ax.set_xlim(0,len(signal))
ax.set_xlabel('Time (turns)')
ax.set_ylabel('Bunch c.\,m. offset (arb. units)')
ax.title.set_text('Bunch \# {:} c.\,m. offset'.format(int(bunchnumber)))
ax.set_title(r'Bunch $\#{}$: Damping time = {:.2f} turns'.format(bunchnumber, damping_time))
return damping_time
def plot_beam_offset_data(ax, filename, folder, smoothing_window_size = 10, N=0):
mean_y_bunch, t0, t1, emitYC02before, emitYC16before, emitYC02after, emitYC16after, emitXC02before, emitXC02after = read_bunch_data(0, filename, folder, N)
size = mean_y_bunch.shape[0]
mean_y_beam = np.empty(shape=(size*H_RF), dtype=np.float64)
turns = np.linspace(0, size, H_RF*size)
for bunchnumber in range(0, H_RF):
mean_y_bunch, t0, t1, emitYC02before, emitYC16before, emitYC02after, emitYC16after, emitXC02before, emitXC02after = read_bunch_data(bunchnumber, filename, folder, N)
mean_y_bunch -= np.mean(mean_y_bunch)
mean_y_beam[bunchnumber::H_RF] = mean_y_bunch
def plot_beam_offset_data(self, ax= None):
if ax == None:
_, ax = plt.subplots()
mean_y_beam= self.mean_y_data.flatten()
turns = np.linspace(0, self.n_turns, self.n_turns*self.n_bunches)
ax.plot(turns, mean_y_beam)
ax.axvline(t0, linestyle='dashed')
ax.axvline(t1, linestyle='dashed')
ax.axvline(self.t0, linestyle='dashed')
ax.axvline(self.t1, linestyle='dashed')
ax.set_xlim(0,)
ax.set_xlabel('Time (turns)')
ax.set_ylabel('Beam c.\,m. offset (arb. units)')
ax.text(mean_y_beam.shape[0]/10/H_RF, np.min(mean_y_beam)/2, '$\\varepsilon_i={:.1e}$,\n $\\varepsilon_f={:.1e}$'.format(emitYC16before, emitYC16after), fontsize=12)
return None
ax.text(len(mean_y_beam)/10/H_RF, np.min(mean_y_beam)/2, '$\\varepsilon_i={:.1e}$,\n $\\varepsilon_f={:.1e}$'.format(self.emitYC16before, self.emitYC16after), fontsize=12)
ax.set_title(r'Beam offset data')
def plot_beam_risetime_fit(self,
avg_window=30,
min_level_select=50,
min_level_factor=1.8,
smoothing_window_size= 2000,
min_points= 3,
min_n_risetimes= 0.5,
show_fitted_signal=False,
ax=None):
if ax==None:
_, ax = plt.subplots()
# Flatten mean_y_data for full beam signal
mean_y_beam= self.mean_y_data.flatten()
def plot_beam_risetime_fit(ax, filename, folder, smoothing_window_size = 2000, N=0):
mean_y_bunch, t0, t1, emitYC02before, emitYC16before, emitYC02after, emitYC16after, emitXC02before, emitXC02after = read_bunch_data(0, filename, folder, N)
size = mean_y_bunch.shape[0]
mean_y_beam = np.empty(shape=(size*H_RF), dtype=np.float64)
turns = np.linspace(0, size, H_RF*size)
for bunchnumber in range(0, H_RF):
mean_y_bunch, t0, t1, emitYC02before, emitYC16before, emitYC02after, emitYC16after, emitXC02before, emitXC02after = read_bunch_data(bunchnumber, filename, folder, N)
mean_y_bunch -= np.mean(mean_y_bunch)
mean_y_beam[bunchnumber::H_RF] = mean_y_bunch
ax.axvline(t0, linestyle='dashed', color='black')
ax.axvline(t1, linestyle='dashdot', color='black')
ax.set_xlim(0,)
ax.set_xlabel('Time (turns)')
ax.set_ylabel('Beam c.\,m. offset (arb. units)')
ax.text(mean_y_beam.shape[0]/10/H_RF, np.min(mean_y_beam)/2, '$\\varepsilon_i={:.1e}$,\n $\\varepsilon_f={:.1e}$'.format(emitYC16before, emitYC16after), fontsize=12)
signal = np.abs(hilbert(mean_y_beam))
min_level = np.min((5*np.mean(signal[:200]), 5*np.mean(signal[-200:])))
if show_fitted_signal:
ax.plot(signal, alpha= 0.5)
signal= np.convolve(signal, np.ones(avg_window)/avg_window, mode='valid')
min_level = np.min((np.mean(signal[:min_level_select]), np.mean(signal[-min_level_select:])))
risetime = fit_risetime(signal,
min_level=min_level,
min_level=min_level*min_level_factor,
smoothing_window_size=smoothing_window_size,
until= None,
start_from_0= False,
min_points=min_points,
min_n_risetimes=min_n_risetimes,
matplotlib_axis=ax)
return risetime/H_RF
def plot_beam_damptime_fit(ax, filename, folder, smoothing_window_size = 2000, N=0):
mean_y_bunch, t0, t1, emitYC02before, emitYC16before, emitYC02after, emitYC16after, emitXC02before, emitXC02after = read_bunch_data(0, filename, folder, N)
size = mean_y_bunch.shape[0]
mean_y_beam = np.empty(shape=(size*H_RF), dtype=np.float64)
turns = np.linspace(0, size, H_RF*size)
for bunchnumber in range(0, H_RF):
mean_y_bunch, t0, t1, emitYC02before, emitYC16before, emitYC02after, emitYC16after, emitXC02before, emitXC02after = read_bunch_data(bunchnumber, filename, folder, N)
mean_y_bunch -= np.mean(mean_y_bunch)
mean_y_beam[bunchnumber::H_RF] = mean_y_bunch
ax.axvline(t0, linestyle='dashed', color='black')
ax.axvline(t1, linestyle='dashdot', color='black')
ax.set_xlim(0,)
ax.set_xlim(0, len(signal))
ax.set_xlabel('Time (turns)')
ax.set_ylabel('Beam c.\,m. offset (arb. units)')
ax.text(mean_y_beam.shape[0]/10/H_RF, np.min(mean_y_beam)/2, '$\\varepsilon_i={:.1e}$,\n $\\varepsilon_f={:.1e}$'.format(emitYC16before, emitYC16after), fontsize=12)
signal = np.flip(np.abs(hilbert(mean_y_beam)))
ax.set_ylabel('Beam c.m. offset (arb. units')
ax.set_title(f'Beam rise time fit: {risetime:.2f} turns')
return risetime
def plot_beam_damptime_fit(self,
avg_window=30,
min_level_select=50,
min_level_factor=1.8,
smoothing_window_size= 2000,
min_points= 3,
min_n_risetimes= 0.5,
show_fitted_signal=False,
ax=None):
if ax==None:
_, ax = plt.subplots()
min_level = np.min((10*np.mean(signal[:200]), 10*np.mean(signal[-200:])))
mean_y_beam= self.mean_y_data.flatten()
mean_y_beam= np.flip(mean_y_beam)
signal = np.abs(hilbert(mean_y_beam))
if show_fitted_signal:
ax.plot(signal, alpha= 0.5)
signal= np.convolve(signal, np.ones(avg_window)/avg_window, mode='valid')
min_level = np.min((np.mean(signal[:min_level_select]), np.mean(signal[-min_level_select:])))
damptime = fit_risetime(signal,
min_level=min_level,
min_level=min_level*min_level_factor,
smoothing_window_size=smoothing_window_size,
until= None,
start_from_0= False,
min_points=min_points,
min_n_risetimes=min_n_risetimes,
matplotlib_axis=ax)
return damptime/H_RF
def plot_beam_spectrum(ax, filename, folder, N=0):
mean_y_bunch, t0, t1, emitYC02before, emitYC16before, emitYC02after, emitYC16after, emitXC02before, emitXC02after = read_bunch_data(0, filename, folder, N)
size = mean_y_bunch.shape[0]
mean_y_beam = np.empty(shape=(size*H_RF), dtype=np.float64)
turns = np.linspace(0, size, H_RF*size)
for bunchnumber in range(0, H_RF):
mean_y_bunch, t0, t1, emitYC02before, emitYC16before, emitYC02after, emitYC16after, emitXC02before, emitXC02after = read_bunch_data(bunchnumber, filename, folder, N)
mean_y_bunch -= np.mean(mean_y_bunch)
mean_y_beam[bunchnumber::H_RF] = mean_y_bunch
fftfreq = H_RF*rfftfreq(mean_y_beam.shape[0])
ax.set_xlim(0, len(signal))
ax.set_xlabel('Time (turns)')
ax.set_ylabel('Beam c.m. offset (arb. units)')
ax.set_title(f'Beam damping time fit: {damptime:.2f} turns')
return damptime
def plot_beam_spectrum(self, ax=None):
if ax is None:
_,ax= plt.subplots()
mean_y_beam= self.mean_y_data.flatten()
turns = np.linspace(0, self.n_turns, self.n_turns*self.n_bunches)
fftfreq = H_RF*rfftfreq(len(mean_y_beam))
absfft = np.abs(rfft(mean_y_beam))
ax.semilogy(fftfreq, absfft/np.max(absfft))
beta_x = CIRCUMFERENCE/(2*pi)/Q_X
beta_y = CIRCUMFERENCE/(2*pi)/Q_Y
sig_x = np.sqrt(1e-9*emitXC02before*beta_x)
sig_y = np.sqrt(1e-12*emitYC02before*beta_y)
sig_x = np.sqrt(1e-9*self.emitXC02before*beta_x)
sig_y = np.sqrt(1e-12*self.emitYC02before*beta_y)
bunch_intensity = 100.6e-3/H_RF*CIRCUMFERENCE/c/e
omega_Hx, omega_Hy = get_omega_i(sig_x, sig_y, ion_mass=2, bunch_intensity=bunch_intensity, bunch_spacing=0.85)
omega_COx, omega_COy = get_omega_i(sig_x, sig_y, ion_mass=28, bunch_intensity=bunch_intensity, bunch_spacing=0.85)
omega_CO2x, omega_CO2y = get_omega_i(sig_x, sig_y, ion_mass=44, bunch_intensity=bunch_intensity, bunch_spacing=0.85)
omega_H = ion_frequency(bunch_intensity, 0.85, sig_x, sig_y, ion="H2", dim="y", express="coupling")
omega_CO = ion_frequency(bunch_intensity, 0.85, sig_x, sig_y, ion="CO", dim="y", express="coupling")
omega_CO2 = ion_frequency(bunch_intensity, 0.85, sig_x, sig_y, ion="CO2", dim="y", express="coupling")
print('Ion frequencies H {:.1f}, {:.1f}, CO {:.1f}, {:.1f}, CO2 {:.1f}, {:.1f}'.format(omega_Hy/OMEGA_REV, omega_H/OMEGA_REV,
omega_COy/OMEGA_REV, omega_CO/OMEGA_REV,
omega_CO2y/OMEGA_REV, omega_CO2/OMEGA_REV))
ax.set_xlim(0, H_RF/2)
ax.set_ylim(1e-3, 1.05)
ax.set_xlabel('Coherent frequency, $\omega/\omega_0$')
ax.set_ylabel('Power spectrum (arb. units)')
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment

The place to be for your source code