SAD でproton synchrotron で lattice の計算と tracking を考えています。
GUIを構築するような manual は多数見かけますが、加速器の計算にかんする manual らしい物が見つからないのですが、存在しないのでしょうか？

加速器計算関連でReference Manual的な物は
http://acc-physics.kek.jp/SAD/SADHelp.HTML

後は、SAD WorkshopのProceedingsやOHOの教科書での
記載ぐらいしか見たことが有りません
##個人的に公開されてるチュートリアル類は有るのかな?

コメントありがとうございます。

> 加速器計算関連でReference Manual的な物は
> http://acc-physics.kek.jp/SAD/SADHelp.HTML

SAD を噂で聞いたことがある程度のレベルの僕には、これだけでは無理そうです。

> 後は、SAD WorkshopのProceedingsやOHOの教科書での
> 記載ぐらいしか見たことが有りません

OHO のテキストを入手してみます。

放射光源の原田です。ＳＡＤの解説を書いてある人にメールしたことがあるん
ですが、それをそのまま引用してみます。お役に立てば幸いです。
あっているかどうか保証しません。（エキスパートの皆様へ：いい加減なこと
書いてあっても、責めないで下さい……）

その１：ＳＡＤの走らせ方と必要なファイル
　
（実行方法と渡すファイル）　
　
　　ＳＡＤの実行形式は sad1.exe で、実行すると、最初に初期化の為の
スクリプト init.n など自動で読み込みます。初期化やＳＡＤ自体の変数
定義、ＧＵＩの為のボタンやスイッチ画像などはＳＡＤフォルダの下のサブ
フォルダにあります。ＳＡＤを実行する為には、どのフォルダに何がある
かの環境変数定義がまず必要で、それを含めたＳＡＤ実行の為のスクリプトが
"gs" です。（普通は gs はゴーストスクリプトを呼んでしまいますので、
path を変えるか、full path で呼んでください。）

＊＊＊＊＊＊＊ｇｓ＊＊＊＊＊＊＊＊＊＊＊＊＊　

#!/bin/csh -f
cd .
setenv SADROOT /SAD
limit stacksize 6250
setenv TCL_LIBRARY ""
setenv TK_LIBRARY ""
setenv UNIX95 ''
setenv SAD_PACKAGES $SADROOT/share/Packages
setenv KBFRAMEDIR $SADROOT/share/KBFrame
ln -fs $1 fort.77
setenv ftn77 $1
setenv FORT77 $ftn77
$SADROOT/bin/`$SADROOT/bin/HostArch`/sad1.exe "OFF CTIME LOG;READ 77; " $2 $3
　（ここで sad1.exe を呼び出し、コマンド（とラティス）を書いたファイルを
　渡します。"HostArch"は CPU (OS) を判定するスクリプトで、ここでは "Linux" が
　返されます。つまり、 $SADROOT/bin/`$SADROOT/bin/HostArch`/sad1.exe は、
　/SAD/bin/Linux/sad1.exe と同じこと。）
exit

＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
　
　　ＳＡＤファイルでラティスとコマンドは完全に分けることが可能です。
（コマンド上からラティスを定義することも可能ですが、普通はまず
ラティスを定義し、その上で計算をします。）　個人的にはラティスと
コマンドを分けたファイルにしていますが、それは単に最初にＳＡＤを
教えてくれた人がそうしていたというだけで、両者を繋げて扱っている
人も大勢います。（２分割する場合、

cat $1 >/scratch/$1_$2.exe
cat $2 >>/scratch/$1_$2.exe
/SAD/bin/gs < /scratch/$1_$2.exe

