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(<INPUT>){
	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(<INPUT>){
	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(<TW>){
	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(<SW>){
       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);
}