print "SUPERFISH has been started.\n";

$file=@ARGV[0];

$file_af   = $file.".af";
$file_T35  = $file.".T35";
$file_in7  = $file.".in7";
$file_SFO  = $file.".SFO";
$file_cell = $file.".dat";
$file_sf7  = "OUTSF7.TXT";
$file_gdf  = $file.".gdf";

$file_name=$current_dir.$file_af;

$it_max=10;
$error=1e-6;

#------- cell shape parameter -----------
$title="S-Band Prebuncher for GPT test";
$mm=1;
$m=1e+3*$mm;
$sec=1;
$MHz=1e+6;

$conv = 0.1;
$dx   = 0.2;
$freq = 2856;
$beta = 0.54822088;

$pi=atan2(1,1)*4.0;
$light=2.99792458e+8*$m/$sec;
$lambda=$light/($freq*$MHz);


&read_cell_data;       # read cell data form input cell geometory file
&tune_b_acc;           # Tune b which is acc structure
&write_result;         # write result cell geometory
&make_gdf;             # make a gdf file


#====================================================
#        tunning the accelerating cavity radius
#====================================================
sub tune_b_acc{

    printf "\n";
    @tune_b_acc[1]=$b2;
    for($i=1; $i<=$it_max; $i++){
	
	write_sf_data(1, 1, $tune_b_acc[$i]);
	
	system "autofish $file_af";
	@sf_frequency[$i]=&read_fr($file_SFO);
	printf "\t%d\t%f\t%f\n",$i,@tune_b_acc[$i],@sf_frequency[$i];
	if(abs((@sf_frequency[$i]-$freq)/$freq)<$error){last};
	
	
	if($i==1){
	    @tune_b_acc[2]=@sf_frequency[1]*@tune_b_acc[1]/$freq;
	}else{
	    @tune_b_acc[$i+1]=
		(@tune_b_acc[$i]-@tune_b_acc[$i-1])
		/(@sf_frequency[$i]-@sf_frequency[$i-1])
		*($freq-@sf_frequency[$i])
		+@tune_b_acc[$i];
	}
    }    
}

#====================================================
#        read data
#====================================================
sub read_cell_data{

    open(CELL_DATA,"<$file_cell") or die "open: $!";

    while(<CELL_DATA>){
	chomp;
	if(/^#/){
	    next;
	}

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

	$a2   = @oneline[0];
	$b2   = @oneline[1];
	$D    = @oneline[2];
	$r    = @oneline[3];
	$gap  = @oneline[4];
	$addL = @oneline[5];
	last;	
    }

    close(CELL_DATA);

    &check_data;

}

#====================================================
#        check data
#====================================================
sub check_data{

    printf "2a    = %f\n",$a2;
    printf "2b    = %f\n",$b2;
    printf "D     = %f\n",$D;
    printf "gap   = %f\n",$gap;
    printf "addL  = %f\n",$addL;
}


#====================================================
#        write result
#====================================================
sub write_result{
    open(CELL_DATA,">>$file_cell") or die "open: $!";

    printf CELL_DATA "\n%f\t%f\t%f\t%f\t%f\t%f",
    $a2, $b2, $D, $r, $gap, $addL;

    close(CELL_DATA);

    printf "\n%f\t%f\t%f\t%f\t%f\t%f",
    $a2, $b2, $D, $r, $gap, $addL;
}

