#!/usr/bin/perl
##################################################################
#                                                                #     
#    This program computes Lick indices for spectral lines       #
#    of a given spectra.                                         #
#                                                                #
##################################################################

use IPC::Open2;
require Math::Spline;
use Cwd;

$dir=cwd;
$version="Version 1.0";

#################
#Initial parameters
$quiet="on";

($help,$greeting)=&DefineInitialTexts;


for(@ARGV){
     if ( /^-q$/ || /^--quiet$/ ){
          $quiet="on";
     }
     elsif( /^-h$/ || /^--help$/ ){
          print $help;
          die "\n";
     }
}

system "clear";
print $greeting;

          
mkdir PACCE_output;
mkdir "$dir/PACCE_output/Equiv_Widths";
mkdir "$dir/PACCE_output/continuum_data_points";
mkdir "$dir/PACCE_output/figures/";
mkdir "$dir/PACCE_output/gnuplot_commands/";


$again=1;
while($again > 0){     

     chdir $dir;
     $fig=0;

     print "\n\tInsert ASCII spectrum file name: ";
     $filename1=<STDIN>;
     chomp($filename1);

     
     # Reading spectrum
     open (FILE, "<$filename1") or die "\n\n\t\*\*\* File: $filename1 not found! \*\*\*\n\n";
          @filetemp=<FILE>;
     close FILE;

     $filename1=$dir."/".$filename1;
     
     print "\n\n\tInsert the file name containing indices definitions: ";
     $filename2=<STDIN>;
     chomp($filename2);


     # Reading index definitions
     open (FILE, "<$filename2") or die "\n\n\t\*\*\* File: $filenema2 not found! \*\*\*\n\n";
          @intervals=<FILE>;
     close FILE;
     $filename2=$dir."/".$filename2;
     
     $save="";
     $show=2;
     while ( $show > 1 ){
          print "\n\n\tDo you want to show figures? (y\/n) ";
          $temp=<STDIN>;
          chomp($temp);
          if ( $temp =~ /^n$/i ){
               $save="\#";               
               $show=0;
          }
          elsif( $temp =~ /^y$/i ){
               $show=1;
          }
     }
     $temp=$filename1;     
     for($temp){s/$dir\///;}
     print "\n\n\tProcessing $again\) $temp\n\n";
     print "\tWait a while ";               



     # Inserting the points with the indices in the spectrum
     @file = &prepareSpectrum(\@intervals,\@filetemp);


     # Estipulates variables precision
     foreach $line (@file){
          chomp($line);
          @temp=split(" ",$line);
          $temp[0]=sprintf "%.10f",$temp[0];
          $temp[1]=sprintf "%.6f",$temp[1];
          push(@spectra,"$temp[0]");
          push(@intensity,"$temp[1]");          
     }
     
     
     
     $INTERV=0;
     chdir "./PACCE_output/continuum_data_points";
     for (@intervals){
          
          if(($_ !~ /^\#/) && ($_ !~ /\|/) ){
               $INTERV++;
               print ".";
               # Reading ion name
               /^(.+)\(/;
               $ion_name=$+;
               for($ion_name){s/\s//g;}
               
               # Searching wavelenght limits to compute absorption lines
               /\((.+)\,(.+)\)/;
               $lambda_1=$1;
               $lambda_2=$2;
               
               ######################################################               
               # Searching wavelenght limits to compute continuum
               ######################################################               

               ##
               ## Case 1: linear adjust of continuum 
               ##
               if (/\[(.+)\:(.+)\]\,\[(.+)\:(.+)\](.*)/){
                    $lambdaL_i=$1;
                    $lambdaL_f=$2;
                    $lambdaR_i=$3;
                    $lambdaR_f=$4;
                    $midpoint=$5;
                    
                    # Taking absorption lines 
                    # @spectra_L -> spectra lines of Left continuum band pass 
                    # @intensity_L -> associated @spectra_L intensity values  
                    # @spectra_R -> spectra lines of Right continuum band pass
                    # @intensity_R -> associated @spectra_R intensity values 
                    # @spectra_C -> spectra lines of Central (hypothetical) continuum band pass 
                    # @intensity_C -> associated @spectra_C intensity values 
                    #      Here @spectra_C is a line fitted using least squares of points 
                    #     extracted from the continuum band pass at Left and Right sides
                    
                    
                    %name_endings=(L_i => '_L', L_f => '_L' ,
                                   R_i => '_R', R_f => '_R' ,
                                   _1 => '_C' , _2 => '_C' );
                                   
                    @name_endings=("L\_i","L\_f","R\_i","R\_f","\_1","\_2");
                    
                    foreach $end (@name_endings){
                         splice(@{"spectra".$name_endings{$end}});
                         splice(@{"intensity".$name_endings{$end}});
                    }
                    splice(@linear_spectra);
                    splice(@linear_intensity);
                    
                    foreach $end (@name_endings){
                         if( ($end =~/i/) || ($end =~ /1/) ){
                              @values=&SearchValues(${"lambda".$end},\@spectra,\@intensity);
                              $i=$values[2];     
                              @{"spectra".$name_endings{$end}}=@spectra;
                              @{"spectra".$name_endings{$end}}=splice(@{"spectra".$name_endings{$end}},$i);
                              @{"intensity".$name_endings{$end}}=@intensity;
                              @{"intensity".$name_endings{$end}}=splice(@{"intensity".$name_endings{$end}},$i);
                         }
                         elsif( ($end =~/f/) || ($end =~ /2/) ){
                              @values=&SearchValues(${"lambda".$end},\@{"spectra".$name_endings{$end}},
                                                       \@{"intensity".$name_endings{$end}});
                              $i=$values[2];
                              $i++;
                              splice(@{"spectra".$name_endings{$end}},$i);
                              splice(@{"intensity".$name_endings{$end}},$i);
                         }                    
                    }     

                    ###############################################################
                    #     Important Comment: Explanation about procedure above
                    ###############################################################
                    #
                    #foreach $end (@name_endings){
                    #     print "\@spectra$name_endings{$end}=";
                    #     for(@{"spectra".$name_endings{$end}}){
                    #     print $_,",";                              
                    #     }
                    #     print "\n\n";                         
                    #}
                    #
                    # Here we mantain some explanations about the above procedure
                    # which is repeated for @spectra_C, @spectra_L and @spectra_R
                    # with respective intensities @intensity_C, @intensity_L and 
                    # @intensity_R. 
                    #@values=&SearchValues($lambdaL_i,\@spectra,\@intensity);
                    #$i=$values[2];
                    #@spectra_L=@spectra;
                    #@spectra_L=splice(@spectra_L,$i);
                    #
                    #@values=&SearchValues($lambdaL_f,\@spectra_L,\@intensity_L);
                    #$i=$values[2];
                    #$i++;
                    #splice(@spectra_L,$i);
                    #
                    ###############################################################


                    # Calculating coefficients for a linear continuum ajust.
                    @spectra_L_R=(@spectra_L,@spectra_R);
                    @intensity_L_R=(@intensity_L,@intensity_R);
                    @a_b= &LeastSquaresMethod(\@spectra_L_R,\@intensity_L_R);


                    # Defining linear and angular coefficients for continuum of the line
                    $a=$a_b[0];
                    $b=$a_b[1];
                    # Signal to noise and number of pixeis in the continuum
                    $SN_ratio=$a_b[2];
                    $NpixFc=$a_b[3];
                    $rms=$a_b[4];

                    if($midpoint=~/\*/){
                         splice(@X);
                         splice(@Y);
                         foreach $end ("_L","_R"){
                              $s=0;
                              $sum_s=0;
                              foreach $s (@{"spectra".$end}){$sum_s+=$s;}
                              $tot=scalar(@{"spectra".$end});
                              $mean_s=$sum_s/$tot;
                              push(@X,$mean_s);

                              $i=0;
                              $sum_i=0;
                              foreach $i (@{"intensity".$end}){$sum_i+=$i;}
                              $mean_i=$sum_i/$tot;
                              push(@Y,$mean_i);
                         }
                         $a=($Y[1]-$Y[0])/($X[1]-$X[0]);
                         $b=$Y[1]-$a*$X[1];
                    }

                     $lim_inf=$spectra_C[0];                                        
                    $lim_sup=$spectra_C[-1];
                    $area_C = &ComputeArea($a,$b,$lim_inf,$lim_sup);
                    for ( $X=$spectra_L[0]; $X<=$spectra_R[-1]; $X+=0.025 ){
                                $Y=$a*$X+$b;
                                $X=sprintf "%.4f",$X;
                                $Y=sprintf "%.10f",$Y;
                         push(@linear_spectra,$X);
                         push(@linear_intensity,$Y);
                    }
                    @spectra_L_C_R=(@spectra_L,@spectra_C,@spectra_R);
                    @intensity_L_C_R=(@intensity_L,@intensity_C,@intensity_R);
                    splice(@l_or_sp);
                    @l_or_sp=("linear","linear\_spectra","linear\_intensity","spectra\_L\_C\_R","intensity\_L\_C\_R");
                                        
               }
               ##
               ## Case 2: Spline adjust of continnum 
               ##
               elsif( /\{(.+)\}.*\<(.+)\:(.+)\>/ ){
                    # Defining the nearest point of ASCII spectra file, compared 
                    # to user's guess for continuum, from Band Pass in Spectroscopic 
                    # Indices Definitions. These points are stored in @lambdas_sp 
                    # and the respective intensities in @intensity_sp to be used
                    @lambdas=split(",",$1);                    
                    @sp_SNcont=($2,$3);
                    
                    splice(@sp_SNcont_lamb);
                    splice(@sp_SNcont_flux);
                    splice(@lambdas_sp);
                    splice(@intensity_sp);
                    splice(@spline_spectra);
                    splice(@spline_intensity);

                    for($tt=0;$tt<=$#spectra;$tt++){
                         if (($spectra[$tt] >= $sp_SNcont[0]) && ($spectra[$tt] <= $sp_SNcont[1]) ){
                              push(@sp_SNcont_lamb,$spectra[$tt]);
                              push(@sp_SNcont_flux,$intensity[$tt]);
                         }
                    }
                    
                    @a_b= &LeastSquaresMethod(\@sp_SNcont_lamb,\@sp_SNcont_flux);
                    $SN_ratio=$a_b[2];
                    $NpixFc=$a_b[3];

                    foreach $lambda (@lambdas){

                         @values=&SearchValues($lambda,\@spectra,\@intensity);
                         push(@lambdas_sp,$values[0]);
                         push(@intensity_sp,$values[1]);

                    }
                    # Using @lambdas_sp and @intensity_sp to compute a spline. 
                    # Spline points are stored in @spline_spectra and
                    # @spline_intensity
                    $spline=new Math::Spline(\@lambdas_sp,\@intensity_sp);
               
                              for ( $X=$lambdas_sp[0]; $X<=$lambdas_sp[-1]; $X+=0.025 ){
                         $sp=$spline->evaluate($X);
                         $X=sprintf "%.4f",$X;
                         $sp=sprintf "%.10f",$sp;
                         push(@spline_spectra,$X);
                         push(@spline_intensity,$sp);
                    }
                    @spline_spectra_temp=@spline_spectra;
                    @spline_intensity_temp=@spline_intensity;
                    
                    # Searching interval indices of original spectra between $lambda_1 and $lambda_2
                    @values=&SearchValues($lambda_1,\@spectra,\@intensity);
                    $i=$values[2];
                    @spectra_C=@spectra;
                    @intensity_C=@intensity;
                    @spectra_C=splice(@spectra_C,$i);
                    @intensity_C=splice(@intensity_C,$i);
                    
                    @values=&SearchValues($lambda_2,\@spectra_C,\@intensity_C);
                    $i=$values[2];
                    $i++;
                    splice(@spectra_C,$i);
                    splice(@intensity_C,$i);

                    $lim_inf=$spectra_C[0];                                        
                    $lim_sup=$spectra_C[-1];           
                    # Defining superior and inferior limits for areas integration, in case 
                    # of a spline continuum. These integration limits are also the nearest
                    # values from guessed values inserted by the user for the Band Pass.
                    @values=&SearchValues($spectra_C[0],\@spline_spectra,\@spline_intensity);
                    $i=$values[2];

                    @spline_spectra=splice(@spline_spectra,$i);
                    @spline_intensity=splice(@spline_intensity,$i);

                    @values=&SearchValues($spectra_C[-1],\@spline_spectra,\@spline_intensity);
                    $i=$values[2];                    
                    $i++;          
                    splice(@spline_spectra,$i);
                    splice(@spline_intensity,$i);
                    
                                        
                    $area_C=0;
                    $spline_spectra_ant = $spline_spectra[0];
                    $spline_intensity_ant = $spline_intensity[0];
                    for ($i=1;$i<=$#spline_spectra;$i++){
                         @x_sp=("$spline_spectra_ant","$spline_spectra[$i]");
                         @y_sp=("$spline_intensity_ant","$spline_intensity[$i]");
                         @a_b_sp=&LeastSquaresMethod(\@x_sp,\@y_sp);
                         $a=$a_b_sp[0];
                         $b=$a_b_sp[1];
                         $trash=$a_b_C[2];
                         $area_C += &ComputeArea($a,$b,$spline_spectra_ant,$spline_spectra[$i]);
                         $spline_spectra_ant = $spline_spectra[$i];
                         $spline_intensity_ant = $spline_intensity[$i];
                    }

                    splice(@l_or_sp);
                    @l_or_sp=("spline","spline\_spectra\_temp","spline\_intensity\_temp","lambdas\_sp","intensity\_C");
                    
               }
                              
               ######################################################
               # Calculating area below line
               # $area_A -> area below spectral line
               # $area_C -> area below continuum line
               
               $area_A = 0;
               $spectra_C_ant = $spectra_C[0];
               $intensity_C_ant = $intensity_C[0];
               for ($i=1;$i<=$#spectra_C;$i++){
                    @x_C=("$spectra_C_ant","$spectra_C[$i]");
                    @y_C=("$intensity_C_ant","$intensity_C[$i]");
                    @a_b_C=&LeastSquaresMethod(\@x_C,\@y_C);
                    $a=$a_b_C[0];
                    $b=$a_b_C[1];
                    $trash=$a_b_C[2];
                    $area_A += &ComputeArea($a,$b,$x_C[0],$x_C[1]);                    
                    $spectra_C_ant=$spectra_C[$i];
                    $intensity_C_ant=$intensity_C[$i];
               }

               # Calculating Equivalent Width
               $eqw=(($area_C-$area_A)/($area_C))*($spectra_C[-1]-$spectra_C[0]);

               # The errors is beeing computed using Eq 7 of Vollmann & Eversberg 
               # (Astron. Nachr., 2006, 327, 9). Note we have assumed that the 
               # ratio between the areas under the continuum and line are to 
               # the normalized one.
               
               $FC_F=$area_C/$area_A;
# Old way erros was definede               
#               $error=(sqrt(1+($FC_F))*((($spectra_C[-1]-$spectra_C[0])-$eqw)/($SN_ratio*sqrt($NpixFc))));

               $error=(sqrt(1+($FC_F))*((($spectra_C[-1]-$spectra_C[0])-$eqw)/($SN_ratio)));

               push (@lick_indices,"$eqw $error");
               
               # Calculating Flux
               $flux=($area_A-$area_C);
               push (@line_flux,$flux);

               # Calculating Line Center and Continuum at Line Center
               $mean_lambda=($lim_inf+$lim_sup)/2;
               @values=&SearchValues($mean_lambda,\@{$l_or_sp[0]."\_spectra"},\@{$l_or_sp[0]."\_intensity"});
               $mean_lambda=sprintf "%.2f",$values[0];
               $mean_intensity=sprintf "%.3f",$values[1];
               

               # Printing data files to be used by gnuplot
               $temp = $filename1;
               $temp =~ /\/.+\/(.+)/;
               $temp = $+;                    
               
               open FILE, ">".$l_or_sp[0]."\_cont\_".$ion_name."\_".$mean_lambda."\_".$temp;
                    $i=0;
                    for(@{$l_or_sp[1]}){
                         print FILE $_," ",${$l_or_sp[2]}[$i],"\n";
                         $i++;
                    }
               close FILE;

               splice(@spline_spectra_temp);
               splice(@spline_intensity_temp);               

               # Preparing commands for gnuplot
               &PrepareGnuplot($l_or_sp[3],$l_or_sp[4]); 
               
               
               # Preparing final output file
               $output_line=sprintf "%17s  %12.2f  %10.2f %10.2f  %13.3f\n",$ion_name,$mean_lambda,$eqw,$error,$mean_intensity;
               push(@output_1,$output_line);     

               
          }
          elsif( $_ !~ /^\#/ ){
               
               # Preparing to compute arithmetic mean
               s/\s//g;
               /^(.+)\|.+\|/; $name=$+;
               /\|(.+)\|/; $line=$+;
               @colors_limits=split(/\,/,$line);
               splice(@colors);
               splice(@mean_colors);
               push(@colors,"\n\# $name\n   ");
               push(@mean_colors," ");
               for (@colors_limits){
                    
                    # Extracting range values
                    /(.+)\:(.+)/;
                    $color_1=$1;
                    $color_2=$2;
                    $sum_intensity=0;

                    @values=&SearchValues($color_1,\@spectra,\@intensity);
                    $lim1=$values[2];
                    
                    @values=&SearchValues($color_2,\@spectra,\@intensity);
                    $lim2=$values[2];
                    
                    # Computes arithmethic mean and store in $temp
                    for( $i=$lim1;$i<=$lim2;$i++ ){
                         $sum_intensity += $intensity[$i];                         
                    }
                    
                    $N=$lim2-$lim1+1;
                    
                    $temp=$sum_intensity/$N;
                    
                    for( $ik=$lim1;$ik<=$lim2;$ik++ ){
                         $sum_RMS += ($intensity[$ik]-$temp)**2;;                         
                    }
                    $temperrColor=(sqrt($sum_RMS/($N-1)));
                    
                    $printing=sprintf "    %4.0f\:%4.0f  %8.2f  %8.2f\n",$color_1,$color_2,$temp,$temperrColor;
                    push(@output_2,$printing);
               }
          }
          
     }
     if($cont_plot!=9){
          &Gnuplot($filename1);
          $cont_plot=0;
     }

     $tempr = $filename1;
     $tempr =~ /\/.+\/(.+)/;
     $tempr = $+;                    


     $temprr = $filename2;
     $temprr =~ /\/.+\/(.+)/;
     $temprr = $+;                    

     chdir "$dir/PACCE_output/Equiv_Widths";
     
     open OUT,">"."EqW_".$tempr;

          print OUT "#" x 80;
          printf OUT "\n\#%79s\n","\#";
          $spaces=49-length($tempr);
          $spaces2=49-length($temprr);          
          printf OUT "\#%27s %.60s %$spaces.s","Spectrum:",$tempr,"\#";
          print OUT "\#\n";
          printf OUT "\#%27s %.60s %$spaces2.s","Index Definitions:",$temprr,"\#";
          print OUT "\#\n";
          printf OUT "\#%79s\n","\#";

          print OUT "#" x 80;
          printf OUT "\n\#%72s%7s\n","","\#";
          printf OUT "\#%16s"," Ion";
          printf OUT "%16s","Line Center";
          printf OUT "%11s","EW (A)";
          printf OUT "%11s","err (A)";          
          printf OUT "%20s","Cont. Line Center\n";          
          printf OUT "\#%16s","----";
          printf OUT "%16s","-----------";
          printf OUT "%11s","-------";
          printf OUT "%11s","-------";
          printf OUT "%21s\n","-------------------","\#";
          
          print OUT @output_1,"\#\n\#\n";
          print OUT "#" x 80;
          printf OUT "\n\#%79s\n","\#";
          printf OUT "\#%56s%23s\n","Mean Continuum at selected intervals","\#";
          printf OUT "\#%79s\n","\#";
          print OUT "\#" x 80,"\n";
          print OUT "\#    Interval      Mean      Error\n";
          print OUT "\#               Continuum\n";
          print OUT @output_2,"\n\#\n\#";

     close OUT;

     splice (@output_1);
     splice (@output_2);
     splice(@commands_body);

     $repeat=1;
     while ( $repeat > 0 ){
          print "\n\n\tRepeat for other spectra? (y\/n) ";
          $temp=<STDIN>;
          chomp($temp);
          if ( $temp =~ /^n$/i ){
               $repeat=0;
               $again=0;
          }
          elsif( $temp =~ /^y$/i ){
               $repeat=0;
               $again++;
          }
     }
     print "\n";
          
     splice(@spectra);
     splice(@intensity);
     
}


## Appendix
     
&LapsedTime;
     
##################################################################
#                    Subroutines
##################################################################


#Prepare the spectrum inserting the lambda Flambda values os the index defintions into the spectrum

sub prepareSpectrum{

     my @tab = @{$_[0]};
     my @spec = @{$_[1]};
     my @preparedSpec;
        my @linedefs;
     my $line;
     my $refPoint;
     my $flx;
     my $lamb;
     my @inicial;
     my @final;
        my @sortedIndices;
     my @newSpecLamb;
     my @newSpecFlux;
     my @newPoints;
     my $tr;

     foreach $line (@tab){

          if($line !~ /^\s*\#/){
               $line =~ /\((.+)\,(.+)\)/;
               push(@linedefs,"$1","$2");
          }
     }


     foreach $refPoint (@linedefs){
          $ref=$cnt=$ref1=$ii=0;
          for($ii=0;$ii<=$#spec;$ii++){($lamb, $flx) = split(" ",$spec[$ii]);
               if($lamb >= $refPoint){
               $ref1=$ii;
               $cnt=$cnt+1;
               }
          }
          $ref=$ref1-$cnt;
          push (@inicial,$spec[$ref]);
          push (@final,$spec[$ref+1]);
          for($inicial[0]) {($x_inf, $y_inf) = split(" ",$inicial[0]);}
          for($final[0]) {($x_sup, $y_sup) = split(" ",$final[0]);}

          $a=($y_sup - $y_inf)/($x_sup - $x_inf);
          $b=($y_inf)-($a*$x_inf);
          $FluxRefPoint=$a*$refPoint+$b;

          push(@newPoints, "$refPoint  $FluxRefPoint \n");
          splice(@inicial);
          splice(@final);
     }


     @newSpectra=(@spec, @newPoints);


     for($nn=0;$nn<=$#newSpectra;$nn++){
     ($ll, $ff) = split(" ",$newSpectra[$nn]);
     push(@newSpecLamb, "$ll ");
     push(@newSpecFlux, "$ff ");
     }
       
     @sortedIndices = sort{ $newSpecLamb[ $a ] <=> $newSpecLamb[ $b ] } 0 .. $#newSpecLamb;

     @newSpecLamb = @newSpecLamb[ @sortedIndices ];
     @newSpecFlux = @newSpecFlux[ @sortedIndices ];

     for($jj=0;$jj<=$#newSpecLamb;$jj++){
          if($newSpecLamb[$jj]!=$newSpecLamb[$jj+1]){
               push(@preparedSpec, "$newSpecLamb[$jj]  $newSpecFlux[$jj]\n");
          }
     }
#     print @preparedSpec;
     return @preparedSpec;
}
##########################################################################


# Adjusts a line for a given set of points     
sub LeastSquaresMethod {

     # Associates external variables with internal
     my @X=@{$_[0]};
     my @Y=@{$_[1]};     
     my $total= scalar(@X); # -> total entries of array @X 
     # Least Squares Method applied to @X,@Y 
     # to find the best line equation
     my $i=0;
     my $j=0;
     my $soma_Xi=0;
     my $soma_Xi2=0;
     my $soma_Yi=0;
     my $soma_XiYi=0;
     my $temp_RMS=0;
     my $N=$total;
     for( $i=0;$i<=$total;$i++){ 
          $soma_Xi+=$X[$i]; 
          $soma_Xi2+=$X[$i]**2; 
          $soma_Yi+=$Y[$i];
          $soma_XiYi+=$X[$i]*$Y[$i];
     }
     #media     
     my $media = $soma_Yi/$total;
     for( $j=0;$j<$total;$j++){ 
     $temp_RMS+=($Y[$j] - $media)**2;
     }
        $rms = (sqrt($temp_RMS/($total)));
     # Inverse Matrix Determinant
     my $det=($N*$soma_Xi2)-($soma_Xi**2);
        # Signal to noise of the line     
     if($rms eq 0){$rms=1e-15;}
     my $SN_ratio=$media/$rms;
     
     # line equation: y = a.x + b
     my $a=(1/$det)*(($N*$soma_XiYi)-($soma_Xi*$soma_Yi));
     my $b=(1/$det)*(-($soma_Xi*$soma_XiYi)+($soma_Xi2*$soma_Yi));

     my @result;     
     $result[0]=sprintf "%.10f",$a;
     $result[1]=sprintf "%.10f",$b;
     $result[2]=sprintf "%.10f",$SN_ratio;
     $result[3]=$total;     
     $result[4]=$rms;

     return @result;
}

##################################################################

sub ComputeArea {
     my $a,$b,$lim_inf,$lim_sup;
    
     my ($a,$b,$lim_inf,$lim_sup)=@_;

     my $base1=$a*$lim_inf+$b;
     my $base2=$a*$lim_sup+$b;
     my $h=$lim_sup-$lim_inf;

     my $area=(($base1+$base2)*$h)/2;
     $area=sprintf "%.10f\n",$area;
     

     return $area;
     
}

##################################################################

sub NearestPoint {
     my $p1=$_[0];
     my $ref=$_[1];
     my $p2=$_[2];
     if ( ($p2-$ref) >= ($ref-$p1) ){
          return $p1;
     }
     else{return $p2;}     
}

##########################################################

sub MidPoint{
     for($j=0;$j<=$#spectra;$j++){
     ($lamb,$flx) = split(" ",$spectra[$j]);
     
          if(($lamb >=$cont_L_1) && ($lamb <=$cont_L_2)) {
          $cntL=$cntL+1;
          $media_L=$flx+$media_L;
          }
             
          if(($lamb >=$cont_R_1) && ($lamb <=$cont_R_2)) {
          $cntR=$cntR+1;
          $media_R=$flx+$media_R;
          }

     
     
     }
     
     $media_L=$media_L/$cntL;
     $media_R=$media_R/$cntR;

}

##################################################################

sub SearchValues {
     my $i=0;
     my $lambda=$_[0];
     my $spectra_array_name=$_[1];     
     my $intensity_array_name=$_[2];
     while ( ${$spectra_array_name}[$i] < $lambda ){$i++;}
     my $k=$i-1;
     my @result;
     $result[0]=&NearestPoint(${$spectra_array_name}[$k],$lambda,${$spectra_array_name}[$i]);
     if ( $result[0] == ${$spectra_array_name}[$k]){
          $result[1]=${$intensity_array_name}[$k];
          $result[2]=$k;
     }
     else{
          $result[1]=${$intensity_array_name}[$i];
          $result[2]=$i;
     }
     return @result;
}
##################################################################

sub LapsedTime {
     $begin = $^T;
     $now = time();
     $lapse=$now - $begin;     
     printf "\n%70s","#" x 66;
     print "\n\n\tTotal elapsed Time: ",$lapse," s\n\n";
}

##################################################################

sub PrepareGnuplot {
     
     $cont_plot++;
     
     my $array1=$_[0];
     my $array2=$_[1];

     my $X_i=${$array1}[0];
     my $X_f=${$array1}[-1];     
     
     my $IonName=$ion_name;
     $IonName =~ ( s/\_/ /g );
     
     my $print1="\#";
     if ($cont_plot == 8){$print1="";} 
     my $print2="\#";
     if($cont_plot == 4){$print2="";}
     my $print3="\#";
     my $title="";
     if($cont_plot == 2){
          $print3="";
          $title="Band Pass ".($fig*9+1)." to ".(($fig+1)*9);          
     }
     
     my @positions_X=("","0.08","0.39","0.7","0.08","0.39","0.7","0.08","0.39","0.7");
     my @positions_Y=("","0.7","0.7","0.7","0.39","0.39","0.39","0.08","0.08","0.08");
     my @Y=&SmallestAndGreatest($array2);
     
     
     $filename1 =~ /\/.+\/(.+)/;
     my $filename=$+;
     
     @commands_body=(@commands_body,
          "\#set key on left width 5.0 height 2.0\n",
          "\n",
          "##----------------\n",
          "##     Plot $cont_plot\n",
          "##----------------\n",
          "set origin $positions_X[$cont_plot],$positions_Y[$cont_plot]\n",
          "unset label\n",
          "unset title\n",
          "unset xlabel\n",
          "unset ylabel\n",
          "unset xrange\n",
          "unset yrange\n",
          "unset arrow\n",
          "set arrow from $lambda_1,graph 0 to $lambda_1,graph 1 nohead linestyle 4\n",
          "set arrow from $lambda_2,graph 0 to $lambda_2,graph 1 nohead linestyle 4\n",
          "$print1 set xlabel \"Lambdas\" font \"Helvetica,10\" offset \"+0,-0.1\"\n",
          "$print2 set ylabel \"Intensities\" font \"Helvetica,10\" offset \"+0.9,+0\"\n",       
          "$print3 set title \"$filename - $title\" font \"Helvetica,14\" \n",
          "\#set format x \" \"\n",
          "set xrange [$X_i:$X_f]\n",
          "#set key on\n",
              "plot \"$filename1\" title \"$IonName\" with lines, ",
              "  \"$dir/PACCE_output/continuum_data_points/$l_or_sp[0]\_cont\_$ion_name\_$mean_lambda\_$filename\" title \"Continuum\" with lines lt 3\n",
              "\n"
     );

     if ($cont_plot==9){          
          &Gnuplot($filename);
          $cont_plot=0;
          splice(@commands_body);
          
     }

}


##################################################################

sub Gnuplot {

     my $filename=$_[0];
     $fig++;
     $filename =~ s/\/.+\///;
     $filename =~ s/\..{2,3}$/_$fig/;
     
     if($cont_plot < 9){
          for(@commands_body){
               if(/^\sset\stitle/){
                    my $temp1=$fig*9;
                    my $temp2=($fig-1)*9+$cont_plot;
                    s/to\s$temp1/to $temp2/;
               }
          }
     }
     
     for (my $i=0;$i<=$show;$i++){
            my @commands_header=(
                 "####################################\n",
               "#          Setting output          #\n",
               "####################################\n",
               "reset\n",
               "\n",
               "set terminal postscript eps color \"Helvetica\" 8\n",
               "set output \"$dir/PACCE_output/figures/$filename.eps\"\n",
               "\n",
               "$save set terminal x11 font \"Helvetica,8\" title \"PACCE \(gnuplot\)\"\n",
               "$save set output \"STDOUT\"\n",
               "\n",
               "####################################\n",
               "#    Setting General Parameters    #\n",
               "####################################\n",
               "\n",
               "set encoding iso_8859_1\n",
               "set tics in\n",
               "set xtics rotate\n",
               "set size 1,1\n",
               "set origin 0,0\n",
               "set bmargin 0\n",
               "set tmargin 0\n",
               "set lmargin 0\n",
               "set rmargin 0\n",
               "set style line 4 linetype 2 linewidth 2\n",
               "set style line 2 linetype 3 linewidth 1\n",
               "\n",
               "set multiplot\n",
               "\n",
               "set size 0.26,0.25\n"
           );
           
          my @commands_tail=(
               "unset multiplot\n",
               "reset"
          );
          
          my @commands=(@commands_header,@commands_body,@commands_tail);
#  Aqui o padrão de busca !~ não era adequado, pois ele só é executado na negação . Substitui ele por =~.
          if ($save =~ "\#" ){
               open OUT,">$dir/PACCE_output/gnuplot_commands/gnuplot_commands_$filename.txt";
                    for (@commands){print OUT $_;}
               close OUT;
          }     

          open2 (\*OUTPUT,\*INPUT,"gnuplot -persist") or
            die ("\n\n\t\*\*\* You must have \"gnuplot\" installed! \*\*\*\n\n");
                 foreach $command (@commands){
                    print INPUT $command;
               }
          close INPUT;
          close OUTPUT;
          $save = "\#";          
     }
     if($show==1){$save = "";}
}

##################################################################

sub SmallestAndGreatest {     
     my @array = sort {$a <=> $b} @{$_[0]};     
     my @smallest_and_greatest=("$array[0]","$array[-1]");     
     return @smallest_and_greatest;
}


##################################################################

sub DefineInitialTexts {

     my $help="
     #########################################################
     PACCE - Flux and Equivalent Widths Program
     ----   ----------------------------------
     $version

     This program computes equivalent widths and fluxes for spectral lines.

     Usage:       perl $0 [options]

     Options:
        -t , --textmode       use textmode, especially if you want to make
                          a batch file with:
                            + spectra filename,
                            + table filename,
                            + y/n (option for images exhibition), and
                            + y/n (option for process repetition).

        -q , --quiet       don't show initial messages.

        -h , --help       show this help.
     #########################################################\n";


     my $greeting="
     #################################################################

          PACCE - Flux and Equivalent Widths Program
          
     #################################################################

          You need the following programs:

               - Perl;
               - Perl's \"Math::Derivative\" package;
               - Perl's \"Math::Spline\" package;
               - Gnuplot (www.gnuplot.info)

          Further details: http://arxiv.org/abs/1105.0318

     #################################################################\n";

     return $help,$greeting;
}
##################################################################

