use strict;
use Class::Struct;
use File::Find;
use File::Basename;
use File::Spec;
use Math::Trig;
our($mm, $sec, $MHz, $MW, $deg);
our($freq, $omega, $light, $lambda, $eta);
our($init_phase, $phase_shift, $dp_ri);
our($set);
our($input_file, $tw_file, $sw_file, $GPT_file, $para_file);
our($inputRF, $inputDIR);
our(@tw_result, @sw_result);
our(@cell);
our($in_SW, $out_SW);
struct ResutSFO =>{
freq => '$',
U => '$',
Q => '$',
ro_v_Q => '$',
Rs => '$',
Ez_S => '$',
Ez_E => '$'
};
struct TW_Result => {
i => '$',
N => '$',
shape => 'CELL_Shape',
vg => '$',
SS => 'ResutSFO',
gdf_SS => '$',
OO => 'ResutSFO',
gdf_OO => '$'
};
struct CELL_Shape =>{
beta => '$',
a2 => '$',
b2 => '$',
t => '$',
D => '$'
};
struct SW_Shape => {
pipe_L => '$',
a2 => '$',
ra => '$',
b2 => '$',
D_half => '$'
};
struct SW_Result => {
i => '$',
shape => 'SW_Shape',
SS => 'ResutSFO',
gdf => '$',
};
struct SET_GPT => {
i => '$',
j => '$',
KC => '$',
TC => '$'
};
struct CELL => {
zs => '$',
ze => '$',
RF => '$',
phase => '$',
a2 => '$',
D => '$',
t => '$',
vg => '$',
alpha => '$',
Q => '$',
ampSS => '$',
ampOO => '$',
gdf_SS => '$',
gdf_OO => '$'
};
struct SW_CELL => {
zs => '$',
ze => '$',
phase => '$',
pipe_L => '$',
a2 => '$',
amp => '$',
gdf => '$'
};
#------- cell shape parameter -----------
&init($ARGV[0]);
&read_input_data;
&read_tw_data;
&read_sw_data;
&set_cell_para;
&set_sw_para;
&write_cell_para;
&write_gpt_SW_txt($in_SW, "in");
&write_gpt_TW_txt;
&write_gpt_SW_txt($out_SW, "out");
#=============================================================================
# initialize
#=============================================================================
sub init{
my $input = $_[0];
my($s,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst);
my($path, $fname, $ext);
my($abspath, $abspath_file);
fileparse_set_fstype("Unix");
($fname, $path, $ext) = fileparse($input, qr/\.[^\.]+$/);
printf "\n";
printf "input file: %s\n", $input;
printf "file name : %s\n", $fname;
printf "path : %s\n", $path;
printf "ext : %s\n", $ext;
printf "\n";
if($ext ne ".dat"){
printf "\nfile extension error!!! file extension should be \".dat\". \n";
exit(0);
}
chdir $path;
$input_file = $fname.".dat";
$tw_file = $fname.".tw";
$sw_file = $fname.".sw";
$GPT_file = $fname.".forGPT";
$para_file = $fname.".para";
unlink($GPT_file, $para_file);
$mm = 1;
$sec = 1;
$MHz = 1e+6;
$MW = 1e+6;
$deg = pi/180.0;
$freq = 2856*$MHz;
$omega = 2*pi*$freq;
$light = 2.99792458e+8;
$lambda = $light/($freq);
$eta = 0.95; # Real Q / SUPERFISH Q
$init_phase = 0.0; # No.1 cell phase [deg]
$phase_shift =-120.0; # phase shift/cell [deg]
$dp_ri = -90; # phase difference between real and imag
($s,$min,$hour,$mday,$mon,$year,$wday,$yday,$isdst) = localtime(time);
$year += 1900;
$mon += 1;
$abspath_file = File::Spec->rel2abs($fname.".dat");
$abspath = (fileparse($abspath_file))[1];
open(GPT,">".$GPT_file) or die "open: $!";
printf GPT "#=========================================================\n";
printf GPT "# TW Acelerating Structure \n";
printf GPT "# $year.$mon.$mday $hour:$min:$s\n";
printf GPT "# data dir : %s\n", $abspath;
printf GPT "# input data : %s\n", $fname.".dat";
printf GPT "#==========================================================\n\n";
close GPT;
$set = SET_GPT->new();
}
#=============================================================================
# read data
#=============================================================================
sub read_input_data{
my($i);
my $find_dir=0;
my $find_RF_pwr = 0;
open(INPUT,"<".$input_file) or die "open: $!";
while(){
chomp;
if(/^#/){next;} # 先頭に'#'があるとコメント文扱い
if(/#/){s/(^.+)#.*/$1/} # '#'以降はコメント文扱い
if($find_dir or /^\$directory/){ # 先頭に $directory があるまでスキップ
if($find_dir == 0){$find_dir = 1;next;}
}else{
next;
}
if(/^\$end/){last;} # データの終わり
s/^\s*(.*?)\s*$/$1/; # 前後の空白を削除
$inputDIR = $_;
}
while(){
chomp;
if(/^#/){next;} # 先頭に'#'があるとコメント文扱い
if(/#/){s/(^.+)#.*/$1/} # '#'以降はコメント文扱い
if($find_RF_pwr or /^\$input_RF/){ # 先頭に$input_RFがあるまでスキップ
if($find_RF_pwr == 0){$find_RF_pwr = 1;next;}
}else{
next;
}
if(/^\$end/){last;} # データの終わり
s/^\s*(.*?)\s*$/$1/; # 前後の空白を削除
$inputRF = $_ * $MW;
}
close(INPUT);
printf "\t directory : $inputDIR\n";
printf "\t input_RF/MW : %.2f\n",$inputRF*1e-6;
if(not ($find_dir and $find_RF_pwr)){
printf "input data error !!\n";
printf "\t file : $input_file\n";
printf "\t directory : $find_dir\n";
printf "\t input_RF : $find_RF_pwr\n";
exit;
}
}
#=============================================================================
# read tw data
#=============================================================================
sub read_tw_data{
my($i);
my(@in_dat);
$i = 1;
open(TW,"<".$tw_file) or die "open: $!";
while(){
chomp;
if(/^#/){next;} # 先頭に'#'があるとコメント文扱い
s/^\s*(.*?)\s*$/$1/; # 前後の空白を削除
@in_dat = split(/\s+/,$_); # 空白でデータを分割
$tw_result[$i] = TW_Result->new();
$tw_result[$i]->shape(CELL_Shape->new());
$tw_result[$i]->SS(ResutSFO->new());
$tw_result[$i]->OO(ResutSFO->new());
$tw_result[$i]->i($in_dat[0]);
$tw_result[$i]->shape->beta($in_dat[1]);
$tw_result[$i]->shape->a2($in_dat[2]);
$tw_result[$i]->shape->b2($in_dat[3]);
$tw_result[$i]->shape->D($in_dat[4]);
$tw_result[$i]->shape->t($in_dat[5]);
$tw_result[$i]->vg($in_dat[6]);
$tw_result[$i]->SS->freq($in_dat[7]);
$tw_result[$i]->SS->U($in_dat[8]);
$tw_result[$i]->SS->Q($in_dat[9]);
$tw_result[$i]->SS->ro_v_Q($in_dat[10]);
$tw_result[$i]->SS->Rs($in_dat[11]);
$tw_result[$i]->SS->Ez_S($in_dat[12]);
$tw_result[$i]->SS->Ez_E($in_dat[13]);
$tw_result[$i]->gdf_SS($in_dat[14]);
$tw_result[$i]->OO->freq($in_dat[15]);
$tw_result[$i]->OO->U($in_dat[16]);
$tw_result[$i]->OO->Q($in_dat[17]);
$tw_result[$i]->OO->ro_v_Q($in_dat[18]);
$tw_result[$i]->OO->Rs($in_dat[19]);
$tw_result[$i]->OO->Ez_S($in_dat[20]);
$tw_result[$i]->OO->Ez_E($in_dat[21]);
$tw_result[$i]->gdf_OO($in_dat[22]);
$i++;
}
close TW;
$set->KC($i-1);
printf "\tkind of tw cell : %d\n", $set->KC;
}
#=============================================================================
# read sw data
#=============================================================================
sub read_sw_data{
my($i);
my(@in_dat);
$i = 1;
open(SW,"<".$sw_file) or die "open: $!";
while(){
chomp;
if(/^#/){next;} # 先頭に'#'があるとコメント文扱い
s/^\s*(.*?)\s*$/$1/; # 前後の空白を削除
@in_dat = split(/\s+/,$_); # 空白でデータを分割
$sw_result[$i] = SW_Result->new();
$sw_result[$i] -> shape(SW_Shape->new());
$sw_result[$i] -> SS(ResutSFO->new());
$sw_result[$i]->i($in_dat[0]);
$sw_result[$i]->shape->pipe_L($in_dat[1]);
$sw_result[$i]->shape->a2($in_dat[2]);
$sw_result[$i]->shape->ra($in_dat[3]);
$sw_result[$i]->shape->b2($in_dat[4]);
$sw_result[$i]->shape->D_half($in_dat[5]);
$sw_result[$i]->SS->freq($in_dat[6]);
$sw_result[$i]->SS->U($in_dat[7]);
$sw_result[$i]->SS->Q($in_dat[8]);
$sw_result[$i]->SS->ro_v_Q($in_dat[9]);
$sw_result[$i]->SS->Rs($in_dat[10]);
$sw_result[$i]->SS->Ez_S($in_dat[11]);
$sw_result[$i]->SS->Ez_E($in_dat[12]);
$sw_result[$i]->gdf($in_dat[13]);
$i++;
}
close SW_DATA;
if($i != 3){
printf "sw regin read error !!\n";
printf "\tnumber of regions : $i\n";
exit;
}
}
#=============================================================================
# set the cell parameters
#=============================================================================
sub set_cell_para {
my($i);
my($zs, $zc, $ze);
my($a2, $D, $t);
my($RFs, $RFc, $RFe, $phase, $L);
my($vg, $Q, $alpha);
my($U_SS, $U_OO, $U); # Stored Energy
my($sign_SS, $sign_OO); # sign of Ez_E +1/-1
my($ampSS, $ampOO);
my($gdf_SS, $gdf_OO);
$ze = $0;
$RFe = $inputRF;
$phase = $init_phase - $phase_shift;
for($i=1; $i<=$set->KC; $i++){
$cell[$i] = CELL->new();
#-------- set phase ---------
$phase += $phase_shift;
$cell[$i]->phase($phase);
#-------- set z coordinate ---------
$zs = $ze;
$L = $tw_result[$i]->shape->D;
$ze = $zs + $L;
$zc = 0.5*($zs + $ze);
$cell[$i]->zs($zs);
$cell[$i]->ze($ze);
#-------- set shape ---------
$a2 = $tw_result[$i]->shape->a2;
$D = $tw_result[$i]->shape->D;
$t = $tw_result[$i]->shape->t;
$cell[$i]->a2($a2);
$cell[$i]->D($D);
$cell[$i]->t($t);
#-------- set RF parameters ---------
$vg = $tw_result[$i]->vg;
$Q = $tw_result[$i]->SS->Q;
$alpha = $omega/(2.0 * $vg*$light * $Q);
$RFs = $RFe;
$RFe = $RFs*exp(-2.0*$alpha*$L*1e-3/$eta);
$RFc = 0.5*($RFs+$RFe);
$cell[$i]->vg($vg);
$cell[$i]->alpha($alpha);
$cell[$i]->Q($Q);
$cell[$i]->RF($RFc);
#-------- set field multiplication factor ---------
$U_SS = $tw_result[$i]->SS->U; # 1.5 cell Stored energy
$U_OO = $tw_result[$i]->OO->U; # 1.5 cell Stored energy
$U = $RFc/($vg*$light)*1.5*$D*1e-3; # 1.5 cell Stored energy
$sign_SS = $tw_result[$i]->SS->Ez_E <=> 0;
$sign_OO = $tw_result[$i]->OO->Ez_E <=> 0;
$ampSS = $sign_SS*sqrt(0.5*$U/$U_SS); # $U_SS and $U_OO are
$ampOO = $sign_OO*sqrt(0.5*$U/$U_OO); # 1.5 cell Stored energy
$cell[$i]->ampSS($ampSS);
$cell[$i]->ampOO(-$ampOO);
#-------- set field map file ---------
$gdf_SS = $inputDIR.$tw_result[$i]->gdf_SS;
$gdf_OO = $inputDIR.$tw_result[$i]->gdf_OO;
$cell[$i]->gdf_SS($gdf_SS);
$cell[$i]->gdf_OO($gdf_OO);
}
$set->TC($i-1);
printf "\tTotal TW cells : %d\n", $set->TC;
}
#=============================================================================
# set the standing reagion parameter
#=============================================================================
sub set_sw_para {
my($zs, $ze, $L, $phase, $amp, $sign);
my($Ez_SS_S, $Ez_OO_S, $Ez_in_TW1_abs, $Ez_zero_TW1);
my($Ez_SS_E, $Ez_OO_E, $Ez_in_TW2_abs, $Ez_zero_TW2);
my($Ez_in_SW1_abs, $Ez_zero_SW1);
my($Ez_in_SW2_abs, $Ez_zero_SW2);
my($K_final, $T_final);
#----------- input region ----------
$in_SW = SW_CELL->new();
$ze = $cell[1]->zs;
$L = $sw_result[1]->shape->pipe_L + $sw_result[1]->shape->D_half;
$zs = $ze - $L;
$in_SW->zs($zs);
$in_SW->ze($ze);
$in_SW->pipe_L($sw_result[1]->shape->pipe_L);
$in_SW->a2($sw_result[1]->shape->a2);
$in_SW->gdf($inputDIR.$sw_result[1]->gdf);
$Ez_SS_S = ($cell[1]->ampSS)*($tw_result[1]->SS->Ez_S);
$Ez_OO_S = ($cell[1]->ampOO)*($tw_result[1]->OO->Ez_S);
$Ez_in_TW1_abs = sqrt($Ez_SS_S**2 + $Ez_OO_S**2);
$Ez_in_SW1_abs = abs($sw_result[1]->SS->Ez_E);
$phase = $cell[1]->phase;
$Ez_zero_TW1 = $Ez_SS_S*cos($phase)+$Ez_OO_S*cos($phase+$dp_ri);
$Ez_zero_SW1 = $sw_result[1]->SS->Ez_E*cos($phase);
$sign = ($Ez_zero_TW1*$Ez_zero_SW1) <=> 0;
$amp = $sign*$Ez_in_TW1_abs/$Ez_in_SW1_abs;
$in_SW->phase($phase);
$in_SW->amp($amp);
#----------- output region ----------
$out_SW = SW_CELL->new();
$K_final=$set->KC;
$T_final=$set->TC;
$zs = $cell[$set->TC]->ze;
$L = $sw_result[2]->shape->pipe_L + $sw_result[2]->shape->D_half;
$ze = $zs + $L;
$out_SW->zs($zs);
$out_SW->ze($ze);
$out_SW->pipe_L($sw_result[2]->shape->pipe_L);
$out_SW->a2($sw_result[2]->shape->a2);
$out_SW->gdf($inputDIR.$sw_result[2]->gdf);
$Ez_SS_E = ($cell[$T_final]->ampSS)*($tw_result[$K_final]->SS->Ez_E);
$Ez_OO_E = ($cell[$T_final]->ampOO)*($tw_result[$K_final]->OO->Ez_E);
$Ez_in_TW2_abs = sqrt($Ez_SS_E**2 + $Ez_OO_E**2);
$Ez_in_SW2_abs = abs($sw_result[2]->SS->Ez_S);
$phase = $cell[$T_final]->phase;
$Ez_zero_TW2 = $Ez_SS_E*cos($phase)+$Ez_OO_E*cos($phase+$dp_ri);
$Ez_zero_SW2 = $sw_result[2]->SS->Ez_S*cos($phase + $phase_shift);
$sign = ($Ez_zero_TW2*$Ez_zero_SW2) <=> 0;
$amp = $sign*$Ez_in_TW2_abs/$Ez_in_SW2_abs;
$out_SW->phase($phase + $phase_shift);
$out_SW->amp(-$amp);
}
#=============================================================================
# set the cell parameters
#=============================================================================
sub write_cell_para {
my($i);
open(PARA,">".$para_file) or die "open: $!";
printf PARA "#";
printf PARA "zs\t";
printf PARA "ze\t";
printf PARA "RF\t";
printf PARA "phase\t";
printf PARA "a2\t";
printf PARA "D\t";
printf PARA "t\t";
printf PARA "vg\t";
printf PARA "alpha\t";
printf PARA "Q\t";
printf PARA "ampSS\t";
printf PARA "ampOO\t";
printf PARA "gdf_SS\t";
printf PARA "gdf_OO\n";
for($i=1; $i<=$set->TC; $i++){
printf PARA "%.3f\t", $cell[$i]->zs;
printf PARA "%.3f\t", $cell[$i]->ze;
printf PARA "%.3f\t", $cell[$i]->RF;
printf PARA "%.1f\t", $cell[$i]->phase;
printf PARA "%.3f\t", $cell[$i]->a2;
printf PARA "%.3f\t", $cell[$i]->D;
printf PARA "%.3f\t", $cell[$i]->t;
printf PARA "%.3f\t", $cell[$i]->vg*100;
printf PARA "%.3f\t", $cell[$i]->alpha;
printf PARA "%.1f\t", $cell[$i]->Q;
printf PARA "%e\t", $cell[$i]->ampSS;
printf PARA "%e\t", $cell[$i]->ampOO;
printf PARA "%s\t", $cell[$i]->gdf_SS;
printf PARA "%s\n", $cell[$i]->gdf_OO;
}
close(PARA);
}
#=============================================================================
# write TW cell parameters for GPT input
#=============================================================================
sub write_gpt_TW_txt{
my($i);
my($zc, $R, $L);
open(GPT,">>".$GPT_file) or die "open: $!";
for($i=1; $i<=$set->TC; $i++){
printf GPT "#----- cell #%02d -------\n", $i;
#------------- Short-Short field ---------------
printf GPT "map25D_TM(\"wcs\", \"z\", ";
printf GPT "(StrZ%+.3f)*mm, ", $cell[$i]->zs;
printf GPT "\"%s\", ", $cell[$i]->gdf_SS;;
printf GPT "\"R\", \"Z\", \"Er\", \"Ez\", \"H\", ";
printf GPT "%e, ", $cell[$i]->ampSS;
printf GPT "0, ";
printf GPT "(StrP%+.2f)/deg, ",$cell[$i]->phase,;
printf GPT "omega);\n";
#------------- Open-Open field ---------------
printf GPT "map25D_TM(\"wcs\", \"z\", ";
printf GPT "(StrZ%+.3f)*mm, ", $cell[$i]->zs;
printf GPT "\"%s\", ", $cell[$i]->gdf_OO;;
printf GPT "\"R\", \"Z\", \"Er\", \"Ez\", \"H\", ";
printf GPT "%e, ", $cell[$i]->ampOO;
printf GPT "0, ";
printf GPT "(StrP%+.2f)/deg, ",$cell[$i]->phase+$dp_ri,;
printf GPT "omega);\n";
#------------- center disk ---------------
$zc = ($cell[$i]->zs + $cell[$i]->ze)/2.0;
$R = $cell[$i]->a2/2.0;
$L = $cell[$i]->t;
printf GPT "rmax(\"wcs\", \"z\", ";
printf GPT "(StrZ%+.3f)*mm, ", $zc;
printf GPT "%.3f*mm, ", $R;
printf GPT "%.3f*mm", $L;
printf GPT ");\n";
}
close(GPT);
}
#=============================================================================
# write SW cell parameters for GPT input
#=============================================================================
sub write_gpt_SW_txt{
my($SW, $type);
my($zc);
$SW = $_[0];
$type = $_[1];
open(GPT,">>".$GPT_file) or die "open: $!";
printf GPT "#----- SW region -------\n";
printf GPT "map25D_TM(\"wcs\", \"z\", ";
printf GPT "(StrZ%+.3f)*mm, ", $SW->zs;
printf GPT "\"%s\", ", $SW->gdf;;
printf GPT "\"R\", \"Z\", \"Er\", \"Ez\", \"H\", ";
printf GPT "%e, ", $SW->amp;
printf GPT "0, ";
printf GPT "(StrP%+.2f)/deg, ",$SW->phase,;
printf GPT "omega);\n";
#------------- input/output pipe ---------------
if($type eq "in"){
$zc = $SW->zs + $SW->pipe_L/2.0;
}elsif($type eq "out"){
$zc = $SW->ze - $SW->pipe_L/2.0;
}else{
printf "SW region type error !!\n";
printf "Type : %s\n", $type;
exit;
}
printf GPT "rmax(\"wcs\", \"z\", ";
printf GPT "(StrZ%+.3f)*mm, ", $zc;
printf GPT "%.3f*mm, ", ($SW->a2)/2.0;
printf GPT "%.3f*mm", $SW->pipe_L;
printf GPT ");\n";
close(GPT);
}