use strict;
use Class::Struct;
use File::Find;
use File::Basename;
use Math::Trig;

our($file);
our($SFPara);
our($mm, $sec, $MHz);
our($freq, $light, $lambda);
our($it_max, $error);
our($z_min, $z_max);
our($cal);
our(@input, @tw_result);
our($seg_Wall_L, $seg_Wall_R);
our($U_SS, @rL_SS, @rR_SS, @HL_SS, @HR_SS, $N_divR);
our($U_OO, @rL_OO, @rR_OO, @EL_OO, @ER_OO);
our($inSW, $outSW, @sw_result);

struct File => {
    input => '$',
    base  => '$',
    af    => '$',
    seg   => '$',
    T35   => '$',
    in7   => '$',
    SFO   => '$',
    tw    => '$',
    GPT   => '$',
    sf7   => '$',
    gdf   => '$',
    sw    => '$'
};

struct forSF => {
    title => '$',
    conv  => '$',
    dx    => '$'
};

struct INPUT => {
    shape => 'CELL_Shape',
};

struct ResutSFO =>{
    freq   => '$',
    U      => '$',
    Q      => '$',
    ro_v_Q => '$',
    Rs     => '$',
    Ez_S   => '$',
    Ez_E   => '$'
};

struct Result => {
    i      => '$',
    N      => '$',
    shape  => 'CELL_Shape',
    vg     => '$',
    SS     => 'ResutSFO',
    gdf_SS => '$',
    OO     => 'ResutSFO',
    gdf_OO => '$'
};

struct CELL_Shape =>{
    beta      => '$',
    a2        => '$',
    b2        => '$',
    t         => '$',
    D         => '$'
};

struct CalMode => {
    i    => '$',
    mode => '$',
    N    => '$'
};

struct SW_Region => {
    pipe_L => '$',
    a2     => '$',
    ra     => '$',
    b2     => '$',
    D_half => '$'
};

struct SW_Result => {
    i     => '$',
    shape => 'SW_Region',
    SS    => 'ResutSFO',
    gdf   => '$',
};

print "SUPERFISH has been started.\n";


#------- cell shape parameter -----------

&init($ARGV[0]);
&read_cell_data;       # read cell data form input cell geometory file

for($cal->i(1); $cal->i <= $cal->N; $cal->i($cal->i+1)){
    printf "\n------------ Cell %d ---------------", $cal->i;
    &tune_b_acc;           # Tune b which is acc structure
    &cal_SF_two;           # Calculate two modes (SO, OS) by SF
    &cal_vg;
    &write_result;
}

&mk_seg(2, 5);
&init_SW_region;
&mk_sw_region;
&rm_files;

#=============================================================================
#        initialize
#=============================================================================
sub init{
    my $input = $_[0];
    my($path, $fname, $ext);

    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;

    $file = File->new();
    
    $file -> input($fname.$ext);
    $file -> base($fname);
    $file -> af($fname.".af");
    $file -> seg($fname.".seg");
    $file -> T35($fname.".T35");
    $file -> in7($fname.".in7");
    $file -> SFO($fname.".SFO");
    $file -> tw( $fname.".tw");
    $file -> GPT($fname.".gpt");
    $file -> sf7("OUTSF7.TXT");
    $file -> gdf($fname.".gdf");
    $file -> sw( $fname.".sw");

    $SFPara = forSF->new();
    $SFPara -> title("S-Band Accelerator Structure");
    $SFPara -> conv(0.1);
    $SFPara -> dx(0.2);

    $cal = CalMode->new();
    $cal->i(1);

    $mm=1;
    $sec=1;
    $MHz=1e+6;

    $freq=2856*$MHz;
    $light=2.99792458e+8;
    $lambda=$light/($freq);

    $it_max=10;                   # 2b調整時の最大イタレーション回数
    $error=1e-6;                  # 2b調整時の周波数エラーの許容値
    $seg_Wall_L = 1;              # 左端の壁のセグメント番号
    $seg_Wall_R = 10;             # 右
    $N_divR = 1000;               # ポインティングベクトル計算の壁の分割数


    unlink($file->af, $file->T35, $file -> in7, $file->SFO);
    unlink($file->tw, $file->sw);

    opendir(DIR,"./") or die "could not open directory: $!";
    while ($fname = readdir(DIR)){
	if($fname =~ /\.gdf$/){unlink($fname);next;}
	if($fname =~ /\.TXT$/){unlink($fname);next;}
	if($fname =~ /\.TBL$/){unlink($fname);next;}
	if($fname =~ /\.log$/){unlink($fname);next;}
	if($fname =~ /\.INF$/){unlink($fname);next;}
	if($fname =~ /\.PMI$/){unlink($fname);next;}
    }

    open(CELL_DATA,">".$file->tw) or die "open: $!";
    printf CELL_DATA "#i\t";
    printf CELL_DATA "beta\t";
    printf CELL_DATA "2a\t";
    printf CELL_DATA "2b\t";
    printf CELL_DATA "D\t";
    printf CELL_DATA "t\t";
    printf CELL_DATA "vg\t";
    printf CELL_DATA "f(SS)\t";
    printf CELL_DATA "U\t";
    printf CELL_DATA "Q\t";
    printf CELL_DATA "r/Q\t";
    printf CELL_DATA "Rs\t";
    printf CELL_DATA "Ez_S\t";
    printf CELL_DATA "Ez_E\t";
    printf CELL_DATA "gdf\t";
    printf CELL_DATA "f(OO)\t";
    printf CELL_DATA "U\t";
    printf CELL_DATA "Q\t";
    printf CELL_DATA "r/Q\t";
    printf CELL_DATA "Rs\t";
    printf CELL_DATA "Ez_S\t";
    printf CELL_DATA "Ez_E\t";
    printf CELL_DATA "gdf\n";

    close CELL_DATA;

    open(SW_DATA,">".$file->sw) or die "open: $!";
    printf SW_DATA "#i\t";
    printf SW_DATA "pipe_L\t";
    printf SW_DATA "2a\t";
    printf SW_DATA "ra\t";
    printf SW_DATA "2b\t";
    printf SW_DATA "Dh\t";
    printf SW_DATA "f(SS)\t";
    printf SW_DATA "U\t";
    printf SW_DATA "Q\t";
    printf SW_DATA "r/Q\t";
    printf SW_DATA "Rs\t";
    printf SW_DATA "Ez_S\t";
    printf SW_DATA "Ez_E\t";
    printf SW_DATA "gdf\n";

    close SW_DATA;
}


