#!/usr/bin/perl
#				Gattari & Manini	v.1.3 - Jun. 2003
# convolve a x-sorted data (from STDIN) 'x f(x)'
# with a Gaussian of width '$sigma'.
# NOTE:    "$comment"-ed lines are passed through.

$comment = "#";

# defaults values:
$sigma=10; # change with -s or -l
$fatt=9;   # integ. between y-$fatt*$sigma and y+$fatt*$sigma ; change with -f
$pt_per_sigma=10;  # y-step   $delta_y = $sigma/$pt_per_sigma ; change with -p

while($p=shift(@ARGV)){
#tutti i parametri della linea di comando che non cominciano col "-" vengono considerati nomi di file di configurazione, a meno che non sono l'argomento di qualche opzione che inizia col meno.
    if($p=~/^[^\-]/){
	if(-e $p){
	    push(@file_di_input,$p);
	}else{
	    print STDERR "File $p inesistemte!!!\n";
	}
	next;
    }
# [-s sigma] assegna la larghezza sigma della gaussiana con cui viene fatta la convoluzione
    if($p eq "-s"){
	$sigma=shift(@ARGV);
	if($sigma=~/^\-/){
	    die "ERRORE Command line: -s deve avere un argomento!\n";
	}	
	$sigma+=0.0;
	if($sigma<=0){
	    die "ERRORE Command line: sigma non puo' essere negativo\n";
	}
    }
# [-l larghezza] assegna la semilarghezza a mezza altezza (HWHM) della gaussian
    elsif($p eq "-l"){
	$sigma=shift(@ARGV);
	if($sigma=~/^\-/){
	    die "ERRORE Command line: -l deve avere un argomento!\n";
	}	
	$sigma+=0.0;
	if($sigma<=0){
	    die "ERRORE Command line: sigma non puo' essere negativo\n";
	}
	$sigma*=0.424660900144009521;	#	sigma=larghezza/(sqrt(2 log 2))
    }
# [-f fatt] setta $fatt che e' il numero di sigma entro cui viene integrato.
# senza argomenti viene settata a 0 e quindi viene fatto l'integrale su tutto. 
    elsif($p eq "-f"){
	if($#ARGV==-1){
	    $fatt=0.0;
	}else{
	    $fatt=shift(@ARGV);
	    if($fatt=~/^\-/){
		unshift(@ARGV,$fatt);
		$fatt=0.0;
	    }else{
		$fatt+=0.0;
		if($fatt<=0){
		    die "ERRORE Command line: fatt non puo' essere negativo\n";
		}
	    }
	}
    }
# [-p punti] e' il numero di punti per sigma che vengono tracciati. Se  sigma=10 e punti=5 
# allora verra' calcolato un valore di F ogni delta_y=2. 
    elsif($p eq "-p"){
	$pt_per_sigma=shift(@ARGV);
	if($pt_per_sigma=~/^\-/){
	    die "ERRORE Command line: -p deve avere un argomento!\n";
	}	
	$pt_per_sigma+=0;
	if($pt_per_sigma<=0){
	    die "ERRORE Command line: punti_per_sigma non puo' essere negativo\n";
	}
    }
# [-y y_min:y_max oppure -y y] setta l'intervallo in cui viene calcolata la convoluzione a 'y_min','y_max' o a '-y','y'
    elsif($p eq "-y"){
	@yy=split(/:/,shift(@ARGV));
	if($#yy==0){
	    $y_max=shift(@yy)+0.0;
	    $y_min=-$y_max;
	}elsif($#yy==1){
	    $y_min=shift(@yy)+0.0;
	    $y_max=shift(@yy)+0.0;
	}else{
	    die "Valore dell'intervallo di energia errato!\n";
	}
	($y_max>$y_min) || die "y_max <= y_min\n";
    }
    else{
	die "ERRORE: $p non e' un opzione valida\n";
    }
}

#se non ci sono file di input aspetta i dati dallo STDIN '-'
if($#file_di_input==-1){
    push(@file_di_input,"-");
}

$fattsigma=$fatt*$sigma;
$sigmared=$sigma/(($fatt+3-abs($fatt-3))/2); # the () means min($fatt,3)
$xtmp=1.e230;

@x,@fx; # store data in here
foreach $nome_file (@file_di_input){
    open(INPUT,$nome_file) || die "Impossibile aprire il file $nome_file\n";
    foreach $riga (<INPUT>){
	if( ($riga=~/^$comment/ || $riga=~/^\s+$/)){
	    print $riga;     # let original comments go through
	}else{
	    $riga=~s/^\s+//;
	    @dati = split(/\s+/,$riga);
#  fill with extra interpolated data if actual data are scarce:
	    while(($xtmp+=$sigmared)<$dati[0]){
		$interpf=($dati[1]-$oldf)/($dati[0]-$oldx)*($xtmp-$oldx)+$oldf;
		push(@x,$xtmp);
		push(@fx,$interpf);
	    }
	    push(@x,$dati[0]);
	    push(@fx,$dati[1]);
	    $oldx=$dati[0];
	    $xtmp=$oldx;
	    $oldf=$dati[1];
	}
    }
    close(INPUT);
}

$y_min=$x[0]-$fattsigma;
$y_max=$x[$#x]+$fattsigma;

print "$comment  Spettro convoluto con una gaussiana di larghezza\n";
print "$comment  sigma: $sigma\n";

$delta_y = $sigma/$pt_per_sigma;

#faccio la convoluzione e printo Fx in stdout
for($y=$y_min;$y<$y_max;$y+=$delta_y){
    while($x[$i_min]<($y-$fattsigma)){
	$i_min++;
    }
    $Fy=0;
#    for($j=$i_min;$j<$#x && $x[$j]<($y+$fattsigma);$j++){
#	$Fy+=gaussiana($y-$x[$j],$sigma)*$fx[$j]*($x[$j+1]-$x[$j]); 
#    }
    $funcold=gaussiana($y-$x[$i_min],$sigma)*$fx[$i_min];
    for($j=$i_min+1;$j<=$#x && $x[$j]<($y+$fattsigma);$j++){
	$func=gaussiana($y-$x[$j],$sigma)*$fx[$j];
	$Fy+=($funcold+$func)*($x[$j]-$x[$j-1]); 
	$funcold=$func;
    }
    $Fy*=0.5;
    printf "%25.16g  %25.16g\n", $y, $Fy;
}



sub gaussiana{	# args: x,sigma
    local($result);
#    local($PI=3.14159265358979);
    local($sqr2PI=2.5066282746310005024);
    $result=$_[0]/$_[1];
    $result=exp(-$result*$result/2);
    $result/=$sqr2PI*$_[1];
}