広島大学の原田寛之です。

粒子の初期分布を
a = {{0},{0},{0},{0},{0},{0.01},{1}};
とし、いくつかの計算より、進行方向の変位ｚを求め、その差を比較しました。
しかしながら、TrackParticlesによる結果を理解することができません。
ソースコードの方を眺めてみましたが、理解するには時間がかかりそうです。
そこで、今回、TrackParticlesによる進行方向の計算に関して、お聞きしたいことがあり、
ご連絡いたしました。

計算手法として以下の（１）〜（３）を行いました。

手法（１）：TrackParticlesを用いた計算
TrackParticlesを用いて、初期分布を1回周回させると、運動量のずれの分、同期位相Z=0から
ずれが生じます。
現在計算されるスリッページファクターは、η＜０であるため、zは増加します。
ハーモニックナンバーを１としているため、位相は２πz/Lとなります。Lは周長です。
＜実際の計算＞
b = TrackParticles[{1,a},1];
z = b[5,1];

手法（２）：スリッページファクターηからの計算
ηとΔp/pから進行方向の変位zは、
z = -η Δp/p L
となりますので、以下のように計算しました。
＜実際の計算＞
slip = MomentumCompaction /. Emittance[Emittance->Flase][6];
dp = 0.01;
L = 348.333;
z = -slip * dp * L;

手法（３）：６次元行列による計算
SADから6次元行列を取り出し、行列計算で変位zを求めました。
＜実際の計算＞
m = TransferMatrices /. Last[Emittance[Matrix->True]];
m = Last[m]; 　　!!!　一周の行列の導出
b = m.Take[a,{1,6}];
z = b[5];

上記の（１）、（２）、（３）を比較すると、（２）と（３）は一致しますので、zは進行方向の
位相だと考えております。しかし、（１）つまりTrackParticlesにおけるzが何を意味し、計算を
行っているのか、わかりません。
何か特別なおまじないが抜けているのか、最初からTrackParticlesにおけるzの定義が異なるのか
教えていただきたいです。

お手数ですが、宜しく御願い致します。

　ｚは設計粒子に対する進行方向の位置のズレを[m]単位で表したものだと思っています。
0.01は、ビーム進行方向に 1cm 後ろにずらして計算を始めることになります。位相に
直す場合、RF の波長で割ってやれば位相になると思います。
　　（思います、と書いているのは確信がないからですが……済みません。）　計算例も
作れると思うのですが、後日。

原田健太郎様、

ご返信いただきありがとうございます。

私もｚは同期粒子（設計粒子）に対する進行方向の位置のずれを[m]単位で表したものだと
思っておりましたが、実際に同期位相(z=0)でΔp/p0が0.01つまり1%のずれを持つ粒子を
TrackParticles[]で１周させた場合の結果が、６次元行列による結果とスリッページ
ファクターからの導いた結果に対して違いが見られているため、ｚの定義に違いが
あるのか、もしくはTrackParticles[]を用いる際に、例えばRFSWのようなおまじないが
別途必要であるのかと思っております。

説明が不得手でわかりにくかったとは思いますが、宜しく御願いいたします。

原田寛之
J-PARC 3GeV RCS

(2)と(3)が一致する件ですが、たしか内部的にはEMITモジュールによって 1周分の
Transfer Matrixを計算した結果を MomentumCompactionに焼き直しているから
当然の結果だと
#<η/ρ>ではなくγ^-2を含んだ Slipage Factorが返ってきたと思います...

不一致の程度が示されていないので、一般論ですが
TrackParticlesに関しては TRACKモジュールを使っているので、計算エンジンが異なる
EMITと一般には結果が一致しないと思います
また、RADや FLUCが ONになっている場合は、違う結果が出るはずです
#もちろん、SRが無く線形近似が成立する範囲では、数値誤差の範囲で一致すべきものですが

森田様、

ご返事いただき、ありがとうございます。

> (2)と(3)が一致する件ですが、たしか内部的にはEMITモジュールによって 1周分の
> Transfer Matrixを計算した結果を MomentumCompactionに焼き直しているから
> 当然の結果だと
> #<η/ρ>ではなくγ^-2を含んだ Slipage Factorが返ってきたと思います...
現在、陽子で計算しているので、γ^-2によってSlippage Factorが変化しており、
出力される数値は手計算とも一致しております。

