Zenn
Open8

加速器について

電子線形加速器
nikuniku

電子を操作するなら、電子顕微鏡か電子線形加速器かな、ということで加速器の勉強して、あわよくば法律に引っかからないレベルの加速器をつくりたい（ブラウン管だってある意味加速器だろうし）。

ということで、大昔に購入して放置した
Particle Accelerator Physics , Helmut Wiedemann (著)
https://amzn.asia/d/3DwTQSk
を読んでいくことにした。

nikuniku

まずは、prebuncherとは何ぞやということが知りたくなったので
15.1.1 Prebuncher
から読み始める。英語は苦手なので時間がかかりそう。

nikuniku

Fig.15.1. をmatplotlibで再現したい

どういうステップを踏めば再現できるか考えてみる。

草案１

初期化＝15.1 aのプロット
・横軸＝位相、縦軸＝エネルギーなので、横軸を有限の個数でサンプリング=＞位相のリストを生成
・サンプリングされた位置は一個の粒子とみなす
・各粒子は一定のエネルギーを持つ＝＞位相のリストに対応したエネルギーのリストを生成

エネルギー変調=15.1 bのプロット
・プリバンチャーで与えられるエネルギーは位相で決まっている＝＞aで作った位相のリストを生成
・位相ごとに与えるエネルギーのリストが生成できる＝＞位相のリストに対応したエネルギーのリストを生成
・aのときに作成した粒子のエネルギーのリスト＋bの与えられるエネルギーのリストを加算する

ドリフトスペースを動いた後の位相変調＝15.1 cのプロット
・bの時点で各粒子を（エネルギー, 位相）というタプルのリストで考えられる。
・ドリストスペースでは粒子のエネルギーは一定、位相はエネルギーによって変わるのでその方程式から、時間経過を与えると、どの程度位相が変わるのかの（エネルギー、位相）＋経過時間を入力して（エネルギー、位相）を出力する関数を作り、経過時間ごとに（エネルギー、位相）のタプルのリストをプロットする

草案１に対する自分のコメント

最初はエネルギーのリストと位相のリストが別々で、最終的に（エネルギー、位相）のタプルのリストになるので、最初から（エネルギー、位相）のタプルにしたほうが良い気がする

nikuniku

とりあえずFig. 15.1 aのプロット

import matplotlib.pyplot as plt
import numpy as np

# 位相のリストを作成
phi_start = 0
phi_end = np.pi
phi_sample = 100
phi_list = np.linspace(phi_start, phi_end, phi_sample)

# それぞれの初期位相のサンプリング位置を粒子とみなしE0のリストを作成
E0 = 1
E0_list = np.full(len(phi_list), E0)

# 初期状態1.5.1 aのプロット
fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(phi_list, E0_list)


Fig 15.1.a をエネルギー1[MeV]として描画

nikuniku

つぎに 15.1 bのプロット
(15.1aのプロットのコードに以下のコードを追加）

# rfの電圧
V0 = 0.1

# rf位相に対し、与えられるeV
def get_eV(phase):
    return V0*np.sin(phase)

# rf位相に対し、与えられるeVのリスト
eV_list = np.apply_along_axis(get_eV, 0, phi_list)

# rfにより粒子は加速されエネルギーが加算される
E_list = E0_list + eV_list
print(E_list)


fig = plt.figure()
ax = fig.add_subplot(111)
ax.plot(phi_list, E_list)
ax.set_xlabel("rf phase")
ax.set_ylabel("energy of particle")


Fig 15.1 b を加速電圧0.1[MeV]として上記コードでプロットしたもの

nikuniku

15.1 cの計算のために位相の変化量を計算するための各関数を書き出してみた

import numpy as np
# γを計算する関数
def get_gamma_from_E_mc2(E, mc2):
    return E/mc2

E = 1.0 # 1MeV
mc2 = 0.511 # 0.511MeV電子の静止エネルギー

gamma = get_gamma_from_E_mc2(E, mc2)
print(f"γ:{gamma}")

# Δγを計算する関数
def get_dgamma_from_dE_mc2(dE, mc2):
    return dE/mc2

dE = 0.1 # 0.1Vのrf電圧を印加
dgamma = get_dgamma_from_dE_mc2(dE, mc2)
print(f"Δγ：{dgamma}")

# βを計算する関数
def get_beta_from_gamma(gamma):
    return(np.sqrt(1-1/(gamma**2)))

beta = get_beta_from_gamma(gamma)
print(f"β：{beta}")

# Δβを計算する関数
def get_dbeta_from_gamma_dgamma_beta(gamma, dgamma, beta):
    return dgamma/(beta*(gamma**3))

dbeta = get_dbeta_from_gamma_dgamma_beta(gamma, dgamma, beta)
print(f"Δβ：{dbeta}")

# Δvを計算する関数
c = 299792458.0 #光速m/s
def get_dv_from_dbeta(dbeta):
    return 299792458.0*dbeta

dv = get_dv_from_dbeta(dbeta)
print(f"Δv：{dv}")

# 位相の変化Δφを計算する関数
def get_dphi_from_dv_dt_wavelength(dv, dt, wavelength):
    return 2*np.pi*dv*dt/wavelength
    
wavelength = c/2856E+6
dphi = get_dphi_from_dv_dt_wavelength(dv, 1E-9, wavelength)
print(f"Δφ：{dphi}")
nikuniku

クソコードだけど一応描画

dE_phi = np.vstack([eV_list, phi_list])

print(dE_phi.shape)
wavelength = c/2856E+6
E = 1.0 # 1MeV
mc2 = 0.511 # 0.511MeV電子の静止エネルギー
dt = 1E-9

def update_dE_phi(E, mc2, dt, wavelength, de_phi):
    dphi = get_dphi_from_E0_dE_mc2_dt_wavelength(E, de_phi[0], mc2, dt, wavelength)
    return np.array([de_phi[0], de_phi[1]+dphi])


fig = plt.figure()
ax = fig.add_subplot(111)

ax.set_xlabel("rf phase")
ax.set_ylabel("energy of particle")

updated_dE_phi = np.apply_along_axis(lambda de_phi: update_dE_phi(E, mc2, 0, wavelength, de_phi), 0, dE_phi)
ax.plot(updated_dE_phi[1], updated_dE_phi[0]+E)
updated_dE_phi = np.apply_along_axis(lambda de_phi: update_dE_phi(E, mc2, 1E-9, wavelength, de_phi), 0, dE_phi)
ax.plot(updated_dE_phi[1], updated_dE_phi[0]+E)
updated_dE_phi = np.apply_along_axis(lambda de_phi: update_dE_phi(E, mc2, 2E-9, wavelength, de_phi), 0, dE_phi)
ax.plot(updated_dE_phi[1], updated_dE_phi[0]+E)
updated_dE_phi = np.apply_along_axis(lambda de_phi: update_dE_phi(E, mc2, 3E-9, wavelength, de_phi), 0, dE_phi)
ax.plot(updated_dE_phi[1], updated_dE_phi[0]+E)
updated_dE_phi = np.apply_along_axis(lambda de_phi: update_dE_phi(E, mc2, 4E-9, wavelength, de_phi), 0, dE_phi)
ax.plot(updated_dE_phi[1], updated_dE_phi[0]+E)