Yamamoto's Laboratory
はじめに
 
設定
 
入力
  座標変換
  初期粒子分布
 
電磁場取り込み
  SUPERFISH
  電子銃
  進行波型加速管
  定在波型加速管
 
User elemetns
  基本事項
  作成方法
  公開elements
 
結果処理
GDF ファイル
 
Python
  EasyGDF
 
 
 
 
 
 
 
 
研究内容 加速器 GPT User elements 作成方法

GPTUser elements  作成方法シミュレーション機能の追加

GPTのユーザーエレメントの追加の方法を示します．

目次

はじめに

GPTには，電磁石やRF空洞，ビームをコントロールする様々なエレメントが用意されています．単純なケースでは，それらでシミュレーションは可能です．しかし，ちょっと特殊なシステムだと，予め用意されているエレメントでは計算できないことがあります．そのときには，ユーザーがプログラムを記述して，エレメントを追加します．これこそ，GPTの最も優れた特徴で，ユーザーが新たにエレメントを定義できるようになっています．

新規にエレメント作成し，それを実行することは極めて簡単です．ちょっとしたC言語とコンピューターの知識があれば，誰でも，安全にエレメントを追加できます．安全というのは，GPTの他の部分を破壊することなく，ユーザーのコードが実行できるということです．これまで，私は他のトラッキングコードを使ったこともありますが，ユーザーの機能追加という点では，GPTは特に優れていると感じています．

エレメントの追加方法については，付属の「Programmer's Reference」に書かれています．あまり親切に書かれていませんが，必要な情報は得られます．GPTを使う人は，このリファレンスを必ず読んで欲しいと思います．また，このような解析プログラムは「こういう風につくらなくてはならない」というようなこと分かります．GPTはコンピューター科学を勉強した人がつくったプログラムでしょう．

これからは，物理の解析プログラムであろうとも，コンピューター科学の知識は必須と考えています．GPTの「Programmer's Reference」を読むと，そのことを強く感じます．大学で，英語が必須になっているのと同じように，コンピューター科学も必須にすべきと思います．高校で，その学習ができればもっと良いのですが．

話が逸れましたが，ここでは簡単な例を用いて，エレメントを追加する方法を示します．また，「User elements 公開エレメント」に，私の作成したエレメントを公開しています．それも参考になるでしょう．また，連絡していただければ，特殊なエレメントをボランティアーで作成します．

準備

エレメントの動作を定義するプログラムは，C言語で書きます．その実行モジュール作成には，コンパイラーが必要です．Ver 2.8の場合，「GPT News」の 21-May-08 や 18-Oct-07 にしたがい，コンパイラー「Visual C++ 2008 Express Edition」をインストールしてください．さらに，21-May-08 に書かれている設定も必要です．

64ビット環境の場合は，少し事情が異なります．「GPT 設定 64ビットのWindows 7 でコンパイル」に書いていますので，参考にしてください．

エレメントの追加方法

ここでは，ユーザーが定義したエレメントを追加する方法を示します．

追加エレメントの動作とコマンド

進行波空洞

追加するエレメントは，電場の成分が z方向のみの進行波の空洞とします．その半径は r，全長は L とします．電場は z 方向のみで，Ezsin(ωt-kz-φ)とします(下図参照)．

EzTwCavityの様子．半径r，全長Lの軸対称空洞中をz方向のみに電場がある波が，軸に沿って移動する．

端の処理

計算速度の向上のためには，ちょっとばかりの工夫が必要です．空洞の入口や出口で電場がステップ関数で急激に変化しないようにします．マクロ粒子にステップ関数的に電場が印可さると，積分の時間ステップの計算が非常にやっかいになり，計算量の増大を招きます．GPTは計算精度を保証するために，積分のステップ間隔を自動的に調節しています．

ステップ関数で電場が変化する入口と出口で，誤差関数を使いなめらかに電場が変化するようにします．これは，GPTのビルトインエレメント「Trwlinac」などでも使われているテクニックです．実際の電場は，下図のようにします．誤差関数を乗じている領域の 3[mm]は，Trwlinacと同一です．

両端でなめらかに電場が変化するように，端の3[mm]は誤差関数を乗じています(赤線)．

