PF 制御計算機が今年夏に更新されたのですが、最近、各種制御の為に SAD を連続で呼び出すと
だんだん計算機が重くなり、仕舞に落とす（ユーザープロセスで制御計算機本体が落ちていいのか
どうかはまあ、置いといて……）という現象が生じ、制御グループの方が SAD を更新しようとして
いるのですが、古い SAD と新しい SAD で計算結果が違います。
　　具体的には、design_example.sad で、emittance の計算結果が違います。関連部分のみ
抜粋すると、

MOMENTUM = 1GEV;
DRIFT L1 =(L =1 )
;
BEND B =(L =2 ANGLE =.2617993877991494 )
;
QUAD QRD =(L =.5 K1 =-.1078621179953663 )
QRF =(L =1 K1 =.19457865291109586 )
QSD =(L =1 K1 =-.21812755206494103 )
QSF =(L =1 K1 =.15519425429347913 )
QD =(L =1 K1 =-.1560624058480439 )
QF =(L =1 K1 =.1269625766602043 )
;
SEXT SF =(L =1 K2 =.12116257080688167 )
SD =(L =1 K2 =-.20401978160297035 )
;
CAVI CA1 =(L =1 VOLT =1e+06 HARM =100 )
;
LINE

sample = (IP1 -QRD -L1 -L1 -L1 -CA1 -L1
-L1 -L1 -QRF -L1 -L1 -L1 -B
-L1 -L1 -L1 -QSD -L1 -L1 -L1
-L1 -L1 -L1 -L1 -L1 -QSF L1
SF L1 B L1 L1 L1 QD
L1 SD L1 B L1 L1 L1
QF L1 SF L1 B L1 L1
L1 QD L1 SD L1 B L1
L1 L1 QF L1 SF L1 B
L1 L1 L1 QD L1 SD L1
B L1 L1 L1 QF L1 SF
L1 B L1 L1 L1 QD L1
SD L1 B L1 L1 L1 QF
L1 SF L1 B L1 L1 L1
QD L1 SD L1 B L1 L1
L1 QSF L1 L1 L1 L1 L1
L1 L1 L1 QSD L1 L1 L1
B L1 L1 L1 QRF L1 L1
L1 CA1 L1 L1 L1 QRD -QRD
-L1 -L1 -L1 -CA1 -L1 -L1 -L1
-QRF -L1 -L1 -L1 -B -L1 -L1
-L1 -QSD -L1 -L1 -L1 -L1 -L1
-L1 -L1 -L1 -QSF -L1 -L1 -L1
-B -L1 -SD -L1 -QD -L1 -L1
-L1 -B -L1 -SF -L1 -QF -L1
-L1 -L1 -B -L1 -SD -L1 -QD
-L1 -L1 -L1 -B -L1 -SF -L1
-QF -L1 -L1 -L1 -B -L1 -SD
-L1 -QD -L1 -L1 -L1 -B -L1
-SF -L1 -QF -L1 -L1 -L1 -B
-L1 -SD -L1 -QD -L1 -L1 -L1
-B -L1 -SF -L1 -QF -L1 -L1
-L1 -B -L1 -SD -L1 -QD -L1
-L1 -L1 -B -L1 -SF -L1 QSF
L1 L1 L1 L1 L1 L1 L1
L1 QSD L1 L1 L1 B L1
L1 L1 QRF L1 L1 L1 CA1
L1 L1 L1 QRD -IP1 )
;

ON LOG;
ON RAD FLUC COD RFSW RADCOD;
FFS;
$FORM = "15.10";
PageWidth = 1999;

USE sample;
cell calc;
emit;
disp $$$;

bpos = LINE["POSITION","B*"];