#=============================================================================
#        make seg_file for calculating Q factor
#=============================================================================
sub mk_seg{
    my($i);
    my($seg_start, $seg_end);

    $seg_start = $_[0];
    $seg_end   = $_[1];

    open(SEG,">".$file->seg) or die "open: $!";

    printf SEG "FieldSegments\n";
    for($i=$seg_start; $i<$seg_end; $i++){
	printf SEG "%d ",$i;
    }
    printf SEG "%d\n",$i;
    printf SEG "EndData\n";
    printf SEG "End\n";

    close(SEG);

}

#=============================================================================
#        make seg_file for calculating vg
#=============================================================================
sub mk_seg2{
    my($segL, $segR);

    $segL = $_[0];
    $segR = $_[1];

    open(SEG,">".$file->seg) or die "open: $!";

    printf SEG "FieldSegments\n";
    printf SEG "%d %d\n", $segL, $segR;
    printf SEG "EndData\n";
    printf SEG "End\n";

    close(SEG);

}

#=============================================================================
#        read data
#=============================================================================
sub read_cell_data{
    my($i);
    my(@in_dat);
    my $find = 0;

    $i=1;

    open(CELL_DATA,"<".$file->input) or die "open: $!";

    while(<CELL_DATA>){
	chomp;
	if(/^#/){next;}               # 先頭に'#'があるとコメント文扱い
	if(/#/){s/(^.+)#.*/$1/}       # '#'以降はコメント文扱い
	if($find or /^\$cell/){       # 先頭に$cellがあるまでスキップ
	    if($find == 0){$find = 1;next;}
	}else{
	    next;
	}
	if(/^\$end/){last;}           # データの終わり
	s/^\s*(.*?)\s*$/$1/;          # 前後の空白を削除

	printf "%d\t%s\n",$i, $_;

	@in_dat = split(/\s+/,$_);

	@input[$i] = INPUT->new();
	@input[$i]->shape(CELL_Shape->new());

	$input[$i] -> shape -> beta($in_dat[0]);
	$input[$i] -> shape -> a2(  $in_dat[1]);
	$input[$i] -> shape -> b2(  $in_dat[2]);
	$input[$i] -> shape -> t(   $in_dat[3]);
	$input[$i] -> shape -> D(   $in_dat[0]*$lambda/3.0*1e3);

	$i++
    }

    close(CELL_DATA);

    $cal->N($i-1);

    &check_data;

}

#=============================================================================
#        tunning the accelerating cavity radius
#=============================================================================
sub tune_b_acc{
    my(@tune_b_acc, @sf_frequency, $i, $j);

    $i = $cal->i;

    printf "\n";
    $tune_b_acc[1]=$input[$i]->shape->b2;    # 初期値
    for($j=1; $j<=$it_max; $j++){
	write_sf_data(0, 0);
	system "autofish ".$file->af;
	$sf_frequency[$j]=&read_fr;
	printf "\t%d\t%f\t%f\n",$j,$tune_b_acc[$j],$sf_frequency[$j]/$MHz;
	if(abs(($sf_frequency[$j]-$freq)/$freq)<$error){last};
	
	
	if($j==1){
	    $tune_b_acc[2]=$sf_frequency[1]*$tune_b_acc[1]/$freq;
	}else{
	    $tune_b_acc[$j+1]=
		($tune_b_acc[$j]-$tune_b_acc[$j-1])
		/($sf_frequency[$j]-$sf_frequency[$j-1])
		*($freq-$sf_frequency[$j])
		+$tune_b_acc[$j];
	}
	    $input[$i]->shape->b2($tune_b_acc[$j+1]);
    }    
}

#=============================================================================
#        calculate OS and SO modes
#=============================================================================
sub cal_SF_two{
    my($i);
    my($freq, $U, $Q, $r, $r_ov_Q, $Rs);
    my($f_error);
    my($Ez_S, $Ez_E);
    my($gdf_file);

    &mk_seg($seg_Wall_L+1, $seg_Wall_R-1);

    $i = $cal->i;
    $tw_result[$i]=Result->new();

    $tw_result[$i]->i($i);

    $tw_result[$i]->shape(CELL_Shape->new());
    $tw_result[$i]->shape($input[$i]->shape);

    #--------------- calculate SS mode ------------------
    write_sf_data(1, 1);
    system "autofish ".$file->af;
    printf "\t ---- Short - Short mode ----\n";
    ($freq, $U, $Q, $r_ov_Q, $Rs)=&read_SFO;
    ($Ez_S, $Ez_E) = &side_Ez(0.0, $tw_result[$i]->shape->D);

    printf "\t Frequency    \t%.5f\n", $freq/$MHz;
    printf "\t Stored Energy\t%e\n", $U;
    printf "\t Q value      \t%.2f\n", $Q;
    printf "\t r/Q          \t%e\n", $r_ov_Q;
    printf "\t Rs           \t%e\n", $Rs;
    printf "\t Ez(start)    \t%e\n", $Ez_S;
    printf "\t Ez(end)      \t%e\n", $Ez_E;

    $tw_result[$i]->SS(ResutSFO->new());
    $tw_result[$i]->SS->freq($freq);
    $tw_result[$i]->SS->U($U);
    $tw_result[$i]->SS->Q($Q);
    $tw_result[$i]->SS->ro_v_Q($r_ov_Q);
    $tw_result[$i]->SS->Rs($Rs);
    $tw_result[$i]->SS->Ez_S($Ez_S);
    $tw_result[$i]->SS->Ez_E($Ez_E);
    $gdf_file = &mk_gdf_file_name("SS");
    $tw_result[$i]->gdf_SS($gdf_file);
    &make_gdf($gdf_file);             # make a gdf file

    #--------------- calculate OO mode ------------------ 
    write_sf_data(0, 0);
    system "autofish ".$file->af;
    printf "\t ---- Open - Open mode ----\n";
    ($freq, $U, $Q, $r_ov_Q, $Rs)=&read_SFO;
    ($Ez_S, $Ez_E)=&side_Ez(0.0, $tw_result[$i]->shape->D);

    printf "\t Frequency\t%.5f\n", $freq/$MHz;
    printf "\t Stored Energy\t%e\n", $U;
    printf "\t Q value\t%.2f\n", $Q;
    printf "\t r/Q\t%e\n", $r_ov_Q;
    printf "\t Rs\t%e\n", $Rs;
    printf "\t Ez(start)    \t%e\n", $Ez_S;
    printf "\t Ez(end)      \t%e\n", $Ez_E;

    $tw_result[$i]->OO(ResutSFO->new());
    $tw_result[$i]->OO->freq($freq);
    $tw_result[$i]->OO->U($U);
    $tw_result[$i]->OO->Q($Q);
    $tw_result[$i]->OO->ro_v_Q($r_ov_Q);
    $tw_result[$i]->OO->Rs($Rs);
    $tw_result[$i]->OO->Ez_S($Ez_S);
    $tw_result[$i]->OO->Ez_E($Ez_E);
    $gdf_file = &mk_gdf_file_name("OO");
    $tw_result[$i]->gdf_OO($gdf_file);
    &make_gdf($gdf_file);             # make a gdf file

    #------------- frequency check ------------
    $f_error = abs($tw_result[$i]->OO->freq - $tw_result[$i]->SS->freq);
    if(100*$error< $f_error/$tw_result[$i]->OO->freq){
	printf "Frequency error !!\n";
	printf "\tSS mode : %f\n",$tw_result[$i]->SS->freq/$MHz;
	printf "\tOO mode : %f\n",$tw_result[$i]->OO->freq/$MHz;
	exit;
    }

}

#=============================================================================
#        calculate group velocity
#=============================================================================
sub cal_vg{
    my($i);
    my($freq_SS, $Q_SS, $r_SS, $r_ov_Q_SS, $Rs_SS);
    my($freq_OO, $Q_OO, $r_OO, $r_ov_Q_OO, $Rs_OO);
    my($f_err, $U_err);                              # error
    my($a, $b, $Lz);
    my($U_tw, $RF_pwr_L, $RF_pwr_R);            # Stored Energy and RF power
    my($vg_L, $vg_R);

    $i = $cal->i;

    $a = $tw_result[$i]->shape->a2/2.0;
    $b = $tw_result[$i]->shape->b2/2.0;
    $Lz =$tw_result[$i]->shape->D*1.5;

    #--------------- calculate SS mode ------------------
    write_sf_data(1, 1);
    system "autofish ".$file->af;
    printf "\t ---- Short - Short mode for vg cal ----\n";
    ($freq_SS, $U_SS, $Q_SS, $r_ov_Q_SS, $Rs_SS)=&read_SFO;
    printf "\t Frequency\t%f.5\n", $freq_SS/$MHz;
    # --- mode check ---------
    $f_err = abs($freq_SS - $tw_result[$i]->SS->freq);
    $U_err = abs($U_SS - $tw_result[$i]->SS->U);
    if(1e-6 < $f_err/$freq_SS || 1e-6 < $U_err/$U_SS){
	printf "mode check error at cal_vg  SS mode !!\n";
	exit(0);
    }

    &set_field("Ht", 0.0, 0.0, $b, \@rL_SS, \@HL_SS);
    &set_field("Ht", $Lz, 0.0, $a, \@rR_SS, \@HR_SS);

    #--------------- calculate OO mode ------------------ 
    write_sf_data(0, 0);
    system "autofish ".$file->af;
    printf "\t ---- Open - Open mode for vg cal ----\n";
    ($freq_OO, $U_OO, $Q_OO, $r_ov_Q_OO, $Rs_OO)=&read_SFO;
    printf "\t Frequency\t%f.5\n", $freq_OO/$MHz;
    # --- mode check ---------
    $f_err = abs($freq_OO - $tw_result[$i]->OO->freq);
    $U_err = abs($U_OO - $tw_result[$i]->OO->U);
    if(1e-6 < $f_err/$freq_OO || 1e-6 < $U_err/$U_OO){
	printf "mode check error at cal_vg  OO mode !!\n";
	exit(0);
    }

    &set_field("Er", 0.0, 0.0, $b, \@rL_OO, \@EL_OO);
    &set_field("Er", $Lz, 0.0, $a, \@rR_OO, \@ER_OO);

    #------------- frequency_check ------------
    $f_err = abs($freq_SS - $freq_OO);
    if(100*$error< $f_err/$freq_SS){
	printf "Frequency error at cal_vg !!\n";
	printf "\tSS mode : %f\n",$freq_SS/$MHz;
	printf "\tOO mode : %f\n",$freq_OO/$MHz;
	exit;
    }

    ($U_tw, $RF_pwr_L, $RF_pwr_R) = &cal_RF_pwr;
    printf("\tStoroed Energy in 1.5 cell : %e [Joules]\n", $U_tw);

    $vg_L = $Lz*1e-3*$RF_pwr_L/$U_tw;
    printf("\tRF power flow(left wall)   : %e [Watt]\n",   $RF_pwr_L);
    printf("\tGrout velocity(left wall)  : %e [m/sec]\n",  $vg_L);
    printf("\t\tvg/c(left wall)          : %e\n",          $vg_L/$light);

    $vg_R = $Lz*1e-3*$RF_pwr_R/$U_tw;
    printf("\tRF power flow(right wall)   : %e [Watt]\n",   $RF_pwr_R);
    printf("\tGrout velocity(right wall)  : %e [m/sec]\n",  $vg_R);
    printf("\t\tvg/c(rith wall)          : %e\n",          $vg_R/$light);

    $tw_result[$i]->vg($vg_L/$light);

}



#=============================================================================
#        check data
#=============================================================================
sub check_data{
    my($i);

    printf "\n=============== input data check ==============\n";

    for($i=1; $i<=$cal->N; $i++){
	printf "Cell Number %d\n", $i;
	printf "beta  \t%f\n",$input[$i]->shape->beta;
	printf "2a    \t%f\n",$input[$i]->shape->a2;
	printf "2b    \t%f\n",$input[$i]->shape->b2;
	printf "D     \t%f\n",$input[$i]->shape->D;
	printf "t     \t%f\n",$input[$i]->shape->t;
	printf "\n";
    }

}


#=============================================================================
#        write result
#=============================================================================
sub write_result{
    my($i);
    $i = $cal->i;

    open(TW_DATA,">>".$file->tw) or die "open: $!";

    printf TW_DATA "%d\t",   $tw_result[$i]->i;
    printf TW_DATA "%.2f\t", $tw_result[$i]->shape->beta;
    printf TW_DATA "%.3f\t", $tw_result[$i]->shape->a2;
    printf TW_DATA "%.3f\t", $tw_result[$i]->shape->b2;
    printf TW_DATA "%.3f\t", $tw_result[$i]->shape->D;
    printf TW_DATA "%.3f\t", $tw_result[$i]->shape->t;
    printf TW_DATA "%e\t",   $tw_result[$i]->vg;
    printf TW_DATA "%.5f\t", $tw_result[$i]->SS->freq/$MHz;
    printf TW_DATA "%e\t",   $tw_result[$i]->SS->U;
    printf TW_DATA "%.2f\t", $tw_result[$i]->SS->Q;
    printf TW_DATA "%e\t",   $tw_result[$i]->SS->ro_v_Q;
    printf TW_DATA "%e\t",   $tw_result[$i]->SS->Rs;
    printf TW_DATA "%e\t",   $tw_result[$i]->SS->Ez_S;
    printf TW_DATA "%e\t",   $tw_result[$i]->SS->Ez_E;
    printf TW_DATA "%s\t",   $tw_result[$i]->gdf_SS;
    printf TW_DATA "%.5f\t", $tw_result[$i]->OO->freq/$MHz;
    printf TW_DATA "%e\t",   $tw_result[$i]->OO->U;
    printf TW_DATA "%.2f\t", $tw_result[$i]->OO->Q;
    printf TW_DATA "%e\t",   $tw_result[$i]->OO->ro_v_Q;
    printf TW_DATA "%e\t",   $tw_result[$i]->OO->Rs;
    printf TW_DATA "%e\t",   $tw_result[$i]->OO->Ez_S;
    printf TW_DATA "%e\t",   $tw_result[$i]->OO->Ez_E;
    printf TW_DATA "%s\n",   $tw_result[$i]->gdf_OO;

    close TW_DATA;
    
}

#=============================================================================
#        write superfish data
#=============================================================================
sub write_sf_data {
    my($nbslf, $nbsrt);
    my($beta, $a, $b, $D, $t);
    my($xdri, $ydri);
    my($i);

    $i = $cal->i;

    $nbslf = $_[0];
    $nbsrt = $_[1];
    $beta  = $input[$i]->shape->beta;
    $a     = ($input[$i]->shape->a2)/2.0;
    $b     = ($input[$i]->shape->b2)/2.0;
    $D     = $input[$i]->shape->D;
    $t     = $input[$i]->shape->t;
    $xdri  = 1.0*$D;
    $ydri  = 0.7*$b;

    open(SF_DATA,">".$file->af) or die "open: $!";

    printf SF_DATA "%s\n", $SFPara -> title;
    printf SF_DATA " \$reg kprob=1, ";
    printf SF_DATA "dx=%f ,", $SFPara->dx;
    printf SF_DATA "conv=%f ,", $SFPara->conv; 
    printf SF_DATA "kprob=1,\n";
    printf SF_DATA " xdri=$xdri, ydri=$ydri,\n";
    printf SF_DATA " nbsup=1, nbslo=0, nbslf=$nbslf, nbsrt=$nbsrt\n";
    printf SF_DATA " freq=%f,\n",$freq/$MHz;
    printf SF_DATA " kmethod=1,\n";
    printf SF_DATA " beta=%f\$\n",$beta;

    # ------ accelerating cell (half half)  -----------
    printf SF_DATA "!---------- half half cell -----------\n";
    printf SF_DATA " \$po x=%f, y=%f \$\n", 0.0, 0.0;
    printf SF_DATA " \$po x=%f, y=%f \$\n", 0.0, $b;
    printf SF_DATA " \$po x=%f, y=%f \$\n", 0.5*$D-0.5*$t, $b;
    printf SF_DATA " \$po x=%f, y=%f \$\n", 0.5*$D-0.5*$t, $a+0.5*$t;
    printf SF_DATA
	" \$po nt=2, r=%f, theta=270, x0=%f, y0=%f \$\n",
	0.5*$t, 0.5*$D, $a+0.5*$t;

    # ------ accelerating cell (full cell)  -----------
    printf SF_DATA "!---------- full cell -----------\n";
    printf SF_DATA
	" \$po nt=2, r=%f, theta=0, x0=%f, y0=%f \$\n",
	0.5*$t, 0.5*$D, $a+0.5*$t;
    printf SF_DATA " \$po x=%f, y=%f \$\n", 0.5*$D+0.5*$t, $b;
    printf SF_DATA " \$po x=%f, y=%f \$\n", 1.5*$D-0.5*$t, $b;
    printf SF_DATA " \$po x=%f, y=%f \$\n", 1.5*$D-0.5*$t, $a+0.5*$t;
    printf SF_DATA
	" \$po nt=2, r=%f, theta=270, x0=%f, y0=%f \$\n",
	0.5*$t, 1.5*$D, $a+0.5*$t;

    printf SF_DATA " \$po x=%f, y=%f \$\n", 1.5*$D, 0.0;
    printf SF_DATA " \$po x=%f, y=%f \$\n", 0.0, 0.0;

    close SF_DATA;

}

#=============================================================================
#        READ frequency (***.SFO)
#=============================================================================
sub read_fr{

    my($resonance_fr);
    my(@oneline);

    open(SFOUT,"<".$file->SFO) or die "open:$file  $!";;
    while(<SFOUT>){
	chomp;
	
	if(/^Frequency\s+=/){
		@oneline=split(/\s+/,$_);
		$resonance_fr=@oneline[2];
		next;
	}
    }

    close SFOUT;
    return $resonance_fr*$MHz;
}

#=============================================================================
#        READ frequency and Rs etc (***.SFO)
#=============================================================================
sub read_SFO{

    my($f, $U, $Q, $r, $r_ov_Q);
    my(@oneline);

    $f      = 0;
    $U      = 0;
    $Q      = 0;
    $r      = 0;
    $r_ov_Q = 0;
    
    open(sfout,"<".$file->SFO) or die "open:$file  $!";

    while(<sfout>){
	chomp;
	
	if(/^All calculated values below refer to the mesh geometry only./){
	    last;
	}
    }

    while(<sfout>){
	chomp;
	
	if(/^Frequency\s+=/){
		@oneline=split(/\s+/,$_);
		$f=$oneline[2]*$MHz;
	}

	if(/^Stored\s+energy/){
		@oneline=split(/\s+/,$_);
		$U=$oneline[3];
	}


	if(/^Q\s+=/){
		@oneline=split(/\s+/,$_);
		$Q=$oneline[2];
	}

	if(/Z\*T\*T/){
		@oneline=split(/\s+/,$_);
		$r=$oneline[6];
	}

	if(/r\/Q/){
		@oneline=split(/\s+/,$_);
		$r_ov_Q=$oneline[2];
	}

	if(/^Wall segments:/){last;}
    }

    close sfout;
    return ($f, $U, $Q, $r_ov_Q, $r);
}


#=============================================================================
#    read Ez at side
#=============================================================================
sub side_Ez{
    my($nz);
    my($d, $min_z, $max_z, $r);
    my($i);
    my($mode);
    my($Ez_S, $Ez_E);
    my($line, @data);

    $min_z = $_[0];
    $max_z = $_[1];
    $r = 0.0;
    $nz=1;

    open(IN7,">".$file->in7) or die "open: $!";

    printf IN7 "line noscreen\n";
    printf IN7 "%.3f\t%.3f\t%.3f\t%.3f\n", $min_z, $r ,$max_z, $r;
    printf IN7 "%d\t%d\n", $nz;
    printf IN7 "end";

    close(IN7);

    system "sf7";

    open(SFOUT,"<".$file->sf7);

    while(<SFOUT>){
	chomp;
	s/^\s*(.*?)\s*$/$1/;
	if(/^\(mm\)\s+\(mm\)/){last;}
    }

    $line = <SFOUT>;
    chomp($line);
    $line =~ s/^\s*(.*?)\s*$/$1/;    # 前後の空白の削除
    @data = split(/\s+/,$line);
    $Ez_S = $data[2];

    $line = <SFOUT>;
    chomp($line);
    $line =~ s/^\s*(.*?)\s*$/$1/;
    @data = split(/\s+/,$line);    # 前後の空白の削除
    $Ez_E = $data[2];

    close(SFOUT);

    return($Ez_S, $Ez_E);

}


#=============================================================================
#    set Er or Ht field
#=============================================================================
sub set_field{

    my($fld, $z, $rmin, $rmax, $r, $EH);
    my($Nr, $dr);
    my($idx, $fac);
    my($i, $line, @data);

    $fld  = $_[0];       # Er or Ht
    $z    = $_[1];       # H_theta を格納する z座標
    $rmin = $_[2];       #                    r座標(最小)
    $rmax = $_[3];       #                    r座標(最大)
    $r    = $_[4];       # r座標の配列の参照渡し
    $EH   = $_[5];       # 磁場を格納する配列の参照渡し

    if($fld eq "Er"){
	$idx = 3;
	$fac = 1e+6;
    }elsif($fld eq "Ht"){
	$idx = 5;
	$fac = 1.0;
    }else{
	printf("Error occured at set_field !!\n");
	printf("1st arg of set_field should be \"Er\" or\" Ht\".\n");
	printf("\t arg:%s\n", $fld);
	exit(0);
    }

    $dr   = ($rmax - $rmin)/$N_divR;

    open(IN7,">".$file->in7) or die "open: $!";
    printf IN7 "line noscreen\n";
    printf IN7 "%.3f\t%.3f\t%.3f\t%.3f\n", $z, $rmin ,$z, $rmax;
    printf IN7 "%d\n", $N_divR;
    printf IN7 "end";
    close(IN7);

    system "sf7";

    open(OUTSF7,"<".$file->sf7) or die "at set_field:$file  $!";

    while(<OUTSF7>){
	chomp;
	s/^\s*(.*?)\s*$/$1/;
	if(/^\(mm\)\s+\(mm\)\s+\(MV\/m\)\s+\(MV\/m\)\s+\(MV\/m\)\s+\(A\/m\)/){last;}
    }

    for($i=0; $i<=$N_divR; $i++){
	$line = <OUTSF7>;
	chomp($line);
	$line =~ s/^\s*(.*?)\s*$/$1/;    # 前後の空白の削除
	if($line=~/^\s*$/){last;}        # 空行ならば読み込み終了
	@data = split(/\s+/,$line);	
	if(0.001 < abs($data[0]-$z) || 0.001 < abs($data[1]-$i*$dr)){
	    printf("Error occured at set_field !!\n");
	    printf("\t\t%d\t\tZ = %f\t\t R = %f\n", $i, $data[0], $data[1]);
	    exit(0);
	}
	$$r[$i]  = $data[1];
	$$EH[$i] = $fac*$data[$idx];
    }

    if(0.001<abs($$r[0]-$rmin) || 0.001<abs($$r[$N_divR]-$rmax)){
	printf("Error occured at set_field !!\n");
	printf("r[0]:%f\t\trmin:%f\n", $$r[0], $rmin);
	printf("r[%d]:%f\t\trmax:%f\n", $N_divR, $$r[$N_divR], $rmax);
	exit(0);
    }

    close OUTSF7;

}

#=============================================================================
#    calculate RF power case of stored energy 2 [Joules] in 1.5 cells
#=============================================================================
sub cal_RF_pwr{
    my($U_unit);
    my($fac_H, $fac_E);
    my($E, $H, $r, $dr);
    my($E, $H);
    my($pwrL, $pwrR);
    my($i);

    $U_unit = 1;                         # " [Joules]"

    $fac_H = sqrt($U_unit/$U_SS);
    $fac_E = sqrt($U_unit/$U_OO);

    #========= left side calculation ========
    $pwrL = 0.0;
    for($i=0; $i<$N_divR; $i++){
	$E = $fac_E*($EL_OO[$i]+$EL_OO[$i+1])/2.0;
	$H = $fac_H*($HL_SS[$i]+$HL_SS[$i+1])/2.0;
	if(0.001<abs($rL_SS[$i]-$rL_OO[$i])){
	    printf("error occured at cal_RF_pwr !!\n");
	    printf("\t left side r coodinate.\n");
	    printf("\t%d\trL_SS:%f\trL_OO:%f\n", $rL_SS[$i], $rL_OO[$i]);
	    exit(0);
	}else{
	    $r  = ($rL_SS[$i]+$rL_SS[$i+1])/2.0*1e-3;   # r  [m]
	    $dr = ($rL_SS[$i+1]-$rL_SS[$i])*1e-3;       # dr [m]
	}
	$pwrL += 2*pi*$r*$dr*$E*$H;	
    }
    	$pwrL /= 2.0;


    #========= Right side calculation ========
    $pwrR = 0.0;
    for($i=0; $i<$N_divR; $i++){
	$E = $fac_E*($ER_OO[$i]+$ER_OO[$i+1])/2.0;
	$H = $fac_H*($HR_SS[$i]+$HR_SS[$i+1])/2.0;
	if(0.001<abs($rR_SS[$i]-$rR_OO[$i])){
	    printf("error occured at cal_RF_pwr !!\n");
	    printf("\t right side r coodinate.\n");
	    printf("\t%d\trR_SS:%f\trR_OO:%f\n", $rR_SS[$i], $rR_OO[$i]);
	    exit(0);
	}else{
	    $r  = ($rR_SS[$i]+$rR_SS[$i+1])/2.0*1e-3;   # r  [m]
	    $dr = ($rR_SS[$i+1]-$rR_SS[$i])*1e-3;       # dr [m]
	}
	
	$pwrR += 2.0*pi*$r*$dr*$E*$H;
	
    }
    	$pwrR /= 2.0;


    return($U_unit*2.0, $pwrL, $pwrR);

}


#=============================================================================
#    make a gdf file
#=============================================================================
sub make_gdf{
    my($nz, $nr);
    my($d, $min_z, $max_z, $min_r, $max_r);
    my($i);
    my($gdf_file);

    $gdf_file = $_[0];

    $i = $cal->i;

    $d=0.1;
    $min_z = 0.0;
    $max_z = $tw_result[$i]->shape->D;
    $min_r = 0.0;
    $max_r = $tw_result[$i]->shape->a2/2.0+1.0;

    $nz=int(($max_z-$min_z)/$d);
    $nr=int(($max_r-$min_r)/$d);

    open(IN7,">".$file->in7) or die "open: $!";

    printf IN7 "rect noscreen\n";
    printf IN7 "%.3f\t%.3f\t%.3f\t%.3f\n", $min_z, $min_r ,$max_z, $max_r;
    printf IN7 "%d\t%d\n", $nz, $nr;
    printf IN7 "end";

    close(IN7);

    system "sf7";
    system "fish2gdf -o ".$gdf_file." ". $file->sf7;
    
}

#=============================================================================
#        initialize standing wave region
#=============================================================================
sub mk_gdf_file_name{
    my($mode, $file_name, $add, $ext, $gdf_file);
    my($i);

    $mode = $_[0];
    $i = $cal->i;

    $file_name = $file->gdf;

    $file_name =~ s/^(.*)[.].*$/$1/;    # 拡張子削除
    $ext = $file->gdf;
    $ext =~ s/^(.*)[.](.*)$/$2/;  # 拡張子を取り出し
    $add = sprintf("%02d", $i);
    $gdf_file = $file_name.'_'.$add.'_'.$mode.'.'.$ext;

    return($gdf_file);

}

#=============================================================================
#        initialize standing wave region
#=============================================================================
sub init_SW_region{
    my($N);
    my(@in_dat);
    my $find = 0;
    my $i = 0;

    $inSW  = SW_Region -> new();
    $outSW = SW_Region -> new();

    $inSW -> ra(($input[1]->shape->t)/2.0);
    $inSW -> b2($input[1]->shape->b2);
    $inSW -> D_half(($input[1]->shape->D - $input[1]->shape->t)/2.0);

    $N = $cal->N;

    $outSW -> ra(($input[$N]->shape->t)/2.0);
    $outSW -> b2($input[$N]->shape->b2);
    $outSW -> D_half(($input[$N]->shape->D - $input[$N]->shape->t)/2.0);

    open(CELL_DATA,"<".$file->input) or die "open: $!";

    while(<CELL_DATA>){
	chomp;
	if(/^#/){next;}               # 先頭に'#'があるとコメント文扱い
	if(/#/){s/(^.+)#.*/$1/}       # '#'以降はコメント文扱い
	if($find or /^\$sw/){       # 先頭に$cellがあるまでスキップ
	    if($find == 0){$find = 1;next;}
	}else{
	    next;
	}
	if(/^\$end/){last;}           # データの終わり
	s/^\s*(.*?)\s*$/$1/;          # 前後の空白を削除

	@in_dat = split(/\s+/,$_);

	if($i==0){
	    $inSW  -> pipe_L($in_dat[0]);
	    $inSW  -> a2($in_dat[1]);
	}else{
	    $outSW -> pipe_L($in_dat[0]);
	    $outSW -> a2($in_dat[1]);
	}

	$i++;
    }

    close(CELL_DATA);

}


#=============================================================================
#        write superfish data for input SW region
#=============================================================================
sub mk_sw_region{
    my($freq, $U, $Q, $r, $r_ov_Q, $Rs);
    my($Ez_S, $Ez_E);
    my($gdf_file);

    printf "\n\n------------ RF input coupler SW region ---------------\n";
    &write_sf_inSW;
    system "autofish ".$file->af;
    ($freq, $U, $Q, $r_ov_Q, $Rs)=&read_SFO;
    ($Ez_S, $Ez_E) = &side_Ez(0, $inSW->pipe_L + $inSW->D_half);

    printf "\t Frequency    \t%f.5\n", $freq/$MHz;
    printf "\t Stored Energy\t%e\n", $U;
    printf "\t Q value      \t%.2f\n", $Q;
    printf "\t r/Q          \t%e\n", $r_ov_Q;
    printf "\t Rs           \t%e\n", $Rs;
    printf "\t Ez(start)    \t%e\n", $Ez_S;
    printf "\t Ez(end)      \t%e\n", $Ez_E;

    $sw_result[1] = SW_Result->new();
    $sw_result[1] -> i(1);

    $sw_result[1] -> shape(SW_Region->new());
    $sw_result[1] -> shape($inSW);

    $sw_result[1]->SS(ResutSFO->new());
    $sw_result[1]->SS->freq($freq);
    $sw_result[1]->SS->U($U);
    $sw_result[1]->SS->Q($Q);
    $sw_result[1]->SS->ro_v_Q($r_ov_Q);
    $sw_result[1]->SS->Rs($Rs);
    $sw_result[1]->SS->Ez_S($Ez_S);
    $sw_result[1]->SS->Ez_E($Ez_E);
    $gdf_file = &mk_gdf_file_name_sw("inSW");
    $sw_result[1]->gdf($gdf_file);
    &make_gdf_sw($sw_result[1]);

    printf "\n\n------------ RF output coupler SW region ---------------\n";
    &write_sf_outSW;
    system "autofish ".$file->af;
    ($freq, $U, $Q, $r_ov_Q, $Rs)=&read_SFO;
    ($Ez_S, $Ez_E) = &side_Ez(0, $outSW->pipe_L + $outSW->D_half);

    printf "\t Frequency    \t%f.5\n", $freq/$MHz;
    printf "\t Stored Energy\t%e\n", $U;
    printf "\t Q value      \t%.2f\n", $Q;
    printf "\t r/Q          \t%e\n", $r_ov_Q;
    printf "\t Rs           \t%e\n", $Rs;
    printf "\t Ez(start)    \t%e\n", $Ez_S;
    printf "\t Ez(end)      \t%e\n", $Ez_E;

    $sw_result[2] = SW_Result->new();
    $sw_result[2] -> i(2);

    $sw_result[2] -> shape(SW_Region->new());
    $sw_result[2] -> shape($outSW);

    $sw_result[2]->SS(ResutSFO->new());
    $sw_result[2]->SS->freq($freq);
    $sw_result[2]->SS->U($U);
    $sw_result[2]->SS->Q($Q);
    $sw_result[2]->SS->ro_v_Q($r_ov_Q);
    $sw_result[2]->SS->Rs($Rs);
    $sw_result[2]->SS->Ez_S($Ez_S);
    $sw_result[2]->SS->Ez_E($Ez_E);
    $gdf_file = &mk_gdf_file_name_sw("outSW");
    $sw_result[2]->gdf($gdf_file);
    &make_gdf_sw($sw_result[2]);

    &write_sw_result;
}


#=============================================================================
#        write superfish data for input SW region
#=============================================================================
sub write_sf_inSW {
    my($nbslf, $nbsrt);
    my($pl, $a, $ra, $b, $D);
    my($xdri, $ydri);

    $pl = $inSW-> pipe_L;
    $a  = ($inSW -> a2)/2;
    $ra = $inSW -> ra;
    $b  = ($inSW -> b2)/2;
    $D  = $inSW -> D_half;


    $nbslf = 1;
    $nbsrt = 1;
    $xdri  = 0.99*($pl+$D);
    $ydri  = 0.99*$b;

    open(SF_DATA,">".$file->af) or die "open: $!";

    printf SF_DATA "%s\n", $SFPara -> title;
    printf SF_DATA " \$reg kprob=1, ";
    printf SF_DATA "dx=%f ,", $SFPara->dx;
    printf SF_DATA "conv=%f ,", $SFPara->conv; 
    printf SF_DATA "kprob=1,\n";
    printf SF_DATA " xdri=$xdri, ydri=$ydri,\n";
    printf SF_DATA " nbsup=1, nbslo=0, nbslf=$nbslf, nbsrt=$nbsrt\n";
    printf SF_DATA " freq=%f,\n",$freq/$MHz;
    printf SF_DATA " kmethod=1,\n";
    printf SF_DATA " beta=%f\$\n",1;

    printf SF_DATA "!---------- input SW region -----------\n";
    printf SF_DATA " \$po x=%f, y=%f \$\n", 0.0, 0.0;
    printf SF_DATA " \$po x=%f, y=%f \$\n", 0.0, $a;
    printf SF_DATA " \$po x=%f, y=%f \$\n", $pl-$ra, $a;
    printf SF_DATA
	" \$po nt=2, r=%f, theta=0, x0=%f, y0=%f \$\n",
	$ra, $pl-$ra, $a+$ra;
    printf SF_DATA " \$po x=%f, y=%f \$\n", $pl, $b;
    printf SF_DATA " \$po x=%f, y=%f \$\n", $pl+$D, $b;
    printf SF_DATA " \$po x=%f, y=%f \$\n", $pl+$D, 0.0;
    printf SF_DATA " \$po x=%f, y=%f \$\n", 0.0, 0.0;

    close SF_DATA;

}

#=============================================================================
#        write superfish data for output SW region
#=============================================================================
sub write_sf_outSW {
    my($nbslf, $nbsrt);
    my($pl, $a, $ra, $b, $D);
    my($xdri, $ydri);

    $pl = $outSW-> pipe_L;
    $a  = ($outSW -> a2)/2;
    $ra = $outSW -> ra;
    $b  = ($outSW -> b2)/2;
    $D  = $outSW -> D_half;


    $nbslf = 1;
    $nbsrt = 1;
    $xdri  = 0.0;
    $ydri  = 0.99*$b;

    open(SF_DATA,">".$file->af) or die "open: $!";

    printf SF_DATA "%s\n", $SFPara -> title;
    printf SF_DATA " \$reg kprob=1, ";
    printf SF_DATA "dx=%f ,", $SFPara->dx;
    printf SF_DATA "conv=%f ,", $SFPara->conv; 
    printf SF_DATA "kprob=1,\n";
    printf SF_DATA " xdri=$xdri, ydri=$ydri,\n";
    printf SF_DATA " nbsup=1, nbslo=0, nbslf=$nbslf, nbsrt=$nbsrt\n";
    printf SF_DATA " freq=%f,\n",$freq/$MHz;
    printf SF_DATA " kmethod=1,\n";
    printf SF_DATA " beta=%f\$\n",1;

    printf SF_DATA "!---------- output SW region -----------\n";
    printf SF_DATA " \$po x=%f, y=%f \$\n", 0.0, 0.0;
    printf SF_DATA " \$po x=%f, y=%f \$\n", 0.0, $b;
    printf SF_DATA " \$po x=%f, y=%f \$\n", $D, $b;
    printf SF_DATA " \$po x=%f, y=%f \$\n", $D, $a+$ra;
    printf SF_DATA
	" \$po nt=2, r=%f, theta=270, x0=%f, y0=%f \$\n",
	$ra, $D+$ra, $a+$ra;
    printf SF_DATA " \$po x=%f, y=%f \$\n", $D+$pl, $a;
    printf SF_DATA " \$po x=%f, y=%f \$\n", $pl+$D, 0.0;
    printf SF_DATA " \$po x=%f, y=%f \$\n", 0.0, 0.0;

    close SF_DATA;

}


#=============================================================================
#        write result
#=============================================================================
sub write_sw_result{
    my($i);

    $i=1;

    open(SW_DATA,">>".$file->sw) or die "open: $!";

    for($i=1; $i<=2; $i++){
	printf SW_DATA "%d\t",   $sw_result[$i]->i;
	printf SW_DATA "%.3f\t", $sw_result[$i]->shape->pipe_L;
	printf SW_DATA "%.3f\t", $sw_result[$i]->shape->a2;
	printf SW_DATA "%.3f\t", $sw_result[$i]->shape->ra;
	printf SW_DATA "%.3f\t", $sw_result[$i]->shape->b2;
	printf SW_DATA "%.3f\t", $sw_result[$i]->shape->D_half;
	printf SW_DATA "%.5f\t", $sw_result[$i]->SS->freq/$MHz;
	printf SW_DATA "%e\t",   $sw_result[$i]->SS->U;
	printf SW_DATA "%.2f\t", $sw_result[$i]->SS->Q;
	printf SW_DATA "%e\t",   $sw_result[$i]->SS->ro_v_Q;
	printf SW_DATA "%e\t",   $sw_result[$i]->SS->Rs;
	printf SW_DATA "%e\t",   $sw_result[$i]->SS->Ez_S;
	printf SW_DATA "%e\t",   $sw_result[$i]->SS->Ez_E;
	printf SW_DATA "%s\n",   $sw_result[$i]->gdf;
    }

    close SW_DATA;
    
}

#=============================================================================
#    make a gdf file
#=============================================================================
sub make_gdf_sw{
    my($nz, $nr);
    my($d, $min_z, $max_z, $min_r, $max_r);
    my($gdf_file);
    my($sw);
    my($command);

    $sw = SW_Result->new();
    $sw = $_[0];

    $d=0.1;
    $min_z = 0.0;
    $max_z = $sw->shape->D_half + $sw->shape->pipe_L;
    $min_r = 0.0;
    $max_r = ($sw->shape->a2)/2.0+1;

    $nz=int(($max_z-$min_z)/$d);
    $nr=int(($max_r-$min_r)/$d);

    # --- make .in7 file -----------
    open(IN7,">".$file->in7) or die "open: $!";
    printf IN7 "rect noscreen\n";
    printf IN7 "%.3f\t%.3f\t%.3f\t%.3f\n", $min_z, $min_r ,$max_z, $max_r;
    printf IN7 "%d\t%d\n", $nz, $nr;
    printf IN7 "end";
    close(IN7);

    # --- make gdf file ---------
    $gdf_file = $sw->gdf;
    system "sf7";
    $command = sprintf("fish2gdf -o %s %s",$gdf_file, $file->sf7); 
    system $command;
    
}

#=============================================================================
#    make a gdf file name for SW region
#=============================================================================
sub mk_gdf_file_name_sw{
    my($file_name, $ext, $gdf_file);
    my($mode);

    $mode = $_[0];
    $file_name = $file->gdf;
    $file_name =~ s/^(.*)[.].*$/$1/;    # 拡張子削除
    $ext = $file->gdf;
    $ext =~ s/^(.*)[.](.*)$/$2/;  # 拡張子を取り出し    
    $gdf_file = $file_name.'_'.$mode.'.'.$ext;

    return($gdf_file);
}


#=============================================================================
#      delete files
#=============================================================================
sub rm_files{
    my($fname);

    opendir(DIR,"./") or die "could not open directory: $!";

    unlink($file->af);
    unlink($file->seg);
    unlink($file->T35);
    unlink($file->in7);
    unlink($file->SFO);
    unlink($file->sf7);
    unlink("OUTAUT.TXT");
    unlink("OUTFIS.TXT");
    unlink("TAPE35.INF");
    unlink($file->base.".PMI");

}