後で述べる検証プログラムでは，次のようになりました．

  • ステップ関数で電場を変化:2分38秒
  • 誤差関数を使って，ステップの部分をなめらかにした場合:2分13秒

思ったほどの効果は無いですが，少し早くなっています．

コマンド

GPTのマニュアルに従うと，この進行波型の空洞のコマンドは，次のように指定することができます．

	trwEzlinac(ECS, Ez,phi, w, k, r, L)

	ECS    Element coodinate system
	Ez     Electric field;Ez. in [V/m].
	phi    Phase factor; φ.in [radian]
	w      Angular frequency:ω.in [s-1].
	k      Longitudinal wavenumber k in [m-1].
	r      Cavity radius in [m].
	L      Cavity length in [m].

エレメントの追加とプログラム

GPT-Element properties の設定

  1. GPTを起動します．
  2. 左下の「Elems」と「Progs」のモード選択は，Elemsでなくてはなりません．Progsの場合には，Elemsをクリックし，モードを変更します．
    3. Elements→Elemsモードに．
  3. メインメニューのElementsのNewを選択します．
    4. Elements→Newを選択．
  4. するとダイアログ「General GPT-Element properties」が現れます．次のように設定します．
    • Name:trwEzlinac
    • Filename:trwEzlinac
    • Description:Travelling wave Ez linac
    • Element typeは，Local fieldsを選択．
    • Parametersは，順にEz, phi, w, k, r, LをAddボタンを使って，追加します．追加する順序がコマンドの引数の並びになりますので，気をつけなくてはなりません．
    この設定が完了したら，「次へ(N)>」ボタンをクリックします．
    5. General GPT-Element properiesの設定．
  5. ダイアログ「Electromagnetic field equations」が現れます．電磁場の設定ですが，個々では，面倒なので，後のプログラム中で設定します．ただし，Length of elementは，Lを選んでおきましょう．
    Electromagnetic field equations の設定．
    設定が終わったら，「完了」ボタンをクリックします．
  6. すると，プログラムのひな形ができあがります．これを，次の「プログラミング」の節で完成させます．
    プログラムのひな形．

プログラミング

先ほど，作成したソースコードのひな形には，次の3つの重要の要素があります．

構造体「trwEzlinac_info」 ユーザー定義関数のデータを保管する構造体です．
関数「trwEzlinac_init」 プログラム実行時，1度に呼び出されます．
関数「trwEzlinac_sim」 計算ステップ毎に呼び出されます．

この構造体と関数の動作を理解できたならば，ひな形のプログラムの修正を行います．下に示すように，不必要なコメント文や実行文を削除し，赤文字の部分を追加すれば，「trwEzlinac(ECS, Ez,phi, w, k, r, L)」(trwEzlinac.c)の完成です．

	  /* trwEzlinac.c: Travelling wave Ez linac */

	  #include <stdio.h>
	  #include <math.h>
	  #include "elem.h"

	  #define SHAPELEN 0.001        /* Shaping over 1 mm */

	  /* Info structure containing all relevant parameters for this element */
	  struct trwEzlinac_info
	  {
	  double Ez ;
	  double phi ;
	  double w ;
	  double k ;
	  double r ;
	  double L ;
	  } ;

	  /* Forward declaration of the routine calculating the electromagnetic fields */
	  static int trwEzlinac_sim(gptpar *par,double t,struct trwEzlinac_info *info) ;

	  /* Initialization routine */
	  void trwEzlinac_init(gptinit *init)
	  {
	  struct trwEzlinac_info *info ;

	  /* Read Element Coordinate System (ECS) from parameter list */
	  gptbuildECS( init ) ;

	  /* Print usage line when the number of parameters is incorrect */
	  if( gptgetargnum(init)!=6 )
	  gpterror( "Syntax: %s(ECS,Ez,phi,w,k,r,L)\n", gptgetname(init) ) ;

	  /* Allocate memory for info structure */
	  info = gptmalloc( sizeof(struct trwEzlinac_info) ) ;

	  /* Read all parameters as doubles and store them in info structure */
	  info->Ez       = gptgetargdouble(init,1) ;
	  info->phi      = gptgetargdouble(init,2) ;
	  info->w        = gptgetargdouble(init,3) ;
	  info->k        = gptgetargdouble(init,4) ;
	  info->r        = gptgetargdouble(init,5) ;
	  info->L        = gptgetargdouble(init,6) ;

	  /* Register the routine calculating the electromagnetic fields to the GPT kernel */
	  gptaddEBelement( init, trwEzlinac_sim, gptfree, GPTELEM_LOCAL, info ) ;
	  }


	  /* The following routine calculates the electromagnetic fields */
	  static int trwEzlinac_sim(gptpar *par,double t,struct trwEzlinac_info *info)
	  {
	  double Ez, phi, w, k, r, L;
	  double shape;               /* smoothing factor for border */
	  double RXY;

	  /* Check element boundaries, return 0 when outside (local) element */
	  if( fabs(Z) > info->L/2+3*SHAPELEN ) return( 0 ) ; 
	  if( X*X+Y*Y > gptSQR(info->r+3*SHAPELEN) ) return( 1 ) ;
	  
	  Ez       = info->Ez ;
	  phi      = info->phi ;
	  w        = info->w ;
	  k        = info->k ;
	  r        = info->r ;
	  L        = info->L ;

	  

	    RXY = sqrt(X*X+Y*Y);
	    shape = erfshape((Z+(L/2.0+3*SHAPELEN))/SHAPELEN) 
            * erfshape(((L/2.0+3*SHAPELEN)-Z)/SHAPELEN) 
            * erfshape(((r+3*SHAPELEN)-RXY)/SHAPELEN); 
	    EZ = shape*Ez*cos(w*t-k*Z+phi);
	  


	  /* Return 1 to notify particle is INSIDE element */
	  return( 1 ) ;
	  }

