放射積分の値を表示する方法はありませんでしょうか？
RADINT というコマンドはありますが、まともに動作している
感じがしません。
どうぞよろしくお願い申しあげます。
放射積分計算のスクリプト
（そのまま実行できる条件）
・B は水平ベンドのみ。
・名前は B* で、それ以外なく、B* すべては蹴り角がゼロではない。
・B の角度は ANGLE で入っている。
・逆ベンドはあっても大丈夫にしたつもり。
・最後の resemit の単位は nmrad です。
（そのまま実行できないとき）
・垂直ベンドあり → 大幅な？スクリプト書き換え。
・蹴り角ゼロのベンドがある → ANGLE = 1e-15 とか、超微小の蹴り角を入れればいい。
ON LOG RAD FLUC COD RFSW RADCOD;
FFS;
!NPARA = 4;
$FORM = "15.10";
PageWidth = 1999;
USE RING;
cell calc emit;
bposi = LINE["POSITION","B*"];
GetIntegralCH[num_]:=(
IH = 0;
! Print["Middle Parameters of ",LINE["NAME",num]];
bex = Twiss["BX",num+0.5];
alx = Twiss["AX",num+0.5];
gax = (1+alx*alx)/bex;
etx = Twiss["EX",num+0.5];
epx = Twiss["EPX",num+0.5];
L = LINE["L",num]/2;
tha = LINE["ANGLE",num]/2;
rho = L/tha;
k1 = -1*LINE["K1",num];
k = k1/L/2;
kef2= k - 1/rho/rho;
! Print[LINE["NAME",num]];
! Print["BX = ", bex," AX = ",alx," GX = ", gax];
! Print["EX = ", etx," EPX= ",epx];
! Print["rho= ", rho," L = ",L, " tha = ",tha ];
! Pring["k1 = ", k1, " k = ", k];
! Print["kef= ", kef2];
! FFS["Abort;"];
If[kef2>0,
kef = Sqrt[kef2];
aaa = 2*(-8*kef^2*rho*epx*Sinh[kef*L]+8*kef^3*(kef^2*etx*rho+1)*rho*epx*L)/4/kef^5/rho^2;
bbb = 2*(kef^2*Sinh[2*kef*L]+4*kef^5*rho^2*epx^2*L-2*kef^3*L)/4/kef^5/rho^2;
ccc = 2*(-8*(1+kef^2*rho*etx)*Sinh[kef*L]+6*kef*L+4*kef^5*rho^2*etx^2*L+Sinh[2*kef*L]+8*kef^3*rho*etx*L)/4/kef^5/rho^2;
! Print["abc = ",aaa,bbb,ccc];
IH = IH + aaa * alx + bbb * bex + ccc * gax;
,
kef = Sqrt[-kef2];
aaa = 2*(8*kef^2*rho*epx*Sin[kef*L]+8*kef^3*(kef^2*etx*rho-1)*rho*epx*L)/4/kef^5/rho^2;
bbb = 2*(-kef^2*Sin[2*kef*L]+4*kef^5*rho^2*epx^2*L+2*kef^3*L)/4/kef^5/rho^2;
ccc = 2*(-8*(1-kef^2*rho*etx)*Sin[kef*L]+6*kef*L+4*kef^5*rho^2*etx^2*L+Sin[2*kef*L]-8*kef^3*rho*etx*L)/4/kef^5/rho^2;
! Print["abc = ",aaa,bbb,ccc];
IH = IH + aaa * alx + bbb * bex + ccc * gax;
! Print[IH];
];
IH = IH /rho/rho/Abs[rho];
IH
);
GetIntegralEX[num_]:=(
IE = 0;
! Print["Middle Parameters of ",LINE["NAME",num]];
etx = Twiss["EX",num+0.5];
epx = Twiss["EPX",num+0.5];
L = LINE["L",num]/2;
tha = LINE["ANGLE",num]/2;
rho = L/tha;
k1 = -1*LINE["K1",num];
k = k1/L/2;
kef2= k - 1/rho/rho;
! Print[LINE["NAME",num]];
! Print["BX = ", bex," AX = ",alx," GX = ", gax];
! Print["EX = ", etx," EPX= ",epx];
! Print["rho= ", rho," L = ",L, " tha = ",tha ];
! Pring["k1 = ", k1, " k = ", k];
! Print["kef= ", kef2];
! FFS["Abort;"];
If[kef2>0,
kef = Sqrt[kef2];
IE = ((etx*kef^2*rho+1)*Sinh[kef*L]+epx*kef*rho*Cosh[kef*L]-kef*L-kef*rho*epx)/kef^3/rho;
etx = Twiss["EX",num+0.5];
epx = -1*Twiss["EPX",num+0.5];
IE = IE + ((etx*kef^2*rho+1)*Sinh[kef*L]+epx*kef*rho*Cosh[kef*L]-kef*L-kef*rho*epx)/kef^3/rho;
,
kef = Sqrt[-kef2];
IE = ((etx*kef^2*rho-1)*Sin[kef*L]-epx*kef*rho*Cos[kef*L]+kef*L+kef*rho*epx)/kef^3/rho;
etx = Twiss["EX",num+0.5];
epx = -1*Twiss["EPX",num+0.5];
IE = IE + ((etx*kef^2*rho-1)*Sin[kef*L]-epx*kef*rho*Cos[kef*L]+kef*L+kef*rho*epx)/kef^3/rho;
! Print[IE];
];
IE
);
Jx = 0;
Jy = 0;
Jz = 0;
DFac = 0;
RadIntI1 = 0;
RadIntI2 = 0;
RadIntI3 = 0;
RadIntI4 = 0;
RadIntI5 = 0;
Do[
! Print[i," ",LINE["NAME",bposi[i]]," ",LINE["K1",bposi[i]]==0];
IH = GetIntegralCH[bposi[i]];
RadIntI5 = RadIntI5 + IH;
! Print[bposi[i]," ",LINE["NAME",bposi[i]]," ",IH,aaa,bbb,ccc,alx,bex,gax,L/tha];
,{i,Length[bposi]}];
Do[
IE = GetIntegralEX[bposi[i]];
L = LINE["L",bposi[i]];
tha = LINE["ANGLE",bposi[i]];
rho = L/tha;
k1 = -LINE["K1",bposi[i]]/L;
RadIntI1 = RadIntI1 + IE / rho;
RadIntI2 = RadIntI2 + tha*tha/L;
RadIntI3 = RadIntI3 + L / rho / rho / Abs[rho];
RadIntI4 = RadIntI4 + IE*(1/rho^3-2*k1/rho);
! Print[bposi[i], IE];
,{i,Length[bposi]}];
Print[RadIntI1,RadIntI2,RadIntI3,RadIntI4,RadIntI5];
DFac = RadIntI4 / RadIntI2;
sademit = Emittance[];
Jx = 1 - DFac;
Jy = 1; !!!!!!! only varid for the case without vertical bendings.
Jz = 2 + DFac;
Print[Jx,Jy,Jz];
cq = 3.832/10000;
gamma = MOMENTUM / 1000000/0.51099906;
resemit = cq * gamma * gamma * RadIntI5 / RadIntI2 / Jx;
Print[resemit];
Print[(Emittances/.sademit)[1]];
stop;
stop;
> 放射積分の値を表示する方法はありませんでしょうか？
> RADINT というコマンドはありますが、まともに動作している
> 感じがしません。
> どうぞよろしくお願い申しあげます。
>
"DETERMINANT=0 MATRIX-INV NOT FOUND"で落ちるのは、src/TCAL.fの R0, R0I配列の未初期化参照バグ
もっとも、BEND/QUAD/SEXT/CAVIしか数えてないので、MULTによる衝突点モデルやSOL+MULTによる変更電磁石モデルなどはまともに計算出来ないと思われる
> 放射積分の値を表示する方法はありませんでしょうか？
> RADINT というコマンドはありますが、まともに動作している
> 感じがしません。
> どうぞよろしくお願い申しあげます。
>
6極の計算まわりがバグバグしてる
* src/SEXTU.fで正しくCOMMONブロックが参照されていない
* src/INTGRL.fでSEXTU呼び出し前にCOMMONブロックの変数を正しく初期化していない
* UNITM4の第二引数がREAL(8)型なのに、単精度のリテラルが積んである
その辺直して、script/design_example.sadを計算すると
tune/emittance/compaction等はEMITと概ね一致する
RADINTのコード解析まとめ
バグ
* BODYをステップ評価する際に2回目以降、TCAL内部でR0/R0Iの未初期化参照が発生する
** 暗黙のSAVE属性を与える Fortranコンパイラでは動いてしまう
* UINTM4のREAL(8)型の第2引数にREAL(4)型のリテラルを渡している
* SEXTUがTMATRコモンブロックを参照していない
* SEXTU呼び出し前に必要なコモン変数(DL/DKL2)を初期化していない
** SEXTに関しては、1995年のCVSによる管理開始時から動かなかったものと思われる
放射積分評価と対象パラメータ
* BEND(L, ANGLE, E1, E2)
** BODYは、シンプソン積分
** EDGEは、thin lens評価
* QUAD(L, K1)
** BODYは、シンプソン積分
** EDGE無し
* SEXT(L, K2)
** BODYは、中点評価
** EDGE無し
* CAVI(VOLT,HARM,FREQ)
** 電圧の積算
** HARM/FREQは、最終値を控えるのみ
"すぐに分かる"制限事項
* 反転設置の非対称BENDの評価が正しくない
* K0/AE1/AE2を持ったBENDの評価は正しくない
* エレメントのDX/DYオフセットは評価されない
* SEXT内の軌道による転送行列の摂動は評価されない
* "LINEAR CHROMATICITY"に ROTATEが反映されない
** 例えば、90° ROTATEかつ K1反転した QUADを並べると"LINEAR CHROMATICITY"が反転する
表示について
* "LINEAR CHROMATICITY"は、無次元化された(2πで割った) natural linear chromaticity