Skip to content
Snippets Groups Projects
Commit a7520315 authored by Vadim Gubaidulin's avatar Vadim Gubaidulin
Browse files

utilities and postprocessing for 

instability growth/damp measurements
parent 12ee4676
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Growth_damp.ipynb 0 → 100644
+ 555
− 0
View file @ a7520315
Source diff could not be displayed: it is too large. Options to address this: view the blob.
growth_damp_utils.py 0 → 100644
+ 158
− 0
View file @ a7520315
import h5py as hp
import numpy as np
from scipy.constants import pi, c, m_p, e
from matplotlib import pyplot as plt
from scipy.signal import hilbert
from aps_figures.aps_one_column import *
from numpy.fft import rfft, rfftfreq
from FITX import fit_risetime
from mbtrack2.instability.ions import *
CIRCUMFERENCE = 354.10
H_RF = 416
Q_X = 18.16
Q_Y = 10.22
OMEGA_REV = 2*pi*c/CIRCUMFERENCE
r_p = e**2/(m_p*c**2)
def read_bunch_data(bunchnumber, filename, folder):
f = hp.File(folder+filename+'.hdf5')
BPM_data = f['raw'][:]
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_C02_before'])
emitYC16after = np.array(f['emitZ_C02_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'])
# print('Horizontal tune before {:.3f} and after {:.3f} the measurement'.format(Qxi, Qxf))
# print('Vertical tune before {:.3f} and after {:.3f} the measurement'.format(Qyi, Qyf))
filling_pattern_before = np.array(f['filling_pattern_before'])
filling_pattern_after = np.array(f['filling_pattern_after'])
f.close()
return BPM_data[bunchnumber,:], t0, t1, emitYC02before, emitYC16before, emitYC02after, emitYC16after, emitXC02before, emitXC02after
def plot_bunch_offset_data(ax, bunchnumber, filename, folder):
mean_y_bunch, t0, t1, emitYC02before, emitYC16before, emitYC02after, emitYC16after, emitXC02before, emitXC02after = read_bunch_data(bunchnumber, filename, folder)
ax.plot(mean_y_bunch-np.mean(mean_y_bunch))
ax.axvline(t0, linestyle='dashed')
ax.axvline(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(ax, bunchnumber, filename, folder):
mean_y_bunch, t0, t1, emitYC02before, emitYC16before, emitYC02after, emitYC16after, emitXC02before, emitXC02after = read_bunch_data(bunchnumber, filename, folder)
mean_y_bunch-=np.mean(mean_y_bunch)
fftfreq = rfftfreq(mean_y_bunch.shape[0])
absfft = np.abs(rfft(mean_y_bunch))
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):
mean_y_bunch, t0, t1, emitYC02before, emitYC16before, emitYC02after, emitYC16after, emitXC02before, emitXC02after = read_bunch_data(bunchnumber, filename, folder)
mean_y_bunch -= np.mean(mean_y_bunch)
signal = np.abs(hilbert(mean_y_bunch))
if show_fitted_signal:
ax.plot(signal)
min_level = np.min((5.*np.mean(signal[:200]), 5.*np.mean(signal[-200:])))
risetime = fit_risetime(signal,
min_level=min_level,
smoothing_window_size=smoothing_window_size,
matplotlib_axis=ax)
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 risetime
def plot_bunch_damping_time_fit(ax, bunchnumber, filename, folder, smoothing_window_size = 2, show_fitted_signal=False):
mean_y_bunch, t0, t1, emitYC02before, emitYC16before, emitYC02after, emitYC16after, emitXC02before, emitXC02after = read_bunch_data(bunchnumber, filename, folder)
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))
if show_fitted_signal:
ax.plot(signal)
min_level = np.min((5*np.mean(signal[:200]), 5*np.mean(signal[-200:])))
damping_time = fit_risetime(signal,
min_level=min_level,
smoothing_window_size=smoothing_window_size,
matplotlib_axis=ax)
ax.set_xlim(0,)
ax.set_xlabel('Reverse time (turns)')
ax.set_ylabel('Bunch c.\,m. offset (arb. units)')
ax.title.set_text('Bunch \# {:} c.\,m. offset'.format(int(bunchnumber)))
return damping_time
def plot_beam_offset_data(ax, filename, folder, smoothing_window_size = 10):
mean_y_bunch, t0, t1, emitYC02before, emitYC16before, emitYC02after, emitYC16after, emitXC02before, emitXC02after = read_bunch_data(0, filename, folder)
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)
mean_y_bunch -= np.mean(mean_y_bunch)
mean_y_beam[bunchnumber::H_RF] = mean_y_bunch
ax.plot(turns, mean_y_beam)
ax.axvline(t0, linestyle='dashed')
ax.axvline(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)
signal = np.abs(hilbert(mean_y_bunch))
signal_reversed = np.flip(signal)
min_level = np.min((5*np.mean(signal[:200]), 5*np.mean(signal[-200:])))
damping_time = fit_risetime(signal_reversed,
min_level=min_level,
smoothing_window_size=smoothing_window_size,
matplotlib_axis=None)
risetime = fit_risetime(signal_reversed,
min_level=min_level,
smoothing_window_size=smoothing_window_size,
matplotlib_axis=None)
return risetime/H_RF, damping_time/H_RF
def plot_beam_spectrum(ax, filename, folder):
mean_y_bunch, t0, t1, emitYC02before, emitYC16before, emitYC02after, emitYC16after, emitXC02before, emitXC02after = read_bunch_data(0, filename, folder)
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)
mean_y_bunch -= np.mean(mean_y_bunch)
mean_y_beam[bunchnumber::H_RF] = mean_y_bunch
fftfreq = rfftfreq(mean_y_bunch.shape[0])
absfft = np.abs(rfft(mean_y_bunch))
ax.plot(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)
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_xlim(0, .5)
ax.set_ylim(0, 1.05)
ax.set_xlabel('Coherent frequency, $\omega/\omega_0$')
ax.set_ylabel('Power spectrum (arb. units)')
# ax.text(mean_y_beam.shape[0]/15/H_RF, np.min(mean_y_beam)/2, '$\\varepsilon_i={:.1e}$,\n $\\varepsilon_f={:.1e}$'.format(emitYC16before, emitYC16after), fontsize=12)
return None
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