などと別にシェルスクリプトを用意して、ファイルを１つにまとめてから
gs に渡す必要があります。
　　いずれにせよ、gs に渡すファイルは１つで、前半（ラティス部分）と
後半（コマンド部分）に分けられます。
　
（ラティスファイル）
　
　　以下は先日のメールで書いたラティスファイルです。これは「トップ
レベル」で読み込まれます。ＳＡＤプログラムには２つの段階があり、
「トップレベル」と「ＦＦＳ」です。コマンドの最初に大抵は ffs と
書きますが、この"ffs"以降はＳＡＤの中で ffs という「スクリプト言語
インタープリタ」が立ち上がった状態になります。「トップレベル直接」と
「ＦＦＳ上」ではラティスもコマンドも書き方（文法）が違います。
（どちらでも同じ内容の計算が出来ることもありますが、ラティスは
トップレベルで渡した方が簡単、コマンドはＦＦＳを使った方が簡単です。）
ＳＡＤのコマンドと言う場合、多くはＦＦＳ上の関数、ＦＦＳ上の
スクリプトになります。
　　なお、ＦＦＳとは Final Focus System （最終収束系）のことで、
生出さんがリニアコライダーの衝突点オプティクス設計を前提に開発した
言語なのでそう呼ばれています。
　
　
（ラティスファイルの例）
　　（トップレベルに読み込まれます。）
!********************************************************************
!*
!* NSLS XRING
!*
!********************************************************************
　　（行の中で！以降はコメントアウトになります。トップレベルでもＦＦＳ
でも同じです。）
MOMENTUM=2.8GeV
;
MARK MKN =()
MKL =()
　　（マーク。計算では "LINE（ラティスのある部分）" を選んで計算をします。
　　　その"LINE"の最初は必ず MARK 要素である必要があります。MARK 要素は
　　　ＳＡＤの為の要素で、中身を表示させる（例えば type MKN;）と、
　　　 MARK MKN =(AX =1.00353875799E-14 BX =1.5971069053692
　　　 AY =-1.045786698E-15 BY =.4347850999073 EX =-3.10233731805E-6
　　　 EPX =-1.2621474416E-15 DP =.01 SIGZ =.0555606415224
　　　 EMITX =1.262173642876E-7 EMITY =1E-12 )
　　　などとなります。つまり、その点のトゥイスパラメータ、中心軌道、運動量の
　　　ずれ、マッチングの時に考える運動量幅、エミッタンスなどが格納されています。）
;
CAVI CAVI=(VOLT= 1.12MV HARM = 30)
　　　（空洞。周波数、もしくはハーモニック数のどちらかを入れるともう片方が
　　　　自動的に決まります。ちなみに58.88MHz だとハーモニクスが整数にならな
　　　　かったのですが、このファイルは周長が違っているでしょうか。）
;
BEND B =(L = 2.7 K1 = -0.0724921758 ANGLE = .39269908 E1 = 0.5 E2=0.5)
　　　（偏向磁石。ステアリングとメインの偏向磁石の区別はありません。変数の定義で
　　　　区別します。ANGLE は「中心軌道」を曲げ、K0 は中心軌道は曲げずに、
　　　　「ＣＯＤ」を作ります。ANGLE も K0 も rad で数字としては同じですが、
　　　　ステアリングは K0 に角度を入れます。勿論、「度」でもよくて、
　　　　ANGLE = 22.5 DEG でもよいです。E1、E2 は端の角度で、E1 が入り口、E2 が
　　　　出口、偏向電磁石の曲げ角に対しての割合で入力します。
　　　　E1=0.5 E2=0.5 が矩形、E1=0 E2=0 はセクタになります。
　　　　ちなみに偏向磁石を２に割った場合、
　　　　　BH1 =(L = 2.7/2 K1 = -0.0724921758/2 ANGLE = .39269908/2 E1=1 E2=0)
　　　　　BH2 =(L = 2.7/2 K1 = -0.0724921758/2 ANGLE = .39269908/2 E1=0 E2=1)
　　　　の様に、前半と後半で E の値を変える必要があります。）
;
QUAD
Q1 =(L = .45 K1 = -0.675839592 )
Q2 =(L = .8 K1 = 1.069849888 )
Q3 =(L = .4132 K1 = -0.57926284872 )
Q4 =(L = .225 K1 = 0.2923738695 )
　　　（４極電磁石。特に難しい話はありません。Ｆは正、Ｄは負。Ｋ１は
　　　　仰るとおり、 (B'L)/(B rho)です。以下、K2もK3も同じです。）
;
SEXT SF =(L =0 K2 = 0 )
SD =(L =0 K2 = 0 )
　　　（長さを０．２にしたら、両脇のドリフトスペースを０．１ずつ削ればいいだけ
　　　　です。実は薄レンズかどうかでＳＡＤ内部で呼び出すサブルーチンが変わるの
　　　　ですが、コマンドやラティス上は何も変わりません。）
;
DRIFT D1 =(L = 2.25 )
D2 =(L = 0.685 )
D3 =(L = 0.3484 )
D4 =(L = 0.70825)
D5 =(L = 0.9 )
D6 =(L = 0.25 )
　　（ドリフトスペースです。）
;
LINE
NCELLH =(MKN D1 Q1 D2 Q2 D3 Q3 D4
B
D5 SD D5 SF D6 Q4
MKL)
NCELL =(NCELLH -NCELLH)
RING =(4*NCELL CAVI 4*NCELL )
　　　（計算に使うラインを定義します。ＳＡＤではビームラインといいます。
　　　　頭に - を付けると、後ろからひっくり返して（偏向磁石の E1 と E2 の
　　　　方向なども含めて）並べます。すなわち、-MARKを見ると、アルファや
　　　　分散ダッシュの符号は MARK の逆になります。同じセルを並べたい場合、
　　　　4*NCELL などとすることができます。（かけ算は頭に付けること。））
）
;
　
　　以上がラティスの定義です。

ラティス例：

!********************************************************************
!*
!* NSLS XRING
!*
!********************************************************************

MOMENTUM=2.8GeV
;
MARK MKN =()
MKL =()
;
CAVI CAVI=(VOLT= 2MV HARM = 288)
;
BEND B =(L = 2.7 K1 = -0.0724921758 ANGLE = .39269908 E1 = 0.5 E2=0.5)
;
QUAD
Q1 =(L = .45 K1 = -0.675839592 )
Q2 =(L = .8 K1 = 1.069849888 )
Q3 =(L = .4132 K1 = -0.57926284872 )
Q4 =(L = .225 K1 = 0.2923738695 )
;
SEXT SF =(L =0 K2 = 0 )
SD =(L =0 K2 = 0 )
;
DRIFT D1 =(L = 2.25 )
D2 =(L = 0.685 )
D3 =(L = 0.3484 )
D4 =(L = 0.70825)
D5 =(L = 0.9 )
D6 =(L = 0.25 )
;
LINE
NCELLH =(MKN D1 Q1 D2 Q2 D3 Q3 D4
B
D5 SD D5 SF D6 Q4
MKL)
NCELL =(NCELLH -NCELLH)
RING =(4*NCELL CAVI 4*NCELL )
;

ON LOG;
ON RAD FLUC COD RFSW RADCOD;
FFS;
USE RING;
cell calc;
emit;
disp;
stop;
stop;

マニュアル１：加速器特有の計算コマンド関連（ＷＥＢのみ）
http://acc-physics.kek.jp/SAD/SADHelp.HTML
マニュアル２：ＧＵＩの作り方、スクリプト言語（ＦＦＳ）の文法基礎
SAD/Tkinter Manual (in Japanese, .pdf 706K 10/29/1997)
http://acc-physics.kek.jp/SAD/SADTkinter.pdf
　
　　その２：基本的なコマンドいくつか
　
　　まずは最も基本的なコマンドについていくつか述べておきます。
　
１、オプティクスや各種パラメータの計算
　
＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
ON LOG;
ON RAD FLUC COD RFSW RADCOD;
　　（これらはおまじないです。LOG はいちいちコマンドを復唱しないという
　　　フラグ、あとはトラッキングにおいて放射あり、縦方向を入れて６次元
　　　計算あり、などの条件を決めます。（FLAGはマニュアル１に載ってます。）
FFS;
　　（ここからいよいよＦＦＳです。）
$FORM = "15.10";
　　（小数を全１５桁にする。すなわち、整数部４桁＋"."＋小数部１０桁）
PageWidth = 1999;
　　（ファイルに書き出す場合、勝手に改行を入れない。ファイル読み書きの
　　　仕方についてはその３で書きます。その際、これを指定しないと200文字
　　　くらいで勝手に改行されて、エクセルなどで読み込む時に泣けます。）

USE RING;
　　　（ラティスファイルの LINE 定義の中で RING を使う。）
cell;
　　　（"cell" は単なるフラグ（条件判断変数）です。CELL でなければ INS に
　　　　なります。
　　　　INS は MARK の中の トゥイスパラメータを初期値にして、それを順に転送
　　　　していくモード、CELL は周期解を求めるモードです。リングの途中だけ
　　　　使いたい場合など、CELL で全体を計算、MARK にパラメータをセーブして、
　　　　その後で直線部を USE で呼び出し、INS で計算、等とします。）
calc;
　　　（"calc" は閉軌道（ＣＯＤ）を求めた後、周期解（トゥイスパラメータの
　　　　計算）を求めるコマンドです。閉軌道はあってもオプティクスは発散する
　　　　場合もありますし、誤差などを入れた場合にはそもそも何をどうやっても
　　　　軌道が閉じない場合もあります^_^;）
emit;
　　　（"emit"はエミッタンス、バンチ長など各種平衡パラメータを計算します。
　　　　calc まではＸＹの４次元ですが、emit は RF が入って６次元になります。
　　　　頭の呪文（ON RAD FLUC COD RFSW RADCOD;）を変えるとＲＦ関連を弾くことも
　　　　できます。）
SAVE ALL;
　　　（パラメータを保存しておきます、一応。）
stop;
stop;
　　　（最初の stop で FFS が終わり、次の stop でトップレベルが終わります。）

＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊＊
　
　　ちなみにＦＦＳ上では後ろに [] がつくかどうかで「コマンド」か「関数」か、
区別があります。cell calc emitは全て「コマンド」で、「関数」ではありません。
関数は大文字小文字を区別しますが、コマンドは区別しません。CELL でも cellでも
同じです。
「コマンド」を関数の中で使う場合には、
FFS["cell calc"];
などと、FFS 関数を使って呼び出す必要があります。例えば、

Do[
USE RING;
cell calc emit;
,{i,1,10}];

とやっても何も起きません。（ＳＡＤは何もせずにひっそり受け流してくれる。）

Do[
FFS["USE RING;"];
FFS["cell calc emit;",6];
,{i,1,10}];

としないといけなません。FFS["中身",結果表示先]
（装置番号６番がディスプレイです。）

２、よく使うコマンド＆関数リファレンス

・disp
　その場所のトゥイスパラメータを表示します。場所を指定しないと LINE 全部を
書き出します。

例：
disp Q4;
結果：
disp Q4;
AX BX NX EX EPX Element Length Value s(m) AY BY NY EY EPY DetR #
-3.0889 10.4015 .56819 1.37415 .41081 Q4.1 .22500 .2923739 10.404850 .76890 2.33132 .37149 .00000 .00000 .0000 15
5.1E-14 1.59711 9.14379 -3.1E-6 -1.E-15 $$$ .00000 0 170.077600 1.6E-14 .43479 6.20232 .00000 .00000 .0000 258
　　（必ず LINE の最後（リング終点）の値もくっついて出力されます。
　　左から、水平アルファ、ベータ、チューン、分散、分散ダッシュ、ラティスの
　　要素名、その要素の長さ、その要素の第１変数（Ｂなら曲げ角 rad、ＱならK1、
　　ＳならK2、ドリフトならＬ）、その要素までの軌道長、垂直アルファ、ベータ、
　　チューン、分散、分散ダッシュ、カップリングを考えた時の変換行列行列式、
　　です。）
・type
　Ｂ、Ｑ、ＳＸなど要素の中身を表示。

例： type Q4;
;
結果：
QUAD Q4 =(L =.225 K1 =.2923738695 )
;

・Twiss[欲しいパラメータ、場所]
　その場所のトゥイスパラメータを表示します。
　欲しいパラメータ： AX, BX, EX, EPX, AY, BY, EY, EPY、DX、DY （などWEB上のHELP参照。）
　場所：要素名（＝ライン上に最初に出てくるその要素）、ライン上の要素名、
　　　　　　　　　　　　ライン上の要素番号（小数で中を指定するのも可）などで指定。
　　　例えばラティスではＢは１回しか定義せず、それを１６回使います。ＳＡＤでは
　　　それを後ろに数字を付けて区別します。B.1、B.2……B.15、B.16など。ただ単に
　　　Twiss["BX","B"] とやるとTwiss["BX","B.1"]と同じです。
　　　また、Twiss の DX、DY はＣＯＤ（中心軌道）です。
例１：
Print[Twiss["BX","Q10A"]];
結果：
10.4014995745
例２：
Print[Twiss["BY","Q10A*"]];
結果：
{ 2.3313239517, 2.1620176780, 2.3313239517, 2.1620176780,......}
　（リングにQ4は１６個。（２分割しているので……Q4Hとかすべきでした。））
例３：
Print[Twiss[{"BX","BY"},{129,129.5,130}]];
　（Q4.1 は 15 番。Q4 の入り口、中央、出口の BX, BY を表示。）
結果：
{{ 10.4014995745, 10.9315382133, 11.1121567045},{ 2.3313239517, 2.2036557040, 2.1620176780}}
　（最初の３つが BX、後ろの３つが BY）
　　
・LINE[欲しい値、要素]
　その要素のパラメータを表示、もしくは代入します。
　欲しい値：POSITION、NAME、ELEMENT、S、L、ANGLE、K0、K1、K2、GEO、DX、DY
　ここでの DX, DY は据え付け誤差または設置のオフセットです。
　
例１：
Print[LINE[{"POSITION","NAME","ELEMENT","S","L","K1"},"Q4"]];
結果；
{ 15.0000000000,"Q4.1","Q4", 10.4048500000, .2250000000, .2923738695}
ライン上の番号（15番目）、名前、ファミリー名、軌道長（入り口で10.40m）、
　　　　　　　Ｑの長さ（.225cm）、B'L/(B rho) = 0.2923
例２：
Print[LINE["GEO",16]];
結果：
（GEOは幾何学上の場所を表示します。座標原点はLINEの始点です。要するに、最初の
　Q4終端（セル中央）は、x（東西方向）=10.38m、y（南北方向）=1.394mの場所にあり
　ます。）
{{ 10.3878134062, 1.3939705855, .0000000000},
{ .3926990800, .0000000000, .0000000000}}
　{{X,Y,Z},{θ、φ、ψ？^_^;}}　（ラティスが水平面内にしかないなら、
　θだけありません。角度３つは、ビームの進行方向を向いた単位ベクトルの
　３次元の角度？……HELPを参照。）
例３：
Print[LINE["LENGTH"]];
結果：
258.0000000000
　　リング（ビームライン）の要素数。
　

****************************************************************************
input
****************************************************************************

!********************************************************************
!*
!* NSLS XRING
!*
!********************************************************************

MOMENTUM=2.8GeV
;
MARK MKN =()
MKL =()
;
MONI PM =()
;
CAVI CAVI=(VOLT= 1.12MV HARM = 30 )!FREQ = 58.88MHZ)
;
BEND B =(L = 2.7 K1 = -0.0724921758 ANGLE = .39269908 E1 = 0.5 E2=0.5)
;
QUAD
Q1 =(L = .45 K1 = -0.675839592 )
Q2 =(L = .8 K1 = 1.069849888 )
Q3 =(L = .4132 K1 = -0.57926284872 )
Q4 =(L = .225 K1 = 0.2923738695 )
;
SEXT SF =(L =0 K2 = 0 )
SD =(L =0 K2 = 0 )
;
DRIFT D1 =(L = 2.25 )
D2 =(L = 0.685 )
D3 =(L = 0.3484 )
D4 =(L = 0.70825)
D5 =(L = 0.9 )
D6 =(L = 0.25 )
;
LINE
NCELLH =(MKN D1 Q1 D2 Q2 D3 Q3 D4
B
D5 SD D5 SF D6 Q4
MKL)
NCELL =(NCELLH -NCELLH)
RING =(4*NCELL CAVI 4*NCELL )
;
!
! display for tk
!
ON LOG;
ON RAD FLUC COD RFSW RADCOD;

