Programma per calcolare la forma media della valanga
Il programma utilizza due moduli esterni:
- Chart::Graph::Xmgrace (see docs) per creare automaticamente il file xmgrace e quindi il flie di immagine .png
- Math::Spline e Math::Derivative per fare uno spline delle valanghe da mediare
#!/usr/bin/perl -w
# Include modules
use Chart::Graph::Xmgrace qw(xmgrace);
#require Math::Spline;
use Math::Spline qw(spline linsearch binsearch);
use Math::Derivative qw(Derivative2);
foreach $file (@ARGV)
{
print "# Processing $file...";
@xv=();
@yv=();
@lenght=();
$lenghtmin=1e32;
$nava=0;
open(INPFILE,"<$file") || die;
while(<INPFILE>){
next if(/^\n/);
# first point of the avalanche
$nava++;
($x,$avax,,)=split(/,/);
@aval=($avax);
@xaval=($x);
$norm=0.0;
while(<INPFILE>){
last if(/^\n/);
($x,$avax,,)=split(/,/);
$norm+=$avax*($x-$xaval[$#xaval]);
push @xaval, $x;
push @aval, $avax;
}
for($i=0;$i<@aval;$i++){
$aval[$i]/=$norm;
}
push @xv, @xaval;
push @yv, @aval;
$lenghtmin=( $#aval+1) if($#aval+1<$lenghtmin);
push @lenght,$#aval+1;
}
close(INPFILE);
$nbins=$lenghtmin;
if($nbins<10 || $nava==0){
print "no data!\n";
next;
}
# compute average
@xave=@yave=();
$k=0;
for($i=0;$i<$nava;$i++)
{
@y=@x=();
for($j=0;$j<$lenght[$i];$j++){
next if($j>0 && $xv[$k]==$x[$#x]);
push @x, $xv[$k];
push @y, $yv[$k];
$k++;
}
@y2=Derivative2(\@x,\@y);
$index=binsearch(\@x,0.0);
$spline=new Math::Spline(\@x,\@y);
for($j=0;$j<=$nbins;$j++){
$index=linsearch(\@x,$j/$nbins,$index);
$y_interp=spline(\@x,\@y,\@y2,$index,$j/$nbins);
$yave[$j]+= $y_interp
}
}
for($i=0;$i<=$nbins;$i++)
{
$yave[$i]/=$nava;
$xave[$i]=($i/$nbins);
}
$#xave=$nbins;
$#yave=$nbins;
# Write out the average and the theoretic prediction (for instance, semicircle)
$pi=3.14159265358979;
@xth=@yth=();
for($i=0;$i<=$nbins;$i++)
{
$xth[$i]=$i/$nbins;
$yth[$i]=8./$pi*sqrt($xth[$i]*(1.-($xth[$i])));
}
# Create xmgrace plot
xmgrace(
############################## Graph options
{ "title" => "Average Avalanche",
"subtitle" =>"$file",
"type of graph" => "XY graph",
"output type" => "png",
"output file" => "$file.png",
"x-axis label" => "Normalized Velocity",
"y-axis label" => "Rescaled Displacement",
"grace output file" => "$file.agr",
},
############################## First set
[ { "title" => "Avalanches",
"set presentation" => "XY",
"options" => {
"line" => {
"type" => "0",
# "color" => "8",
# "linewidth" => "1",
"linestyle" => "0",
},
"symbol" => {
"symbol type" => "1",
"size" => "0.1",
"color" => "1",
"fill pattern" => "1",
"fill color" => "1",
},
"fill" => {
"type" => "0",
},
},
# "data format" => "file",
"data format" => "columns",
},
############################## Second set
[\@xv,\@yv]],
# "data.dat"],
[ { "title" => "Average",
"set presentation" => "XY",
"options" => {
"line" => {
"type" => "1",
"color" => "8",
"linewidth" => "4",
"linestyle" => "1",
},
"symbol" => {
"symbol type" => "0",
"size" => "0.1",
"color" => "1",
"fill pattern" => "1",
"fill color" => "1",
},
"fill" => {
"type" => "0",
},
},
# "data format" => "file",
"data format" => "columns",
},
[\@xave,\@yave]],
# "avadata.dat"],
############################## Third set
[ { "title" => "Semicircle",
"set presentation" => "XY",
"options" => {
"line" => {
"type" => "1",
"color" => "6",
"linewidth" => "3",
"linestyle" => "1",
},
"symbol" => {
"symbol type" => "0",
"size" => "0.1",
"color" => "1",
"fill pattern" => "1",
"fill color" => "1",
},
"fill" => {
"type" => "0",
},
},
# "data format" => "file",
"data format" => "columns",
},
[\@xth,\@yth]]
);
print " ...done\n";
}
--
TWikiGuest - 21 Sep 2004