GetIntegralH[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;
! Print["BX = ", bex," AX = ",alx," GX = ", gax];
! Print["EX = ", etx," EPX= ",epx];
! Print["rho= ", rho," L = ",L];
aaa = 2*etx*epx*tha+2*epx*rho*(Sin[tha]-tha)+2*(etx-rho)*(1-Cos[tha])+rho/2*(1-Cos[2*tha]);
bbb = (tha+2*epx*epx*tha-Sin[2*tha]/2+4*epx*(1-Cos[tha]))/2;
ccc = (2*etx*etx*tha+4*etx*rho*(Sin[tha]-tha)+rho*rho*(3*tha-4*Sin[tha]+Sin[2*tha]/2))/2;
IH = IH + aaa * alx + bbb * bex + ccc * gax;
bex = Twiss["BX",num+0.5];
alx = -1*Twiss["AX",num+0.5];
gax = (1+alx*alx)/bex;
etx = Twiss["EX",num+0.5];
epx = -1*Twiss["EPX",num+0.5];
aaa = 2*etx*epx*tha+2*epx*rho*(Sin[tha]-tha)+2*(etx-rho)*(1-Cos[tha])+rho/2*(1-Cos[2*tha]);
bbb = (tha+2*epx*epx*tha-Sin[2*tha]/2+4*epx*(1-Cos[tha]))/2;
ccc = (2*etx*etx*tha+4*etx*rho*(Sin[tha]-tha)+rho*rho*(3*tha-4*Sin[tha]+Sin[2*tha]/2))/2;
IH = IH + aaa * alx + bbb * bex + ccc * gax;
! Print[IH];
IH
);


IntH = 0;
Do[
IntH = IntH + GetIntegralH[bpos[i]];
,{i,Length[bpos]}];

Print[IntH];
cq = 3.832/10000;
gamma = 1000/0.51099906;
emittance = cq * gamma * gamma * IntH /2 /Pi /rho;
Print[emittance];

stop;
stop;

です。下側はシンクロトロン積分の解析式でエミッタンスを計算する積分式です。
（BMの中央から、上流側、下流側の両方向へ積分式を計算し、和を求めます。）

　　古い SAD では、emit の結果は、
*** Welcome to SAD Ver.1.0.10.2b1 built at 2007-09-10 11:46:34 JST ***
Emittance X = 2.45855E-7 m
ですが、afsad、amsad、PF でコンパイルし直した新しい sad では、
*** Welcome to SAD Revision 3995 built at 2012-11-22 20:42:37 +0900 ***
Emittance X = 4.58933E-7 m
です。上記の計算式では、244.9820119666 nmrad （古い方とほぼ同じ）と出ます。
afsad、amsadでは、新しい方とほぼ同じ結果になります。
　
　　ちなみに、PF のラティスで計算すると、emittance の計算結果にこんなに
大きな差は出ません。
旧： Emittance X = 3.69803E-8 m
新： Emittance X = 3.71388E-8 m
解析式：37.1172052922 nmrad
　
　　design example のラティスは、強いて言えば若干和共鳴に近いですが、
他に何か不運なケースなのでしょうか？　何かノウハウがありましたらお願い
致します。

試しに実行される場合は、
MARK IP1 = ();
をラティスに加えてください。

> 　　PF 制御計算機が今年夏に更新されたのですが、最近、各種制御の為に SAD を連続で呼び出すと
> だんだん計算機が重くなり、仕舞に落とす（ユーザープロセスで制御計算機本体が落ちていいのか
> どうかはまあ、置いといて……）という現象が生じ、制御グループの方が SAD を更新しようとして
> いるのですが、古い SAD と新しい SAD で計算結果が違います。
> 　　具体的には、design_example.sad で、emittance の計算結果が違います。関連部分のみ
> 抜粋すると、
(snip)
> 　　古い SAD では、emit の結果は、
> *** Welcome to SAD Ver.1.0.10.2b1 built at 2007-09-10 11:46:34 JST ***
> Emittance X = 2.45855E-7 m
> ですが、afsad、amsad、PF でコンパイルし直した新しい sad では、
> *** Welcome to SAD Revision 3995 built at 2012-11-22 20:42:37 +0900 ***
> Emittance X = 4.58933E-7 m
> です。上記の計算式では、244.9820119666 nmrad （古い方とほぼ同じ）と出ます。
> afsad、amsadでは、新しい方とほぼ同じ結果になります。
> 　
> 　　ちなみに、PF のラティスで計算すると、emittance の計算結果にこんなに
> 大きな差は出ません。
> 旧： Emittance X = 3.69803E-8 m
> 新： Emittance X = 3.71388E-8 m
> 解析式：37.1172052922 nmrad
> 　
> 　　design example のラティスは、強いて言えば若干和共鳴に近いですが、
> 他に何か不運なケースなのでしょうか？　何かノウハウがありましたらお願い
> 致します。
>
仕様変更で無いのに、値が変わるのはバグである可能性があります(今回は、解析解とも異なる値です)
で、原因となるコードの変更を特定し、物理モデルとしての正しさを検証する必要があると思われます
CVS MAIN trunkに関しては、V1.0.8.19.3b以降は定期的に CVS tagを打ち込んであるので
tagをたどりながら2分法を行なえばエミッタンスが変化した commitを絞り込めます
SAD_1_0_10で始まるCVS tagを探してみてください
# V1.0.10.4.19a10から V1.0.10.5aの間は、変更アナウンスに追従しきれず tagが打てていません

