use strict;
use Class::Struct;
use Math::Trig;
use File::Find;
use File::Basename;
use FileHandle;

struct cell_struct=>{
    beta   => '$',
    a2     => '$',
    b2     => '$',
    t      => '$',
    zs     => '$',
    zc     => '$',
    ze     => '$',
    D      => '$',
    p0     => '$'
};

struct structure_struct=>{
    start_str => '$',          # z-position of 1st cell center [mm]   *.in
    BaPhase   => '$',          # initial phase of buncher      [deg]  *.bat
    N         => '$'           # number of cells                      *.dat
};


our($structure, @cell);
our($structure_file, $GDF_in, $GDF_out, $temp_time, $temp_phase);
our($frequency, $omega, $omega_deg, $lambda, $k_deg, $c);

$structure = structure_struct->new();


&initialize;
&read_cell_data_file;
&set_cell_pzD;
&print_cell_zp;
&print_structure_data;

&gdf2a;
&time2phase;
&asci2gdf;
&delete_temp;

#=============================================================================
#        read structure data
#=============================================================================
sub initialize{

    $structure_file = "../../TW/TW_structure.dat";
    $GDF_in         = "../TW_Buncher.gdf";
    $GDF_out        = "phase.gdf";
    $temp_time      = "temp_time.txt";
    $temp_phase     = "temp_phase.txt";

    $frequency = 2.856e9;
    $omega     = 2*pi*$frequency;
    $omega_deg = $frequency*360;

    $c      = 2.99792458e+8;   # speed of light [m/sec]
    $lambda = $c/$frequency;   # wavelength [m]
    $k_deg  = 360.0/$lambda;

    $structure->start_str(100);    # start point: center of 1st cell(input cpl)
    $structure->BaPhase(30);       # initial phase [deg]
    
}

