...
function [isok, BDES0, BDES1, BDES2] = E300_calc_QS_3(z_ob, z_im, QS, m12_req, m34_req)keep_KQS0_eq_KQS2=1; % set to one if QS0=QS2
E0 = 10;
isok = 1;
if(nargin < 4); m12_req = 0; end;
if(nargin < 5); m34_req = 0; end;% initial guesses (2012 values)
KQS0_0 = -0.3;
KQS1_0 = 0.23;
KQS2_0 = -0.3;
mytol = (0.01^2 + 0.01^2 );% linac z locations of QS1 and QS2
z_QS0 = 1996.98249; % [m], middle of quad
z_QS1 = 1999.206615; % [m], middle of quad
z_QS2 = 2001.431049; % [m], middle of quadLEFF_QS0 = 1; % [m]
LEFF_QS1 = 1; % [m]
LEFF_QS2 = 1; % [m]OO = zeros(2,2);
d1 = (z_QS0-LEFF_QS0/2) - z_ob ;
d2 = (z_QS1-LEFF_QS1/2) - (z_QS0+LEFF_QS0/2);
d3 = (z_QS2-LEFF_QS2/2) - (z_QS1+LEFF_QS1/2);
d4 = z_im - (z_QS2+LEFF_QS2/2);M_01 = [1 d1; 0 1];
M4_01 = [M_01 OO; OO M_01];
M_02 = [1 d2; 0 1];
M4_02 = [M_02 OO; OO M_02];
M_03 = [1 d3; 0 1];
M4_03 = [M_03 OO; OO M_03];
M_04 = [1 d4; 0 1];
M4_04 = [M_04 OO; OO M_04];
options = optimset('TolX',1e-8);
if keep_KQS0_eq_KQS2
[fit_result, chi2] = fminsearch(@transportError_2, [KQS0_0 KQS1_0],options);
BDES0 = fit_result(1) * (E0+QS) * LEFF_QS0 / 0.0299792458;
BDES1 = fit_result(2) * (E0+QS) * LEFF_QS1 / 0.0299792458;
BDES2 = fit_result(1) * (E0+QS) * LEFF_QS2 / 0.0299792458;
else
[fit_result, chi2] = fminsearch(@transportError, [KQS0_0 KQS1_0 KQS2_0],options);
BDES0 = fit_result(1) * (E0+QS) * LEFF_QS0 / 0.0299792458;
BDES1 = fit_result(2) * (E0+QS) * LEFF_QS1 / 0.0299792458;
BDES2 = fit_result(3) * (E0+QS) * LEFF_QS2 / 0.0299792458;
end
if(chi2 > mytol)
isok = 0;
warning('could not converge to solution');
BDES0 = NaN;
BDES1 = NaN;
BDES2 = NaN;
endBMAX = 385; % max value, from SCP
if(abs(BDES1) > BMAX || abs(BDES2) > BMAX || abs(BDES0) > BMAX)
isok = 0;
warning('solution is outside QS range');
end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function chi2 = transportError(K)
% QS0 transport matrix
k = abs(K(1));
phi = LEFF_QS0*sqrt(k);
M_F = [cos(phi) (1/sqrt(k))*sin(phi)
-sqrt(k)*sin(phi) cos(phi)];
M_D = [cosh(phi) (1/sqrt(k))*sinh(phi)
sqrt(k)*sinh(phi) cosh(phi)];
M4_D1 = [M_D OO; OO M_F];
% QS1 transport matrix
k = abs(K(2));
phi = LEFF_QS1*sqrt(k);
M_F = [cos(phi) (1/sqrt(k))*sin(phi)
-sqrt(k)*sin(phi) cos(phi)];
M_D = [cosh(phi) (1/sqrt(k))*sinh(phi)
sqrt(k)*sinh(phi) cosh(phi)];
M4_F = [M_F OO; OO M_D];
% QS2 transport matrix
k = abs(K(3));
phi = LEFF_QS2*sqrt(k);
M_F = [cos(phi) (1/sqrt(k))*sin(phi)
-sqrt(k)*sin(phi) cos(phi)];
M_D = [cosh(phi) (1/sqrt(k))*sinh(phi)
sqrt(k)*sinh(phi) cosh(phi)];
M4_D2 = [M_D OO; OO M_F];% dump line optics
M4 = M4_04*M4_D2*M4_03*M4_F*M4_02*M4_D1*M4_01;
chi2 = (M4(1,2)-m12_req)^2 + (M4(3,4)-m34_req)^2;
endfunction chi2 = transportError_2(K)
% QS0 transport matrix
k = abs(K(1));
phi = LEFF_QS0*sqrt(k);
M_F = [cos(phi) (1/sqrt(k))*sin(phi)
-sqrt(k)*sin(phi) cos(phi)];
M_D = [cosh(phi) (1/sqrt(k))*sinh(phi)
sqrt(k)*sinh(phi) cosh(phi)];
M4_D1 = [M_D OO; OO M_F];
% QS1 transport matrix
k = abs(K(2));
phi = LEFF_QS1*sqrt(k);
M_F = [cos(phi) (1/sqrt(k))*sin(phi)
-sqrt(k)*sin(phi) cos(phi)];
M_D = [cosh(phi) (1/sqrt(k))*sinh(phi)
sqrt(k)*sinh(phi) cosh(phi)];
M4_F = [M_F OO; OO M_D];
% QS2 transport matrix
k = abs(K(1));
phi = LEFF_QS2*sqrt(k);
M_F = [cos(phi) (1/sqrt(k))*sin(phi)
-sqrt(k)*sin(phi) cos(phi)];
M_D = [cosh(phi) (1/sqrt(k))*sinh(phi)
sqrt(k)*sinh(phi) cosh(phi)];
M4_D2 = [M_D OO; OO M_F];% dump line optics
M4 = M4_04*M4_D2*M4_03*M4_F*M4_02*M4_D1*M4_01;
chi2 = (M4(1,2)-m12_req)^2 + (M4(3,4)-m34_req)^2;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%end
Python code to reproduce MATLAB settings above (Nov 13, 2023)
##########################################################################
# setup
###############################################################################
#% linac z locations of QS1 and QS2
#z_QS0 = 1996.98249; % [m], middle of quad
#z_QS1 = 1999.206615; % [m], middle of quad
#z_QS2 = 2001.431049; % [m], middle of quadxstart = 1993.27370 # [m] FILS
xstart = 1992.82000 # [m] PIC_CENTq0d_pos = 1996.98249 # [m]
q0d_length = 1.00 # [m]
q0d_G = 1.0 # [T/m]q1d_pos = 1999.206615 # [m]
q1d_length = 1.00 # [m]
q1d_G = 1.0 # [T/m]q2d_pos = 2001.431049 # [m]
q2d_length = 1.00 # [m]
q2d_G = 1.0 # [T/m]
xend = 2017.529977792559 # [m] PRDMP
xend = 2015.259850655461 # [m] DTOTR
xend = 2015.629843995481 # [m] LFOV
def get_free_matrix(length):
return np.array([[1.0, length], [0.0, 1.0]])
def get_quad_matrix(length, G, eps):
# length: [m]
# G: [T/m]
# eps: [GeV]
# omega^2 = qGc^3/(\gamma mc^2)
# phi = \omega L /c
# sqrt((Tesla / meter) * (elementary charge ) * (speed of light)^3 / GeV)/ (speed of light)
# \approx 0.5475 (1/meter)
omegaoverc = np.sqrt(G/eps) * 0.5475 #[1/m]
phi = omegaoverc * length
return np.array([[np.cos(phi), (1.0/omegaoverc) * np.sin(phi)], [-omegaoverc * np.sin(phi), np.cos(phi)]])def get_quad_matrix_h(length, G, eps):
# length: [m]
# G: [T/m]
# eps: [GeV]
omegaoverc = np.sqrt(G/eps) * 0.5475 #[1/m]
phi = omegaoverc * length
return np.array([[np.cosh(phi), (1.0/omegaoverc) * np.sinh(phi)], [omegaoverc * np.sinh(phi), np.cosh(phi)]])
eps = 8.0 #[GeV]#d1 = (z_QS0-LEFF_QS0/2) - z_ob ;
#d2 = (z_QS1-LEFF_QS1/2) - (z_QS0+LEFF_QS0/2);
#d3 = (z_QS2-LEFF_QS2/2) - (z_QS1+LEFF_QS1/2);
#d4 = z_im - (z_QS2+LEFF_QS2/2);
def get_transport(G_in, G_out, eps, getx = True):
G_in = np.abs(G_in)
G_out = np.abs(G_out)current = get_free_matrix(q0d_pos-q0d_length/2.0-xstart)
next = None
if(getx == True):
next = get_quad_matrix_h(q0d_length, G_out, eps)
else:
next = get_quad_matrix(q0d_length, G_out, eps)
current = np.matmul(next, current)
next = get_free_matrix((q1d_pos-q1d_length/2.0)-(q0d_pos+q0d_length/2.0))
current = np.matmul(next, current)
if(getx == True):
next = get_quad_matrix(q1d_length, G_in, eps)
else:
next = get_quad_matrix_h(q1d_length, G_in, eps)
current = np.matmul(next, current)
next = get_free_matrix((q2d_pos-q2d_length/2.0)-(q1d_pos+q1d_length/2.0))
current = np.matmul(next, current)
if(getx == True):
next = get_quad_matrix_h(q2d_length, G_out, eps)
else:
next = get_quad_matrix(q2d_length, G_out, eps)
current = np.matmul(next, current)
next = get_free_matrix(xend-(q2d_pos+q2d_length/2.0))
current = np.matmul(next, current)
return current
#####
# To understand Matlab code:
# (omega/c)^2 = k = qG/(\gamma mc)
# G = (\gamma mc^2) k / (q c) -> 1e9 energy (Volt) * k / c
#
# 299792458 / 1e9 = 0.299792458
# 1 kilogauss = 0.1 Tesla
#
# Tesla / Volt = second / meter^2 -> Tesla = Volt * second / meter^2
# BDES0 = k * (E0+QS) * LEFF_QS0 / 0.0299792458;
#from scipy.optimize import fmin
def minimize(x):
eps = 7.0
Gin = x[0]
Gout = x[1]
tmpx = get_transport(Gin, Gout, eps, getx = True)
tmpy = get_transport(Gin, Gout, eps, getx = False)
return (tmpx[0,1]-0.0)**2 + (tmpy[0,1]-0.0)**2res = fmin(minimize, [0.23,-0.3], xtol=1e-12, ftol=1e-12)
print(res)
print(res * 10.0) # KiloGauss instead of Tesla