PS.
amorita branchを使っているようですが、SAD Homepage移管直前あたりから source archiveは更新されていません
といか、移管後は、更新手段が未整備＆公式ページからリンク圏外なので無保守＆無保証です
# 現在 revision 4069で、TRPTモードの SOLバグパッチは未マージです

CVS repository mirrorに関しては、移管前と同じサーバーで稼働中ですが、
サーバー管理を引き継ぐ人がいないと維持出来なくなります

> > 　　PF 制御計算機が今年夏に更新されたのですが、最近、各種制御の為に SAD を連続で呼び出すと
> > だんだん計算機が重くなり、仕舞に落とす（ユーザープロセスで制御計算機本体が落ちていいのか
> > どうかはまあ、置いといて……）という現象が生じ、制御グループの方が SAD を更新しようとして
> > いるのですが、古い SAD と新しい SAD で計算結果が違います。
> > 　　具体的には、design_example.sad で、emittance の計算結果が違います。関連部分のみ
> > 抜粋すると、
> (snip)
> > 　　古い SAD では、emit の結果は、
> > *** Welcome to SAD Ver.1.0.10.2b1 built at 2007-09-10 11:46:34 JST ***
> > Emittance X = 2.45855E-7 m
> > ですが、afsad、amsad、PF でコンパイルし直した新しい sad では、
> > *** Welcome to SAD Revision 3995 built at 2012-11-22 20:42:37 +0900 ***
> > Emittance X = 4.58933E-7 m
> > です。上記の計算式では、244.9820119666 nmrad （古い方とほぼ同じ）と出ます。
> > afsad、amsadでは、新しい方とほぼ同じ結果になります。
> > 　
> > 　　ちなみに、PF のラティスで計算すると、emittance の計算結果にこんなに
> > 大きな差は出ません。
> > 旧： Emittance X = 3.69803E-8 m
> > 新： Emittance X = 3.71388E-8 m
> > 解析式：37.1172052922 nmrad
> > 　
> > 　　design example のラティスは、強いて言えば若干和共鳴に近いですが、
> > 他に何か不運なケースなのでしょうか？　何かノウハウがありましたらお願い
> > 致します。
> >
> 仕様変更で無いのに、値が変わるのはバグである可能性があります(今回は、解析解とも異なる値です)
> で、原因となるコードの変更を特定し、物理モデルとしての正しさを検証する必要があると思われます
> CVS MAIN trunkに関しては、V1.0.8.19.3b以降は定期的に CVS tagを打ち込んであるので
> tagをたどりながら2分法を行なえばエミッタンスが変化した commitを絞り込めます
> SAD_1_0_10で始まるCVS tagを探してみてください
> # V1.0.10.4.19a10から V1.0.10.5aの間は、変更アナウンスに追従しきれず tagが打てていません
>
CVS tag binary searchと file revision binary searchにて、design_example.sadのEMITXが変わった変更を特定しました

MAIN trunk V1.0.10.4a05 - V1.0.10.4a07間の変更に属する
src/trade.f revision 1.19 - 1.20間の変更で EMITXの評価値が変わっています

oide氏による version number V1.0.10.4a06(src/trade.f revision 1.20)から現在のEMITXの値になっています

cvs log src/trade.f 抜粋
----------------------------
revision 1.20
date: 2010/01/18 13:08:29;  author: oide;  state: Exp;  lines: +5 -4
V1.0.10.4a06
----------------------------