鋭い読者ならば，関数「trwEzlinac_sim」のr方向の領域の選択の文，

if( X*X+Y*Y > gptSQR(info->r+3*SHAPELEN) ) return( 1 ) ;

の「return(1)」は，「return(0)」とすべきと考えるかもしれません．こうすると途中で「gpt: Step size too small at T=?????」というようなエラーメッセージを出し，計算が終了します．あるいは，長時間計算が終了しないことがあります．理由は分かりませんが，「return(1)」とすれば問題なく動作します．

プログラムが完成したら，メインメニューのFileのSaveを選択し，ソースコードを保存します．

ソースコードの保存．

コンパイル

  1. メインメニューのElementsのCompile GPT-Elemsを選択し，ソースコードをコンパイルします．
    2. ソースコードのコンパイル．

追加したエレメントの使用方法と確認

使用方法

使用方法は簡単です．以下に示すGPTのインプットファイルを参考にしてください．これを，バッチファイル(trwEzlinac_test.in)を使って，実行させれば図のような結果が得られます．

Energy=10e+6;         # beam energy [MeV]
G=1-qe*Energy/(me*c^2);
setparticles("beam",10000,me,qe,-1e-9);
setrxydist("beam","u", 0.03, 0.06);
setphidist("beam","u", 0, 2*pi);
setzdist("beam","u", 0.3, 0.6);
setGBrxydist("beam","u",0, 0);
setGBphidist("beam","u",0, 0);
setGdist("beam","u",G, 0);

trwEzlinac("wcs", "z", 1.0, 10.0e6 , pi/2, 2*pi*2856e6, 2*pi/0.105, 0.05, 0.4);

tout(0,5e-9,0.02e-9);
電場Ezの強弱を色(赤→緑&arr;青)で表しています．緑色はEz=0です．また，アニメーション(trwEzlinac.mpg)もあります．

動作の確認

進行波になっていることは，先ほどのアニメーションで分かるでしょう．また，空洞の領域もきちんと設定されていることも理解できます．電場の強度は，次のようなプロットを描くと確認ができます．この図から，追加したエレメントの動作は問題ないことが分かります．

z座標と電場の関係．わかりにくいですが，色はマクロ粒子の半径方向の位置に依存しています．下の方に，最大電場と最小電場が書かれています．

ページ作成情報

参考資料

  1. Gneral Particle Tracer — Programmer's Reference Version 2.80 —

更新履歴

2010年頃 ページの新規作成

no counter

last update:2013/11/16 15:09:38