FFS;
Get["modules_general.n"];
$FORM = "15.10";
PageWidth = 1999;

USE RING;
cell calc emit;
type MKN;
disp Q4;
type Q4;
Print[Twiss["BX","Q4"]];
Print[Twiss["BY","Q4*"]];
Print[Twiss[{"BX","BY"},{15,15.5,16}]];
Print[LINE[{"POSITION","NAME","ELEMENT","S","L","K1"},"Q4"]];
Print[LINE["GEO",16]];
Print[LINE["LENGTH"]];
disp all;

stop;
stop;

　　その３：リング基礎パラメータの計算（色収差の計算）
　
　　とりあえず６極ゼロなので、線形＆非線形の色収差を求めてみます。
　　
ON LOG;
ON RAD FLUC COD RFSW RADCOD;
FFS;
$FORM = "15.10";
PageWidth = 1999;
USE RING;
cell calc emit;
disp $$$;
tunex = Table[{0,0},{41}];
tuney = Table[{0,0},{41}];
　　　（リストの定義。-2% から +2% まで、0.1% 区切りでベータトロンチューンを
　　計算してみます。点の数は 41 点です。Table[] はリストを作る為のコマンドです。
　　これで{{0,0},{0,0},{0,0}.....,{0.0}}という多重のリストが作れます。多重で
　　なければ、定義は必要なしでいきなり代入して使えます。（色々そうです。この
　　曖昧さは好みの分かれるところです。））
Do[
DP0 = (i - 21) * 0.001;
　　　　　　（DP0 はビームの重心の運動量を変えます。FFS の「変数」です。必ず
　　　　　大文字です。同じく DP というのもありますが、こちらはフィッティングする
　　　　　際のエネルギーの幅を表します。小文字の dp は勝手に定義して使って大丈夫
　　　　　ですが、大文字の DP や DP0 は予約変数なので使ってはまずいです。）
FFS["cell calc"];
　　　　　　（Do[] の中なので FFS[]で囲う必要があります。FFS["cell calc",6]と
　　　　　すると、いちいち Matched..... と計算結果を出します。）
tunex[i] = {DP0,Twiss["NX","$$$"]/2/Pi};
tuney[i] = {DP0,Twiss["NY","$$$"]/2/Pi};
　　　　　　（リストに値を代入します。リスト名[n] はそのリストのｎ番目の要素を
　　　　　指します。なお、disp や オプティクスのマッチングで使う NX、NY は２πで
　　　　　割ったチューンですが、Twiss[]の出力は割ってありません。）
Print[i,DP0,tunex[i,2],tuney[i,2]];
　　　　　　（リスト名[n,m]は、多重のリストにおいて、ｎ番目リストのｍ番目要素を
　　　　　表します。）
,{i,1,41,1}];

Print[PolynomialFit[tunex,3]];
Print[PolynomialFit[tuney,3]];

　（PolynomialFit[リスト、次数] で、リストを多項式フィットします。
結果は　y = C0 + C1 x + C2 x2 + C3 x3... で、{{C0,C1,C2,C3..},{残差}}です。
今回の場合、
水平
{{ 9.1435987631, -22.3461488358, 116.3539887101, -861.430745192},{(Resi...)}}
　horizontal tune = 9.1435 - 22.3461 DP/P + 116.3540 (DP/P)2 ......
垂直
{{ 6.2021978838, -16.5611445238, 47.7018126339, -430.677007587},{(Resi...)}}
　vertical tune = 6.2022 - 16.5611 DP/P + 47.7018 (DP/P)2 ......
なので、水平線形色収差は -22、垂直は -16 となります。（このフィットの結果では
必ずしも DP/P=0 で近似曲線と実際のチューンが重なりません。））

　（この後で計算を続ける様な場合、DP0 = 0; と戻し忘れて悩むことがあります。）

stop;
stop;
　

!********************************************************************
!*
!* NSLS XRING
!*
!********************************************************************

MOMENTUM=2.5GeV
;
MARK MKN =()
MKL =()
;
MONI PM =()
;
CAVI CAVI=(VOLT= 1.12MV HARM = 30 )
;
BEND B =(L = 2.7 K1 = -0.0724921758 ANGLE = .39269908 E1 = 0.5 E2=0.5)
;
QUAD
Q1 =(L = .45 K1 = -0.675839592 )
Q2 =(L = .8 K1 = 1.069849888 )
Q3 =(L = .4132 K1 = -0.57926284872 )
Q4 =(L = .225 K1 = 0.2923738695 )
;
SEXT SF =(L =0.2 K2 = 0 )
SD =(L =0.2 K2 = 0 )
;
DRIFT D1 =(L = 2.25 )
D2 =(L = 0.685 )
D3 =(L = 0.3484 )
D4 =(L = 0.70825)
D5 =(L = 0.8 )
D7 =(L = 0.7 )
D6 =(L = 0.15 )
;
LINE
NCELLH =(MKN D1 Q1 D2 Q2 D3 Q3 D4
B
D5 SD D7 SF D6 Q4
MKL)
NCELL =(NCELLH -NCELLH)
RING =(4*NCELL CAVI 4*NCELL )
;