#=============================================================================
#        read input file
#=============================================================================
sub read_cell_data_file{

    my($fh);
    my(@in_dat);

    $fh = new FileHandle;
    
    printf "read input data form %s\n", $structure_file;

    $fh->open("<".$structure_file) or die "open: $!";

    while(<$fh>){
	chomp;
	if(/^#/){next;}               # 先頭に'#'があるとコメント文
	if(/#/){s/(^.+)#.*/$1/}       # '#' 以降はコメント
	s/^\s*(.*?)\s*$/$1/;          # 前後の空白を削除

	if(/^\$directory/){&read_dir($fh);}
	if(/^\$input_RF/){&read_RF($fh);}
	if(/^\$cell/){&read_cell($fh);}
    }

    $fh->close;

}

#=============================================================================
#        read directory
#=============================================================================
sub read_dir{
    my $fh;

    $fh = $_[0];

    while(<$fh>){
	chomp;
	if(/^#/){next;}               # 先頭に'#'があるとコメント文
	if(/#/){s/(^.+)#.*/$1/}       # '#' 以降はコメント
	s/^\s*(.*?)\s*$/$1/;          # 前後の空白を削除
	if(/^\$end/){last;};

	printf "\tdirectrory:%s\n", $_;

	last;
    }    
}

#=============================================================================
#        read RF
#=============================================================================
sub read_RF{
    my $fh;

    $fh = $_[0];

    while(<$fh>){
	chomp;
	if(/^#/){next;}               # 先頭に'#'があるとコメント文
	if(/#/){s/(^.+)#.*/$1/}       # '#' 以降はコメント
	s/^\s*(.*?)\s*$/$1/;          # 前後の空白を削除
	if(/^\$end/){last;};

	printf "\tRF Power:%.3f [MW]\n", $_;

	last;
    }    
}

#=============================================================================
#        read cell
#=============================================================================
sub read_cell{
    my $fh;
    my @input;
    my $i;

    $fh = $_[0];

    printf "\tbeta\t2a\t2b\tt\trb\tn\n";

    $i=0;
    while(<$fh>){
	chomp;
	if(/^#/){next;}               # 先頭に'#'があるとコメント文
	if(/#/){s/(^.+)#.*/$1/}       # '#' 以降はコメント
	s/^\s*(.*?)\s*$/$1/;          # 前後の空白を削除
	if(/^\$end/){last;};
	@input = split(/\s+/,$_);

	our $one_cell = cell_struct->new();

	$one_cell->beta($input[0]);
	$one_cell->a2($input[1]);
	$one_cell->b2($input[2]);
	$one_cell->t($input[3]);

	$i++;
	$cell[$i]=$one_cell;

	printf "\t%.3f\t%.3f\t%.3f\t%.3f\n", 
	       $cell[$i]->beta, $cell[$i]->a2, $cell[$i]->b2, $cell[$i]->t;
    }

    $structure->N($i);
}

#=============================================================================
# check every cell data 
#=============================================================================
sub set_cell_pzD{
    my($D, $i, $N);
    
    #------------ sw region ---------
    our $one_cell = cell_struct->new();

    $D=$cell[1]->beta * $lambda /3.0;
    $one_cell->D($D);
    $one_cell->beta($cell[1]->beta);
    $one_cell->p0($structure->BaPhase);
    $one_cell->ze($structure->start_str*1e-3);
    $one_cell->zs($one_cell->ze-$D/2.0);
    $one_cell->zc(($one_cell->zs + $one_cell->ze)/2.0);

    $cell[0] = $one_cell;

    #------------  in tw region ---------

    for($i=1; $i<=$structure->N; $i++){
	$D=$cell[$i]->beta * $lambda /3.0;

	$cell[$i]->D($D);
	$cell[$i]->p0($structure->BaPhase-($i-1)*120);
	$cell[$i]->zs($cell[$i-1]->ze);
	$cell[$i]->ze($cell[$i]->zs+$D);
	$cell[$i]->zc(($cell[$i]->zs + $cell[$i]->ze)/2.0);
    }

    #------------- sw region ---------
    $N=$structure->N;
    our $one_cell = cell_struct->new();

    $D=$cell[$N]->beta * $lambda /3.0;
    $one_cell->D($D);
    $one_cell->beta($cell[$N]->beta);
    $one_cell->p0($cell[$N]->p0-120);
    $one_cell->zs($cell[$N]->ze);
    $one_cell->ze($cell[$N]->ze+$D/2.0);
    $one_cell->zc(($one_cell->zs + $one_cell->ze)/2.0);

    $cell[$N+1] = $one_cell;

}

#=============================================================================
# print cell posiotion initial phase
#=============================================================================
sub print_cell_zp{
    my($i);

    printf "\n---------------------------------------------------------------------\n";
    printf "\t#\tzs\t\tzc\t\tze\t\tphase\n";
    printf "=====================================================================\n";
    

    for($i=0; $i<=$structure->N+1; $i++){
	printf "\t%d\t", $i;
	printf "%10.3f\t", $cell[$i]->zs * 1e3;
	printf "%10.3f\t", $cell[$i]->zc * 1e3;
	printf "%10.3f\t", $cell[$i]->ze * 1e3;
	printf "%7.1f\n", $cell[$i]->p0;
    }

    printf "---------------------------------------------------------------------\n\n";

}

#=============================================================================
# print structure data 
#=============================================================================
sub print_structure_data{

    printf "\n-------------------------------------------------------\n";
    printf "\t 1st cell center position : %.3f[mm]\n",  $structure->start_str;
    printf "\t Buncher initial phase    : %.3f[deg]\n", $structure->BaPhase;
    printf "\t Numbe of cells           : %d\n",        $structure->N;
    printf "-------------------------------------------------------\n\n";

}




#=============================================================================
# convert gdf file to ascii file
#=============================================================================
sub gdf2a{

    system("gdf2a -w 15 -o $temp_time $GDF_in z t");

}


#=============================================================================
# convert to phase from time
#=============================================================================
sub time2phase{

    my($z, $time);
    my(@input);
    my($n, $dp);
    my $print_zp = 0;

    open(TIME, "<$temp_time") or die "file:$temp_time\n $!";
    open(PHAS, ">$temp_phase") or die "file:$temp_phase\n $!";

    foreach(<TIME>){
	if(/^\s*\d/){
	    chomp;
	    s/^\s+//;                        # 行頭の空白を削除
	    @input = split(/\s+/,$_);
	    $z    = $input[0];
	    $time = $input[1];
	    if( ($z< $cell[0]->zs) || ($cell[$structure->N+1]->ze < $z)){next;}
	    ($n, $dp) = &z2phase($z, $time);
	    printf PHAS "   %f  %f\n", $z , $dp;
	}elsif(/^\s*z\s+t/){
	    if($print_zp !=1){
		print PHAS "                 z                phase\n";
	    }
	    $print_zp =1;
	}else{
	    print PHAS "#$_";
	}
    }

    close TIME;
    close PHAS;

}
#=============================================================================
# find cell 
#=============================================================================
sub z2phase{

    my $z = $_[0];
    my $t = $_[1];
    my($dz);
    my($n, $dp, $p);

    $n = &find_cell($z);
    $dz = $z-$cell[$n]->zs;

    if($n==0 || $n==$structure->N+1){
	$p = $cell[$n]->p0 + 360*$frequency*$t;
    }else{
	$p = $cell[$n]->p0 + 360*($frequency*$t - $dz/($cell[$n]->beta * $lambda));
    }

    $dp = $p-int($p/360.0)*360;

    if($dp<-180){$dp += 360.0;}
    if(180<$dp) {$dp -= 360.0;}
 
    return ($n, $dp);

}

#=============================================================================
# find cell 
#=============================================================================
sub find_cell{
    my $z = $_[0];
    my($min, $max, $tst);

    $min = 0;
    $max = $structure->N+1;

     do{
	$tst = int(($max + $min)/2);

	if($cell[$tst]->ze < $z){
	    $min = $tst+1;
	}else{
	    $max = $tst;
	}

    }while($max != $min);

    return $max;
}


#=============================================================================
# convert gdf file to ascii file
#=============================================================================
sub asci2gdf{

    system("asci2gdf -o $GDF_out $temp_phase");

}

#=============================================================================
# delete temp file
#=============================================================================
sub delete_temp{

    system("rm -f $temp_time");
    system("rm -f $temp_phase");

}