#====================================================
#        write superfish data
#====================================================
sub write_sf_data {
    local($nbslf, $nbsrt);

    $nbslf = $_[0];
    $nbsrt = $_[1];
    $b2    = $_[2];
    $a     = $a2/2.0;
    $b     = $b2/2.0;
    $xdri  = 0.0;
    $ydri  = 0.99*$b;

    open(SF_DATA,">$file_af") or die "open: $!";

    printf SF_DATA "$title\n";
    printf SF_DATA " \$reg kprob=1, dx=$dx, conv=$conv, 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;
    printf SF_DATA " kmethod=1,\n";
    printf SF_DATA " beta=%f\$\n",$beta;

    # ------ PB cell (lefr half)  -----------
    printf SF_DATA "!--- PB cell (left half)  ---\n";
    printf SF_DATA " \$po x=%f, y=%f \$\n", 0.0, 0.0;
    printf SF_DATA " \$po x=%f, y=%f \$\n", -($D/2+$addL), 0.0;
    printf SF_DATA " \$po x=%f, y=%f \$\n", -($D/2+$addL), $a;
    printf SF_DATA " \$po x=%f, y=%f \$\n", -($gap/2+$r),  $a;
    printf SF_DATA
	" \$po nt=2, r=%f, theta=0, x0=%f, y0=%f \$\n",
	$r,  -($gap/2+$r),  $a+$r;
    printf SF_DATA
	" \$po nt=2, r=%f, theta=90, x0=%f, y0=%f \$\n",
	$r,  -($gap/2+$r),  $a+$r;
    printf SF_DATA " \$po x=%f, y=%f \$\n", -$D/2+$r, $a+2*$r;
    printf SF_DATA
	" \$po nt=2, r=%f, theta=180, x0=%f, y0=%f \$\n",
	$r, -$D/2+$r, $a+3*$r;
    printf SF_DATA " \$po x=%f, y=%f \$\n", -$D/2, $b-$r;
    printf SF_DATA
	" \$po nt=2, r=%f, theta=90, x0=%f, y0=%f \$\n",
	$r, -$D/2+$r, $b-$r;
    printf SF_DATA " \$po x=%f, y=%f \$\n", 0.0, $b;

    # ------PB cell (right half)  -----------
    printf SF_DATA "!--- PB cell (right half)  ---\n";
    printf SF_DATA " \$po x=%f, y=%f \$\n", $D/2-$r, $b;
    printf SF_DATA
	" \$po nt=2, r=%f, theta=0, x0=%f, y0=%f \$\n",
	$r, $D/2-$r, $b-$r;
    printf SF_DATA " \$po x=%f, y=%f \$\n", $D/2, $a+3*$r;
    printf SF_DATA
	" \$po nt=2, r=%f, theta=270, x0=%f, y0=%f \$\n",
	$r, $D/2-$r, $a+3*$r;
    printf SF_DATA " \$po x=%f, y=%f \$\n", 0.5*$gap+$r, $a+2*$r;
    printf SF_DATA
	" \$po nt=2, r=%f, theta=180, x0=%f, y0=%f \$\n",
	$r,  $gap/2+$r,  $a+$r;
    printf SF_DATA
	" \$po nt=2, r=%f, theta=270, x0=%f, y0=%f \$\n",
	$r,  $gap/2+$r,  $a+$r;
    printf SF_DATA " \$po x=%f, y=%f \$\n", 0.5*$D+$addL, $a;
    printf SF_DATA " \$po x=%f, y=%f \$\n", 0.5*$D+$addL, 0.0;
    printf SF_DATA " \$po x=%f, y=%f \$\n", 0.0, 0.0;
    close SF_DATA;

}

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

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

    close sfout;
    return $resonance_fr;
}

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

    my($file);
    my($f, $U, $Q, $r, $r_ov_Q);

    $file   = $_[0];
    $f      = 0;
    $U      = 0;
    $Q      = 0;
    $r      = 0;
    $r_ov_Q = 0;
    
    open(sfout,"<$file") 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];
	}

	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);
}

#====================================================
#    make a gdf file
#====================================================
sub make_gdf{
    my($NZP, $NZM);
    my($d, $max_r);
    
    $d=0.1;
    $min_r = 0;
    $max_r = $a2/2;
    $min_z = -$D/2-$add_L;
    $max_z =  $D/2+$add_L;

    open(IN7,">$file_in7") or die "open: $!";

    printf IN7 "rect noscreen\n";
    printf IN7 "%.2f\t%.2f\t%.2f\t%.2f\n", $min_z, $min_r ,$max_z, $max_r;
    printf IN7 "%d\t%d\n", ($max_z-$min_z)/0.1, ($max_r-$min_r)/0.1;
    printf IN7 "end";

    close(IN7);

    system "sf7";
    system "fish2gdf -o $file_gdf $file_sf7";

    
}