!
! display chromaticity
!
ON LOG;
ON RAD FLUC COD RFSW RADCOD;

FFS;
$FORM = "15.10";
PageWidth = 1999;

USE RING;
cell calc emit;
disp $$$;

tunex = Table[{0,0},{41}];
tuney = Table[{0,0},{41}];

Do[
DP0 = (i - 21) * 0.001;
FFS["cell calc"];
tunex[i] = {DP0,Twiss["NX","$$$"]/2/Pi};
tuney[i] = {DP0,Twiss["NY","$$$"]/2/Pi};
Print[i,DP0,tunex[i,2],tuney[i,2]];
,{i,1,41,1}];

Print[PolynomialFit[tunex,3]];
Print[PolynomialFit[tuney,3]];

stop;
stop;

　　その４：ＳＡＤによるオプティクスのグラフの書き方

（準備：既知の未解決エラー）
　　
１、topdrawer形式
　
（main branch の SAD では必ず起きるエラーがあります。森田氏はもしか
すると、直してくれているかも知れません。）
　
・縦軸凡例の文字（√βｘとか）が必ず文字化けします。（DEC の
　true unix 上の topdrawer では大丈夫なんですが……）
　原因（未解決）
　topdrawer形式のファイルの中の
　 TITLE 1.10 7.33 CENTER ANGLE 90. SPACE= 7.1 '2B06O0X1 (2M062O1)'
　 CASE 'MGUUDXLX MLUUUDU '
　が軸の凡例指定です。上の行の'2B06....'が文字指定、CASE はその文字の
　役割（大文字、小文字、数学記号……）です。これが、Linux 版の
　topdrawer で解読できないみたいです。下のエラーも、true unix の
　topdrawer では起きません。
　
・軸のラベル大きさが「負」だというエラーメッセージが出ます。
　Ｘで使っていればいいんですが、teraterm でグラフを書かせると、
　上書きされてグラフが見えなくなります。でも、これは解決できます。
　twsdraw.f と tdlat.f の軸ラベルの"-"を取ってしまえばいいです。
　例えば、以下の如し。
sed s/"sngl(-2.4"/"sngl(2.4"/g oldsad/src/twsdrw.f
sed s/"sngl(-2.0"/"sngl(2.0"/g oldsad/src/twsdrw.f
sed s/"TITLE SIZE -1.6"/"TITLE SIZE 1.6"/g oldsad/src/twsdrw.f
sed s/"tsize=-3d0"/"tsize=3d0"/g oldsad/src/twsdrw_new.f
sed s/"TITLE SIZE -1.6"/"TITLE SIZE 1.6"/g oldsad/src/tdlat.f
　
２、TclTk 形式
　
・やっぱり軸ラベルが文字化けすることがあります。ＯＳとバージョンに
　依って、起きたり起きなかったり。（何とかして欲しいです……）
　
（グラフの書き方）
　
１、topdrawer形式（個人的には普段はこちらを使っています。）
　
　draw "書き出したいパラメータ"　"ラティス先頭の文字"
です。
　draw bx & by q*
だと、βx と βy（のルート）を書き出し、下にラティスの絵を加えて、
Ｑで始まる要素（種類（ＢＱＳ…）に依らず）の名前を表示します。
q* をやめて、draw bx & by だけだと、下にラティスの絵を書きません。
名前表示なしでラティスの絵だけ、はできません。
　draw bx by q*
だと、βx と βy が同じところに上書きされます。

コマンド例
$DisplayFunction=TopDrawer;
　　（これをやっておけば、tdr が見つからなくても tdr 形式のファイルを
　　書き出せます。）
draw bx q*;
　　（画面に tdr 形式のテキストデータがずらずらと並びます。これを個別に
　　ファイルに保存する為に、下のようにします。）
out 15 draw bx & by q* end;
　　（普段の出力（ディスプレイ）は６番です。それを out 15 で装置番号１５番に
　　変更します。これは普通は "fort.15" というファイルになります。そのファイルに
　　向かって "draw bx&by q*" の結果を書き出し、end で出力を６番に戻します。
　　fort.15 にβxとβyのグラフを tdr 形式で書き出す為の１行です。）
FFS["draw bx by & ex b*",12];
　　（個人的に好きなのはこちら。FFS["コマンド"、出力装置番号]ですので、
　　これで "fort.12" に上と同じβxとβyのグラフを書き出せます。）
FFS["draw nx & ny & epx b*",13];
　　（"fort.13"に水平分散関数、垂直分散関数を書き出し。）
FFS["draw geo",14];
　　（これだけは TclTk で同じグラフが書けません。リングのラティスを
　　書き出します。"fort.14"）
　
２、TclTk 形式

　　例えば以下のように書きます。WEB の Tkinter のマニュアルに中身については
詳しく書いてあるので参照してください。（TclTkはややこしいです。）

w = KBMainFrame["sample", namae, Title->"Optics Plot!"];
c1=Canvas[namae,Width->640,Height->400];
$DisplayFunction=CanvasDrawer;
Canvas$Widget=c1;
pl=OpticsPlot[{{"BX"},{"BY"},{"EX"}}];
TkWait[];
　
　　コマンドサンプルと出力されたグラフのサンプルを添付します。
td -d postscr fort.12 などと .ps 形式にして、それを Acrobat distiller で
pdf にしました。tcltk 形式は画面キャプチャです。

!********************************************************************
!*
!* NSLS XRING
!*
!********************************************************************

MOMENTUM=2.5GeV
;
MARK MKN =()
MKL =()
;
MONI PM =()
;
CAVI CAVI=(VOLT= 1.12MV HARM = 30 )
;
BEND B =(L = 2.7 K1 = -0.0724921758 ANGLE = .39269908 E1 = 0.5 E2=0.5)
;
QUAD
Q1 =(L = .45 K1 = -0.675839592 )
Q2 =(L = .8 K1 = 1.069849888 )
Q3 =(L = .4132 K1 = -0.57926284872 )
Q4 =(L = .225 K1 = 0.2923738695 )
;
SEXT SF =(L =0.2 K2 = 0 )
SD =(L =0.2 K2 = 0 )
;
DRIFT D1 =(L = 2.25 )
D2 =(L = 0.685 )
D3 =(L = 0.3484 )
D4 =(L = 0.70825)
D5 =(L = 0.8 )
D7 =(L = 0.7 )
D6 =(L = 0.15 )
;
LINE
NCELLH =(MKN D1 Q1 D2 Q2 D3 Q3 D4
B
D5 SD D7 SF D6 Q4
MKL)
NCELL =(NCELLH -NCELLH)
RING =(4*NCELL CAVI 4*NCELL )
;
!
! to write graph!
!
ON LOG;
ON RAD FLUC COD RFSW RADCOD;

FFS;
$FORM = "15.10";
PageWidth = 1999;

USE RING;
cell calc emit;

$DisplayFunction=TopDrawer;
draw bx q*;
out 15 draw bx & by q* end;
FFS["draw bx by & ex b*",12];
FFS["draw nx & ny & epx b*",13];
FFS["draw geo",14];

w = KBMainFrame["sample", namae, Title->"Optics Plot!"];
c1=Canvas[namae,Width->640,Height->400];
$DisplayFunction=CanvasDrawer;
Canvas$Widget=c1;
pl=OpticsPlot[{{"BX"},{"BY"},{"EX"}}];
TkWait[];
stop;
stop;