>
> 不一致の程度が示されていないので、一般論ですが
> TrackParticlesに関しては TRACKモジュールを使っているので、計算エンジンが異なる
> EMITと一般には結果が一致しないと思います
> また、RADや FLUCが ONになっている場合は、違う結果が出るはずです
> #もちろん、SRが無く線形近似が成立する範囲では、数値誤差の範囲で一致すべきものですが
一昔前に、別のことでお伺いした際に、計算エンジンが異なることはお聞きしましたので、
私もその違いだと思っております。
現在、RADやFLUCはOFFにしております。
参考までに不一致の程度を以下に示します。
町田さんのようにリンクせずにすみません。

synchroous momentum[eV]　　 (1)Z[m] 　　(2)Z[m]　　 (3)Z[m]　　 (1)/(2)
610254963.10968　　 2.38986　　 2.40630　　 2.40630　　 0.99317
694211799.05375 　　2.19205 　　2.20951 　　2.20951 　　0.99210
899417856.20906 　　1.75522 　　1.77374 　　1.77374 　　0.98956
1165299128.23321 　　1.31059 　　1.32851 　　1.32851 　　0.98651
1462640381.25931 　　0.95775 　　0.97399 　　0.97399 　　0.98332
1776864907.58976 　　0.70360 　　0.71797 　　0.71797 　　0.97999
2097421271.47762 　　0.52658 　　0.53931 　　0.53931 　　0.97640
2414624275.10440 　　0.40400 　　0.41543 　　0.41543 　　0.97248
2718990723.23688 　　0.31868 　　0.32913 　　0.32913 　　0.96824
3001306749.98293 　　0.25887 　　0.26860 　　0.26860 　　0.96378
3252898246.42824 　　0.21682 　　0.22602 　　0.22602 　　0.95929
3465939332.85819 　　0.18749 　　0.19632 　　0.19632 　　0.95505
3633741981.58432 　　0.16764 　　0.17620 　　0.17620 　　0.95139
3751003986.03732 　　0.15519 　　0.16359 　　0.16359 　　0.94864
3814004148.65103 　　0.14894 　　0.15725 　　0.15725 　　0.94710

森田様、

お手数をおかけしてしまい、申し訳ないです。

> (2)と(3)が一致する件ですが、たしか内部的にはEMITモジュールによって 1周分の
> Transfer Matrixを計算した結果を MomentumCompactionに焼き直しているから
> 当然の結果だと
> #<η/ρ>ではなくγ^-2を含んだ Slipage Factorが返ってきたと思います...
>
> 不一致の程度が示されていないので、一般論ですが
> TrackParticlesに関しては TRACKモジュールを使っているので、計算エンジンが異なる
> EMITと一般には結果が一致しないと思います
> また、RADや FLUCが ONになっている場合は、違う結果が出るはずです
> #もちろん、SRが無く線形近似が成立する範囲では、数値誤差の範囲で一致すべきものですが

森田さんにお聞きするのが正しいかどうかはわかりませんが、TRACKモジュールに関する計算
手法がどこかに記録として残っていないのでしょうか？もしくは、ソースコード群のどの
ファイルがTRACKモジュールに相当するのかを教えていただけないでしょうか？
（例えば、tracke.fなど)
参考までに、現在使用しているSADは、amoritaバージョンです。

大変お手数ですが、宜しく御願い致します。

原田　寛之

> 森田さんにお聞きするのが正しいかどうかはわかりませんが、TRACKモジュールに関する計算
> 手法がどこかに記録として残っていないのでしょうか？
>
たぶん間違い
歴史的には最古に属するコードなので、私は知らない
生出さんに聞いて見るべし！

#ひょっとすると、最初期のコードに関する論文が出てくるかもしれませんが、
#CVSのcommitlogにも改訂履歴がきちんと残ってないので、どの様な意図で
#修正されてきたかはたぶん分りません(Bugfixなのか、近似の改良なのかの区別も...)
#変更理由に関しては、直接コードを書いた人が忘れている場合は、全ては闇の中 Orz