cvs diff -r1.19 -r1.20 src/trade.f
Index: src/trade.f
===================================================================
RCS file: /SAD/cvsroot/oldsad/src/trade.f,v
retrieving revision 1.19
retrieving revision 1.20
diff -d -u -r1.19 -r1.20
--- src/trade.f 17 Jan 2010 07:19:11 -0000      1.19
+++ src/trade.f 18 Jan 2010 13:08:29 -0000      1.20
@@ -1,4 +1,4 @@
-c$Header: /SAD/cvsroot/oldsad/src/trade.f,v 1.19 2010/01/17 07:19:11 oide Exp $
+c$Header: /SAD/cvsroot/oldsad/src/trade.f,v 1.20 2010/01/18 13:08:29 oide Exp $
       subroutine trade(trans,beam,cod,bx,by,bz,br,
      $     bxx,bxy,byy,dldx,
      $     f1,s,ala,al)
@@ -43,7 +43,7 @@
         call tclr(radi,36)
         radi(6,1)=-h1sq*tlc*2.d0*
      $       (pzi*(-vx*bxy+vy*bxx)+vz*(pxi*bxy-pyi*bxx))
-     1       -h1sq*brad*dtdsc*al*dldx
+     1       -h1sq*brad*dtdsc*dldx
         radi(6,2)=(dp*dtds**2*pxi-h1sq*tlc*2.d0*
      $       ((vx*by-vy*bx)*pxi/pzi-vy*bz+vz*by))/pr
         radi(6,3)=-h1sq*tlc*2.d0*
@@ -54,6 +54,7 @@
         radi(6,6)=-2.d0*p1*p0*tlc*brad
      1       -dp*dtds**2*sq/pr
      1       +2.d0*h1sq*(qx*vx+qy*vy+vz**2)/pr*tlc
+c        write(*,'(1p6g15.7)')(radi(6,i),i=1,6)
         call tradp(radi,pxi,pyi,dp,pr,ar)
         call tmuld(trans,radi)
         do i=1,6
@@ -142,8 +143,8 @@
       subroutine tradp(trans,pxi,pyi,dp,pr,ar)
       implicit none
       real*8 trans(6,6),pxi,pyi,dp,pr,ar
-      trans(6,1)=trans(6,1)-ar*trans(6,4)/pr
-      trans(6,3)=trans(6,3)+ar*trans(6,2)/pr
+      trans(6,1)=trans(6,1)-ar*trans(6,4)
+      trans(6,3)=trans(6,3)+ar*trans(6,2)
       trans(2,1)=pxi*trans(6,1)
       trans(2,2)=pxi*trans(6,2)
       trans(2,3)=pxi*trans(6,3)

design_example.sadのEMITXを変えたのは、この部分の変更でした。

変数 alは、おそらく長さの次元を持つので、radi(6,1)の第2項の次元を変える変更
(dldxは、dl/dxのことで無次元量と思われる)

> Index: src/trade.f
> ===================================================================
> RCS file: /SAD/cvsroot/oldsad/src/trade.f,v
> retrieving revision 1.19
> retrieving revision 1.20
> diff -d -u -r1.19 -r1.20
> --- src/trade.f 17 Jan 2010 07:19:11 -0000 1.19
> +++ src/trade.f 18 Jan 2010 13:08:29 -0000 1.20
> @@ -43,7 +43,7 @@
> call tclr(radi,36)
> radi(6,1)=-h1sq*tlc*2.d0*
> $ (pzi*(-vx*bxy+vy*bxx)+vz*(pxi*bxy-pyi*bxx))
> - 1 -h1sq*brad*dtdsc*al*dldx
> + 1 -h1sq*brad*dtdsc*dldx
> radi(6,2)=(dp*dtds**2*pxi-h1sq*tlc*2.d0*
> $ ((vx*by-vy*bx)*pxi/pzi-vy*bz+vz*by))/pr
> radi(6,3)=-h1sq*tlc*2.d0*