　　その５：色収差補正
　
　　線形の色収差が消えるように６極の値を決めます。
具体的には、運動量のずれの２点 +DP と -DP において
チューンが理想的なチューンと等しくなるように６極の
値をフィッティングで決めます。
　
{nx0, ny0}=Twiss[{"NX","NY"},"$$$"]/2/Pi;
　　（設計運動量におけるチューンを求め、水平を nx0、
　　垂直を ny0 という変数に入れます。）
fit nx nx0 2;
　　（fit （場所）（フィットしたい変数）（行き先の値）（点数）
　　で、フィッティングの計算内容を指定します。
　　場所はリングの最後の点ですが、ここでは省略しています。
　　fit $$$ nx nx0 2; または fit MKN.16 nx nx0 2;
　　としても同じです。フィットしたい変数はチューンなので、
　　nx、ny とします。bx by ax ay ex epx なども使えます。
　　行き先の値はこの値に持って行きなさい、という値です。
　　ここでは、設計運動量におけるチューンです。
　　点数は（下で書く） DP で指定する運動量幅 -DP から +DP までの
　　間の何点を取ってフィッティングの計算をするかの指定です。
　　２だと、-DP と +DP の２点、３だと -DP, 0, +DP となります。
　　自由度がなくても数を増やすことはできますが、当然として、
　　結果は完全にはフィットできていない状態になります。なお、
　　６極のファミリー数を増やしてから DP を 0.01 (=1%) くらいに
　　して同じフィッティングをすると、高次の色収差を補正する
　　ことができます。（普通はハーモニック６極は運動量でなくて
　　振幅依存のチューンシフトを補正しますが。））
fit ny ny0 2;
DP=0.0001;
　　（運動量の幅を指定します。デフォルトは 0.01 （=1%)です。
　　MARK 要素に格納されている DP はこの DP です。）
free S*;
　　（フィッティングで何を振るか（何が自由変数か）の指定です。
　　S で始まる要素（種類は問わない）を変化させます。ステアリングを
　　STRH、STRV 等という名前で作ると、このときにステアリングも
　　振られます。勿論、free SF, SD; とすればよいです。また、
　　変数を指定しない限り、その要素の「第１変数」を振ります。
　　free B*; だと ANGLE を変えます。結合Ｂの収束力を変えたい
　　場合は free B* K1; とします。）
go;
　　（フィッティングの計算をします。go で始めて計算を行います。）
FITP 1;
　　（フィッティングの際の運動量方向の点数を１に戻します。これを
　　忘れてチューンやオプティクスのフィッティングをすると、fit の
　　（点数）を書かないと全て（上で指定した）２にされて泣きます。）
type S*;
　　（求まった６極の値を表示します。）
SAVE S*;
　　（６極の値を保存します。これをせずに USE で新しく LINE を
　　展開すると、値が消えます。SAVE だけだと、その時点の全ての
　　要素の「第１変数」を保存します。例えば結合ＢのＫ１を変えて
　　あったとしても、SAVE だけでは保存されません。（angle が
　　保存されます。）　第１変数以外を保存するには、SAVE B* ALL; と
　　しないといけません。）　ちなみに K2 は６極の第１変数なので、
　　SAVE だけか、もしくは SAVE S*; で大丈夫です。）
　
　　以上で色収差補正ができました。ＦＦＳのコマンドは大文字と
小文字の区別がありません。実はその時の気分によって、適当です。
fitp も save も小文字でいいし、FIT も GO も大文字でもいいです。

!********************************************************************
!*
!* NSLS XRING
!*
!********************************************************************

MOMENTUM=2.5GeV
;
MARK MKN =()
MKL =()
;
MONI PM =()
;
CAVI CAVI=(VOLT= 1.12MV HARM = 30 )
;
BEND B =(L = 2.7 K1 = -0.0724921758 ANGLE = .39269908 E1 = 0.5 E2=0.5)
;
QUAD
Q1 =(L = .45 K1 = -0.675839592 )
Q2 =(L = .8 K1 = 1.069849888 )
Q3 =(L = .4132 K1 = -0.57926284872 )
Q4 =(L = .225 K1 = 0.2923738695 )
;
SEXT SF =(L =0.2 K2 = 0 )
SD =(L =0.2 K2 = 0 )
;
DRIFT D1 =(L = 2.25 )
D2 =(L = 0.685 )
D3 =(L = 0.3484 )
D4 =(L = 0.70825)
D5 =(L = 0.8 )
D7 =(L = 0.7 )
D6 =(L = 0.15 )
;
LINE
NCELLH =(MKN D1 Q1 D2 Q2 D3 Q3 D4
B
D5 SD D7 SF D6 Q4
MKL)
NCELL =(NCELLH -NCELLH)
RING =(4*NCELL CAVI 4*NCELL )
;
!
! chromaticity correction
!
ON LOG;
ON RAD FLUC COD RFSW RADCOD;

FFS;
$FORM = "15.10";
PageWidth = 1999;

USE RING;
cell calc;
{nx0, ny0}=Twiss[{"NX","NY"},"$$$"]/2/Pi;

fit nx nx0 2;
fit ny ny0 2;
DP=0.0001;
free S*;
go;
FITP 1;
type S*;
SAVE S*;

tunex = Table[{0,0},{41}];
tuney = Table[{0,0},{41}];

Do[
DP0 = (i - 21) * 0.001;
FFS["cell calc"];
tunex[i] = {DP0,Twiss["NX","$$$"]/2/Pi};
tuney[i] = {DP0,Twiss["NY","$$$"]/2/Pi};
! Print[i,DP0,tunex[i,2],tuney[i,2]];
,{i,1,41,1}];

Print[PolynomialFit[tunex,3]];
Print[PolynomialFit[tuney,3]];

stop;
stop;

　　その６：チューンのマッチングとモジュール化
　
　　チューンをマッチングする。

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

（色収差補正とノーマルセルのチューンのマッチングをモジュール化する。
ChromCorrect[]:=(
Print["Chromaticity Correction"];
FFS["cell calc emit;"];
{nx0, ny0}=Twiss[{"NX","NY"},"***"]/2/Pi;
FFS["FIT;"];
FFS["NX nx0 2;"];
FFS["NY ny0 2;"];
FFS["DP=0.0001;"];
FFS["FREE S*;"];
FFS["GO;"];
FFS["FITP 1;"];
);
　（前回のスクリプトをモジュール化しただけです。）
（モジュールの定義の仕方は以下の通り。
　　モジュール名[受け渡し変数]:=(
　　　　　　　　　内容
　　　　　　　　　);
　　　受け渡し変数がない場合、変数はグローバルとローカルで完全に区別が
　　なくなります。モジュール外の変数をそのまま呼び出せ、値を変えて、
　　ということが可能です。何か受け渡すと、少し複雑になります。
　　　また、";"の意味ですが、関数やモジュールの中か外で変わります。
１、直接 FFS の上で";"がある場合、コマンドや関数を復唱しないという意味に
　なります。
　例えば、
　 a=1
　 a=2;
　と言うのを FFS 上で実行すると、
　 a=1
　Out[4]:= 1
　 a=2;　
　など、";"を付けないと出力を復唱します。実は FFS ではコマンドを打つと
　まずは全て IN[] という関数に引き渡されます。その後、それが変数なのか、
　コマンドなのか、関数なのか判断され、それぞれの実行がなされます。
　IN[] が読みとった内容を画面に書き出す（OUT[]として）かどうかが";"に
　よって決まります。
２、関数やモジュールの中で";"をつけると、それはコマンド同士の行を繋げると
　いう意味になります。実を言えば、関数やモジュール内では全てを繋げないと
　いけません。それが";"です。（どこかの行に";"を付け忘れると、エラーに
　なります。ただし、モジュールの場合、最後の行に値を持ってきて、";"を
　付けないと、それは返値になります。例えば、
ModulePlus[a_,b_]:=(
Print[a, "plus", b];
a+b
);
Print[ModulePlus[1,2]];
　等とすると、"a + b"が返値となり、すなわち、ModulePlus[]を一つの数として
　扱うことができるようになります。（ModulePlus[1,2] + 2 とかできる。）
　　今回のモジュール（色収差補正とチューンのマッチング）では返値はありません。）