> もしくは、ソースコード群のどのファイルがTRACKモジュールに相当するのかを
> 教えていただけないでしょうか？（例えば、tracke.fなど)
>
単一のファイルではないし、ファイル名の規則性も完全ではないので見分けるのは簡単では無いし
末端だけ見ても何も分らないので subroutineコールの連鎖を grepしながらズルズルと
手繰り寄せるのが早いと思います
#同系列のsubroutineでも微妙に渡される引数が違うとか、引数の解説コメントが無いので
#呼び出し元も含めて読まないと意味が分からない

TrackParticle[]に渡す座標系と内部変数の対応関係は、TrackParticleの実装を読むべし

一応、ファイル名の良くある規則
tf*.f -> FFSの subroutine
itf*.f -> FFSの function
t*e.f -> EMITモジュール関連
t*.f -> TRACKモジュール関連

Debuggerでトラップ仕込んでコールトレースを取る手も有りますが、原始的な grepでたどる
(global等の ソースコードタグ付けツールを併用すると効果的です)のに、いくら何でも
src/MAIN.fから読み始めるのは骨なのでコードをたどる上でのヒント

o SADScript側の関数からたどる場合
src/tfinitn.fのitfunallocから関数IDを拾って src/tfefun*.fの table gotoから実体を見つける

o FFSコマンドからたどる場合
src/tffsa.fのコマンドラインパーサーから実体を見つける

実際のビームラインエレメントごとにトラッキングを行ってる部分ですと、tturn1()@src/ttrun.fが
エントリポイントです。内部のtable gotoでビームラインエレメント毎の subroutineへ分岐して
行きます。例えば、QUADの実装の代表的なエントリポイント tquad()@src/tquad.f
エレメント毎のパラメータテーブルを参照するマジックナンバーは、src/initb1.fで kytbl(*, *)配列に
設定されてますが、kytbl(*,*)を使わずに即値でかかれているコードが多いため可読性・保守性が
損なわれています(いつも読んでる人は、頭にテーブルが入ってるんでしょうが...) Orz

森田様、

ご丁寧に説明いただき、ありがとうございます。

森田さんのヒントから追っていくコードをしぼることができそうです。
読んでもわからない、もしくは見つからない場合は、直接生出さんに
コンタクトを取って、ご教授を御願いしてみます。

原田　寛之

すでに原田さんへはメールで伝えてありますが、一応ポストしておきます。図の添付をどうするのか不明なので、図はありません。


この問題は、高次のmomentum compaction (slippage factor)
によるものと思われます。CALC(TWiss["DZ","$$$")で求めた
dzと、dp/p=0のときのm[[5,6]]にdp/pをかけた
ものとの差を、dp/pに対してプロットすると、添付図のよう
になります。このカーブはdp/p=0では傾きがゼロになってい
るので、m[[5,6]]の計算は正しいと思われます。

高次のslippage factorはlatticeの非線形と1/
gamma^2のdp/pに対する非線形の両方から生じます。比較のために
300 GeVの場合のプロットを二枚目に示します。この場合は1/
gamma^2の項は小さいので、latticeの非線形によるものだけ
が見えていると思われます。180 MeVの場合は1/gamma^2
の項の方が大きいようです。

ちなみに作図のスクリプトは以下の通りです。

e=Emittance[OneTurnInformation->True];
m=OneTurnTransferMatrix/.e;
z={1,{{0},{0},{0},{0},{0},{0},{1}}};

l=Table[dp=i*0.01/10;z[[2,6,1]]=dp;z1=TrackParticles[z,1];
{dp,z1[[2,5,1]]-m[[5,6]]*dp},
{i,-10,10}];
ListPlot[l,GridLines->{Automatic,Automatic},PlotJoined->True,
FrameLabel->{"`fD`np/p","`fD`nz - `fD`nz`dLinear`n","P = 180 MeV/c"}]

Update[];

なお、TrackParticlesとCALC(Twiss["DZ"])ではわず
かな差がありますが、これは例えば六極の計算にTrackParticles
の方にだけより高次の項が入っているためだと思われます。