From de57412741847f4184408d7de6f071d2e53f4816 Mon Sep 17 00:00:00 2001 From: giadarol Date: Thu, 17 Nov 2022 13:33:25 +0100 Subject: [PATCH 01/20] Renormalize eigenvectors after propagation --- examples/tapering/001b_fccee_taper.py | 14 ++++---- xtrack/twiss.py | 52 +++++++++++++++++++++++++++ 2 files changed, 60 insertions(+), 6 deletions(-) diff --git a/examples/tapering/001b_fccee_taper.py b/examples/tapering/001b_fccee_taper.py index 08c7f2ea6..d08832361 100644 --- a/examples/tapering/001b_fccee_taper.py +++ b/examples/tapering/001b_fccee_taper.py @@ -5,11 +5,13 @@ with open('line_no_radiation.json', 'r') as f: line = xt.Line.from_dict(json.load(f)) + + tracker = line.build_tracker() # Introduce some closed orbit -line['qc1l1.1..1'].knl[0] += 1e-6 -line['qc1l1.1..1'].ksl[0] += 1e-6 +#line['qc1l1.1..1'].knl[0] += 1e-6 +#line['qc1l1.1..1'].ksl[0] += 1e-6 # Initial twiss (no radiation) tracker.configure_radiation(model=None) @@ -25,6 +27,8 @@ # Twiss(es) with radiation tw_real_tracking = tracker.twiss(method='6d', matrix_stability_tol=3., eneloss_and_damping=True) +tw_real_tracking_4d = tracker.twiss(method='4d', matrix_stability_tol=3., + eneloss_and_damping=True) tw_sympl = tracker.twiss(radiation_method='kick_as_co', method='6d') tw_scale_as_co = tracker.twiss( radiation_method='scale_as_co', @@ -46,15 +50,13 @@ plt.subplot(2,1,1) plt.plot(tw_no_rad.s, tw_sympl.betx/tw_no_rad.betx - 1) plt.plot(tw_no_rad.s, tw_scale_as_co.betx/tw_no_rad.betx - 1) -#tw.betx *= (1 + delta_beta_corr) -#plt.plot(tw_no_rad.s, tw.betx/tw_no_rad.betx - 1) +plt.plot(tw_no_rad.s, tw_real_tracking.betx/tw_no_rad.betx - 1) plt.ylabel(r'$\Delta \beta_x / \beta_x$') plt.subplot(2,1,2) plt.plot(tw_no_rad.s, tw_sympl.bety/tw_no_rad.bety - 1) plt.plot(tw_no_rad.s, tw_scale_as_co.bety/tw_no_rad.bety - 1) -#tw.bety *= (1 + delta_beta_corr) -#plt.plot(tw_no_rad.s, tw.bety/tw_no_rad.bety - 1) +plt.plot(tw_no_rad.s, tw_real_tracking.bety/tw_no_rad.bety - 1) plt.ylabel(r'$\Delta \beta_y / \beta_y$') plt.figure(10) diff --git a/xtrack/twiss.py b/xtrack/twiss.py index 1e2f7551f..d49606094 100644 --- a/xtrack/twiss.py +++ b/xtrack/twiss.py @@ -374,6 +374,58 @@ def _propagate_optics(tracker, W_matrix, particle_on_co, Ws[:, 4, :] = (tracker.record_last_track.zeta[:6, i_start:i_stop+1] - zeta_co).T / scale_eigen Ws[:, 5, :] = (tracker.record_last_track.ptau[:6, i_start:i_stop+1] - ptau_co).T / particle_on_co._xobject.beta0[0] / scale_eigen + # Re normalize eigenvalues + print('Re-normalizing eigenvalues...') + v1 = Ws[:, :, 0] + 1j * Ws[:, :, 1] + v2 = Ws[:, :, 2] + 1j * Ws[:, :, 3] + v3 = Ws[:, :, 4] + 1j * Ws[:, :, 5] + + S = lnf.S + S_v1_imag = v1 * 0.0 + S_v2_imag = v2 * 0.0 + S_v3_imag = v3 * 0.0 + for ii in range(6): + for jj in range(6): + if S[ii, jj] !=0: + S_v1_imag[:, ii] += S[ii, jj] * v1.imag[:, jj] + S_v2_imag[:, ii] += S[ii, jj] * v2.imag[:, jj] + S_v3_imag[:, ii] += S[ii, jj] * v3.imag[:, jj] + + nux = x_co * 0.0 + 0j + nuy = x_co * 0.0 + 0j + nuzeta = x_co * 0.0 + 0j + + for ii in range(6): + nux += v1.real[:, ii] * S_v1_imag[:, ii] + nuy += v2.real[:, ii] * S_v2_imag[:, ii] + nuzeta += v3.real[:, ii] * S_v3_imag[:, ii] + + nux = np.sqrt(nux) + nuy = np.sqrt(nuy) + nuzeta = np.sqrt(nuzeta) + + for ii in range(6): + v1[:, ii] /= nux + v2[:, ii] /= nuy + v3[:, ii] /= nuzeta + + print('Done re-normalizing eigenvalues.') + + #for ii in range(n_elem): + # print(f'renormalizing eigenvalues for element {ii}/{n_elem}', end='\r', flush=True) + # nux[ii] = np.sqrt(np.matmul(np.matmul(v1[ii, :].real, lnf.S[:, :]), v1[ii, :].imag)) + # v1[ii, :] = v1[ii, :] / nux[ii] + # nuy[ii] = np.sqrt(np.matmul(np.matmul(v2[ii, :].real, lnf.S[:, :]), v2[ii, :].imag)) + # v2[ii, :] = v2[ii, :] / nuy[ii] + # nuzeta[ii] = np.sqrt(np.matmul(np.matmul(v3[ii, :].real, lnf.S), v3[ii, :].imag)) + # v3[ii, :] = v3[ii, :] / nuzeta[ii] + Ws[:, :, 0] = np.real(v1) + Ws[:, :, 1] = np.imag(v1) + Ws[:, :, 2] = np.real(v2) + Ws[:, :, 3] = np.imag(v2) + Ws[:, :, 4] = np.real(v3) + Ws[:, :, 5] = np.imag(v3) + # Rotate eigenvectors to the Courant-Snyder basis phix = np.arctan2(Ws[:, 0, 1], Ws[:, 0, 0]) phiy = np.arctan2(Ws[:, 2, 3], Ws[:, 2, 2]) From 366080c9133049da588c7966683eb3ca2ab25995 Mon Sep 17 00:00:00 2001 From: giadarol Date: Thu, 17 Nov 2022 14:20:33 +0100 Subject: [PATCH 02/20] Working on examples --- examples/tapering/001b_fccee_taper.py | 12 +-- examples/tapering/001d_compare_twiss_modes.py | 80 +++++++++++++++++++ 2 files changed, 86 insertions(+), 6 deletions(-) create mode 100644 examples/tapering/001d_compare_twiss_modes.py diff --git a/examples/tapering/001b_fccee_taper.py b/examples/tapering/001b_fccee_taper.py index d08832361..785210427 100644 --- a/examples/tapering/001b_fccee_taper.py +++ b/examples/tapering/001b_fccee_taper.py @@ -48,26 +48,26 @@ plt.figure(2) plt.subplot(2,1,1) -plt.plot(tw_no_rad.s, tw_sympl.betx/tw_no_rad.betx - 1) -plt.plot(tw_no_rad.s, tw_scale_as_co.betx/tw_no_rad.betx - 1) +#plt.plot(tw_no_rad.s, tw_sympl.betx/tw_no_rad.betx - 1) +#plt.plot(tw_no_rad.s, tw_scale_as_co.betx/tw_no_rad.betx - 1) plt.plot(tw_no_rad.s, tw_real_tracking.betx/tw_no_rad.betx - 1) plt.ylabel(r'$\Delta \beta_x / \beta_x$') plt.subplot(2,1,2) -plt.plot(tw_no_rad.s, tw_sympl.bety/tw_no_rad.bety - 1) -plt.plot(tw_no_rad.s, tw_scale_as_co.bety/tw_no_rad.bety - 1) +#plt.plot(tw_no_rad.s, tw_sympl.bety/tw_no_rad.bety - 1) +#plt.plot(tw_no_rad.s, tw_scale_as_co.bety/tw_no_rad.bety - 1) plt.plot(tw_no_rad.s, tw_real_tracking.bety/tw_no_rad.bety - 1) plt.ylabel(r'$\Delta \beta_y / \beta_y$') plt.figure(10) plt.subplot(2,1,1) -plt.plot(tw_no_rad.s, tw_no_rad.x, 'k') +#plt.plot(tw_no_rad.s, tw_no_rad.x, 'k') #plt.plot(tw_no_rad.s, tw_real_tracking.x, 'b') plt.plot(tw_no_rad.s, tw_sympl.x, 'r') plt.plot(tw_no_rad.s, tw_scale_as_co.x, 'g') plt.subplot(2,1,2) -plt.plot(tw_no_rad.s, tw_no_rad.y, 'k') +#plt.plot(tw_no_rad.s, tw_no_rad.y, 'k') #plt.plot(tw_no_rad.s, tw_real_tracking.y, 'b') plt.plot(tw_no_rad.s, tw_sympl.y, 'r') plt.plot(tw_no_rad.s, tw_scale_as_co.y, 'g') diff --git a/examples/tapering/001d_compare_twiss_modes.py b/examples/tapering/001d_compare_twiss_modes.py new file mode 100644 index 000000000..1a691f8c7 --- /dev/null +++ b/examples/tapering/001d_compare_twiss_modes.py @@ -0,0 +1,80 @@ +import json +import xtrack as xt + +# Build line and tracker +with open('line_no_radiation.json', 'r') as f: + line = xt.Line.from_dict(json.load(f)) + + + +tracker = line.build_tracker() + +# Introduce some closed orbit +#line['qc1l1.1..1'].knl[0] += 1e-6 +#line['qc1l1.1..1'].ksl[0] += 1e-6 + +# Initial twiss (no radiation) +tracker.configure_radiation(model=None) +tw_no_rad = tracker.twiss(method='4d', freeze_longitudinal=True) + +# Enable radiation +tracker.configure_radiation(model='mean') + +# - Set cavity lags to compensate energy loss +# - Taper magnet strengths +tracker.compensate_radiation_energy_loss() + +# Twiss(es) with radiation +tw_real_tracking = tracker.twiss(method='6d', matrix_stability_tol=3., + eneloss_and_damping=True) +tw_real_tracking_4d = tracker.twiss(method='4d', matrix_stability_tol=3., + eneloss_and_damping=True) +tw_sympl = tracker.twiss(radiation_method='kick_as_co', method='6d') +tw_scale_as_co = tracker.twiss( + radiation_method='scale_as_co', + method='6d', + matrix_stability_tol=0.5) + + +import matplotlib.pyplot as plt +plt.close('all') + +print('Non sympltectic tracker:') +print(f'Tune error = error_qx: {abs(tw_real_tracking.qx - tw_no_rad.qx):.3e} error_qy: {abs(tw_real_tracking.qy - tw_no_rad.qy):.3e}') +print('Sympltectic tracker:') +print(f'Tune error = error_qx: {abs(tw_sympl.qx - tw_no_rad.qx):.3e} error_qy: {abs(tw_sympl.qy - tw_no_rad.qy):.3e}') +print ('Preserve angles:') +print(f'Tune error = error_qx: {abs(tw_scale_as_co.qx - tw_no_rad.qx):.3e} error_qy: {abs(tw_scale_as_co.qy - tw_no_rad.qy):.3e}') +plt.figure(2) + +plt.subplot(2,1,1) +plt.plot(tw_no_rad.s, tw_sympl.betx/tw_no_rad.betx - 1) +plt.plot(tw_no_rad.s, tw_scale_as_co.betx/tw_no_rad.betx - 1) +plt.plot(tw_no_rad.s, tw_real_tracking.betx/tw_no_rad.betx - 1) +plt.ylabel(r'$\Delta \beta_x / \beta_x$') + +plt.subplot(2,1,2) +plt.plot(tw_no_rad.s, tw_sympl.bety/tw_no_rad.bety - 1) +plt.plot(tw_no_rad.s, tw_scale_as_co.bety/tw_no_rad.bety - 1) +plt.plot(tw_no_rad.s, tw_real_tracking.bety/tw_no_rad.bety - 1) +plt.ylabel(r'$\Delta \beta_y / \beta_y$') + +plt.figure(10) +plt.subplot(2,1,1) +plt.plot(tw_no_rad.s, tw_no_rad.x, 'k') +plt.plot(tw_no_rad.s, tw_real_tracking.x, 'b') +plt.plot(tw_no_rad.s, tw_sympl.x, 'r') +plt.plot(tw_no_rad.s, tw_scale_as_co.x, 'g') + +plt.subplot(2,1,2) +plt.plot(tw_no_rad.s, tw_no_rad.y, 'k') +plt.plot(tw_no_rad.s, tw_real_tracking.y, 'b') +plt.plot(tw_no_rad.s, tw_sympl.y, 'r') +plt.plot(tw_no_rad.s, tw_scale_as_co.y, 'g') + +plt.figure(3) +plt.subplot() +plt.plot(tw_no_rad.s, tracker.delta_taper) +plt.plot(tw_real_tracking.s, tw_real_tracking.delta) + +plt.show() \ No newline at end of file From 4eb0b580eb3712fa94b9eeaa0b213ade3d6f2341 Mon Sep 17 00:00:00 2001 From: giadarol Date: Thu, 17 Nov 2022 14:27:06 +0100 Subject: [PATCH 03/20] Clean up example --- examples/tapering/001b_fccee_taper.py | 47 +++++++-------------------- xtrack/twiss.py | 8 ----- 2 files changed, 11 insertions(+), 44 deletions(-) diff --git a/examples/tapering/001b_fccee_taper.py b/examples/tapering/001b_fccee_taper.py index 785210427..3ed3ce6e2 100644 --- a/examples/tapering/001b_fccee_taper.py +++ b/examples/tapering/001b_fccee_taper.py @@ -5,14 +5,8 @@ with open('line_no_radiation.json', 'r') as f: line = xt.Line.from_dict(json.load(f)) - - tracker = line.build_tracker() -# Introduce some closed orbit -#line['qc1l1.1..1'].knl[0] += 1e-6 -#line['qc1l1.1..1'].ksl[0] += 1e-6 - # Initial twiss (no radiation) tracker.configure_radiation(model=None) tw_no_rad = tracker.twiss(method='4d', freeze_longitudinal=True) @@ -25,56 +19,37 @@ tracker.compensate_radiation_energy_loss() # Twiss(es) with radiation -tw_real_tracking = tracker.twiss(method='6d', matrix_stability_tol=3., +tw = tracker.twiss(method='6d', matrix_stability_tol=3., eneloss_and_damping=True) -tw_real_tracking_4d = tracker.twiss(method='4d', matrix_stability_tol=3., - eneloss_and_damping=True) -tw_sympl = tracker.twiss(radiation_method='kick_as_co', method='6d') -tw_scale_as_co = tracker.twiss( - radiation_method='scale_as_co', - method='6d', - matrix_stability_tol=0.5) - import matplotlib.pyplot as plt plt.close('all') -print('Non sympltectic tracker:') -print(f'Tune error = error_qx: {abs(tw_real_tracking.qx - tw_no_rad.qx):.3e} error_qy: {abs(tw_real_tracking.qy - tw_no_rad.qy):.3e}') -print('Sympltectic tracker:') -print(f'Tune error = error_qx: {abs(tw_sympl.qx - tw_no_rad.qx):.3e} error_qy: {abs(tw_sympl.qy - tw_no_rad.qy):.3e}') -print ('Preserve angles:') -print(f'Tune error = error_qx: {abs(tw_scale_as_co.qx - tw_no_rad.qx):.3e} error_qy: {abs(tw_scale_as_co.qy - tw_no_rad.qy):.3e}') +print(f'Tune error - error_qx: ' + f'{abs(tw.qx - tw_no_rad.qx):.3e}' + ' error_qy: ' + f'{abs(tw.qy - tw_no_rad.qy):.3e}') + plt.figure(2) plt.subplot(2,1,1) -#plt.plot(tw_no_rad.s, tw_sympl.betx/tw_no_rad.betx - 1) -#plt.plot(tw_no_rad.s, tw_scale_as_co.betx/tw_no_rad.betx - 1) -plt.plot(tw_no_rad.s, tw_real_tracking.betx/tw_no_rad.betx - 1) +plt.plot(tw_no_rad.s, tw.betx/tw_no_rad.betx - 1) plt.ylabel(r'$\Delta \beta_x / \beta_x$') plt.subplot(2,1,2) -#plt.plot(tw_no_rad.s, tw_sympl.bety/tw_no_rad.bety - 1) -#plt.plot(tw_no_rad.s, tw_scale_as_co.bety/tw_no_rad.bety - 1) -plt.plot(tw_no_rad.s, tw_real_tracking.bety/tw_no_rad.bety - 1) +plt.plot(tw_no_rad.s, tw.bety/tw_no_rad.bety - 1) plt.ylabel(r'$\Delta \beta_y / \beta_y$') plt.figure(10) plt.subplot(2,1,1) -#plt.plot(tw_no_rad.s, tw_no_rad.x, 'k') -#plt.plot(tw_no_rad.s, tw_real_tracking.x, 'b') -plt.plot(tw_no_rad.s, tw_sympl.x, 'r') -plt.plot(tw_no_rad.s, tw_scale_as_co.x, 'g') +plt.plot(tw_no_rad.s, tw.x, 'g') plt.subplot(2,1,2) -#plt.plot(tw_no_rad.s, tw_no_rad.y, 'k') -#plt.plot(tw_no_rad.s, tw_real_tracking.y, 'b') -plt.plot(tw_no_rad.s, tw_sympl.y, 'r') -plt.plot(tw_no_rad.s, tw_scale_as_co.y, 'g') +plt.plot(tw_no_rad.s, tw.y, 'g') plt.figure(3) plt.subplot() plt.plot(tw_no_rad.s, tracker.delta_taper) -plt.plot(tw_real_tracking.s, tw_real_tracking.delta) +plt.plot(tw.s, tw.delta) plt.show() \ No newline at end of file diff --git a/xtrack/twiss.py b/xtrack/twiss.py index d49606094..b28bb2ee4 100644 --- a/xtrack/twiss.py +++ b/xtrack/twiss.py @@ -411,14 +411,6 @@ def _propagate_optics(tracker, W_matrix, particle_on_co, print('Done re-normalizing eigenvalues.') - #for ii in range(n_elem): - # print(f'renormalizing eigenvalues for element {ii}/{n_elem}', end='\r', flush=True) - # nux[ii] = np.sqrt(np.matmul(np.matmul(v1[ii, :].real, lnf.S[:, :]), v1[ii, :].imag)) - # v1[ii, :] = v1[ii, :] / nux[ii] - # nuy[ii] = np.sqrt(np.matmul(np.matmul(v2[ii, :].real, lnf.S[:, :]), v2[ii, :].imag)) - # v2[ii, :] = v2[ii, :] / nuy[ii] - # nuzeta[ii] = np.sqrt(np.matmul(np.matmul(v3[ii, :].real, lnf.S), v3[ii, :].imag)) - # v3[ii, :] = v3[ii, :] / nuzeta[ii] Ws[:, :, 0] = np.real(v1) Ws[:, :, 1] = np.imag(v1) Ws[:, :, 2] = np.real(v2) From 3a03612b686783cf6542544db19bb824467173fc Mon Sep 17 00:00:00 2001 From: giadarol Date: Thu, 17 Nov 2022 15:07:37 +0100 Subject: [PATCH 04/20] Fix in tracker.py --- xtrack/tracker.py | 2 +- xtrack/twiss.py | 5 +---- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/xtrack/tracker.py b/xtrack/tracker.py index 7891d007f..f36108ff2 100644 --- a/xtrack/tracker.py +++ b/xtrack/tracker.py @@ -58,7 +58,7 @@ def __init__( _element_ref_data=None, ): self.config = TrackerConfig() - self.XTRACK_MULTIPOLE_NO_SYNRAD=True + self.config.XTRACK_MULTIPOLE_NO_SYNRAD=True if sequence is not None: raise ValueError( diff --git a/xtrack/twiss.py b/xtrack/twiss.py index b28bb2ee4..9b2f77d30 100644 --- a/xtrack/twiss.py +++ b/xtrack/twiss.py @@ -374,8 +374,7 @@ def _propagate_optics(tracker, W_matrix, particle_on_co, Ws[:, 4, :] = (tracker.record_last_track.zeta[:6, i_start:i_stop+1] - zeta_co).T / scale_eigen Ws[:, 5, :] = (tracker.record_last_track.ptau[:6, i_start:i_stop+1] - ptau_co).T / particle_on_co._xobject.beta0[0] / scale_eigen - # Re normalize eigenvalues - print('Re-normalizing eigenvalues...') + # Re normalize eigenvectors v1 = Ws[:, :, 0] + 1j * Ws[:, :, 1] v2 = Ws[:, :, 2] + 1j * Ws[:, :, 3] v3 = Ws[:, :, 4] + 1j * Ws[:, :, 5] @@ -409,8 +408,6 @@ def _propagate_optics(tracker, W_matrix, particle_on_co, v2[:, ii] /= nuy v3[:, ii] /= nuzeta - print('Done re-normalizing eigenvalues.') - Ws[:, :, 0] = np.real(v1) Ws[:, :, 1] = np.imag(v1) Ws[:, :, 2] = np.real(v2) From 8870fbf1f2bbd385e711dffd04ef00e484224b46 Mon Sep 17 00:00:00 2001 From: giadarol Date: Thu, 17 Nov 2022 17:12:07 +0100 Subject: [PATCH 05/20] All asserts in example are passing --- examples/tapering/000_taper_clic_dr.py | 50 ++++++++++++++------------ 1 file changed, 27 insertions(+), 23 deletions(-) diff --git a/examples/tapering/000_taper_clic_dr.py b/examples/tapering/000_taper_clic_dr.py index a2c16c54f..2acf6f6c6 100644 --- a/examples/tapering/000_taper_clic_dr.py +++ b/examples/tapering/000_taper_clic_dr.py @@ -20,7 +20,7 @@ tracker.compensate_radiation_energy_loss() # Twiss(es) with radiation -tw_real_tracking = tracker.twiss(method='6d', matrix_stability_tol=3., +tw = tracker.twiss(method='6d', matrix_stability_tol=3., eneloss_and_damping=True) tw_sympl = tracker.twiss(radiation_method='kick_as_co', method='6d') tw_scale_as_co = tracker.twiss( @@ -28,11 +28,13 @@ method='6d', matrix_stability_tol=0.5) +p0corr = 1 + tracker.delta_taper + import matplotlib.pyplot as plt plt.close('all') print('Non sympltectic tracker:') -print(f'Tune error = error_qx: {abs(tw_real_tracking.qx - tw_no_rad.qx):.3e} error_qy: {abs(tw_real_tracking.qy - tw_no_rad.qy):.3e}') +print(f'Tune error = error_qx: {abs(tw.qx - tw_no_rad.qx):.3e} error_qy: {abs(tw.qy - tw_no_rad.qy):.3e}') print('Sympltectic tracker:') print(f'Tune error = error_qx: {abs(tw_sympl.qx - tw_no_rad.qx):.3e} error_qy: {abs(tw_sympl.qy - tw_no_rad.qy):.3e}') print ('Preserve angles:') @@ -40,17 +42,15 @@ plt.figure(2) plt.subplot(2,1,1) -plt.plot(tw_no_rad.s, tw_sympl.betx/tw_no_rad.betx - 1) -plt.plot(tw_no_rad.s, tw_scale_as_co.betx/tw_no_rad.betx - 1) -#tw.betx *= (1 + delta_beta_corr) -#plt.plot(tw_no_rad.s, tw.betx/tw_no_rad.betx - 1) +plt.plot(tw.s, tw.betx*p0corr/tw_no_rad.betx-1) +plt.plot(tw_no_rad.s, tw_sympl.betx*p0corr/tw_no_rad.betx - 1) +plt.plot(tw_no_rad.s, tw_scale_as_co.betx*p0corr/tw_no_rad.betx - 1) plt.ylabel(r'$\Delta \beta_x / \beta_x$') plt.subplot(2,1,2) -plt.plot(tw_no_rad.s, tw_sympl.bety/tw_no_rad.bety - 1) -plt.plot(tw_no_rad.s, tw_scale_as_co.bety/tw_no_rad.bety - 1) -#tw.bety *= (1 + delta_beta_corr) -#plt.plot(tw_no_rad.s, tw.bety/tw_no_rad.bety - 1) +plt.plot(tw.s, tw.bety*p0corr/tw_no_rad.bety-1) +plt.plot(tw_no_rad.s, tw_sympl.bety*p0corr/tw_no_rad.bety - 1) +plt.plot(tw_no_rad.s, tw_scale_as_co.bety*p0corr/tw_no_rad.bety - 1) plt.ylabel(r'$\Delta \beta_y / \beta_y$') plt.figure(10) @@ -69,54 +69,58 @@ plt.figure(3) plt.subplot() plt.plot(tw_no_rad.s, tracker.delta_taper) -plt.plot(tw_real_tracking.s, tw_real_tracking.delta) +plt.plot(tw.s, tw.delta) assert np.isclose(tracker.delta_taper[0], 0, rtol=0, atol=1e-10) assert np.isclose(tracker.delta_taper[-1], 0, rtol=0, atol=1e-10) assert np.isclose(np.max(tracker.delta_taper), 0.00568948, rtol=1e-4, atol=0) assert np.isclose(np.min(tracker.delta_taper), -0.00556288, rtol=1e-4, atol=0) -assert np.allclose(tw_real_tracking.delta, tracker.delta_taper, rtol=0, atol=1e-6) +assert np.allclose(tw.delta, tracker.delta_taper, rtol=0, atol=1e-6) assert np.allclose(tw_sympl.delta, tracker.delta_taper, rtol=0, atol=1e-6) assert np.allclose(tw_scale_as_co.delta, tracker.delta_taper, rtol=0, atol=1e-6) -assert np.isclose(tw_real_tracking.qx, tw_no_rad.qx, rtol=0, atol=5e-4) +assert np.isclose(tw.qx, tw_no_rad.qx, rtol=0, atol=5e-4) assert np.isclose(tw_sympl.qx, tw_no_rad.qx, rtol=0, atol=5e-4) assert np.isclose(tw_scale_as_co.qx, tw_no_rad.qx, rtol=0, atol=5e-4) -assert np.isclose(tw_real_tracking.qy, tw_no_rad.qy, rtol=0, atol=5e-4) +assert np.isclose(tw.qy, tw_no_rad.qy, rtol=0, atol=5e-4) assert np.isclose(tw_sympl.qy, tw_no_rad.qy, rtol=0, atol=5e-4) assert np.isclose(tw_scale_as_co.qy, tw_no_rad.qy, rtol=0, atol=5e-4) -assert np.isclose(tw_real_tracking.dqx, tw_no_rad.dqx, rtol=0, atol=0.2) +assert np.isclose(tw.dqx, tw_no_rad.dqx, rtol=0, atol=0.2) assert np.isclose(tw_sympl.dqx, tw_no_rad.dqx, rtol=0, atol=0.2) assert np.isclose(tw_scale_as_co.dqx, tw_no_rad.dqx, rtol=0, atol=0.2) -assert np.isclose(tw_real_tracking.dqy, tw_no_rad.dqy, rtol=0, atol=0.2) +assert np.isclose(tw.dqy, tw_no_rad.dqy, rtol=0, atol=0.2) assert np.isclose(tw_sympl.dqy, tw_no_rad.dqy, rtol=0, atol=0.2) assert np.isclose(tw_scale_as_co.dqy, tw_no_rad.dqy, rtol=0, atol=0.2) -assert np.allclose(tw_real_tracking.x, tw_no_rad.x, rtol=0, atol=1e-7) +assert np.allclose(tw.x, tw_no_rad.x, rtol=0, atol=1e-7) assert np.allclose(tw_sympl.x, tw_no_rad.x, rtol=0, atol=1e-7) assert np.allclose(tw_scale_as_co.x, tw_no_rad.x, rtol=0, atol=1e-7) -assert np.allclose(tw_real_tracking.y, tw_no_rad.y, rtol=0, atol=1e-7) +assert np.allclose(tw.y, tw_no_rad.y, rtol=0, atol=1e-7) assert np.allclose(tw_sympl.y, tw_no_rad.y, rtol=0, atol=1e-7) assert np.allclose(tw_scale_as_co.y, tw_no_rad.y, rtol=0, atol=1e-7) -assert np.allclose(tw_sympl.betx, tw_no_rad.betx, rtol=0.02, atol=0) -assert np.allclose(tw_scale_as_co.betx, tw_no_rad.betx, rtol=0.003, atol=0) +assert np.allclose(tw.betx*p0corr, tw_no_rad.betx, rtol=0.02, atol=0) +assert np.allclose(tw_sympl.betx*p0corr, tw_no_rad.betx, rtol=0.02, atol=0) +assert np.allclose(tw_scale_as_co.betx, tw_no_rad.betx, rtol=0.02, atol=0) -assert np.allclose(tw_sympl.bety, tw_no_rad.bety, rtol=0.04, atol=0) -assert np.allclose(tw_scale_as_co.bety, tw_no_rad.bety, rtol=0.003, atol=0) +assert np.allclose(tw.bety*p0corr, tw_no_rad.bety, rtol=0.02, atol=0) +assert np.allclose(tw_sympl.bety*p0corr, tw_no_rad.bety, rtol=0.02, atol=0) +assert np.allclose(tw_scale_as_co.bety, tw_no_rad.bety, rtol=0.02, atol=0) +assert np.allclose(tw.dx, tw.dx, rtol=0.00, atol=0.1e-3) assert np.allclose(tw_sympl.dx, tw_no_rad.dx, rtol=0.00, atol=0.1e-3) assert np.allclose(tw_scale_as_co.dx, tw_no_rad.dx, rtol=0.00, atol=0.1e-3) +assert np.allclose(tw.dy, tw.dy, rtol=0.00, atol=0.1e-3) assert np.allclose(tw_sympl.dy, tw_no_rad.dy, rtol=0.00, atol=0.1e-3) assert np.allclose(tw_scale_as_co.dy, tw_no_rad.dy, rtol=0.00, atol=0.1e-3) -eneloss_real_tracking = tw_real_tracking.eneloss_turn +eneloss_real_tracking = tw.eneloss_turn assert np.isclose(line['rf'].voltage*np.sin(line['rf'].lag/180*np.pi), eneloss_real_tracking/4, rtol=1e-5) assert np.isclose(line['rf1'].voltage*np.sin(line['rf1'].lag/180*np.pi), eneloss_real_tracking/4, rtol=1e-5) assert np.isclose(line['rf2a'].voltage*np.sin(line['rf2a'].lag/180*np.pi), eneloss_real_tracking/4*0.6, rtol=1e-5) From f45d1defafb79d4a2265ec3657ae202144b6f0af Mon Sep 17 00:00:00 2001 From: giadarol Date: Thu, 17 Nov 2022 17:24:42 +0100 Subject: [PATCH 06/20] Fixed test tapering --- tests/test_tapering.py | 32 +++++++++++++++++++------------- 1 file changed, 19 insertions(+), 13 deletions(-) diff --git a/tests/test_tapering.py b/tests/test_tapering.py index b65b30299..45ac9a223 100644 --- a/tests/test_tapering.py +++ b/tests/test_tapering.py @@ -24,7 +24,7 @@ def test_tapering(): tracker.compensate_radiation_energy_loss() # Twiss(es) with radiation - tw_real_tracking = tracker.twiss(method='6d', matrix_stability_tol=3., + tw = tracker.twiss(method='6d', matrix_stability_tol=3., eneloss_and_damping=True) tw_sympl = tracker.twiss(radiation_method='kick_as_co', method='6d') tw_scale_as_co = tracker.twiss( @@ -32,52 +32,58 @@ def test_tapering(): method='6d', matrix_stability_tol=0.5) + p0corr = 1 + tracker.delta_taper + assert np.isclose(tracker.delta_taper[0], 0, rtol=0, atol=1e-10) assert np.isclose(tracker.delta_taper[-1], 0, rtol=0, atol=1e-10) assert np.isclose(np.max(tracker.delta_taper), 0.00568948, rtol=1e-4, atol=0) assert np.isclose(np.min(tracker.delta_taper), -0.00556288, rtol=1e-4, atol=0) - assert np.allclose(tw_real_tracking.delta, tracker.delta_taper, rtol=0, atol=1e-6) + assert np.allclose(tw.delta, tracker.delta_taper, rtol=0, atol=1e-6) assert np.allclose(tw_sympl.delta, tracker.delta_taper, rtol=0, atol=1e-6) assert np.allclose(tw_scale_as_co.delta, tracker.delta_taper, rtol=0, atol=1e-6) - assert np.isclose(tw_real_tracking.qx, tw_no_rad.qx, rtol=0, atol=5e-4) + assert np.isclose(tw.qx, tw_no_rad.qx, rtol=0, atol=5e-4) assert np.isclose(tw_sympl.qx, tw_no_rad.qx, rtol=0, atol=5e-4) assert np.isclose(tw_scale_as_co.qx, tw_no_rad.qx, rtol=0, atol=5e-4) - assert np.isclose(tw_real_tracking.qy, tw_no_rad.qy, rtol=0, atol=5e-4) + assert np.isclose(tw.qy, tw_no_rad.qy, rtol=0, atol=5e-4) assert np.isclose(tw_sympl.qy, tw_no_rad.qy, rtol=0, atol=5e-4) assert np.isclose(tw_scale_as_co.qy, tw_no_rad.qy, rtol=0, atol=5e-4) - assert np.isclose(tw_real_tracking.dqx, tw_no_rad.dqx, rtol=0, atol=0.2) + assert np.isclose(tw.dqx, tw_no_rad.dqx, rtol=0, atol=0.2) assert np.isclose(tw_sympl.dqx, tw_no_rad.dqx, rtol=0, atol=0.2) assert np.isclose(tw_scale_as_co.dqx, tw_no_rad.dqx, rtol=0, atol=0.2) - assert np.isclose(tw_real_tracking.dqy, tw_no_rad.dqy, rtol=0, atol=0.2) + assert np.isclose(tw.dqy, tw_no_rad.dqy, rtol=0, atol=0.2) assert np.isclose(tw_sympl.dqy, tw_no_rad.dqy, rtol=0, atol=0.2) assert np.isclose(tw_scale_as_co.dqy, tw_no_rad.dqy, rtol=0, atol=0.2) - assert np.allclose(tw_real_tracking.x, tw_no_rad.x, rtol=0, atol=1e-7) + assert np.allclose(tw.x, tw_no_rad.x, rtol=0, atol=1e-7) assert np.allclose(tw_sympl.x, tw_no_rad.x, rtol=0, atol=1e-7) assert np.allclose(tw_scale_as_co.x, tw_no_rad.x, rtol=0, atol=1e-7) - assert np.allclose(tw_real_tracking.y, tw_no_rad.y, rtol=0, atol=1e-7) + assert np.allclose(tw.y, tw_no_rad.y, rtol=0, atol=1e-7) assert np.allclose(tw_sympl.y, tw_no_rad.y, rtol=0, atol=1e-7) assert np.allclose(tw_scale_as_co.y, tw_no_rad.y, rtol=0, atol=1e-7) - assert np.allclose(tw_sympl.betx, tw_no_rad.betx, rtol=0.02, atol=0) - assert np.allclose(tw_scale_as_co.betx, tw_no_rad.betx, rtol=0.003, atol=0) + assert np.allclose(tw.betx*p0corr, tw_no_rad.betx, rtol=0.02, atol=0) + assert np.allclose(tw_sympl.betx*p0corr, tw_no_rad.betx, rtol=0.02, atol=0) + assert np.allclose(tw_scale_as_co.betx, tw_no_rad.betx, rtol=0.02, atol=0) - assert np.allclose(tw_sympl.bety, tw_no_rad.bety, rtol=0.04, atol=0) - assert np.allclose(tw_scale_as_co.bety, tw_no_rad.bety, rtol=0.003, atol=0) + assert np.allclose(tw.bety*p0corr, tw_no_rad.bety, rtol=0.02, atol=0) + assert np.allclose(tw_sympl.bety*p0corr, tw_no_rad.bety, rtol=0.02, atol=0) + assert np.allclose(tw_scale_as_co.bety, tw_no_rad.bety, rtol=0.02, atol=0) + assert np.allclose(tw.dx, tw.dx, rtol=0.00, atol=0.1e-3) assert np.allclose(tw_sympl.dx, tw_no_rad.dx, rtol=0.00, atol=0.1e-3) assert np.allclose(tw_scale_as_co.dx, tw_no_rad.dx, rtol=0.00, atol=0.1e-3) + assert np.allclose(tw.dy, tw.dy, rtol=0.00, atol=0.1e-3) assert np.allclose(tw_sympl.dy, tw_no_rad.dy, rtol=0.00, atol=0.1e-3) assert np.allclose(tw_scale_as_co.dy, tw_no_rad.dy, rtol=0.00, atol=0.1e-3) - eneloss_real_tracking = tw_real_tracking.eneloss_turn + eneloss_real_tracking = tw.eneloss_turn assert np.isclose(line['rf'].voltage*np.sin(line['rf'].lag/180*np.pi), eneloss_real_tracking/4, rtol=1e-5) assert np.isclose(line['rf1'].voltage*np.sin(line['rf1'].lag/180*np.pi), eneloss_real_tracking/4, rtol=1e-5) assert np.isclose(line['rf2a'].voltage*np.sin(line['rf2a'].lag/180*np.pi), eneloss_real_tracking/4*0.6, rtol=1e-5) From deb2da20ddee95afdc4bb02bde86b7c3620274db Mon Sep 17 00:00:00 2001 From: giadarol Date: Thu, 17 Nov 2022 22:31:27 +0100 Subject: [PATCH 07/20] Refactored example --- examples/tapering/000_taper_clic_dr.py | 187 +++++++++++-------------- xtrack/linear_normal_form.py | 6 +- xtrack/twiss.py | 1 - 3 files changed, 86 insertions(+), 108 deletions(-) diff --git a/examples/tapering/000_taper_clic_dr.py b/examples/tapering/000_taper_clic_dr.py index 2acf6f6c6..fe791de06 100644 --- a/examples/tapering/000_taper_clic_dr.py +++ b/examples/tapering/000_taper_clic_dr.py @@ -2,8 +2,21 @@ import numpy as np import xtrack as xt +case_name = 'clic_dr' +filename = '../../test_data/clic_dr/line_for_taper.json' +configs = [ + {'radiation_method': 'full', 'p0_correction': False, 'cavity_preserve_angle': False, 'beta_rtol': 2e-2}, + {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': False, 'beta_rtol': 2e-2}, + {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 2e-5}, + {'radiation_method': 'kick_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-3}, + {'radiation_method': 'scale_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-5}, +] -with open('../../test_data/clic_dr/line_for_taper.json', 'r') as f: +# case_name = 'fcc-ee' +# filename = 'line_no_radiation.json' + + +with open(filename, 'r') as f: line = xt.Line.from_dict(json.load(f)) tracker = line.build_tracker() @@ -19,111 +32,77 @@ # - Taper magnet strengths tracker.compensate_radiation_energy_loss() -# Twiss(es) with radiation -tw = tracker.twiss(method='6d', matrix_stability_tol=3., - eneloss_and_damping=True) -tw_sympl = tracker.twiss(radiation_method='kick_as_co', method='6d') -tw_scale_as_co = tracker.twiss( - radiation_method='scale_as_co', - method='6d', - matrix_stability_tol=0.5) -p0corr = 1 + tracker.delta_taper + import matplotlib.pyplot as plt plt.close('all') +ifig = 0 +for conf in configs: + + ifig += 1 + + # Twiss(es) with radiation + tracker.config.XTRACK_CAVITY_PRESERVE_ANGLE = conf['cavity_preserve_angle'] + tw = tracker.twiss(radiation_method=conf['radiation_method'], + eneloss_and_damping=(conf['radiation_method'] != 'kick_as_co')) + tracker.config.XTRACK_CAVITY_PRESERVE_ANGLE = False + + if conf['p0_correction']: + p0corr = 1 + tracker.delta_taper + else: + p0corr = 1 + + plt.figure(ifig, figsize=(6.4*1.3, 4.8)) + plt.suptitle(f"Radiation method: {conf['radiation_method']}, " + f"p0 correction: {conf['p0_correction']}, " + f"cavity preserve angle: {conf['cavity_preserve_angle']}") + + plt.subplot(2,1,1) + plt.title(f'error on Qx: {abs(tw.qx - tw_no_rad.qx):.2e} ' + r'$(\Delta \beta_x / \beta_x)_{max}$ = ' + f'{np.max(tw.betx*p0corr/tw_no_rad.betx-1):.2e}') + plt.plot(tw.s, tw.betx*p0corr/tw_no_rad.betx-1) + plt.ylabel(r'$\Delta \beta_x / \beta_x$') + + plt.subplot(2,1,2) + plt.title(f'error on Qy: {abs(tw.qy - tw_no_rad.qy):.2e} ' + r'$(\Delta \beta_y / \beta_y)_{max}$ = ' + f'{np.max(tw.bety*p0corr/tw_no_rad.bety-1):.2e}') + plt.plot(tw.s, tw.bety*p0corr/tw_no_rad.bety-1) + plt.ylabel(r'$\Delta \beta_y / \beta_y$') + + plt.subplots_adjust(hspace=0.35, top=.85) + + assert np.isclose(tracker.delta_taper[0], 0, rtol=0, atol=1e-10) + assert np.isclose(tracker.delta_taper[-1], 0, rtol=0, atol=1e-10) + + assert np.allclose(tw.delta, tracker.delta_taper, rtol=0, atol=1e-6) + + assert np.isclose(tw.qx, tw_no_rad.qx, rtol=0, atol=5e-4) + assert np.isclose(tw.qy, tw_no_rad.qy, rtol=0, atol=5e-4) + + assert np.isclose(tw.dqx, tw_no_rad.dqx, rtol=0, atol=0.2) + assert np.isclose(tw.dqy, tw_no_rad.dqy, rtol=0, atol=0.2) + + assert np.allclose(tw.x, tw_no_rad.x, rtol=0, atol=1e-7) + assert np.allclose(tw.y, tw_no_rad.y, rtol=0, atol=1e-7) + + assert np.allclose(tw.betx*p0corr, tw_no_rad.betx, rtol=conf['beta_rtol'], atol=0) + assert np.allclose(tw.bety*p0corr, tw_no_rad.bety, rtol=conf['beta_rtol'], atol=0) + + assert np.allclose(tw.dx, tw.dx, rtol=0.00, atol=0.1e-3) + + assert np.allclose(tw.dy, tw.dy, rtol=0.00, atol=0.1e-3) + + if conf['radiation_method'] != 'kick_as_co': + eneloss = tw.eneloss_turn + assert eneloss/line.particle_ref.energy0 > 0.01 + assert np.isclose(line['rf'].voltage*np.sin(line['rf'].lag/180*np.pi), eneloss/4, rtol=1e-5) + assert np.isclose(line['rf1'].voltage*np.sin(line['rf1'].lag/180*np.pi), eneloss/4, rtol=1e-5) + assert np.isclose(line['rf2a'].voltage*np.sin(line['rf2a'].lag/180*np.pi), eneloss/4*0.6, rtol=1e-5) + assert np.isclose(line['rf2b'].voltage*np.sin(line['rf2b'].lag/180*np.pi), eneloss/4*0.4, rtol=1e-5) + assert np.isclose(line['rf3'].voltage*np.sin(line['rf3'].lag/180*np.pi), eneloss/4, rtol=1e-5) + -print('Non sympltectic tracker:') -print(f'Tune error = error_qx: {abs(tw.qx - tw_no_rad.qx):.3e} error_qy: {abs(tw.qy - tw_no_rad.qy):.3e}') -print('Sympltectic tracker:') -print(f'Tune error = error_qx: {abs(tw_sympl.qx - tw_no_rad.qx):.3e} error_qy: {abs(tw_sympl.qy - tw_no_rad.qy):.3e}') -print ('Preserve angles:') -print(f'Tune error = error_qx: {abs(tw_scale_as_co.qx - tw_no_rad.qx):.3e} error_qy: {abs(tw_scale_as_co.qy - tw_no_rad.qy):.3e}') -plt.figure(2) - -plt.subplot(2,1,1) -plt.plot(tw.s, tw.betx*p0corr/tw_no_rad.betx-1) -plt.plot(tw_no_rad.s, tw_sympl.betx*p0corr/tw_no_rad.betx - 1) -plt.plot(tw_no_rad.s, tw_scale_as_co.betx*p0corr/tw_no_rad.betx - 1) -plt.ylabel(r'$\Delta \beta_x / \beta_x$') - -plt.subplot(2,1,2) -plt.plot(tw.s, tw.bety*p0corr/tw_no_rad.bety-1) -plt.plot(tw_no_rad.s, tw_sympl.bety*p0corr/tw_no_rad.bety - 1) -plt.plot(tw_no_rad.s, tw_scale_as_co.bety*p0corr/tw_no_rad.bety - 1) -plt.ylabel(r'$\Delta \beta_y / \beta_y$') - -plt.figure(10) -plt.subplot(2,1,1) -plt.plot(tw_no_rad.s, tw_no_rad.x, 'k') -#plt.plot(tw_no_rad.s, tw_real_tracking.x, 'b') -plt.plot(tw_no_rad.s, tw_sympl.x, 'r') -plt.plot(tw_no_rad.s, tw_scale_as_co.x, 'g') - -plt.subplot(2,1,2) -plt.plot(tw_no_rad.s, tw_no_rad.y, 'k') -#plt.plot(tw_no_rad.s, tw_real_tracking.y, 'b') -plt.plot(tw_no_rad.s, tw_sympl.y, 'r') -plt.plot(tw_no_rad.s, tw_scale_as_co.y, 'g') - -plt.figure(3) -plt.subplot() -plt.plot(tw_no_rad.s, tracker.delta_taper) -plt.plot(tw.s, tw.delta) - -assert np.isclose(tracker.delta_taper[0], 0, rtol=0, atol=1e-10) -assert np.isclose(tracker.delta_taper[-1], 0, rtol=0, atol=1e-10) -assert np.isclose(np.max(tracker.delta_taper), 0.00568948, rtol=1e-4, atol=0) -assert np.isclose(np.min(tracker.delta_taper), -0.00556288, rtol=1e-4, atol=0) - -assert np.allclose(tw.delta, tracker.delta_taper, rtol=0, atol=1e-6) -assert np.allclose(tw_sympl.delta, tracker.delta_taper, rtol=0, atol=1e-6) -assert np.allclose(tw_scale_as_co.delta, tracker.delta_taper, rtol=0, atol=1e-6) - -assert np.isclose(tw.qx, tw_no_rad.qx, rtol=0, atol=5e-4) -assert np.isclose(tw_sympl.qx, tw_no_rad.qx, rtol=0, atol=5e-4) -assert np.isclose(tw_scale_as_co.qx, tw_no_rad.qx, rtol=0, atol=5e-4) - -assert np.isclose(tw.qy, tw_no_rad.qy, rtol=0, atol=5e-4) -assert np.isclose(tw_sympl.qy, tw_no_rad.qy, rtol=0, atol=5e-4) -assert np.isclose(tw_scale_as_co.qy, tw_no_rad.qy, rtol=0, atol=5e-4) - -assert np.isclose(tw.dqx, tw_no_rad.dqx, rtol=0, atol=0.2) -assert np.isclose(tw_sympl.dqx, tw_no_rad.dqx, rtol=0, atol=0.2) -assert np.isclose(tw_scale_as_co.dqx, tw_no_rad.dqx, rtol=0, atol=0.2) - -assert np.isclose(tw.dqy, tw_no_rad.dqy, rtol=0, atol=0.2) -assert np.isclose(tw_sympl.dqy, tw_no_rad.dqy, rtol=0, atol=0.2) -assert np.isclose(tw_scale_as_co.dqy, tw_no_rad.dqy, rtol=0, atol=0.2) - -assert np.allclose(tw.x, tw_no_rad.x, rtol=0, atol=1e-7) -assert np.allclose(tw_sympl.x, tw_no_rad.x, rtol=0, atol=1e-7) -assert np.allclose(tw_scale_as_co.x, tw_no_rad.x, rtol=0, atol=1e-7) - -assert np.allclose(tw.y, tw_no_rad.y, rtol=0, atol=1e-7) -assert np.allclose(tw_sympl.y, tw_no_rad.y, rtol=0, atol=1e-7) -assert np.allclose(tw_scale_as_co.y, tw_no_rad.y, rtol=0, atol=1e-7) - -assert np.allclose(tw.betx*p0corr, tw_no_rad.betx, rtol=0.02, atol=0) -assert np.allclose(tw_sympl.betx*p0corr, tw_no_rad.betx, rtol=0.02, atol=0) -assert np.allclose(tw_scale_as_co.betx, tw_no_rad.betx, rtol=0.02, atol=0) - -assert np.allclose(tw.bety*p0corr, tw_no_rad.bety, rtol=0.02, atol=0) -assert np.allclose(tw_sympl.bety*p0corr, tw_no_rad.bety, rtol=0.02, atol=0) -assert np.allclose(tw_scale_as_co.bety, tw_no_rad.bety, rtol=0.02, atol=0) - -assert np.allclose(tw.dx, tw.dx, rtol=0.00, atol=0.1e-3) -assert np.allclose(tw_sympl.dx, tw_no_rad.dx, rtol=0.00, atol=0.1e-3) -assert np.allclose(tw_scale_as_co.dx, tw_no_rad.dx, rtol=0.00, atol=0.1e-3) - -assert np.allclose(tw.dy, tw.dy, rtol=0.00, atol=0.1e-3) -assert np.allclose(tw_sympl.dy, tw_no_rad.dy, rtol=0.00, atol=0.1e-3) -assert np.allclose(tw_scale_as_co.dy, tw_no_rad.dy, rtol=0.00, atol=0.1e-3) - -eneloss_real_tracking = tw.eneloss_turn -assert np.isclose(line['rf'].voltage*np.sin(line['rf'].lag/180*np.pi), eneloss_real_tracking/4, rtol=1e-5) -assert np.isclose(line['rf1'].voltage*np.sin(line['rf1'].lag/180*np.pi), eneloss_real_tracking/4, rtol=1e-5) -assert np.isclose(line['rf2a'].voltage*np.sin(line['rf2a'].lag/180*np.pi), eneloss_real_tracking/4*0.6, rtol=1e-5) -assert np.isclose(line['rf2b'].voltage*np.sin(line['rf2b'].lag/180*np.pi), eneloss_real_tracking/4*0.4, rtol=1e-5) -assert np.isclose(line['rf3'].voltage*np.sin(line['rf3'].lag/180*np.pi), eneloss_real_tracking/4, rtol=1e-5) -plt.show() \ No newline at end of file +plt.show() diff --git a/xtrack/linear_normal_form.py b/xtrack/linear_normal_form.py index 4a5167f0a..120f658ee 100644 --- a/xtrack/linear_normal_form.py +++ b/xtrack/linear_normal_form.py @@ -6,7 +6,7 @@ import numpy as np DEFAULT_MATRIX_RESPONSIVENESS_TOL = 1e-15 -DEFAULT_MATRIX_STABILITY_TOL = 1e-3 +DEFAULT_MATRIX_STABILITY_TOL = None def healy_symplectify(M): # https://accelconf.web.cern.ch/e06/PAPERS/WEPCH152.PDF @@ -69,7 +69,7 @@ def compute_linear_normal_form(M, symplectify=False, only_4d_block=False, M[4:, 4:] = np.array([[np.cos(muz_dummy), np.sin(muz_dummy)], [-np.sin(muz_dummy), np.cos(muz_dummy)]]) - if np.abs(np.linalg.det(M)-1) > stability_tol: + if stability_tol is not None and np.abs(np.linalg.det(M)-1) > stability_tol: raise ValueError( f'The determinant of M is out tolerance. det={np.linalg.det(M)}') @@ -85,7 +85,7 @@ def compute_linear_normal_form(M, symplectify=False, only_4d_block=False, M = healy_symplectify(M) w0, v0 = np.linalg.eig(M) - if np.any(np.abs(w0) > 1. + stability_tol): + if stability_tol is not None and np.any(np.abs(w0) > 1. + stability_tol): raise ValueError('One-turn matrix is unstable. ' f'Magnitudes of eigenvalues are:\n{repr(np.abs(w0))}') diff --git a/xtrack/twiss.py b/xtrack/twiss.py index 9b2f77d30..7451272db 100644 --- a/xtrack/twiss.py +++ b/xtrack/twiss.py @@ -94,7 +94,6 @@ def twiss_from_tracker(tracker, particle_ref=None, method='6d', elif radiation_method == 'scale_as_co': assert isinstance(tracker._context, xo.ContextCpu) # needs to be serial tracker.config.XTRACK_SYNRAD_SCALE_SAME_AS_FIRST = True - tracker.config.XTRACK_CAVITY_PRESERVE_ANGLE = True res = twiss_from_tracker(**kwargs) return res From d0b028a7e75092ccbe99a5d654cbac1a043ed6c0 Mon Sep 17 00:00:00 2001 From: giadarol Date: Thu, 17 Nov 2022 22:38:09 +0100 Subject: [PATCH 08/20] Adapting example to work also for fcc-ee --- examples/tapering/000_taper_clic_dr.py | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/examples/tapering/000_taper_clic_dr.py b/examples/tapering/000_taper_clic_dr.py index fe791de06..412598e59 100644 --- a/examples/tapering/000_taper_clic_dr.py +++ b/examples/tapering/000_taper_clic_dr.py @@ -12,9 +12,15 @@ {'radiation_method': 'scale_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-5}, ] -# case_name = 'fcc-ee' -# filename = 'line_no_radiation.json' - +case_name = 'fcc-ee' +filename = 'line_no_radiation.json' +configs = [ + {'radiation_method': 'full', 'p0_correction': False, 'cavity_preserve_angle': False, 'beta_rtol': 1e-2}, + {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': False, 'beta_rtol': 3e-3}, + {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 2e-5}, + {'radiation_method': 'kick_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-3}, + {'radiation_method': 'scale_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-5}, +] with open(filename, 'r') as f: line = xt.Line.from_dict(json.load(f)) @@ -32,9 +38,6 @@ # - Taper magnet strengths tracker.compensate_radiation_energy_loss() - - - import matplotlib.pyplot as plt plt.close('all') ifig = 0 @@ -74,6 +77,8 @@ plt.subplots_adjust(hspace=0.35, top=.85) + plt.savefig(f'./{case_name}_fig{ifig}.png', dpi=200) + assert np.isclose(tracker.delta_taper[0], 0, rtol=0, atol=1e-10) assert np.isclose(tracker.delta_taper[-1], 0, rtol=0, atol=1e-10) @@ -82,8 +87,8 @@ assert np.isclose(tw.qx, tw_no_rad.qx, rtol=0, atol=5e-4) assert np.isclose(tw.qy, tw_no_rad.qy, rtol=0, atol=5e-4) - assert np.isclose(tw.dqx, tw_no_rad.dqx, rtol=0, atol=0.2) - assert np.isclose(tw.dqy, tw_no_rad.dqy, rtol=0, atol=0.2) + assert np.isclose(tw.dqx, tw_no_rad.dqx, rtol=0, atol=1e-2*tw.qx) + assert np.isclose(tw.dqy, tw_no_rad.dqy, rtol=0, atol=1e-2*tw.qy) assert np.allclose(tw.x, tw_no_rad.x, rtol=0, atol=1e-7) assert np.allclose(tw.y, tw_no_rad.y, rtol=0, atol=1e-7) From 8c2e71147a0290ef55cf5addedf729f8cd2a8025 Mon Sep 17 00:00:00 2001 From: giadarol Date: Thu, 17 Nov 2022 23:11:45 +0100 Subject: [PATCH 09/20] Adapted all thresholds --- examples/tapering/000_taper_clic_dr.py | 26 +++++++++++++------------- 1 file changed, 13 insertions(+), 13 deletions(-) diff --git a/examples/tapering/000_taper_clic_dr.py b/examples/tapering/000_taper_clic_dr.py index 412598e59..1ac03624d 100644 --- a/examples/tapering/000_taper_clic_dr.py +++ b/examples/tapering/000_taper_clic_dr.py @@ -5,21 +5,21 @@ case_name = 'clic_dr' filename = '../../test_data/clic_dr/line_for_taper.json' configs = [ - {'radiation_method': 'full', 'p0_correction': False, 'cavity_preserve_angle': False, 'beta_rtol': 2e-2}, - {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': False, 'beta_rtol': 2e-2}, - {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 2e-5}, - {'radiation_method': 'kick_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-3}, - {'radiation_method': 'scale_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-5}, + {'radiation_method': 'full', 'p0_correction': False, 'cavity_preserve_angle': False, 'beta_rtol': 2e-2, 'q_atol': 5e-4}, + {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': False, 'beta_rtol': 2e-2, 'q_atol': 5e-4}, + {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 2e-5, 'q_atol': 5e-4}, + {'radiation_method': 'kick_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-3, 'q_atol': 5e-4}, + {'radiation_method': 'scale_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-5}, 'q_atol': 5e-4, ] case_name = 'fcc-ee' filename = 'line_no_radiation.json' configs = [ - {'radiation_method': 'full', 'p0_correction': False, 'cavity_preserve_angle': False, 'beta_rtol': 1e-2}, - {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': False, 'beta_rtol': 3e-3}, - {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 2e-5}, - {'radiation_method': 'kick_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-3}, - {'radiation_method': 'scale_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-5}, + {'radiation_method': 'full', 'p0_correction': False, 'cavity_preserve_angle': False, 'beta_rtol': 1e-2, 'q_atol': 5e-4}, + {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': False, 'beta_rtol': 5e-3, 'q_atol': 5e-4}, + {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 3e-4, 'q_atol': 5e-4}, + {'radiation_method': 'kick_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 3e-3, 'q_atol': 7e-4}, + {'radiation_method': 'scale_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-4, 'q_atol': 1e-4}, ] with open(filename, 'r') as f: @@ -84,8 +84,8 @@ assert np.allclose(tw.delta, tracker.delta_taper, rtol=0, atol=1e-6) - assert np.isclose(tw.qx, tw_no_rad.qx, rtol=0, atol=5e-4) - assert np.isclose(tw.qy, tw_no_rad.qy, rtol=0, atol=5e-4) + assert np.isclose(tw.qx, tw_no_rad.qx, rtol=0, atol=conf['q_atol']) + assert np.isclose(tw.qy, tw_no_rad.qy, rtol=0, atol=conf['q_atol']) assert np.isclose(tw.dqx, tw_no_rad.dqx, rtol=0, atol=1e-2*tw.qx) assert np.isclose(tw.dqy, tw_no_rad.dqy, rtol=0, atol=1e-2*tw.qy) @@ -100,7 +100,7 @@ assert np.allclose(tw.dy, tw.dy, rtol=0.00, atol=0.1e-3) - if conf['radiation_method'] != 'kick_as_co': + if case_name == 'clic_dr' and conf['radiation_method'] != 'kick_as_co': eneloss = tw.eneloss_turn assert eneloss/line.particle_ref.energy0 > 0.01 assert np.isclose(line['rf'].voltage*np.sin(line['rf'].lag/180*np.pi), eneloss/4, rtol=1e-5) From 8daaa8f7093f1d82c96957c2d1a899817a06b614 Mon Sep 17 00:00:00 2001 From: giadarol Date: Thu, 17 Nov 2022 23:16:05 +0100 Subject: [PATCH 10/20] Works also for clic-dr --- examples/tapering/000_taper_clic_dr.py | 28 +++++++++++++------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/examples/tapering/000_taper_clic_dr.py b/examples/tapering/000_taper_clic_dr.py index 1ac03624d..d64c8ad34 100644 --- a/examples/tapering/000_taper_clic_dr.py +++ b/examples/tapering/000_taper_clic_dr.py @@ -9,18 +9,18 @@ {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': False, 'beta_rtol': 2e-2, 'q_atol': 5e-4}, {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 2e-5, 'q_atol': 5e-4}, {'radiation_method': 'kick_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-3, 'q_atol': 5e-4}, - {'radiation_method': 'scale_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-5}, 'q_atol': 5e-4, + {'radiation_method': 'scale_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-5, 'q_atol': 5e-4}, ] -case_name = 'fcc-ee' -filename = 'line_no_radiation.json' -configs = [ - {'radiation_method': 'full', 'p0_correction': False, 'cavity_preserve_angle': False, 'beta_rtol': 1e-2, 'q_atol': 5e-4}, - {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': False, 'beta_rtol': 5e-3, 'q_atol': 5e-4}, - {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 3e-4, 'q_atol': 5e-4}, - {'radiation_method': 'kick_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 3e-3, 'q_atol': 7e-4}, - {'radiation_method': 'scale_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-4, 'q_atol': 1e-4}, -] +# case_name = 'fcc-ee' +# filename = 'line_no_radiation.json' +# configs = [ +# {'radiation_method': 'full', 'p0_correction': False, 'cavity_preserve_angle': False, 'beta_rtol': 1e-2, 'q_atol': 5e-4}, +# {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': False, 'beta_rtol': 5e-3, 'q_atol': 5e-4}, +# {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 3e-4, 'q_atol': 5e-4}, +# {'radiation_method': 'kick_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 3e-3, 'q_atol': 7e-4}, +# {'radiation_method': 'scale_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-4, 'q_atol': 1e-4}, +# ] with open(filename, 'r') as f: line = xt.Line.from_dict(json.load(f)) @@ -87,8 +87,8 @@ assert np.isclose(tw.qx, tw_no_rad.qx, rtol=0, atol=conf['q_atol']) assert np.isclose(tw.qy, tw_no_rad.qy, rtol=0, atol=conf['q_atol']) - assert np.isclose(tw.dqx, tw_no_rad.dqx, rtol=0, atol=1e-2*tw.qx) - assert np.isclose(tw.dqy, tw_no_rad.dqy, rtol=0, atol=1e-2*tw.qy) + assert np.isclose(tw.dqx, tw_no_rad.dqx, rtol=0, atol=1.5e-2*tw.qx) + assert np.isclose(tw.dqy, tw_no_rad.dqy, rtol=0, atol=1.5e-2*tw.qy) assert np.allclose(tw.x, tw_no_rad.x, rtol=0, atol=1e-7) assert np.allclose(tw.y, tw_no_rad.y, rtol=0, atol=1e-7) @@ -96,9 +96,9 @@ assert np.allclose(tw.betx*p0corr, tw_no_rad.betx, rtol=conf['beta_rtol'], atol=0) assert np.allclose(tw.bety*p0corr, tw_no_rad.bety, rtol=conf['beta_rtol'], atol=0) - assert np.allclose(tw.dx, tw.dx, rtol=0.00, atol=0.1e-3) + assert np.allclose(tw.dx, tw.dx, rtol=0.0, atol=0.1e-3) - assert np.allclose(tw.dy, tw.dy, rtol=0.00, atol=0.1e-3) + assert np.allclose(tw.dy, tw.dy, rtol=0.0, atol=0.1e-3) if case_name == 'clic_dr' and conf['radiation_method'] != 'kick_as_co': eneloss = tw.eneloss_turn From e50100796cb043f1d2c34f4e40fe413aaa55cecd Mon Sep 17 00:00:00 2001 From: giadarol Date: Thu, 17 Nov 2022 23:19:48 +0100 Subject: [PATCH 11/20] Some renaming --- ...=> 001_taper_and_compare_twiss_methods.py} | 0 .../tapering/001a_fccee_madx_no_radiation.py | 65 --------------- examples/tapering/001b_fccee_taper.py | 55 ------------- examples/tapering/001d_compare_twiss_modes.py | 80 ------------------- ...k.py => 101_fccee_madx_radiation_thick.py} | 0 5 files changed, 200 deletions(-) rename examples/tapering/{000_taper_clic_dr.py => 001_taper_and_compare_twiss_methods.py} (100%) delete mode 100644 examples/tapering/001a_fccee_madx_no_radiation.py delete mode 100644 examples/tapering/001b_fccee_taper.py delete mode 100644 examples/tapering/001d_compare_twiss_modes.py rename examples/tapering/{001c_fccee_madx_radiation_thick.py => 101_fccee_madx_radiation_thick.py} (100%) diff --git a/examples/tapering/000_taper_clic_dr.py b/examples/tapering/001_taper_and_compare_twiss_methods.py similarity index 100% rename from examples/tapering/000_taper_clic_dr.py rename to examples/tapering/001_taper_and_compare_twiss_methods.py diff --git a/examples/tapering/001a_fccee_madx_no_radiation.py b/examples/tapering/001a_fccee_madx_no_radiation.py deleted file mode 100644 index 178414433..000000000 --- a/examples/tapering/001a_fccee_madx_no_radiation.py +++ /dev/null @@ -1,65 +0,0 @@ -from cpymad.madx import Madx - -mad = Madx() - -mad.input(''' -set, format="19.15f"; -option,update_from_parent=true; - -radiationandrf = 1; - -call,file="fccee_t.seq"; - -pbeam = 182.5; -exbeam = 1.46e-9; -eybeam = 2.9e-12; -nbun = 48; -npar = 2.3e11; -halfcrossingangle = 0.015; -ebeam = sqrt( pbeam^2 + emass^2 ); -beam, particle=positron, npart=npar, kbunch=nbun, energy=ebeam, radiate=false, bv=+1, ex=exbeam, ey=eybeam; - -use, sequence = fccee_p_ring; - -voltca1save = voltca1; -voltca2save = voltca2; - - -!rf and radiation off (for 4d twiss) -voltca1 = 0; -voltca2 = 0; -''') - -tw_thick = mad.twiss().dframe() - -mad.input(''' -select, flag=makethin, class=rfcavity, slice = 1, thick = false; -select, flag=makethin, class=rbend, slice = 4, thick = false; -select, flag=makethin, class=quadrupole, slice = 4, thick = false; -select, flag=makethin, class=sextupole, slice = 4; -select, flag=makethin, pattern="qc*", slice=20; -select, flag=makethin, pattern="sy*", slice=20; - -makethin, sequence=fccee_p_ring, style=teapot, makedipedge=false; -use, sequence = fccee_p_ring; -''') - -tw_thin = mad.twiss().dframe() - -# Restore rf_voltage -mad.input(''' - voltca1 = voltca1save; - voltca2 = voltca2save; -''' -) - -import xtrack as xt -import xpart as xp -line = xt.Line.from_madx_sequence(mad.sequence.fccee_p_ring) -line.particle_ref = xp.Particles(mass0=xp.ELECTRON_MASS_EV, - gamma0=mad.sequence.fccee_p_ring.beam.gamma) - -import json -import xobjects as xo -with open('line_no_radiation.json', 'w') as f: - json.dump(line.to_dict(), f, cls=xo.JEncoder) \ No newline at end of file diff --git a/examples/tapering/001b_fccee_taper.py b/examples/tapering/001b_fccee_taper.py deleted file mode 100644 index 3ed3ce6e2..000000000 --- a/examples/tapering/001b_fccee_taper.py +++ /dev/null @@ -1,55 +0,0 @@ -import json -import xtrack as xt - -# Build line and tracker -with open('line_no_radiation.json', 'r') as f: - line = xt.Line.from_dict(json.load(f)) - -tracker = line.build_tracker() - -# Initial twiss (no radiation) -tracker.configure_radiation(model=None) -tw_no_rad = tracker.twiss(method='4d', freeze_longitudinal=True) - -# Enable radiation -tracker.configure_radiation(model='mean') - -# - Set cavity lags to compensate energy loss -# - Taper magnet strengths -tracker.compensate_radiation_energy_loss() - -# Twiss(es) with radiation -tw = tracker.twiss(method='6d', matrix_stability_tol=3., - eneloss_and_damping=True) - -import matplotlib.pyplot as plt -plt.close('all') - -print(f'Tune error - error_qx: ' - f'{abs(tw.qx - tw_no_rad.qx):.3e}' - ' error_qy: ' - f'{abs(tw.qy - tw_no_rad.qy):.3e}') - -plt.figure(2) - -plt.subplot(2,1,1) -plt.plot(tw_no_rad.s, tw.betx/tw_no_rad.betx - 1) -plt.ylabel(r'$\Delta \beta_x / \beta_x$') - -plt.subplot(2,1,2) -plt.plot(tw_no_rad.s, tw.bety/tw_no_rad.bety - 1) -plt.ylabel(r'$\Delta \beta_y / \beta_y$') - -plt.figure(10) -plt.subplot(2,1,1) -plt.plot(tw_no_rad.s, tw.x, 'g') - -plt.subplot(2,1,2) -plt.plot(tw_no_rad.s, tw.y, 'g') - -plt.figure(3) -plt.subplot() -plt.plot(tw_no_rad.s, tracker.delta_taper) -plt.plot(tw.s, tw.delta) - -plt.show() \ No newline at end of file diff --git a/examples/tapering/001d_compare_twiss_modes.py b/examples/tapering/001d_compare_twiss_modes.py deleted file mode 100644 index 1a691f8c7..000000000 --- a/examples/tapering/001d_compare_twiss_modes.py +++ /dev/null @@ -1,80 +0,0 @@ -import json -import xtrack as xt - -# Build line and tracker -with open('line_no_radiation.json', 'r') as f: - line = xt.Line.from_dict(json.load(f)) - - - -tracker = line.build_tracker() - -# Introduce some closed orbit -#line['qc1l1.1..1'].knl[0] += 1e-6 -#line['qc1l1.1..1'].ksl[0] += 1e-6 - -# Initial twiss (no radiation) -tracker.configure_radiation(model=None) -tw_no_rad = tracker.twiss(method='4d', freeze_longitudinal=True) - -# Enable radiation -tracker.configure_radiation(model='mean') - -# - Set cavity lags to compensate energy loss -# - Taper magnet strengths -tracker.compensate_radiation_energy_loss() - -# Twiss(es) with radiation -tw_real_tracking = tracker.twiss(method='6d', matrix_stability_tol=3., - eneloss_and_damping=True) -tw_real_tracking_4d = tracker.twiss(method='4d', matrix_stability_tol=3., - eneloss_and_damping=True) -tw_sympl = tracker.twiss(radiation_method='kick_as_co', method='6d') -tw_scale_as_co = tracker.twiss( - radiation_method='scale_as_co', - method='6d', - matrix_stability_tol=0.5) - - -import matplotlib.pyplot as plt -plt.close('all') - -print('Non sympltectic tracker:') -print(f'Tune error = error_qx: {abs(tw_real_tracking.qx - tw_no_rad.qx):.3e} error_qy: {abs(tw_real_tracking.qy - tw_no_rad.qy):.3e}') -print('Sympltectic tracker:') -print(f'Tune error = error_qx: {abs(tw_sympl.qx - tw_no_rad.qx):.3e} error_qy: {abs(tw_sympl.qy - tw_no_rad.qy):.3e}') -print ('Preserve angles:') -print(f'Tune error = error_qx: {abs(tw_scale_as_co.qx - tw_no_rad.qx):.3e} error_qy: {abs(tw_scale_as_co.qy - tw_no_rad.qy):.3e}') -plt.figure(2) - -plt.subplot(2,1,1) -plt.plot(tw_no_rad.s, tw_sympl.betx/tw_no_rad.betx - 1) -plt.plot(tw_no_rad.s, tw_scale_as_co.betx/tw_no_rad.betx - 1) -plt.plot(tw_no_rad.s, tw_real_tracking.betx/tw_no_rad.betx - 1) -plt.ylabel(r'$\Delta \beta_x / \beta_x$') - -plt.subplot(2,1,2) -plt.plot(tw_no_rad.s, tw_sympl.bety/tw_no_rad.bety - 1) -plt.plot(tw_no_rad.s, tw_scale_as_co.bety/tw_no_rad.bety - 1) -plt.plot(tw_no_rad.s, tw_real_tracking.bety/tw_no_rad.bety - 1) -plt.ylabel(r'$\Delta \beta_y / \beta_y$') - -plt.figure(10) -plt.subplot(2,1,1) -plt.plot(tw_no_rad.s, tw_no_rad.x, 'k') -plt.plot(tw_no_rad.s, tw_real_tracking.x, 'b') -plt.plot(tw_no_rad.s, tw_sympl.x, 'r') -plt.plot(tw_no_rad.s, tw_scale_as_co.x, 'g') - -plt.subplot(2,1,2) -plt.plot(tw_no_rad.s, tw_no_rad.y, 'k') -plt.plot(tw_no_rad.s, tw_real_tracking.y, 'b') -plt.plot(tw_no_rad.s, tw_sympl.y, 'r') -plt.plot(tw_no_rad.s, tw_scale_as_co.y, 'g') - -plt.figure(3) -plt.subplot() -plt.plot(tw_no_rad.s, tracker.delta_taper) -plt.plot(tw_real_tracking.s, tw_real_tracking.delta) - -plt.show() \ No newline at end of file diff --git a/examples/tapering/001c_fccee_madx_radiation_thick.py b/examples/tapering/101_fccee_madx_radiation_thick.py similarity index 100% rename from examples/tapering/001c_fccee_madx_radiation_thick.py rename to examples/tapering/101_fccee_madx_radiation_thick.py From 7943b4e96708166de5e4c01ff862b88245001bbd Mon Sep 17 00:00:00 2001 From: giadarol Date: Fri, 18 Nov 2022 08:40:34 +0100 Subject: [PATCH 12/20] Cleaning up examples --- .../001_taper_and_compare_twiss_methods.py | 20 ++++-- .../tapering/100_generate_fcc_line_no_rad.py | 65 +++++++++++++++++++ 2 files changed, 79 insertions(+), 6 deletions(-) create mode 100644 examples/tapering/100_generate_fcc_line_no_rad.py diff --git a/examples/tapering/001_taper_and_compare_twiss_methods.py b/examples/tapering/001_taper_and_compare_twiss_methods.py index d64c8ad34..035213418 100644 --- a/examples/tapering/001_taper_and_compare_twiss_methods.py +++ b/examples/tapering/001_taper_and_compare_twiss_methods.py @@ -61,19 +61,27 @@ f"p0 correction: {conf['p0_correction']}, " f"cavity preserve angle: {conf['cavity_preserve_angle']}") - plt.subplot(2,1,1) + betx_beat = tw.betx*p0corr/tw_no_rad.betx-1 + bety_beat = tw.bety*p0corr/tw_no_rad.bety-1 + max_betx_beat = np.max(np.abs(betx_beat)) + max_bety_beat = np.max(np.abs(bety_beat)) + spx = plt.subplot(2,1,1) plt.title(f'error on Qx: {abs(tw.qx - tw_no_rad.qx):.2e} ' r'$(\Delta \beta_x / \beta_x)_{max}$ = ' - f'{np.max(tw.betx*p0corr/tw_no_rad.betx-1):.2e}') - plt.plot(tw.s, tw.betx*p0corr/tw_no_rad.betx-1) + f'{max_betx_beat:.2e}') + plt.plot(tw.s, betx_beat) plt.ylabel(r'$\Delta \beta_x / \beta_x$') + plt.ylim(np.max([0.01, 1.1 * max_betx_beat])*np.array([-1, 1])) + plt.xlim([0, tw.s[-1]]) - plt.subplot(2,1,2) + plt.subplot(2,1,2, sharex=spx) plt.title(f'error on Qy: {abs(tw.qy - tw_no_rad.qy):.2e} ' r'$(\Delta \beta_y / \beta_y)_{max}$ = ' - f'{np.max(tw.bety*p0corr/tw_no_rad.bety-1):.2e}') - plt.plot(tw.s, tw.bety*p0corr/tw_no_rad.bety-1) + f'{max_bety_beat:.2e}') + plt.plot(tw.s, bety_beat) plt.ylabel(r'$\Delta \beta_y / \beta_y$') + plt.ylim(np.max([0.01, 1.1 * max_bety_beat])*np.array([-1, 1])) + plt.xlabel('s [m]') plt.subplots_adjust(hspace=0.35, top=.85) diff --git a/examples/tapering/100_generate_fcc_line_no_rad.py b/examples/tapering/100_generate_fcc_line_no_rad.py new file mode 100644 index 000000000..178414433 --- /dev/null +++ b/examples/tapering/100_generate_fcc_line_no_rad.py @@ -0,0 +1,65 @@ +from cpymad.madx import Madx + +mad = Madx() + +mad.input(''' +set, format="19.15f"; +option,update_from_parent=true; + +radiationandrf = 1; + +call,file="fccee_t.seq"; + +pbeam = 182.5; +exbeam = 1.46e-9; +eybeam = 2.9e-12; +nbun = 48; +npar = 2.3e11; +halfcrossingangle = 0.015; +ebeam = sqrt( pbeam^2 + emass^2 ); +beam, particle=positron, npart=npar, kbunch=nbun, energy=ebeam, radiate=false, bv=+1, ex=exbeam, ey=eybeam; + +use, sequence = fccee_p_ring; + +voltca1save = voltca1; +voltca2save = voltca2; + + +!rf and radiation off (for 4d twiss) +voltca1 = 0; +voltca2 = 0; +''') + +tw_thick = mad.twiss().dframe() + +mad.input(''' +select, flag=makethin, class=rfcavity, slice = 1, thick = false; +select, flag=makethin, class=rbend, slice = 4, thick = false; +select, flag=makethin, class=quadrupole, slice = 4, thick = false; +select, flag=makethin, class=sextupole, slice = 4; +select, flag=makethin, pattern="qc*", slice=20; +select, flag=makethin, pattern="sy*", slice=20; + +makethin, sequence=fccee_p_ring, style=teapot, makedipedge=false; +use, sequence = fccee_p_ring; +''') + +tw_thin = mad.twiss().dframe() + +# Restore rf_voltage +mad.input(''' + voltca1 = voltca1save; + voltca2 = voltca2save; +''' +) + +import xtrack as xt +import xpart as xp +line = xt.Line.from_madx_sequence(mad.sequence.fccee_p_ring) +line.particle_ref = xp.Particles(mass0=xp.ELECTRON_MASS_EV, + gamma0=mad.sequence.fccee_p_ring.beam.gamma) + +import json +import xobjects as xo +with open('line_no_radiation.json', 'w') as f: + json.dump(line.to_dict(), f, cls=xo.JEncoder) \ No newline at end of file From e1ebe47cfb5f3a39103fabb14da29452add1901a Mon Sep 17 00:00:00 2001 From: giadarol Date: Fri, 18 Nov 2022 09:39:35 +0100 Subject: [PATCH 13/20] Move renormalization to function --- examples/tapering/000_taper.py | 114 +++++++++++++++++++++++++++++++++ xtrack/twiss.py | 84 ++++++++++++------------ 2 files changed, 158 insertions(+), 40 deletions(-) create mode 100644 examples/tapering/000_taper.py diff --git a/examples/tapering/000_taper.py b/examples/tapering/000_taper.py new file mode 100644 index 000000000..75eec098c --- /dev/null +++ b/examples/tapering/000_taper.py @@ -0,0 +1,114 @@ +import json +import numpy as np +import xtrack as xt + +######################################### +# Load line and twiss with no radiation # +######################################### + +filename = '../../test_data/clic_dr/line_for_taper.json' + +with open(filename, 'r') as f: + line = xt.Line.from_dict(json.load(f)) +tracker = line.build_tracker() + +tracker.configure_radiation(model=None) +tw_no_rad = tracker.twiss(method='4d', freeze_longitudinal=True) + + +############################################### +# Enable radiation and compensate energy loss # +############################################### +tracker.configure_radiation(model='mean') + +# - Set cavity lags to compensate energy loss +# - Taper magnet strengths to avoid optis and orbit distortions +tracker.compensate_radiation_energy_loss() + +############################## +# Twiss to check the results # +############################## + +tw = tracker.twiss(method='6d') + +import matplotlib.pyplot as plt +plt.close('all') +ifig = 0 +for conf in configs: + + ifig += 1 + + # Twiss(es) with radiation + tracker.config.XTRACK_CAVITY_PRESERVE_ANGLE = conf['cavity_preserve_angle'] + tw = tracker.twiss(radiation_method=conf['radiation_method'], + eneloss_and_damping=(conf['radiation_method'] != 'kick_as_co')) + tracker.config.XTRACK_CAVITY_PRESERVE_ANGLE = False + + if conf['p0_correction']: + p0corr = 1 + tracker.delta_taper + else: + p0corr = 1 + + plt.figure(ifig, figsize=(6.4*1.3, 4.8)) + plt.suptitle(f"Radiation method: {conf['radiation_method']}, " + f"p0 correction: {conf['p0_correction']}, " + f"cavity preserve angle: {conf['cavity_preserve_angle']}") + + betx_beat = tw.betx*p0corr/tw_no_rad.betx-1 + bety_beat = tw.bety*p0corr/tw_no_rad.bety-1 + max_betx_beat = np.max(np.abs(betx_beat)) + max_bety_beat = np.max(np.abs(bety_beat)) + spx = plt.subplot(2,1,1) + plt.title(f'error on Qx: {abs(tw.qx - tw_no_rad.qx):.2e} ' + r'$(\Delta \beta_x / \beta_x)_{max}$ = ' + f'{max_betx_beat:.2e}') + plt.plot(tw.s, betx_beat) + plt.ylabel(r'$\Delta \beta_x / \beta_x$') + plt.ylim(np.max([0.01, 1.1 * max_betx_beat])*np.array([-1, 1])) + plt.xlim([0, tw.s[-1]]) + + plt.subplot(2,1,2, sharex=spx) + plt.title(f'error on Qy: {abs(tw.qy - tw_no_rad.qy):.2e} ' + r'$(\Delta \beta_y / \beta_y)_{max}$ = ' + f'{max_bety_beat:.2e}') + plt.plot(tw.s, bety_beat) + plt.ylabel(r'$\Delta \beta_y / \beta_y$') + plt.ylim(np.max([0.01, 1.1 * max_bety_beat])*np.array([-1, 1])) + plt.xlabel('s [m]') + + plt.subplots_adjust(hspace=0.35, top=.85) + + plt.savefig(f'./{case_name}_fig{ifig}.png', dpi=200) + + assert np.isclose(tracker.delta_taper[0], 0, rtol=0, atol=1e-10) + assert np.isclose(tracker.delta_taper[-1], 0, rtol=0, atol=1e-10) + + assert np.allclose(tw.delta, tracker.delta_taper, rtol=0, atol=1e-6) + + assert np.isclose(tw.qx, tw_no_rad.qx, rtol=0, atol=conf['q_atol']) + assert np.isclose(tw.qy, tw_no_rad.qy, rtol=0, atol=conf['q_atol']) + + assert np.isclose(tw.dqx, tw_no_rad.dqx, rtol=0, atol=1.5e-2*tw.qx) + assert np.isclose(tw.dqy, tw_no_rad.dqy, rtol=0, atol=1.5e-2*tw.qy) + + assert np.allclose(tw.x, tw_no_rad.x, rtol=0, atol=1e-7) + assert np.allclose(tw.y, tw_no_rad.y, rtol=0, atol=1e-7) + + assert np.allclose(tw.betx*p0corr, tw_no_rad.betx, rtol=conf['beta_rtol'], atol=0) + assert np.allclose(tw.bety*p0corr, tw_no_rad.bety, rtol=conf['beta_rtol'], atol=0) + + assert np.allclose(tw.dx, tw.dx, rtol=0.0, atol=0.1e-3) + + assert np.allclose(tw.dy, tw.dy, rtol=0.0, atol=0.1e-3) + + if case_name == 'clic_dr' and conf['radiation_method'] != 'kick_as_co': + eneloss = tw.eneloss_turn + assert eneloss/line.particle_ref.energy0 > 0.01 + assert np.isclose(line['rf'].voltage*np.sin(line['rf'].lag/180*np.pi), eneloss/4, rtol=1e-5) + assert np.isclose(line['rf1'].voltage*np.sin(line['rf1'].lag/180*np.pi), eneloss/4, rtol=1e-5) + assert np.isclose(line['rf2a'].voltage*np.sin(line['rf2a'].lag/180*np.pi), eneloss/4*0.6, rtol=1e-5) + assert np.isclose(line['rf2b'].voltage*np.sin(line['rf2b'].lag/180*np.pi), eneloss/4*0.4, rtol=1e-5) + assert np.isclose(line['rf3'].voltage*np.sin(line['rf3'].lag/180*np.pi), eneloss/4, rtol=1e-5) + + +plt.show() diff --git a/xtrack/twiss.py b/xtrack/twiss.py index 7451272db..4b329b3db 100644 --- a/xtrack/twiss.py +++ b/xtrack/twiss.py @@ -373,46 +373,8 @@ def _propagate_optics(tracker, W_matrix, particle_on_co, Ws[:, 4, :] = (tracker.record_last_track.zeta[:6, i_start:i_stop+1] - zeta_co).T / scale_eigen Ws[:, 5, :] = (tracker.record_last_track.ptau[:6, i_start:i_stop+1] - ptau_co).T / particle_on_co._xobject.beta0[0] / scale_eigen - # Re normalize eigenvectors - v1 = Ws[:, :, 0] + 1j * Ws[:, :, 1] - v2 = Ws[:, :, 2] + 1j * Ws[:, :, 3] - v3 = Ws[:, :, 4] + 1j * Ws[:, :, 5] - - S = lnf.S - S_v1_imag = v1 * 0.0 - S_v2_imag = v2 * 0.0 - S_v3_imag = v3 * 0.0 - for ii in range(6): - for jj in range(6): - if S[ii, jj] !=0: - S_v1_imag[:, ii] += S[ii, jj] * v1.imag[:, jj] - S_v2_imag[:, ii] += S[ii, jj] * v2.imag[:, jj] - S_v3_imag[:, ii] += S[ii, jj] * v3.imag[:, jj] - - nux = x_co * 0.0 + 0j - nuy = x_co * 0.0 + 0j - nuzeta = x_co * 0.0 + 0j - - for ii in range(6): - nux += v1.real[:, ii] * S_v1_imag[:, ii] - nuy += v2.real[:, ii] * S_v2_imag[:, ii] - nuzeta += v3.real[:, ii] * S_v3_imag[:, ii] - - nux = np.sqrt(nux) - nuy = np.sqrt(nuy) - nuzeta = np.sqrt(nuzeta) - - for ii in range(6): - v1[:, ii] /= nux - v2[:, ii] /= nuy - v3[:, ii] /= nuzeta - - Ws[:, :, 0] = np.real(v1) - Ws[:, :, 1] = np.imag(v1) - Ws[:, :, 2] = np.real(v2) - Ws[:, :, 3] = np.imag(v2) - Ws[:, :, 4] = np.real(v3) - Ws[:, :, 5] = np.imag(v3) + # Re normalize eigenvectors (needed when radiation is present) + _renormalize_eigenvectors(Ws) # Rotate eigenvectors to the Courant-Snyder basis phix = np.arctan2(Ws[:, 0, 1], Ws[:, 0, 0]) @@ -1058,3 +1020,45 @@ def match_tracker(tracker, vary, targets, **kwargs): tracker.vars[vv] = x0[ii] raise err return fsolve_info + +def _renormalize_eigenvectors(Ws): + # Re normalize eigenvectors + v1 = Ws[:, :, 0] + 1j * Ws[:, :, 1] + v2 = Ws[:, :, 2] + 1j * Ws[:, :, 3] + v3 = Ws[:, :, 4] + 1j * Ws[:, :, 5] + + S = lnf.S + S_v1_imag = v1 * 0.0 + S_v2_imag = v2 * 0.0 + S_v3_imag = v3 * 0.0 + for ii in range(6): + for jj in range(6): + if S[ii, jj] !=0: + S_v1_imag[:, ii] += S[ii, jj] * v1.imag[:, jj] + S_v2_imag[:, ii] += S[ii, jj] * v2.imag[:, jj] + S_v3_imag[:, ii] += S[ii, jj] * v3.imag[:, jj] + + nux = np.squeeze(Ws[:, 0, 0]) * (0.0 + 0j) + nuy = nux * 0.0 + nuzeta = nux * 0.0 + + for ii in range(6): + nux += v1.real[:, ii] * S_v1_imag[:, ii] + nuy += v2.real[:, ii] * S_v2_imag[:, ii] + nuzeta += v3.real[:, ii] * S_v3_imag[:, ii] + + nux = np.sqrt(nux) + nuy = np.sqrt(nuy) + nuzeta = np.sqrt(nuzeta) + + for ii in range(6): + v1[:, ii] /= nux + v2[:, ii] /= nuy + v3[:, ii] /= nuzeta + + Ws[:, :, 0] = np.real(v1) + Ws[:, :, 1] = np.imag(v1) + Ws[:, :, 2] = np.real(v2) + Ws[:, :, 3] = np.imag(v2) + Ws[:, :, 4] = np.real(v3) + Ws[:, :, 5] = np.imag(v3) From eca1473ddbeddcb92d7a5847a87184336f429ffe Mon Sep 17 00:00:00 2001 From: giadarol Date: Fri, 18 Nov 2022 13:40:39 +0100 Subject: [PATCH 14/20] Testing forest method --- examples/tapering/000_taper.py | 6 +++++- xtrack/twiss.py | 34 ++++++++++++++++++++++++++++++++++ 2 files changed, 39 insertions(+), 1 deletion(-) diff --git a/examples/tapering/000_taper.py b/examples/tapering/000_taper.py index 75eec098c..ef7517c0a 100644 --- a/examples/tapering/000_taper.py +++ b/examples/tapering/000_taper.py @@ -6,12 +6,16 @@ # Load line and twiss with no radiation # ######################################### -filename = '../../test_data/clic_dr/line_for_taper.json' +#filename = '../../test_data/clic_dr/line_for_taper.json' +filename = 'line_no_radiation.json' + with open(filename, 'r') as f: line = xt.Line.from_dict(json.load(f)) tracker = line.build_tracker() +line['qc1l1.1..1'].ksl[0] += 1e-6 + tracker.configure_radiation(model=None) tw_no_rad = tracker.twiss(method='4d', freeze_longitudinal=True) diff --git a/xtrack/twiss.py b/xtrack/twiss.py index 4b329b3db..03eab9a10 100644 --- a/xtrack/twiss.py +++ b/xtrack/twiss.py @@ -380,6 +380,7 @@ def _propagate_optics(tracker, W_matrix, particle_on_co, phix = np.arctan2(Ws[:, 0, 1], Ws[:, 0, 0]) phiy = np.arctan2(Ws[:, 2, 3], Ws[:, 2, 2]) phizeta = np.arctan2(Ws[:, 4, 5], Ws[:, 4, 4]) + v1 = Ws[:, :, 0] + 1j * Ws[:, :, 1] v2 = Ws[:, :, 2] + 1j * Ws[:, :, 3] v3 = Ws[:, :, 4] + 1j * Ws[:, :, 5] @@ -411,6 +412,32 @@ def _propagate_optics(tracker, W_matrix, particle_on_co, muy = muy - muy[0] + muy0 muzeta = muzeta - muzeta[0] + muzeta0 + #### + # Forrest method + + BB = np.zeros(shape=(3, len(s_co), 6, 6), dtype=np.float64) + + + for ii in range(3): + Iii = np.zeros(shape=(6, 6)) + Iii[2*ii, 2*ii] = 1 + Iii[2*ii+1, 2*ii+1] = 1 + Sii = lnf.S @ Iii + + Ws_inv = np.linalg.inv(Ws) + + #BB = Ws@ Sii @ Ws_inv + BB[ii, :, :, :] = Ws @ Sii @ Ws_inv + + #### + + betx_forest = BB[0, :, 0, 1] + bety_forest = BB[1, :, 2, 3] + alfx_forest = BB[0, :, 0, 0] + alfy_forest = BB[1, :, 2, 2] + gamx_forest = -BB[0, :, 1, 1] + gamy_forest = -BB[1, :, 3, 3] + W_matrix = [Ws[ii, :, :] for ii in range(len(s_co))] twiss_res_element_by_element = { @@ -440,6 +467,13 @@ def _propagate_optics(tracker, W_matrix, particle_on_co, 'W_matrix': W_matrix, #'delta_disp_minus': delta_disp_minus, # for debug #'delta_disp_plus': delta_disp_plus, # for debug + 'betx_forest': betx_forest, + 'bety_forest': bety_forest, + 'alfx_forest': alfx_forest, + 'alfy_forest': alfy_forest, + 'gamx_forest': gamx_forest, + 'gamy_forest': gamy_forest, + 'BB': BB, } return twiss_res_element_by_element From bb4031b6d7861ab0aaaff3c120296b3485ef0e68 Mon Sep 17 00:00:00 2001 From: giadarol Date: Fri, 18 Nov 2022 13:58:37 +0100 Subject: [PATCH 15/20] Forest method works --- xtrack/twiss.py | 21 ++++++++++++++------- 1 file changed, 14 insertions(+), 7 deletions(-) diff --git a/xtrack/twiss.py b/xtrack/twiss.py index 03eab9a10..b1d56fa7c 100644 --- a/xtrack/twiss.py +++ b/xtrack/twiss.py @@ -413,11 +413,10 @@ def _propagate_optics(tracker, W_matrix, particle_on_co, muzeta = muzeta - muzeta[0] + muzeta0 #### - # Forrest method + # Forest method BB = np.zeros(shape=(3, len(s_co), 6, 6), dtype=np.float64) - for ii in range(3): Iii = np.zeros(shape=(6, 6)) Iii[2*ii, 2*ii] = 1 @@ -426,17 +425,25 @@ def _propagate_optics(tracker, W_matrix, particle_on_co, Ws_inv = np.linalg.inv(Ws) - #BB = Ws@ Sii @ Ws_inv BB[ii, :, :, :] = Ws @ Sii @ Ws_inv - #### - betx_forest = BB[0, :, 0, 1] bety_forest = BB[1, :, 2, 3] alfx_forest = BB[0, :, 0, 0] alfy_forest = BB[1, :, 2, 2] - gamx_forest = -BB[0, :, 1, 1] - gamy_forest = -BB[1, :, 3, 3] + gamx_forest = -BB[0, :, 1, 0] + gamy_forest = -BB[1, :, 3, 2] + + sign_x = np.sign(betx_forest) + sign_y = np.sign(bety_forest) + betx_forest *= sign_x + alfx_forest *= sign_x + gamx_forest *= sign_x + bety_forest *= sign_y + alfy_forest *= sign_y + gamy_forest *= sign_y + + #### W_matrix = [Ws[ii, :, :] for ii in range(len(s_co))] From 11fa220f73de2e6363302a56d20e979891836a72 Mon Sep 17 00:00:00 2001 From: giadarol Date: Fri, 18 Nov 2022 16:02:03 +0100 Subject: [PATCH 16/20] Working in API --- xtrack/linear_normal_form.py | 2 +- xtrack/tracker.py | 1 + xtrack/twiss.py | 101 ++++++++++++++++++----------------- 3 files changed, 53 insertions(+), 51 deletions(-) diff --git a/xtrack/linear_normal_form.py b/xtrack/linear_normal_form.py index 120f658ee..f7f7d83f7 100644 --- a/xtrack/linear_normal_form.py +++ b/xtrack/linear_normal_form.py @@ -6,7 +6,7 @@ import numpy as np DEFAULT_MATRIX_RESPONSIVENESS_TOL = 1e-15 -DEFAULT_MATRIX_STABILITY_TOL = None +DEFAULT_MATRIX_STABILITY_TOL = 1e-3 def healy_symplectify(M): # https://accelconf.web.cern.ch/e06/PAPERS/WEPCH152.PDF diff --git a/xtrack/tracker.py b/xtrack/tracker.py index f36108ff2..1121a40e3 100644 --- a/xtrack/tracker.py +++ b/xtrack/tracker.py @@ -444,6 +444,7 @@ def twiss(self, particle_ref=None, delta0=None, method='6d', matrix_stability_tol=None, symplectify=False, reverse=False, + use_full_inverse=None ): self._check_invalidated() diff --git a/xtrack/twiss.py b/xtrack/twiss.py index b1d56fa7c..96b464be7 100644 --- a/xtrack/twiss.py +++ b/xtrack/twiss.py @@ -47,6 +47,7 @@ def twiss_from_tracker(tracker, particle_ref=None, method='6d', matrix_stability_tol=None, symplectify=False, reverse=False, + use_full_inverse=None ): assert method in ['6d', '4d'], 'Method must be `6d` or `4d`' @@ -206,7 +207,6 @@ def twiss_from_tracker(tracker, particle_ref=None, method='6d', W[2, 5] = dy_dpzeta W[3, 5] = dpy_dpzeta - twiss_res_element_by_element = _propagate_optics( tracker=tracker, W_matrix=W, @@ -216,7 +216,8 @@ def twiss_from_tracker(tracker, particle_ref=None, method='6d', nemitt_x=nemitt_x, nemitt_y=nemitt_y, r_sigma=r_sigma, - delta_disp=delta_disp) + delta_disp=delta_disp, + use_full_inverse=use_full_inverse) twiss_res.update(twiss_res_element_by_element) twiss_res._ebe_fields = twiss_res_element_by_element.keys() @@ -291,7 +292,8 @@ def twiss_from_tracker(tracker, particle_ref=None, method='6d', def _propagate_optics(tracker, W_matrix, particle_on_co, mux0, muy0, muzeta0, ele_start, ele_stop, - nemitt_x, nemitt_y, r_sigma, delta_disp): + nemitt_x, nemitt_y, r_sigma, delta_disp, + use_full_inverse): ctx2np = tracker._context.nparray_from_context_array @@ -395,14 +397,19 @@ def _propagate_optics(tracker, W_matrix, particle_on_co, Ws[:, :, 4] = np.real(v3) Ws[:, :, 5] = np.imag(v3) - betx = Ws[:, 0, 0]**2 + Ws[:, 0, 1]**2 - bety = Ws[:, 2, 2]**2 + Ws[:, 2, 3]**2 + # Computation of twiss parameters + + if use_full_inverse: + betx, alfx, gamx, bety, alfy, gamy = _extract_twiss_parameters_with_inverse(Ws) + else: + betx = Ws[:, 0, 0]**2 + Ws[:, 0, 1]**2 + bety = Ws[:, 2, 2]**2 + Ws[:, 2, 3]**2 - gamx = Ws[:, 1, 0]**2 + Ws[:, 1, 1]**2 - gamy = Ws[:, 3, 2]**2 + Ws[:, 3, 3]**2 + gamx = Ws[:, 1, 0]**2 + Ws[:, 1, 1]**2 + gamy = Ws[:, 3, 2]**2 + Ws[:, 3, 3]**2 - alfx = - Ws[:, 0, 0] * Ws[:, 1, 0] - Ws[:, 0, 1] * Ws[:, 1, 1] - alfy = - Ws[:, 2, 2] * Ws[:, 3, 2] - Ws[:, 2, 3] * Ws[:, 3, 3] + alfx = - Ws[:, 0, 0] * Ws[:, 1, 0] - Ws[:, 0, 1] * Ws[:, 1, 1] + alfy = - Ws[:, 2, 2] * Ws[:, 3, 2] - Ws[:, 2, 3] * Ws[:, 3, 3] mux = np.unwrap(phix)/2/np.pi muy = np.unwrap(phiy)/2/np.pi @@ -412,38 +419,8 @@ def _propagate_optics(tracker, W_matrix, particle_on_co, muy = muy - muy[0] + muy0 muzeta = muzeta - muzeta[0] + muzeta0 - #### - # Forest method - - BB = np.zeros(shape=(3, len(s_co), 6, 6), dtype=np.float64) - - for ii in range(3): - Iii = np.zeros(shape=(6, 6)) - Iii[2*ii, 2*ii] = 1 - Iii[2*ii+1, 2*ii+1] = 1 - Sii = lnf.S @ Iii - - Ws_inv = np.linalg.inv(Ws) - - BB[ii, :, :, :] = Ws @ Sii @ Ws_inv - - betx_forest = BB[0, :, 0, 1] - bety_forest = BB[1, :, 2, 3] - alfx_forest = BB[0, :, 0, 0] - alfy_forest = BB[1, :, 2, 2] - gamx_forest = -BB[0, :, 1, 0] - gamy_forest = -BB[1, :, 3, 2] - - sign_x = np.sign(betx_forest) - sign_y = np.sign(bety_forest) - betx_forest *= sign_x - alfx_forest *= sign_x - gamx_forest *= sign_x - bety_forest *= sign_y - alfy_forest *= sign_y - gamy_forest *= sign_y - - #### + mux = np.abs(mux) + muy = np.abs(muy) W_matrix = [Ws[ii, :, :] for ii in range(len(s_co))] @@ -472,15 +449,6 @@ def _propagate_optics(tracker, W_matrix, particle_on_co, 'muy': muy, 'muzeta': muzeta, 'W_matrix': W_matrix, - #'delta_disp_minus': delta_disp_minus, # for debug - #'delta_disp_plus': delta_disp_plus, # for debug - 'betx_forest': betx_forest, - 'bety_forest': bety_forest, - 'alfx_forest': alfx_forest, - 'alfy_forest': alfy_forest, - 'gamx_forest': gamx_forest, - 'gamy_forest': gamy_forest, - 'BB': BB, } return twiss_res_element_by_element @@ -1103,3 +1071,36 @@ def _renormalize_eigenvectors(Ws): Ws[:, :, 3] = np.imag(v2) Ws[:, :, 4] = np.real(v3) Ws[:, :, 5] = np.imag(v3) + + +def _extract_twiss_parameters_with_inverse(Ws): + + BB = np.zeros(shape=(3, Ws.shape[0], 6, 6), dtype=np.float64) + + for ii in range(3): + Iii = np.zeros(shape=(6, 6)) + Iii[2*ii, 2*ii] = 1 + Iii[2*ii+1, 2*ii+1] = 1 + Sii = lnf.S @ Iii + + Ws_inv = np.linalg.inv(Ws) + + BB[ii, :, :, :] = Ws @ Sii @ Ws_inv + + betx = BB[0, :, 0, 1] + bety = BB[1, :, 2, 3] + alfx = BB[0, :, 0, 0] + alfy = BB[1, :, 2, 2] + gamx = -BB[0, :, 1, 0] + gamy = -BB[1, :, 3, 2] + + sign_x = np.sign(betx) + sign_y = np.sign(bety) + betx *= sign_x + alfx *= sign_x + gamx *= sign_x + bety *= sign_y + alfy *= sign_y + gamy *= sign_y + + return betx, alfx, gamx, bety, alfy, gamy \ No newline at end of file From 9e29f0ea49b7020b7a122d1a7ef76b6b16a38916 Mon Sep 17 00:00:00 2001 From: giadarol Date: Fri, 18 Nov 2022 16:19:46 +0100 Subject: [PATCH 17/20] Automatically handle matrix stability tol and use_full_inverse if radiation is present --- xtrack/tracker.py | 5 +++++ xtrack/twiss.py | 5 +++++ 2 files changed, 10 insertions(+) diff --git a/xtrack/tracker.py b/xtrack/tracker.py index 1121a40e3..74ca05c74 100644 --- a/xtrack/tracker.py +++ b/xtrack/tracker.py @@ -272,6 +272,7 @@ def _init_track_with_collective( self.particles_monitor_class = supertracker.particles_monitor_class self._element_part = _element_part self._element_index_in_part = _element_index_in_part + self._radiation_model = None def _init_track_no_collective( self, @@ -361,6 +362,7 @@ def _init_track_no_collective( self.track_kernel = track_kernel self.track = self._track_no_collective + self._radiation_model = None if compile: _ = self._current_track_kernel # This triggers compilation @@ -485,10 +487,13 @@ def configure_radiation(self, model=None, mode=None): assert model in [None, 'mean', 'quantum'] if model == 'mean': radiation_flag = 1 + self._radiation_model = 'mean' elif model == 'quantum': radiation_flag = 2 + self._radiation_model = 'quantum' else: radiation_flag = 0 + self._radiation_model = None for kk, ee in self.line.element_dict.items(): if hasattr(ee, 'radiation_flag'): diff --git a/xtrack/twiss.py b/xtrack/twiss.py index 96b464be7..54a39e7f1 100644 --- a/xtrack/twiss.py +++ b/xtrack/twiss.py @@ -57,6 +57,11 @@ def twiss_from_tracker(tracker, particle_ref=None, method='6d', if matrix_stability_tol is None: matrix_stability_tol = tracker.matrix_stability_tol + if tracker._radiation_model is not None: + matrix_stability_tol = None + if use_full_inverse is None: + use_full_inverse = True + if particle_ref is None: if particle_co_guess is None and hasattr(tracker, 'particle_ref'): particle_ref = tracker.particle_ref From 1926fc855a6f09fedcce0d92e85a7b9fd0c9839b Mon Sep 17 00:00:00 2001 From: giadarol Date: Fri, 18 Nov 2022 17:03:53 +0100 Subject: [PATCH 18/20] Example for tapering --- examples/tapering/000_taper.py | 101 ++++++++------------------------- 1 file changed, 23 insertions(+), 78 deletions(-) diff --git a/examples/tapering/000_taper.py b/examples/tapering/000_taper.py index ef7517c0a..92a23f628 100644 --- a/examples/tapering/000_taper.py +++ b/examples/tapering/000_taper.py @@ -6,23 +6,19 @@ # Load line and twiss with no radiation # ######################################### -#filename = '../../test_data/clic_dr/line_for_taper.json' -filename = 'line_no_radiation.json' - +filename = '../../test_data/clic_dr/line_for_taper.json' with open(filename, 'r') as f: line = xt.Line.from_dict(json.load(f)) tracker = line.build_tracker() -line['qc1l1.1..1'].ksl[0] += 1e-6 - tracker.configure_radiation(model=None) tw_no_rad = tracker.twiss(method='4d', freeze_longitudinal=True) - ############################################### # Enable radiation and compensate energy loss # ############################################### + tracker.configure_radiation(model='mean') # - Set cavity lags to compensate energy loss @@ -37,82 +33,31 @@ import matplotlib.pyplot as plt plt.close('all') -ifig = 0 -for conf in configs: - - ifig += 1 - - # Twiss(es) with radiation - tracker.config.XTRACK_CAVITY_PRESERVE_ANGLE = conf['cavity_preserve_angle'] - tw = tracker.twiss(radiation_method=conf['radiation_method'], - eneloss_and_damping=(conf['radiation_method'] != 'kick_as_co')) - tracker.config.XTRACK_CAVITY_PRESERVE_ANGLE = False - - if conf['p0_correction']: - p0corr = 1 + tracker.delta_taper - else: - p0corr = 1 - - plt.figure(ifig, figsize=(6.4*1.3, 4.8)) - plt.suptitle(f"Radiation method: {conf['radiation_method']}, " - f"p0 correction: {conf['p0_correction']}, " - f"cavity preserve angle: {conf['cavity_preserve_angle']}") - - betx_beat = tw.betx*p0corr/tw_no_rad.betx-1 - bety_beat = tw.bety*p0corr/tw_no_rad.bety-1 - max_betx_beat = np.max(np.abs(betx_beat)) - max_bety_beat = np.max(np.abs(bety_beat)) - spx = plt.subplot(2,1,1) - plt.title(f'error on Qx: {abs(tw.qx - tw_no_rad.qx):.2e} ' - r'$(\Delta \beta_x / \beta_x)_{max}$ = ' - f'{max_betx_beat:.2e}') - plt.plot(tw.s, betx_beat) - plt.ylabel(r'$\Delta \beta_x / \beta_x$') - plt.ylim(np.max([0.01, 1.1 * max_betx_beat])*np.array([-1, 1])) - plt.xlim([0, tw.s[-1]]) - - plt.subplot(2,1,2, sharex=spx) - plt.title(f'error on Qy: {abs(tw.qy - tw_no_rad.qy):.2e} ' - r'$(\Delta \beta_y / \beta_y)_{max}$ = ' - f'{max_bety_beat:.2e}') - plt.plot(tw.s, bety_beat) - plt.ylabel(r'$\Delta \beta_y / \beta_y$') - plt.ylim(np.max([0.01, 1.1 * max_bety_beat])*np.array([-1, 1])) - plt.xlabel('s [m]') - - plt.subplots_adjust(hspace=0.35, top=.85) - - plt.savefig(f'./{case_name}_fig{ifig}.png', dpi=200) - - assert np.isclose(tracker.delta_taper[0], 0, rtol=0, atol=1e-10) - assert np.isclose(tracker.delta_taper[-1], 0, rtol=0, atol=1e-10) - - assert np.allclose(tw.delta, tracker.delta_taper, rtol=0, atol=1e-6) - - assert np.isclose(tw.qx, tw_no_rad.qx, rtol=0, atol=conf['q_atol']) - assert np.isclose(tw.qy, tw_no_rad.qy, rtol=0, atol=conf['q_atol']) - - assert np.isclose(tw.dqx, tw_no_rad.dqx, rtol=0, atol=1.5e-2*tw.qx) - assert np.isclose(tw.dqy, tw_no_rad.dqy, rtol=0, atol=1.5e-2*tw.qy) - assert np.allclose(tw.x, tw_no_rad.x, rtol=0, atol=1e-7) - assert np.allclose(tw.y, tw_no_rad.y, rtol=0, atol=1e-7) +fig1 = plt.figure(1, figsize=(6.4, 4.8*1.5)) +spbet = plt.subplot(3,1,1) +spco = plt.subplot(3,1,2, sharex=spbet) +spdisp = plt.subplot(3,1,3, sharex=spbet) - assert np.allclose(tw.betx*p0corr, tw_no_rad.betx, rtol=conf['beta_rtol'], atol=0) - assert np.allclose(tw.bety*p0corr, tw_no_rad.bety, rtol=conf['beta_rtol'], atol=0) +spbet.plot(tw['s'], tw['betx']) +spbet.plot(tw['s'], tw['bety']) - assert np.allclose(tw.dx, tw.dx, rtol=0.0, atol=0.1e-3) +spco.plot(tw['s'], tw['x']) +spco.plot(tw['s'], tw['y']) - assert np.allclose(tw.dy, tw.dy, rtol=0.0, atol=0.1e-3) +spdisp.plot(tw['s'], tw['dx']) +spdisp.plot(tw['s'], tw['dy']) - if case_name == 'clic_dr' and conf['radiation_method'] != 'kick_as_co': - eneloss = tw.eneloss_turn - assert eneloss/line.particle_ref.energy0 > 0.01 - assert np.isclose(line['rf'].voltage*np.sin(line['rf'].lag/180*np.pi), eneloss/4, rtol=1e-5) - assert np.isclose(line['rf1'].voltage*np.sin(line['rf1'].lag/180*np.pi), eneloss/4, rtol=1e-5) - assert np.isclose(line['rf2a'].voltage*np.sin(line['rf2a'].lag/180*np.pi), eneloss/4*0.6, rtol=1e-5) - assert np.isclose(line['rf2b'].voltage*np.sin(line['rf2b'].lag/180*np.pi), eneloss/4*0.4, rtol=1e-5) - assert np.isclose(line['rf3'].voltage*np.sin(line['rf3'].lag/180*np.pi), eneloss/4, rtol=1e-5) +spbet.set_ylabel(r'$\beta_{x,y}$ [m]') +spco.set_ylabel(r'(Closed orbit)$_{x,y}$ [m]') +spdisp.set_ylabel(r'$D_{x,y}$ [m]') +spdisp.set_xlabel('s [m]') +fig1.suptitle( + r'$q_x$ = ' f'{tw["qx"]:.5f}' r' $q_y$ = ' f'{tw["qy"]:.5f}' '\n' + r"$Q'_x$ = " f'{tw["dqx"]:.2f}' r" $Q'_y$ = " f'{tw["dqy"]:.2f}' + r' $\gamma_{tr}$ = ' f'{1/np.sqrt(tw["momentum_compaction_factor"]):.2f}' +) -plt.show() +fig1.subplots_adjust(left=.15, right=.92, hspace=.27) +plt.show() \ No newline at end of file From 48f5d3ff2e7f9f750823ca2f252f4a097a1725be Mon Sep 17 00:00:00 2001 From: giadarol Date: Fri, 18 Nov 2022 17:27:31 +0100 Subject: [PATCH 19/20] Example and test --- .../001_taper_and_compare_twiss_methods.py | 33 ++--- tests/test_tapering.py | 123 ++++++++---------- 2 files changed, 71 insertions(+), 85 deletions(-) diff --git a/examples/tapering/001_taper_and_compare_twiss_methods.py b/examples/tapering/001_taper_and_compare_twiss_methods.py index 035213418..92f50c3ed 100644 --- a/examples/tapering/001_taper_and_compare_twiss_methods.py +++ b/examples/tapering/001_taper_and_compare_twiss_methods.py @@ -2,25 +2,26 @@ import numpy as np import xtrack as xt -case_name = 'clic_dr' -filename = '../../test_data/clic_dr/line_for_taper.json' +# case_name = 'clic_dr' +# filename = '../../test_data/clic_dr/line_for_taper.json' +# configs = [ +# {'radiation_method': 'full', 'p0_correction': False, 'cavity_preserve_angle': False, 'beta_rtol': 2e-2, 'q_atol': 5e-4}, +# {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': False, 'beta_rtol': 2e-2, 'q_atol': 5e-4}, +# {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 2e-5, 'q_atol': 5e-4}, +# {'radiation_method': 'kick_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-3, 'q_atol': 5e-4}, +# {'radiation_method': 'scale_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-5, 'q_atol': 5e-4}, +# ] + +case_name = 'fcc-ee' +filename = 'line_no_radiation.json' configs = [ - {'radiation_method': 'full', 'p0_correction': False, 'cavity_preserve_angle': False, 'beta_rtol': 2e-2, 'q_atol': 5e-4}, - {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': False, 'beta_rtol': 2e-2, 'q_atol': 5e-4}, - {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 2e-5, 'q_atol': 5e-4}, - {'radiation_method': 'kick_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-3, 'q_atol': 5e-4}, - {'radiation_method': 'scale_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-5, 'q_atol': 5e-4}, + {'radiation_method': 'full', 'p0_correction': False, 'cavity_preserve_angle': False, 'beta_rtol': 1e-2, 'q_atol': 5e-4}, + {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': False, 'beta_rtol': 5e-3, 'q_atol': 5e-4}, + {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 3e-4, 'q_atol': 5e-4}, + {'radiation_method': 'kick_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 3e-3, 'q_atol': 7e-4}, + {'radiation_method': 'scale_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-4, 'q_atol': 1e-4}, ] -# case_name = 'fcc-ee' -# filename = 'line_no_radiation.json' -# configs = [ -# {'radiation_method': 'full', 'p0_correction': False, 'cavity_preserve_angle': False, 'beta_rtol': 1e-2, 'q_atol': 5e-4}, -# {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': False, 'beta_rtol': 5e-3, 'q_atol': 5e-4}, -# {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 3e-4, 'q_atol': 5e-4}, -# {'radiation_method': 'kick_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 3e-3, 'q_atol': 7e-4}, -# {'radiation_method': 'scale_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-4, 'q_atol': 1e-4}, -# ] with open(filename, 'r') as f: line = xt.Line.from_dict(json.load(f)) diff --git a/tests/test_tapering.py b/tests/test_tapering.py index 45ac9a223..6322a7291 100644 --- a/tests/test_tapering.py +++ b/tests/test_tapering.py @@ -6,13 +6,22 @@ test_data_folder = pathlib.Path( __file__).parent.joinpath('../test_data').absolute() -def test_tapering(): - with open(test_data_folder / 'clic_dr/line_for_taper.json', 'r') as f: +def test_tapering_and_twiss_with_radiation(): + + filename = test_data_folder / 'clic_dr/line_for_taper.json' + configs = [ + {'radiation_method': 'full', 'p0_correction': False, 'cavity_preserve_angle': False, 'beta_rtol': 2e-2, 'q_atol': 5e-4}, + {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': False, 'beta_rtol': 2e-2, 'q_atol': 5e-4}, + {'radiation_method': 'full', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 2e-5, 'q_atol': 5e-4}, + {'radiation_method': 'kick_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-3, 'q_atol': 5e-4}, + {'radiation_method': 'scale_as_co', 'p0_correction': True, 'cavity_preserve_angle': True, 'beta_rtol': 1e-5, 'q_atol': 5e-4}, + ] + + with open(filename, 'r') as f: line = xt.Line.from_dict(json.load(f)) tracker = line.build_tracker() - # Initial twiss (no radiation) tracker.configure_radiation(model=None) tw_no_rad = tracker.twiss(method='4d', freeze_longitudinal=True) @@ -23,69 +32,45 @@ def test_tapering(): # - Taper magnet strengths tracker.compensate_radiation_energy_loss() - # Twiss(es) with radiation - tw = tracker.twiss(method='6d', matrix_stability_tol=3., - eneloss_and_damping=True) - tw_sympl = tracker.twiss(radiation_method='kick_as_co', method='6d') - tw_scale_as_co = tracker.twiss( - radiation_method='scale_as_co', - method='6d', - matrix_stability_tol=0.5) - - p0corr = 1 + tracker.delta_taper - - assert np.isclose(tracker.delta_taper[0], 0, rtol=0, atol=1e-10) - assert np.isclose(tracker.delta_taper[-1], 0, rtol=0, atol=1e-10) - assert np.isclose(np.max(tracker.delta_taper), 0.00568948, rtol=1e-4, atol=0) - assert np.isclose(np.min(tracker.delta_taper), -0.00556288, rtol=1e-4, atol=0) - - assert np.allclose(tw.delta, tracker.delta_taper, rtol=0, atol=1e-6) - assert np.allclose(tw_sympl.delta, tracker.delta_taper, rtol=0, atol=1e-6) - assert np.allclose(tw_scale_as_co.delta, tracker.delta_taper, rtol=0, atol=1e-6) - - assert np.isclose(tw.qx, tw_no_rad.qx, rtol=0, atol=5e-4) - assert np.isclose(tw_sympl.qx, tw_no_rad.qx, rtol=0, atol=5e-4) - assert np.isclose(tw_scale_as_co.qx, tw_no_rad.qx, rtol=0, atol=5e-4) - - assert np.isclose(tw.qy, tw_no_rad.qy, rtol=0, atol=5e-4) - assert np.isclose(tw_sympl.qy, tw_no_rad.qy, rtol=0, atol=5e-4) - assert np.isclose(tw_scale_as_co.qy, tw_no_rad.qy, rtol=0, atol=5e-4) - - assert np.isclose(tw.dqx, tw_no_rad.dqx, rtol=0, atol=0.2) - assert np.isclose(tw_sympl.dqx, tw_no_rad.dqx, rtol=0, atol=0.2) - assert np.isclose(tw_scale_as_co.dqx, tw_no_rad.dqx, rtol=0, atol=0.2) - - assert np.isclose(tw.dqy, tw_no_rad.dqy, rtol=0, atol=0.2) - assert np.isclose(tw_sympl.dqy, tw_no_rad.dqy, rtol=0, atol=0.2) - assert np.isclose(tw_scale_as_co.dqy, tw_no_rad.dqy, rtol=0, atol=0.2) - - assert np.allclose(tw.x, tw_no_rad.x, rtol=0, atol=1e-7) - assert np.allclose(tw_sympl.x, tw_no_rad.x, rtol=0, atol=1e-7) - assert np.allclose(tw_scale_as_co.x, tw_no_rad.x, rtol=0, atol=1e-7) - - assert np.allclose(tw.y, tw_no_rad.y, rtol=0, atol=1e-7) - assert np.allclose(tw_sympl.y, tw_no_rad.y, rtol=0, atol=1e-7) - assert np.allclose(tw_scale_as_co.y, tw_no_rad.y, rtol=0, atol=1e-7) - - assert np.allclose(tw.betx*p0corr, tw_no_rad.betx, rtol=0.02, atol=0) - assert np.allclose(tw_sympl.betx*p0corr, tw_no_rad.betx, rtol=0.02, atol=0) - assert np.allclose(tw_scale_as_co.betx, tw_no_rad.betx, rtol=0.02, atol=0) - - assert np.allclose(tw.bety*p0corr, tw_no_rad.bety, rtol=0.02, atol=0) - assert np.allclose(tw_sympl.bety*p0corr, tw_no_rad.bety, rtol=0.02, atol=0) - assert np.allclose(tw_scale_as_co.bety, tw_no_rad.bety, rtol=0.02, atol=0) - - assert np.allclose(tw.dx, tw.dx, rtol=0.00, atol=0.1e-3) - assert np.allclose(tw_sympl.dx, tw_no_rad.dx, rtol=0.00, atol=0.1e-3) - assert np.allclose(tw_scale_as_co.dx, tw_no_rad.dx, rtol=0.00, atol=0.1e-3) - - assert np.allclose(tw.dy, tw.dy, rtol=0.00, atol=0.1e-3) - assert np.allclose(tw_sympl.dy, tw_no_rad.dy, rtol=0.00, atol=0.1e-3) - assert np.allclose(tw_scale_as_co.dy, tw_no_rad.dy, rtol=0.00, atol=0.1e-3) - - eneloss_real_tracking = tw.eneloss_turn - assert np.isclose(line['rf'].voltage*np.sin(line['rf'].lag/180*np.pi), eneloss_real_tracking/4, rtol=1e-5) - assert np.isclose(line['rf1'].voltage*np.sin(line['rf1'].lag/180*np.pi), eneloss_real_tracking/4, rtol=1e-5) - assert np.isclose(line['rf2a'].voltage*np.sin(line['rf2a'].lag/180*np.pi), eneloss_real_tracking/4*0.6, rtol=1e-5) - assert np.isclose(line['rf2b'].voltage*np.sin(line['rf2b'].lag/180*np.pi), eneloss_real_tracking/4*0.4, rtol=1e-5) - assert np.isclose(line['rf3'].voltage*np.sin(line['rf3'].lag/180*np.pi), eneloss_real_tracking/4, rtol=1e-5) + for conf in configs: + + # Twiss(es) with radiation + tracker.config.XTRACK_CAVITY_PRESERVE_ANGLE = conf['cavity_preserve_angle'] + tw = tracker.twiss(radiation_method=conf['radiation_method'], + eneloss_and_damping=(conf['radiation_method'] != 'kick_as_co')) + tracker.config.XTRACK_CAVITY_PRESERVE_ANGLE = False + + if conf['p0_correction']: + p0corr = 1 + tracker.delta_taper + else: + p0corr = 1 + + assert np.isclose(tracker.delta_taper[0], 0, rtol=0, atol=1e-10) + assert np.isclose(tracker.delta_taper[-1], 0, rtol=0, atol=1e-10) + + assert np.allclose(tw.delta, tracker.delta_taper, rtol=0, atol=1e-6) + + assert np.isclose(tw.qx, tw_no_rad.qx, rtol=0, atol=conf['q_atol']) + assert np.isclose(tw.qy, tw_no_rad.qy, rtol=0, atol=conf['q_atol']) + + assert np.isclose(tw.dqx, tw_no_rad.dqx, rtol=0, atol=1.5e-2*tw.qx) + assert np.isclose(tw.dqy, tw_no_rad.dqy, rtol=0, atol=1.5e-2*tw.qy) + + assert np.allclose(tw.x, tw_no_rad.x, rtol=0, atol=1e-7) + assert np.allclose(tw.y, tw_no_rad.y, rtol=0, atol=1e-7) + + assert np.allclose(tw.betx*p0corr, tw_no_rad.betx, rtol=conf['beta_rtol'], atol=0) + assert np.allclose(tw.bety*p0corr, tw_no_rad.bety, rtol=conf['beta_rtol'], atol=0) + + assert np.allclose(tw.dx, tw.dx, rtol=0.0, atol=0.1e-3) + + assert np.allclose(tw.dy, tw.dy, rtol=0.0, atol=0.1e-3) + + if conf['radiation_method'] != 'kick_as_co': + eneloss = tw.eneloss_turn + assert eneloss/line.particle_ref.energy0 > 0.01 + assert np.isclose(line['rf'].voltage*np.sin(line['rf'].lag/180*np.pi), eneloss/4, rtol=1e-5) + assert np.isclose(line['rf1'].voltage*np.sin(line['rf1'].lag/180*np.pi), eneloss/4, rtol=1e-5) + assert np.isclose(line['rf2a'].voltage*np.sin(line['rf2a'].lag/180*np.pi), eneloss/4*0.6, rtol=1e-5) + assert np.isclose(line['rf2b'].voltage*np.sin(line['rf2b'].lag/180*np.pi), eneloss/4*0.4, rtol=1e-5) + assert np.isclose(line['rf3'].voltage*np.sin(line['rf3'].lag/180*np.pi), eneloss/4, rtol=1e-5) From b2e1b10d0e41c57073be451a4f48bb54de1b9d6f Mon Sep 17 00:00:00 2001 From: giadarol Date: Fri, 18 Nov 2022 18:34:20 +0100 Subject: [PATCH 20/20] Adapted test --- tests/test_tracker.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/test_tracker.py b/tests/test_tracker.py index e168698a0..582f8f529 100644 --- a/tests/test_tracker.py +++ b/tests/test_tracker.py @@ -392,6 +392,7 @@ def test_tracker_hashable_config(): ('AAA', 'ipsum'), ('TEST_FLAG_BOOL', True), ('TEST_FLAG_INT', 42), + ('XTRACK_MULTIPOLE_NO_SYNRAD', True), # active by default ('ZZZ', 'lorem'), )