MatchingNcellTune[]:=(
Print["Matching Tune to ",tunex,tuney,"!"];
FFS["cell calc emit;"];
FFS["fit MKD EX 0;"];
FFS["fit MKD EPX 0;"];
　　　　（MKD （偏向磁石入り口に付けたマーク要素）で分散を消す。）
FFS["fit NX tunex;"];
FFS["fit NY tuney;"];

! FFS["fit Q2 BX 15"];
! FitValue["Q2","BX",dp_,vgoal_,vnow_]:=If[vnow<5,5,If[vnow>30,30,Null]];
　　　　（２行１組。値を「ある値からある値までの間であればいい」というマッチングを
　　　　する場合、こんな風に書きます。上の行で「フィットに使う」という宣言をして、
　　　　下の行で上限と下限を指定します。上の行の行き先値（bx=15）はこの場合、
　　　　使われません。下の行はほとんど呪文で、
　　　　FitValue["場所","パラメータ",dp_,vgoal_,vnow_]:=If[vnow<下限,下限,If[vnow>上限,上限,Null]];
　　　　と書きます。どうして２回ずつ書くのかは、訊かないでください。^_^;この場合、
　　　　Ｑ２の場所で水平βを５ｍから３０ｍの間で抑えなさい、というコマンドに
　　　　なります。）
! FFS["fit Q4 BX 15"];
! FitValue["Q4","BX",dp_,vgoal_,vnow_]:=If[vnow<5,5,If[vnow>25,25,Null]];
! FFS["fit Q1 BY 10"];
! FitValue["Q1","BY",dp_,vgoal_,vnow_]:=If[vnow<5,5,If[vnow>25,25,Null]];
! FFS["fit Q2 BY 15"];
! FitValue["Q2","BY",dp_,vgoal_,vnow_]:=If[vnow<5,5,If[vnow>30,30,Null]];
　　　　（条件付けるとマッチングがうまくないので、外してしまいました。）
FFS["free Q*;"];
FFS["go;"];
! FFS["SAVE"];
FFS["cell calc"];
);

Do[
Do[
FFS["USE NCELL;"];
　　　（マッチング毎に４極の値はリセット。（SAVE せずに USE をやる。））
MatchingNcellTune[];
　　　（モジュール呼び出して、マッチング。）
ChromCorrect[];
　　　（モジュール呼び出して、色収差補正）
Print[tunex, tuney, Twiss["NX","$$$"]/2/Pi,Twiss["NY","$$$"]/2/Pi,EMITX, LINE["K1","Q1"],LINE["K1","Q2"],LINE["K1","Q3"],LINE["K1","Q4"],LINE["K2","SF"],LINE["K2","SD"]];
　　　（画面に結果を書き出します。）
datafile = OpenAppend["result.txt"];
　　　（ファイルに結果を書き出します。OpenAppend[]は繋げ書き、OpenWrite[]はファイルを
　　　消して新しく作ります。）
Write[datafile,tunex, tuney, Twiss["NX","$$$"]/2/Pi,Twiss["NY","$$$"]/2/Pi,EMITX, LINE["K1","Q1"],LINE["K1","Q2"],LINE["K1","Q3"],LINE["K1","Q4"],LINE["K2","SF"],LINE["K2","SD"]];
Close[datafile];
　　　（いちいち開いて閉じて、をしないと、どうもうまく書けない時があります。本当は
　　　ループの外で開いて、ループが終わったら閉じたいのですが……）
,{tuney, 0.6, 0.91,0.05}];
,{tunex, 1.1, 1.41, 0.05}];
　　　（チューンでループします。）
stop;
stop;

!********************************************************************
!*
!* NSLS XRING
!*
!********************************************************************

MOMENTUM=2.8GeV
;
MARK MKN = ()
MKL = ()
MKD = ()
;
MONI PM =()
;
CAVI CAVI=(VOLT= 1.12MV HARM = 30 )
;
BEND B =(L = 2.7 K1 = -0.0724921758 ANGLE = .39269908 E1 = 0.5 E2=0.5)
;
QUAD
Q1 =(L = .45 K1 = -0.675839592 )
Q2 =(L = .8 K1 = 1.069849888 )
Q3 =(L = .4132 K1 = -0.57926284872 )
Q4 =(L = .225 K1 = 0.2923738695 )
;
QUAD Q1 =(L =.45 K1 =-.7822641438882 )
Q2 =(L =.8 K1 =1.077212593357 )
Q3 =(L =.4132 K1 =-.4741679817474 )
Q4 =(L =.225 K1 =.2923740995656 )
;

SEXT SD =(L =.2 K2 =-5.7352107632318 )
SF =(L =.2 K2 =3.6624554033778 )
;
DRIFT D1 =(L = 2.25 )
D2 =(L = 0.685 )
D3 =(L = 0.3484 )
D4 =(L = 0.70825)
D5 =(L = 0.8 )
D7 =(L = 0.7 )
D6 =(L = 0.15 )
;
LINE
NCELLH =(MKN D1 Q1 D2 Q2 D3 Q3 D4
MKD
B
D5 SD D7 SF D6 Q4
MKL)
NCELL =(NCELLH -NCELLH)
RING =(4*NCELL CAVI 4*NCELL )
;
!
! Matching
!

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

ChromCorrect[]:=(
Print["Chromaticity Correction"];
FFS["cell calc emit;"];
{nx0, ny0}=Twiss[{"NX","NY"},"***"]/2/Pi;
FFS["FIT;"];
FFS["NX nx0 2;"];
FFS["NY ny0 2;"];
FFS["DP=0.0001;"];
FFS["FREE S*;"];
FFS["GO;"];
FFS["FITP 1;"];
);

MatchingNcellTune[]:=(
Print["Matching Tune to ",tunex,tuney,"!"];
FFS["cell calc emit;"];
FFS["fit MKD EX 0;"];
FFS["fit MKD EPX 0;"];
FFS["fit NX tunex;"];
FFS["fit NY tuney;"];
! FFS["fit Q2 BX 15"];
! FitValue["Q2","BX",dp_,vgoal_,vnow_]:=If[vnow<5,5,If[vnow>30,30,Null]];
! FFS["fit Q4 BX 15"];
! FitValue["Q4","BX",dp_,vgoal_,vnow_]:=If[vnow<5,5,If[vnow>25,25,Null]];
! FFS["fit Q1 BY 10"];
! FitValue["Q1","BY",dp_,vgoal_,vnow_]:=If[vnow<5,5,If[vnow>25,25,Null]];
! FFS["fit Q2 BY 15"];
! FitValue["Q2","BY",dp_,vgoal_,vnow_]:=If[vnow<5,5,If[vnow>30,30,Null]];
FFS["free Q*;"];
FFS["go;"];
! FFS["SAVE"];
FFS["cell calc"];
);

Do[
Do[
FFS["USE NCELL;"];
MatchingNcellTune[];
ChromCorrect[];
Print[tunex, tuney, Twiss["NX","$$$"]/2/Pi,Twiss["NY","$$$"]/2/Pi,EMITX, LINE["K1","Q1"],LINE["K1","Q2"],LINE["K1","Q3"],LINE["K1","Q4"],LINE["K2","SF"],LINE["K2","SD"]];
datafile = OpenAppend["result.txt"];
Write[datafile,tunex, tuney, Twiss["NX","$$$"]/2/Pi,Twiss["NY","$$$"]/2/Pi,EMITX, LINE["K1","Q1"],LINE["K1","Q2"],LINE["K1","Q3"],LINE["K1","Q4"],LINE["K2","SF"],LINE["K2","SD"]];
Close[datafile];
,{tuney, 0.6, 0.91,0.05}];
,{tunex, 1.1, 1.41, 0.05}];

stop;
stop;

　　その７：ダイナミックアパーチャサーベイ（方法その１）
　
　　関数　DynamicApertureSurvey[] を用いたダイナミックアパーチャ
サーベイの方法について述べる。
　（トラッキングには RF の DPHI が反映されます。入れたままで計算
すると、加減速されてビームが落ちます。また、Ｘは水平、Ｙは垂直、
Ｚは縦（エネルギー・RF位相方向）というノーテーションです。）
　
（スクリプト）
（初期値（ビームの振幅）は全て「シグマ」単位で指定します。シグマは
　自然エミッタンスを考慮に入れた時のビームサイズですが、分散関数と
　エネルギー広がりの分は考えません。以下が単位換算です。）
w = EMITX;
sigx = Sqrt[Twiss["BX",1] * w];
sigy = Sqrt[Twiss["BY",1] * w];
dp = 0.01 / SIGE;
　（普通は垂直エミッタンスゼロなので、水平エミッタンスにβをかけて
　ルートを取ればビームサイズです。Ｚ方向については、「エネルギー広がり」が
　Ｚ方向１シグマの時のエネルギーのずれ具合なので、0.01 をエネルギー
　広がりで割ると、エネルギーのずれが１％の時に何シグマか、が分かります。
　sigx、sigy は１シグマと言った時の振幅　[m] 、dp はエネルギーのずれが
　１％の時が何シグマに対応するか、です。（考えてみれば、変数名 dp は DP が
　あるので避けるべきですね。））

