# coding: utf-8 # ## Generate field map # # \$z\$ is the "vertical" direction, corresponding to the normal vector of the cyclotron plane import h5py import numpy as np from matplotlib import pyplot as plt from mpl_toolkits.mplot3d import Axes3D def makeConstantFieldmap(filename, Bmag=1.0, rmin=0.0, rmax=300.0, theta_min=0.0, theta_max=np.pi/6, zmin=0.0, zmax=1000.0): """ Create a constant field on a rectangular grid, and save to file. """ r = np.linspace(rmin, rmax, 20) theta = np.linspace(theta_min, theta_max, 40) z = np.linspace(zmin, zmax, 20) R, T, Z = np.meshgrid(r, theta, z) X = R*np.cos(T) Y = R*np.sin(T) def fieldmap(x, y, z): def Bx(x, y, z): return 0.0*np.ones(shape=x.shape) def By(x, y, z): return 0.0*np.ones(shape=x.shape) def Bz(x, y, z): return Bmag*np.ones(shape=x.shape) return (Bx(x, y, z), By(x, y, z), Bz(x, y, z)) Bx, By, Bz = fieldmap(X, Y, Z) cols = np.vstack((X.flatten(), Z.flatten(), Y.flatten(), Bx.flatten(), Bz.flatten(), By.flatten()) ).T with open(filename, 'w') as f: f.write('\n'.join(['', 'X [mm]', 'Z [mm]', 'Y [mm]', 'Bx [T]', 'Bz [T]', 'By [T]', '', ''])) np.savetxt(f, cols) # plt.figure(figsize=(6,6)) # plt.plot(X.flatten(), Y.flatten(), 'b+') # plt.title('Points specified in fieldmap') # plt.xlabel('x (cm)') # plt.ylabel('y (cm)') # plt.xlim([-rmin,rmax]) # plt.ylim([-rmin,rmax]) # # plt.figure() # ax = plt.subplot(111, projection='3d') # ax.scatter(X[::], Y[::], zs=Z[::]) # ax.set_xlabel('x') # ax.set_ylabel('y') # ax.set_zlabel('z') return X, Y, Z E = 0.1 # in GeV PMASS = 0.938 gamma = (E+PMASS)/PMASS beta = np.sqrt(1. - 1./gamma**2) P0 = gamma*beta*PMASS Brho = 3.3356 * P0 # in T/m R = 1.0 X, Y, Z = makeConstantFieldmap("testbend.bmap", Bmag=Brho*R, rmin=300, rmax=600)