> 　　PF 制御計算機が今年夏に更新されたのですが、最近、各種制御の為に SAD を連続で呼び出すと
> だんだん計算機が重くなり、仕舞に落とす（ユーザープロセスで制御計算機本体が落ちていいのか
> どうかはまあ、置いといて……）という現象が生じ、制御グループの方が SAD を更新しようとして
> いるのですが、古い SAD と新しい SAD で計算結果が違います。
> 　　具体的には、design_example.sad で、emittance の計算結果が違います。関連部分のみ
> 抜粋すると、
(snip)
> 　　古い SAD では、emit の結果は、
> *** Welcome to SAD Ver.1.0.10.2b1 built at 2007-09-10 11:46:34 JST ***
> Emittance X = 2.45855E-7 m
> ですが、afsad、amsad、PF でコンパイルし直した新しい sad では、
> *** Welcome to SAD Revision 3995 built at 2012-11-22 20:42:37 +0900 ***
> Emittance X = 4.58933E-7 m
> です。上記の計算式では、244.9820119666 nmrad （古い方とほぼ同じ）と出ます。
> afsad、amsadでは、新しい方とほぼ同じ結果になります。
> 　
(snip)
> 　　design example のラティスは、強いて言えば若干和共鳴に近いですが、
> 他に何か不運なケースなのでしょうか？　何かノウハウがありましたらお願い
> 致します。
>
和共鳴という話が有ったので、チューンを [4.2, 4.8]付近で振って調べましたが
NX = 4.2にて EMITX 〜 6.3e-7から NXの増加に対して減少傾向にあり
NX 〜 4.47付近で極小 EMITX 〜 3.9e-7となり増加に転じる模様
(所々、NYに依存して EMITXが下がる点が有りますが EMITX > 3e-7です)

なので、近傍の動作点に EMITX 〜 2.5e-7な点は無い模様

　　森田氏から連絡を頂きました。どうもありがとうございました。
　
・スクリプトの「解析式」では、ダンピングパーティションナンバーを考慮していない。
・古い SAD では、セクタベンドに対する扱いが正しくない。
　
・新しい SAD では、それが正しく考慮されレいる。
・ダンピングパーティションナンバーを考慮して、正しくエミッタンスを計算すると、
　新しい SAD と正しい解析式はほぼ一致する。
　
とのことです。
　
結論：今の新しい SAD は正しい。修正前の SAD でコンバインドやセクタのベンドに
　　　対してエミッタンスを計算すると、正しい値にならない。
　
　　どうもありがとうございました！

> 　　森田氏から連絡を頂きました。どうもありがとうございました。
> 　
> ・スクリプトの「解析式」では、ダンピングパーティションナンバーを考慮していない。
> ・古い SAD では、セクタベンドに対する扱いが正しくない。
> 　
> ・新しい SAD では、それが正しく考慮されレいる。
> ・ダンピングパーティションナンバーを考慮して、正しくエミッタンスを計算すると、
> 　新しい SAD と正しい解析式はほぼ一致する。
> 　
> とのことです。
> 　
> 結論：今の新しい SAD は正しい。修正前の SAD でコンバインドやセクタのベンドに
> 　　　対してエミッタンスを計算すると、正しい値にならない。
> 　
> 　　どうもありがとうございました！
>

関与してる修正箇所は、以下のdldx依存項で、事実関係は以下の通り
* 例題の解析式では、Damping Partition Number Jx 〜 0.56が考量されていない
* radi(6,1)を構成する未修正の第一項と修正された第二項の次元が一致する(次元解析による)
* SOL区間外では、tbende()からの呼び出しのみ dldx != 0となる(BEND)
* dldxパラメターは、水平偏差cod(1)に応じた軌道長の伸びを表すパラメータ

したがって、Rectangler BENDな加速器では、dldx 〜 0/Jx 〜 1となり r1.19と r1.20の差は表面化しない

> Index: src/trade.f
> ===================================================================
> RCS file: /SAD/cvsroot/oldsad/src/trade.f,v
> retrieving revision 1.19
> retrieving revision 1.20
> diff -d -u -r1.19 -r1.20
> --- src/trade.f 17 Jan 2010 07:19:11 -0000 1.19
> +++ src/trade.f 18 Jan 2010 13:08:29 -0000 1.20
> @@ -43,7 +43,7 @@
> call tclr(radi,36)
> radi(6,1)=-h1sq*tlc*2.d0*
> $ (pzi*(-vx*bxy+vy*bxx)+vz*(pxi*bxy-pyi*bxx))
> - 1 -h1sq*brad*dtdsc*al*dldx
> + 1 -h1sq*brad*dtdsc*dldx
> radi(6,2)=(dp*dtds**2*pxi-h1sq*tlc*2.d0*

> * 例題の解析式では、Damping Partition Number Jx 〜 0.56が考量されていない
値を写し間違えました
解析式の評価: Jx 〜 .5287 EMITX 〜 4.8986e-7
SAD EMITの評価: Jx 〜 .5325 EMITX 〜 4.61865e-7