DynamicApertureSurvey[サーベイ範囲、ターン数、結果出力先] です。
　　サーベイ結果出力先は画面（６番）にし、ターン数は１０００ターンまで
計算します。
　　範囲は、{{xmin, xmax}, {ymin, ymax},{Ｚ方向}}で指定します。{Ｚ方向}の
値１つにつき、１行結果を出します。つまり、-1%, 0%, +1% と指定すれば、
３行結果が出ます。今回はバケットが 1.3% なので、-1.5% から +1.5%
まで、0.5% おきに５点計算（５行表示）します。振幅方向ですが、
(xmax,ymax)、(xmin, ymin) で指定される２点間を５０分割して、両端が
あるので５１点を計算します。必ず５０分割します。そして、各点について
何ターンまわったか（振幅方向に 3m を越えたら死亡と認定）計算します。
結果は100ターン毎に 1, 2, 3,..., 8, 9, A, a, B, b... と表示します。
100ターンはまわった（200ターンまでの間で死亡）なら1、1000ターンまわり
きれば A と表示します。なお、結果を自動の画面表示だけでなくて数字で
表示すると、何ターン目で死んだか具体的に分かります。計算は外側
(xmax, ymax)の点から始めて、順に内側へ進みます。Ｚ方向１点につき、
最大で５１点計算しますが、普通は振幅方向外側から順に何点か連続して
1000 ターンまわれば、それよりも小さい振幅では全部まわるだろうという
推測をします。何点連続でまわりきったらその内側の点の計算を省くか、が
DAPWIDTH です。51 だと推測せずに全部計算、10 にすると、10点連続で
A なら、内側は省いて * にします。外側何シグマの振幅から内側でビームが
生き残ったか（ダイナミックアパーチャ）を「スコア」として表示します。）

DAPWIDTH = 51;

FFS["cell calc emit"];
zrange = Range[-1.5 * dp, 1.6 * dp, dp /2];
range = {{0,100},{1,1},zrange};
resx = DynamicApertureSurvey[range,1000,Output->6];

（出力は以下の通り）
Range Xmin: 0.000 Xmax: 100.000
(Ymin: 1.000 Ymax: 1.000)
　　　　　（(x,y) = (0, 1)と(100, 1)シグマ間を５０分割。ここでは水平
　　　方向の値を求めたいのですが、y = 0 としてしまうと、完全に水平
　　　面内に運動が制限され、非常によくまわる ^_^; 結果が出ます。現実
　　　的にする為に、必ず１シグマくらいずらしておきます。）
Zmin: -16.346 Zmax: 16.346
　　　　　　（DP/P 1% は 10.9 シグマです。-1.5% - +1.5% まで。）
Display: 100 turns/character
　　　　　　　　（水平垂直振幅方向）→
NZ 0----!----1----!----2----!----3----!----4----!----5
　　　　　　　　　　　（←50分割（51点）を計算→）
↓（Ｚ方向）
-16.35 0 333333333333222222000000000000000000000000000000000
-10.90 17 AAAAAAAAAAAAAAAAA0000000000000000000000000000000000
-5.45 16 AAAAAAAAAAAAAAAA00000000000000000000000000000000000
0.00 16 AAAAAAAAAAAAAAAA00000000000000000000000000000000000
5.45 16 AAAAAAAAAAAAAAAA10000000000000000000000000000000000
10.90 17 AAAAAAAAAAAAAAAAA1000000000000000000000000000000000
16.35 7 AAAAAAA8A6AAA55651100000000000000000000000000000000
NZ 0----!----1----!----2----!----3----!----4----!----5
Score: 89

　　（-1.5%（バケットの外） では、小振幅でも 300ターンくらいで
落ちてしまいます。バケット内では、およそ運動量に依らずに 17-16 の
「スコア」になっています。(score - 1) * max / 50 がまわった振幅の
シグマの値です。すなわち、スコアが 17 なら、32 シグマまで "A" で
まわったということです。)

Do[
Print[zrange[i],resx[2,2,i,2],(resx[2,2,i,2] - 1) * sigx * 1000 * 2];
　　　　（一応、運動量（シグマ）、スコア、実振幅 [mm] を書き出します。1000を
　　　　かけるのは mm にしたいから、２は100シグマを50分割したスコアだからです。）
,{i,1,Length[zrange]}];
range = {{1,1},{0,100},zrange};
resy = DynamicApertureSurvey[range,1000,Output->6];
Do[
Print[zrange[i],resy[2,2,i,2],(resy[2,2,i,2] - 1) * sigy * 1000 * 2];
,{i,1,Length[zrange]}];
　
　　LINE の始点がダイナミックアパーチャの始点です。βx 1.1m、βy 0.36m なので
水平７ｍｍ、垂直５ｍｍとか、小さいです。入射点のβはもっと大きいですか？
　
　　なお、お気づきかも知れませんが、これは「ＸＹ面内」でアパーチャサーベイを
するには非常に不便にできています。（それをしたい時は別途スクリプトを書いてDoで
まわします。）

!********************************************************************
!*
!* NSLS Xray-RING
!*
!********************************************************************

MOMENTUM=2.8GeV
;
MARK MKN =()
MKL =()
;
MONI PM =()
;
CAVI CAVI=(L=0.8 VOLT= 1.0MV HARM=30)
;
BEND B =(L = 2.7 ANGLE = .392699082 E1 = 0.5 E2=0.5)
;
QUAD
Q1 =(L = .45 K1 = -0.45*1.4650455999)
Q2 =(L = .8 K1 = 0.8*1.3717947451)
Q3 =(L = .45 K1 = -0.45*1.4109287101)
Q4 =(L = .225 K1 = 0.225*1.381258)
;
SEXT SF =(L =0 K2 = 0 )
SD =(L =0 K2 = 0 )
;
SEXT SD =(L =.2 K2 =-8.8021195157993 )
SF =(L =.2 K2 =5.5762274523103 )
;
DRIFT D1 =(L = 2.25 )
D2 =(L = 0.685 )
D3 =(L = 0.33 )
D4 =(L = 0.69)
D5 =(L = 0.8 )
D7 =(L = 0.7 )
D6 =(L = 0.15 )
D1R =(L = 1.85 )
;
LINE
NCELLH =(MKN D1 Q1 D2 Q2 D3 Q3 D4
B
D5 SD D7 SF D6 Q4
MKL)
NCELLR =(MKN D1R Q1 D2 Q2 D3 Q3 D4
B
D5 SD D7 SF D6 Q4
MKL)
NCELL =(NCELLH -NCELLH)
RING =(6*NCELL NCELLH -NCELLR CAVI NCELLR -NCELLH)
;
!
! dynamic aperture survey command 1
!
ON LOG;
ON RAD FLUC COD RFSW RADCOD;

FFS;
$FORM = "15.10";
PageWidth = 1999;

USE RING;
cell calc emit;
disp $$$;
w = EMITX;
sigx = Sqrt[Twiss["BX",1] * w];
sigy = Sqrt[Twiss["BY",1] * w];
Print[sigx * 1000, sigy * 1000];
SAVE ALL;
DAPWIDTH = 51;
FFS["cell calc emit"];
dp = 0.01 / SIGE;
cell calc;
zrange = Range[-1.5 * dp, 1.6 * dp, dp /2];
! Print[zrange];
range = {{0,100},{1,1},zrange};
resx = DynamicApertureSurvey[range,1000,Output->6];
Do[
Print[zrange[i],resx[2,2,i,2],(resx[2,2,i,2] - 1) * sigx * 1000 * 2];
,{i,1,Length[zrange]}];
range = {{1,1},{0,100},zrange};
resy = DynamicApertureSurvey[range,1000,Output->6];
Do[
Print[zrange[i],resy[2,2,i,2],(resy[2,2,i,2] - 1) * sigy * 1000 * 2];
,{i,1,Length[zrange]}];

stop;
stop;

