Conference Room SAD

[Thread display] [Reply] [New arrival display] [Word search] [Notes] [管理用]

No.	： 727
Title	： Re: 放射積分
Date	： 2018/03/26(Mon) 13:01:44
Name	： K.Harada <kentaro.harada@kek.jp>
URL	：

放射積分計算のスクリプト

（そのまま実行できる条件）
　・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;

- Related List tree (Click ▼ to display all articles in a thread.）

▼

放射積分

K.Harada

No.726

Re: 放射積分 - K.Harada 18/03/26-13:01 No.727
Re: 放射積分 - Akio Morita 18/03/27-14:57 No.729
Re: 放射積分 - Akio Morita 18/03/28-00:55 No.730
Re: 放射積分 - Akio Morita 18/03/28-10:47 No.731
- Re^2: 放射積分 - K.Harada 18/03/28-11:41 No.732

- Reply Form (Reply to the above using the following.)

Name
E-Mail
Subject
Message
> 放射積分計算のスクリプト > > （そのまま実行できる条件） > 　・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;
URL
Password (alphanumeric characters within 8)
Subscribe Key (Enter the number of the right image)

- WebForum -

Name
E-Mail
Subject
Message > 放射積分計算のスクリプト > > （そのまま実行できる条件） > 　・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+alxalx)/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 = -1LINE["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(-8kef^2rhoepxSinh[kefL]+8kef^3(kef^2etxrho+1)rhoepxL)/4/kef^5/rho^2; > bbb = 2(kef^2Sinh[2kefL]+4kef^5rho^2epx^2L-2kef^3L)/4/kef^5/rho^2; > ccc = 2(-8(1+kef^2rhoetx)Sinh[kefL]+6kefL+4kef^5rho^2etx^2L+Sinh[2kefL]+8kef^3rhoetxL)/4/kef^5/rho^2; > ! Print["abc = ",aaa,bbb,ccc]; > IH = IH + aaa * alx + bbb * bex + ccc * gax; > , > kef = Sqrt[-kef2]; > aaa = 2(8kef^2rhoepxSin[kefL]+8kef^3(kef^2etxrho-1)rhoepxL)/4/kef^5/rho^2; > bbb = 2(-kef^2Sin[2kefL]+4kef^5rho^2epx^2L+2kef^3L)/4/kef^5/rho^2; > ccc = 2(-8(1-kef^2rhoetx)Sin[kefL]+6kefL+4kef^5rho^2etx^2L+Sin[2kefL]-8kef^3rhoetxL)/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 = -1LINE["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 = ((etxkef^2rho+1)Sinh[kefL]+epxkefrhoCosh[kefL]-kefL-kefrhoepx)/kef^3/rho; > etx = Twiss["EX",num+0.5]; > epx = -1Twiss["EPX",num+0.5]; > IE = IE + ((etxkef^2rho+1)Sinh[kefL]+epxkefrhoCosh[kefL]-kefL-kefrhoepx)/kef^3/rho; > , > kef = Sqrt[-kef2]; > IE = ((etxkef^2rho-1)Sin[kefL]-epxkefrhoCos[kefL]+kefL+kefrhoepx)/kef^3/rho; > etx = Twiss["EX",num+0.5]; > epx = -1Twiss["EPX",num+0.5]; > IE = IE + ((etxkef^2rho-1)Sin[kefL]-epxkefrhoCos[kefL]+kefL+kefrhoepx)/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 + thatha/L; > RadIntI3 = RadIntI3 + L / rho / rho / Abs[rho]; > RadIntI4 = RadIntI4 + IE(1/rho^3-2k1/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;
URL
Password	(alphanumeric characters within 8)
Subscribe Key	(Enter the number of the right image)