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

Topic attachments
I Attachment Action Size Date Who Comment
plpl graph.pl manage 4.0 K 21 Sep 2004 - 16:30 TWikiGuest Compute av. shape of aval. (velocity vs depl.)
elserpm perl-Chart-GRACE-0.95-8.i386.rpm manage 15.3 K 01 Oct 2004 - 14:37 AndreaBaldassarri  
elserpm perl-Math-Derivative-0.01-1.rhfc1.bio.i386.rpm manage 6.9 K 01 Oct 2004 - 14:37 AndreaBaldassarri  
elserpm perl-Math-Spline-0.01-1.rhfc1.bio.i386.rpm manage 9.3 K 01 Oct 2004 - 14:37 AndreaBaldassarri  
Topic revision: r3 - 27 Feb 2006 - 13:56:54 - DmitryNPopov
 
This site is powered by the TWiki collaboration platformCopyright © by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding TWiki? Send feedback