　　その８：トラッキングの方法（その１）
　
　　これでほぼゴールです。（というか、はじめからこのコマンドの使い方
だけで十分だったという話も……^_^;）　挿入光源の高次の磁場を手で入れる
場合、挿入光源直前まで転送してから、挿入光源部分はスクリプト内で手入力
した数式で転送して、またリング最後まで転送、とすればできます。
　
　　実際に粒子の座標を入れて、トラッキングします。サンプルとして、
４極磁石１個分だけトラックしてみて、転送行列が出てくることを確かめます。
　
TrackParticles[{始点要素番号、転送前座標}、終点要素番号]
で、結果は{終点要素番号、転送後座標}です。結果のリストも
TrackParticles[]にそのまま投げられる形です。
座標は
{{x1,x2,x3...},{x'1,x'2,x'3..},{y1,y2,y3...},{y'1,y'2,y'3...},
{z1,z2,z3..},{dP/P1,dP/P2,dP/P3...},{生死１,生死２,生死３...}}
という７つのリストからなります。粒子を１つだけ転送するなら、
{{x1},{x'1},{y1},{y'1},{z1},{dP/P1},{1}}
となります。粒子の運動がある範囲内（確か、水平垂直 3m 以内？）に
留まっていて安定に生きている場合、生死は１、それを飛び出して
お亡くなりになった場合は生死にゼロが入ります。ゼロが入った
粒子に対しては計算をしません。（転送後も転送前そのままです。）
始点と終点が同じ番号の場合、リングを１周して戻ってきます。また、
始点の方が終点よりも後にある場合も、リング最後まで転送し、
最初に戻って終点まで転送します。粒子１個を１周転送したい場合は、
TrackParticles[{1,{{x1},{x'1},{y1},{y'1},{z1},{dP/P1},{1}}},1]
です。
　
（スクリプト）
rqk = Sqrt[LINE["K1","Q2"]/LINE["L", "Q2"]];
　　（粒子を通す４極として、Ｑ２を選びました。√Ｋを計算しておきます。
　　SAD では K = B'L/(B rho) なので L で割ってルートを取ります。）
ql = LINE["L", "Q2"];
　　（長さです。）
Print[rqk, ql];

qkmath = {{0,0},{0,0}};
qkmatv = {{0,0},{0,0}};
qkmath[1] = {Cos[rqk * ql],Sin[rqk * ql]/rqk};
qkmath[2] = {-rqk * Sin[rqk * ql], Cos[rqk * ql]};
qkmatv[1] = {Cosh[rqk * ql], Sinh[rqk * ql]/rqk};
qkmatv[2] = {rqk * Sinh[rqk * ql], Cosh[rqk * ql]};
　　（言うまでもなく、ごく普通の４極電磁石の転送行列です。）

Print[qkmath[1,1],qkmath[1,2],0,0];
Print[qkmath[2,1],qkmath[2,2],0,0];
Print[0,0,qkmatv[1,1],qkmatv[1,2]];
Print[0,0,qkmatv[2,1],qkmatv[2,2]];

Print[Det[qkmath], Det[qkmatv]];
　　（念のため、行列式が１になることを確認。）

before = {5,{{0.001,0,0,0},{0,0.001,0,0},{0,0,0.001,0},{0,0,0,0.001},{0,0,0,0},{0,0,0,0},{1,1,1,1}}};
　　（LINE の中で Q2 の要素番号は５番です。before は５番の「入り口」での
　　値ということです。初振幅は 1mm、初角度発散も 1mrad とします。（実は
　　これでも振幅依存項が効いてきてしまう値です。）　x,,x',y,y' のどれか
　　１つだけに値を入れて転送（それを４粒子分）することで、転送行列を取り
　　出せます。）
after = TrackParticles[before,6];
　　（６番（要は Q2 の出口まで転送します。）

after[2] = after[2] * 1000;
　　（1mm、1mrad だったので、1000 倍して表示します。）

Print[after[2,1,1],after[2,1,2],after[2,1,3],after[2,1,4]];
Print[after[2,2,1],after[2,2,2],after[2,2,3],after[2,2,4]];
Print[after[2,3,1],after[2,3,2],after[2,3,3],after[2,3,4]];
Print[after[2,4,1],after[2,4,2],after[2,4,3],after[2,4,4]];
　　（これらは振幅依存の高次項の影響を除けば、手計算した転送行列と
　　等しくなっています。）

なお、トラッキングは６次元計算です。理想的にカップリングがない場合、
初座標の指定で水平と垂直は独立のままに保てますが、縦方向は必ず
混じってしまいます。（と言っても桁の問題ですが。）

!********************************************************************
!*
!* NSLS Xray-RING
!*
!********************************************************************

MOMENTUM=2.8GeV
;
MARK MKN =()
MKL =()
;
MONI PM =()
;
CAVI CAVI=(L=0.8 VOLT= 1.0MV HARM=30)
;
BEND B =(L = 2.7 ANGLE = .392699082 E1 = 0.5 E2=0.5)
;
QUAD
Q1 =(L = .45 K1 = -0.45*1.4650455999)
Q2 =(L = .8 K1 = 0.8*1.3717947451)
Q3 =(L = .45 K1 = -0.45*1.4109287101)
Q4 =(L = .225 K1 = 0.225*1.381258)
;
SEXT SF =(L =0 K2 = 0 )
SD =(L =0 K2 = 0 )
;
SEXT SD =(L =.2 K2 =-8.8021195157993 )
SF =(L =.2 K2 =5.5762274523103 )
;
DRIFT D1 =(L = 2.25 )
D2 =(L = 0.685 )
D3 =(L = 0.33 )
D4 =(L = 0.69)
D5 =(L = 0.8 )
D7 =(L = 0.7 )
D6 =(L = 0.15 )
D1R =(L = 1.85 )
;
LINE
NCELLH =(MKN D1 Q1 D2 Q2 D3 Q3 D4
B
D5 SD D7 SF D6 Q4
MKL)
NCELLR =(MKN D1R Q1 D2 Q2 D3 Q3 D4
B
D5 SD D7 SF D6 Q4
MKL)
NCELL =(NCELLH -NCELLH)
RING =(6*NCELL NCELLH -NCELLR CAVI NCELLR -NCELLH)
;
!
! trackparticles sample
!
ON LOG;
OFF ECHO;
ON RAD FLUC COD RFSW RADCOD;

FFS;
$FORM = "20.10";
PageWidth = 1999;

USE NCELL;
cell calc emit;
disp;
SAVE ALL;

rqk = Sqrt[LINE["K1","Q2"]/LINE["L", "Q2"]];
ql = LINE["L", "Q2"];
Print[rqk, ql];

qkmath = {{0,0},{0,0}};
qkmatv = {{0,0},{0,0}};
qkmath[1] = {Cos[rqk * ql],Sin[rqk * ql]/rqk};
qkmath[2] = {-rqk * Sin[rqk * ql], Cos[rqk * ql]};
qkmatv[1] = {Cosh[rqk * ql], Sinh[rqk * ql]/rqk};
qkmatv[2] = {rqk * Sinh[rqk * ql], Cosh[rqk * ql]};

Print[qkmath[1,1],qkmath[1,2],0,0];
Print[qkmath[2,1],qkmath[2,2],0,0];
Print[0,0,qkmatv[1,1],qkmatv[1,2]];
Print[0,0,qkmatv[2,1],qkmatv[2,2]];

Print[Det[qkmath], Det[qkmatv]];

before = {5,{{0.001,0,0,0},{0,0.001,0,0},{0,0,0.001,0},{0,0,0,0.001},{0,0,0,0},{0,0,0,0},{1,1,1,1}}};
after = TrackParticles[before,6];

after[2] = after[2] * 1000;

Print[after[2,1,1],after[2,1,2],after[2,1,3],after[2,1,4]];
Print[after[2,2,1],after[2,2,2],after[2,2,3],after[2,2,4]];
Print[after[2,3,1],after[2,3,2],after[2,3,3],after[2,3,4]];
Print[after[2,4,1],after[2,4,2],after[2,4,3],after[2,4,4]];

stop;
stop;

雑談と冗談交じりに、不正確で適当な情報を書いてしまっただけなので、
間違い、スマートでない部分まもちろん、お怒りを感じる部分や、何を
言ってるか！と思われる部分が多々あると存じます。茶化した部分や
差し出がましい部分など、皆様どうかご容赦下さいますよう、どうぞ
よろしくお願い申し上げます。

前回のようにサーバーのディスクが飛んで掲示板が
蒸発する前に、誰かまとめページ作らない?
#FAQとかチュートリアルとか...

現在、掲示板の定期的なバックアップは
行われているのでしょうか? >>鯖管理者様

前回の掲示板蒸発ではご迷惑をお掛けしました。
今は毎深夜に定期的バックアップを取っています。

A.Morita氏の提案「#FAQとかチュートリアルとか...」は魅力的です。
私は今の所手が回りませんが、どなたか考えて手を動かしてくださいませんか？