Revise paste a2b49998df8be

###################################
# THIS PROGRAMM IS DOCUMENTED WITH INLINE PODs
# There are different comments for users and programmers
# Extract pod with
# USERs:       sed '/=PROG/,/=cut/d' p.pl | sed 's/=USER/=/g' | pod2html
# PROGRAMMERs: sed '/=USER/,/=cut/d' p.pl | sed 's/=PROG/=/g' | pod2html
###################################

=PROGhead1 NAME model.pm

model - containing objects to represent a rayinvr model used by p.pl

=cut

{
########################################################################
package model;
use strict;
use warnings;
use Carp;
use lib '/opt/pray/lib';
use Carp qw(cluck);
use feature qw(switch say);
use POSIX ('ceil', 'floor');
use Data::Dumper;
use Graphics::ColorUtils;
use Math::Round;
use File::Copy;

use layer;
use floatingReflector;
use station;
use codes;

#use Tk::Animation;
#use Tk::Splashscreen;

=PROGhead1 C<model>

Draws rayinvr-model for seismic refraction data, keeps data storage, do
administration

=PROGhead2 Tomo2D

=begin html

'model' can also be used for viewing data from tomo2D. Reading tomo-data
is started with button 'Tomo2D' from main PRay GUI which invokes 
<code>
model->tomo [\&model::tomo, $model, "tomorays", \"1"]
</code>

=end html

=PROGhead2 Model keys

contours 
 d2 
 km1 
 debug 
 profilelength 
 stationpositions 
 raycode 
 balloon 
 layers 
 d1 
 tmin 
 km2 
 tmax 
 stations 
 vredstate 
 layer 
 vgrid 
 zmax 
 tomotimes 
 ytscale 
 vred 
 foo 
 nodes 
 box 
 phasecodes r
 ays 
 PicksCal
 readStr
 time
 t1
 t2
 PicksMan
 zmin
 drawnStations
 canvasheigth
 xscale
 canvaswidth
 yscale
 ShowRays
 version
 vnodes
 xmin
 space
 drawnPhases
 xmax
 config
 tomoGrid
 tomorays
 phasecolors

DEBUGING

All station objects are stored in the hash $model->{stations}
x keys(%{$model->{stations}}) with position-km as hashkey.
e.g.:

DB<17> x keys(%{$model->{stations}})                                                                                   
 0  250.000
 1  89.300
 2  175.000
 3  25.000
 4  225.000
 5  50.000
 6  275.000
 7  200.000
 8  100.000
 9  125.000
 10  150.000

DB<6> x $model->getStation(127) Gives you the stationreference

=cut

our ($tree, $debug, $dev, $verbose, $quiet);  # Debugging and verbosity variables

sub new {

=PROGhead2 C<model::new(%args)>

A new model is born, blessed and issued with configuration values. Only argements,
which are also defined in %defaults are allowed.

my $model = new model('space'=> $cns, "time" => $lzd,
			"balloon" => $balloon,
			"zmin" => $config{zmin}, "zmax" => $config{zmax}, 
			"yscale" => $yscale, "ytscale" => $ytscale, "xscale" => $xscale,
			);

In Future it shall be:

my $model = new model(
			["reverseStation" => (1/0)],	stationfile => statxz,
			"working dir" => ... , tmin, tmax
			"zmin" => $config{zmin}, "zmax" => $config{zmax}, 
			);

=cut

#print "model::new()\n";
	my ($class, %args) = @_;
	my $model = bless ({}, $class);
	$model->{rays}=undef;
	#$model->init(%args);
    
    # old $model->init function:
    # Define default parameters for a model object
    my %defaults = ( "space"		=> undef,	# canvas for model
					 "time"			=> undef,	# canvas for ttdiagramm
                     "statusline"    => undef,  # GUI status line 
                     "mainwindow"   => undef,   # MAIN WINDOW
                     "splash"       => undef,
                     "icons"        => undef,
					 "balloon"		=> undef,	# window, needed for balloons
                     "init"         => -1,      # Is this model initialized? -1 no, 0 yes, 1 currently initializing
					 "layer"  		=> 0,		# ??? obsolete??
					 "layers"		=> undef,	# array for storing layer objects
				     "profilelength"=>  undef,
				     "box"			=> undef,
				     "foo"	  		=> "bar",
				     "zmin" 	=> undef,
				     "zmax" 	=> undef,
				     "xscale"		=> undef,
				     "yscale"		=> undef,
				     "ytscale"		=> undef,
				     "vred"			=> 3,
				     "vredstate"	=> 1,
				     "canvaswidth" 	=> 1600,	# Width of cavas for time and model
					 "canvasheigth" => 400,		# Heigth of cavas for time and model
					 "tmin"			=> 0,
					 "tmax"			=> 12,
					 "stations"		=> undef,	# refernce hash with links to stations
					 "stationpositions" => undef, # hash with names and positions of stations
                     "stationkm" => undef,        # hash with positions and names of stations 
					 "drawnStations" => [],
                     "phasecolors"		=> { "00" => "darkgrey",
										 "12" => "red", 
										 "21" => "green",
										 "22" => "blue",
										 "31" => "orange",
										 "32" => "black",
										 "41" => "yellow",
										 "42" => "springgreen",
										 "51" => "green",
										 "52" => "blue",
										 "61" => "orange",
										 "62" => "black",
										 "71" => "yellow",
										 "72" => "springgreen",                                        
									 },
                    # colors for phases and layer boundarys.
					 "nodes"		=> 1, 		# Show and edit nodes
                     "vnodes"		=> 0,       # Show and edit vnodes
                     "annotvnodes"	=> 1,
					 "drawnPhases"	=> undef,	# which phases should be drawn
					 "drawnRays"	=> undef,	# which rays should be drawn
					 "PicksMan"		=> 1,		# Show manual picks?
					 "PicksCal"		=> 1,		# Show calculated picks?
					 "vgrid"		=> 0,
					 "config"		=> undef,	# Reference to configurationhash
                     "rin"          => undef,   # Reference to rin
					 "contours"		=> 0,
                     "contourcolor" => 0,
                     "blocks"		=> 0,
					 "debug"		=> undef,
					 "ShowRays" 	=> 1,		# Show rays? Usefull when too many raypathes block view on model
					 "phasecodes"   => { 21 => 2.1},	#conversion of phasecodes between ZP and RI
					 "raycode"  => { 2.1 => 21},	#conversion of phasecodes between ZP and RI
					 "tomorays" => 0,	# TEMPORARY
					 "tomoGrid" => 0,	# Inital value for displaying tomogrid
					 "tomotimes" => 0, # Flag for read tomo2D traveltime data
					 "version" => '',  # Save model version
                     "glueNodes" => 1, # Moving pinched nodes together?
                     "xz" => 0,         # Show xz overlay files?
                     "xt" => 0,         # Show xt overlay files?
                     "igmasCoords" => {},   # Hash to store igmas coordinates  (see layer->getIgmasCoordintes()
                     "densities" => 0,  # Annotate vnodes with densities after Barton
                     "resolution" => 0,
                     "codes" => undef, # Store codes object
    );
	
	foreach (keys(%defaults)){
		#print "Setting default for $_ -> $defaults{$_}\n";
		$model->{$_}=$defaults{$_};
	}

foreach (keys(%args)){
		say "initModel:  $_  $args{$_}" if $debug; 
		croak "Unknown attribute <$_> for station initialization" unless exists $defaults{$_};
		$model->{$_} = $args{$_};
	}

#foreach (keys(%{$model->{config}})) {
		#say "initModel: apply config $_ $model->{config}->{$_}";
		#$model->{$_} = $model->{config}->{$_} if exists $defaults{$_};
	#}
    
    $model->{vred} = $model->{config}->{vred};
    
	$model->{box} = [0, 0, $model->{canvaswidth}, $model->{canvasheigth}];	# Size of model and runtimediagram in px

# Setup scale parameters
	$model->{km1} = $model->{config}->{xmin};
	$model->{km2} = $model->{config}->{xmax} ;
	$model->{d1} = $model->{config}->{zmin};
	$model->{d2} = $model->{config}->{zmax};
	$model->{t1} = $model->{config}->{tmin};
	$model->{t2} = $model->{config}->{tmax};
	
	$model->{xmin} = $model->{config}->{xmin};
	$model->{xmax} = $model->{config}->{xmax};
    $model->{profilelength} = $model->{xmax} - $model->{xmin};
	
	my $cns = $model->{space};

#say "Exit initalization of your model";
    #print Dumper \$model->{color};
    
    # Setup debugging variables
    $tree  = $model->{debug}->{tree};
    $debug = $model->{debug}->{debug};
    $dev   = $model->{debug}->{dev};
    $verbose = $model->{debug}->{verbose};
    $quiet = $model->{debug}->{quiet};

return $model;

}
sub init {

=PROGhead2 C<model::init(%args)>

Creates startup-splash (if configured) and reads model data. Currently 
started from C<p.pl>. Cannot be run by model->new because stations need 
to be defined in the model.

=cut

printf "(T) %s(@_)\n", commons::whoami() if $tree;

my ($model, %args) = @_;

my $splash;
    my $mw;
    
    # Display splash screen if switched on in configuration file
    if ( $model->{config}->{splash} == 1) {
        print "Display splash screen\n";
        $model->_set("init", 1);   
        $mw = $model->_get("mainwindow");
        
        # Make splash screen
        $splash = $mw->Splashscreen;
    
        $mw->withdraw;
        $splash->Splash;                     # show Splashscreen
        
        $splash->Label(-text => 'Building your Model ...', -bg => 'blue')->
            pack(qw/-fill both -expand 1/);
        
            
        my $animate;
        if (@ARGV) {
            $animate = $mw->Animation;
            foreach (@ARGV) {
                $animate->add_frame($mw->Photo(-file => $_));
            }
        } else {
            #my $gif89 = Tk->findINC('anim.gif');
            my $icons = $model->_get("icons");
            my $gif89 = "$icons/Praysplash.gif";
            #my $gif89 = '/projects/nam2011/bin/pray/icons/Praysplash.gif';
            print "Take gif $gif89\n";
            $animate = $mw->Animation(-format => 'gif', -file => $gif89);
        }
        $animate->set_image(0);
        $animate->start_animation;
    
        
        $splash->Label(-text => "Button me", 
            #-image => $mw->Photo(-file => "$icons/reload.gif")
            -image => $animate
            )-> pack(qw/-fill both -expand 1/);
        
        $model->_set("splash" => $splash);
        $mw->update;
    }
    
    ## READ DATA
    $model->read("vin", "rays", "times");
    $model->_set("init", 0);
    $model->drawAxes();
    
    # Draw stations
    foreach (keys(%{$model->{stations}})){
        if  (ref $model->{stations}{$_} eq "station") {
            $model->{stations}{$_}->drawStation;
        } else {
            print "WARNING !! There's an unblessed station in your model!! Call the priest !!\n";
        }
    }
    
    # Draw user defined xz if configured
    if ($model->{config}->{xz}) {
        foreach my $file (@{$model->{config}->{xz}}) {
            $model->_readxzfile("$file");
        }
    }
    
    # Draw user defined xz if configured
    if ($model->{config}->{xt}) {
        foreach my $file (@{$model->{config}->{xt}}) {
            $model->_readxzfile("$file", 'xt');
        }
    }
    
    if ( $args{splash} == 1) {
        $splash->Destroy;                    # tear down Splashscreen
        $mw->deiconify;                      # show main window
    }
}

sub _get{

=PROGhead2 C<model::_get($key, ..)>

function returns any value stored in model.
if a key is not found, undef is returned
THIS FUNCTION IS INTENDED FOR INTERNAL USE
external is C<model::get>

arguments are a list of keys to retrieve

=cut

#printf "(T) %s(@_)\n", commons::whoami() if $tree;

my ($model, @args) = @_;
	return undef unless @args;

my ($key, $value);
	my @return;
	
	foreach $key (@args){
				
		if (exists $model->{$key}){
			$value = $model->{$key};
			#print "Get ".$key." = ".$value."\n";
		} else {
			print "ERROR: no such key as $key in model!!\n";
			$value = undef;
		}

push @return, $value;
	}
	
	#D Return scalar, if only one argument 
	#D and referece to array if more
	if ( @return == 1) {
		return $value;
	} else {
		return \@return;
	}
} #D model->_get END

sub _set{

=PROGhead2 C<< model::_set( $key => $value, .. ) >>

model->_set can set a reference to data for any existing key in
a model
THIS FUNCTION IS INTENDED FOR INTERNAL USE
external is set

=cut

printf "(T) %s(@_)\n", commons::whoami() if $tree;

my ($model, %args) = @_;
	my ($key, $value);
		
	while (($key, $value) = each(%args)){
		
		if (exists $model->{$key}){
			#print "Set ".$key.", ".$value."\n";	
			$model->{$key} = $value;
		} else {
			print "ERROR: no such key as $key in model!!\n";
            die;
		}
	}	
}

sub status {

=PROGhead2 C<model::status()>

Prints number of layers, average, minimum and maximum depth and velocities, number of stations with phases.

$model->status;

=cut

printf "(T) %s(@_)\n", commons::whoami() if $tree;

say "--------------------------";
	say "This is your status-report";
	
	my ($model, %args) = @_; 
	
    my $s =" B  min  max  mean variance stddev\n"; # Return string
    
    my $n = keys(%{$model->{stations}});
	
    
    say "Your model has <$n> stations\n";
	#foreach (sort {$b <=> $a}(keys(%{$model->{stations}}))) {		
		#$model->{stations}->{$_}->statusstation;
	#}
	my @layers = @{$model->{layers}};
	say "Your model has <".@layers."> layers";
	for (my $i=0; $i<=$#layers; $i++){
		#say "Status layer $i";
		$s .= $model->{layers}[$i]->status(%args);
	}
	say "--------------------------";
    print "Modelparameter : ";
    my @par = keys (%$model);
	print "@par\n";
	#print "ISA: \n";
	#$model->DESTROY;\
    return $s;
}

sub printStatusMessage{

=PROGhead2 C<model::printStatusMessage("\nmsg")>

Prints msg to GUI. Same like main prog in p.pl

=cut

printf "(T) %s(@_)\n", commons::whoami() if $tree;

my ($model, $msg) = @_; 
   
    my $stline = $model->{"statusline"};
    $stline->insert ('end', "$msg");
    $stline->see('end');
    $model->_get("mainwindow")->update();
}

sub getClosestStation{

=PROGhead2 C<model::getClosestStation($km)>

Returns position of closest station for given km. Currently only used when 
C<model::tomoReadRays>.

=cut

printf "(T) %s(@_)\n", commons::whoami() if $tree;

my ($model, @args) = @_;
	# Used if you want to calculate moveout velocity.
	# coding not complete
	my $x = $args[0];
	my @positions = sort{$a <=> $b}(keys(%{$model->{stations}}));
	
	my $dx;
	my $min =  900;
	my $st;
	foreach (@positions) {
		$dx = abs($_ - $x);
		print "pos: $_, dx $dx, min $min\n" if ( $model->debug());
		if ($dx > $min) {
			print "Closest Position:dx $dx min $min pos $st\n" if ( $model->debug());
			last;
		} else {
			#print "Closest Position: $_\n";
			$min = $dx;
			$st = $_;
			print "Closest Position: $min, $dx, $st\n" if ( $model->debug());
		}
		
	}
	print "Returning $st\n" if ( $model->debug());
	return $st;
}

sub addStation{

=PROGhead2 C<model::addStation()>

Initializes new station with C<< new station(%args, "model" => $model) >>
and fills information in model-object:

# Storage place for stations: $model->stations{position}
 $model->{stations}{$st->{position}} = $st;
 
 # Information to find station
 # name <-> position
 $model->{stationpositions}->{$st->{name}} = $st->{position};
 $model->{stationkm}->{$st->{position}} = $st->{name};

Stations are stored at hash key C<< $model->{stations}{$st->{position}} = $st >> !!
Station position is an integer with three numbers behind comma.

Check stations in model in perls debuging console:
 x keys(%{$model->{stations}})

=cut

printf "(T) %s(@_)\n", commons::whoami() if $tree;

my ($model, %args) = @_;
	my $st = new station(%args, "model" => $model);
	#say "at ".$st->{position};
    
    # Storing station!!
	$model->{stations}{$st->{position}} = $st;
    
    # Add information about station. Enables finding the station with 
    # stationname or stationposition
	$model->{stationpositions}->{$st->{name}} = $st->{position};
    $model->{stationkm}->{$st->{position}} = $st->{name};
}
sub edit{

=PROGhead2 C<model::edit()>

used for all kind of model editing like adding, deleting, moving nodes
C<< $model->edit("op" => "edit", "tags" => \@tags, "value" => [$km, $v]); >>

=over

=item 'editPhase'

# For SPrange selection:  tags = [phase, ["tagsofPick1"], ["tagsofPick2"]
		# For changing all picks: tags = [phase, st]
        
=item 'allV'

=item other

calls

layer->edit(%args);

=item return message

=back

=cut

printf "(T) %s(@_)\n", commons::whoami() if $tree;

my ($model, %args) = @_; 
    
    my @tags = @{$args{"tags"}};
    my $op = $args{"op"};
    my $msg = ""; # Return message to model
    
    #print "model::edit tags: @tags\n";
    #print Dumper(\%args);
    
    my $type = $tags[0];
    print "My type is $type\n" if $debug;
    if ($op eq "editPhase") {
        $model->editPhase(%args);
    
    }
    
    if ($op =~ m/.*Shot/) {
        #Find out, which station was clicked       
        my $station = $model->getStation($tags[2]);
        $station->set(%args);
    }
    
    #####
    # THIS HANDLES ALL OTHER OPERATIONS !!!
    if ($type ne "allV" && $op ne "editPhase" && $op !~ m/.*Shot/){
        #print "model::edit($type) args = \n";
        #print Dumper \%args;
    
        my $layer = $tags[1];
        my $node = $tags[2];
        # TODO: Bug: Throws error for floating refl. Seems to work with
        # just one refl, but probably causes error in the future
        $model->getLayer($layer)->edit(%args);
    }
    
    if ($type eq "allV"){
        
        for(my $i = 0; $i < $#{$model->{layers}}; $i++){
            #print "Save v nodes into layer $i\n";
            my $km = $args{value}->[0]->[$i];
            my $vu = $args{value}->[1]->[$i];
            my $vl = $args{value}->[2]->[$i];
            my $vupar = $args{value}->[3]->[$i];
            my $vlpar = $args{value}->[4]->[$i];			
            
            #print "Editging km @$vl\n";
            # Tags are not necessary, but $layer->edit needs some values there
            $model->{layers}->[$i]
                ->edit( "op" => "editAllV", "tags" => ["TYPE", "LAYER", "NODE"], "value" => [$km, $vu, $vl, $vupar, $vlpar]);
        }
    }
    
    $model->order;
}

sub editAllParDerivs {

=PROGhead2 C<model::editAllParDerivs()>

Set all partial derivatives to given value. This sub is used, when calculating
the resolution for a model. So far no other use.

=cut

printf "(T) %s(@_)\n", commons::whoami() if $tree;

my ($model, $switch) = @_;
    
    my $op;
    $op = 'setPar' if ($switch == 1) ;
    $op = 'unsetPar' if ($switch == 0);
    print "model:: editAllParDerivs: $op\n" if $debug;
    
    my @coords = (0,0); # Fake, not needed
    foreach my $layer (@{$model->{layers}}) {
        my @tags = ('LAYER', "B$layer->{number}");
        $layer->edit("tags", \@tags, "value", \@coords, "op" , $op); 
    }
    $model->order();
}

sub editPhase {

=PROGhead2 C<model::editPhase()>

Edit phasecodes in zp-Files. Called with right click on a phase in the 
traveltime canvas.

=cut

printf "(T) %s(@_)\n", commons::whoami() if $tree;

my ($model, %args) = @_;
    my @tags = @{$args{"tags"}};
 
    my $mode = $args{mode};
    (my $oldphase = $args{value}[0]) =~ s/\.//;
    (my $newphase = $args{value}[1]) =~ s/\.//;
    
    # For SPrange selection:  tags = [phase, ["tagsofPick1"], ["tagsofPick2"]
    # For changing all picks: tags = [phase, st]
    print Dumper(\@tags) if $debug;
    
    # Environment
    my $config = $model->{config};
    my $zpdir = $config->{zpdir};
    my ($file, $st, $sp1, $sp2);
    
    # If only selected picks should change, SPs have to be found for them
    if (@tags == 3){	# If there are not 3 elements in tags array, then there's no selection of SPrange
        # Get variables for first pick
        (my $st1 = $tags[1][1]) =~ s/RAYS//;
        (my $ph1 = $tags[1][2]) =~ s/Ph//;
        my ($pkm1, $t1, $unc) = split (/ /, $tags[1][3]);
        $pkm1 =~ s/km//;
        $t1  =~ s/t//;
        $unc =~ s/unc//;
        $ph1 =~ s/\.//;	# Phases in Pickfile are without decimal ### TODO CHANGE THIS TO USE RIN MATCHES

# Calculate offset of this pick
        my $off1 = sprintf ("%.3f", $pkm1 - $model->{stationpositions}->{$st1});
        print "Pickoffset is $off1\n" if $debug;
        
        # Get variables for second pick
        (my $st2 = $tags[2][1]) =~ s/RAYS//;
        (my $ph2 = $tags[2][2]) =~ s/Ph//;
        my ($pkm2, $t2, $unc2) = split (/ /, $tags[2][3]);
        $pkm2 =~ s/km//;
        $t2  =~ s/t//;
        $unc2 =~ s/unc//;
        $ph2 =~ s/\.//;	# Phases in Pickfile are without decimal ### TODO CHANGE THIS TO USE RIN MATCHES

# Calculate offset of this pick
        my $off2 = sprintf ("%.3f", $pkm2 - $model->{stationpositions}->{$st2});
        print "Pickoffset is $off2\n" if $debug;

$file = $model->getStation($st1)->{zpfile};
        $file =~ s/head/unc/;
        
        # Find picks in unc file to get shotpoints
        my $line;
        open (FILE ,"<$zpdir/$file") or die "Can't open myfile: $_!\n";
        while ($line = <FILE>) {

if ($line =~ /$off1\s+$ph1/i) { 
                my @line = split ( /\s+/ , $line);
                $sp1 = $line[1];
                print $line };

if ($line =~ /$off2\s+$ph2/i) { 
                my @line = split ( /\s+/ , $line);
                $sp2 = $line[1];
                print $line };
        }
        print "SP=$sp1, SP = $sp2\n" if $debug;
        
        if (!defined $sp1 || !defined $sp2) {
            print "COULDN'T FIND BOTH SHOTPOINTS\n";
            return undef;}
    
        # Sort shotpoints so sp2 is larger than sp1
        if ($sp2 < $sp1 ) {
            print "Swap sp1 $sp1 with sp2 $sp2\n" if $debug;
            my $t = $sp2;
            $sp2 = $sp1;
            $sp1 = $t;
            print "Now: sp1 $sp1, sp2 $sp2\n" if $debug;}
        
        #print "My stationnumber is $st1, ph $ph1 km $pkm1, t $t1, unc $unc\nOpen file $file\n";
        $st = $st1;
    } # Now SPs for picks are found
    else {
        $st = $tags[1];
        print "CHange all picks for station $st\n" if $debug;
    }
    
    $file = $model->getStation($st)->{zpfile};
    $file =~ s/head/pick/;
    #$file = "100st135.h.rloc.pickbak";
    my ($sp, $t) = commons::readPicks("$zpdir/$file");
    die "No file $file found\n" if (!defined $sp );

if ($sp1) {
        print "CALL editPhase($sp, $t,$sp, $t, 'sp1' => $sp1  'sp2' => $sp2,".
              "'oldphase' => $oldphase, 'newphase' => $newphase, 'mode' => $mode \n";
        ($sp, $t) = commons::editPhase($sp, $t, 'sp1' => $sp1, 'sp2' => $sp2,
            "oldphase" => $oldphase, "newphase" => $newphase, "mode" => $mode);
    } else {
        print "CALL editPhase($sp, $t, 'oldphase' $oldphase 'newphase' $newphase 'mode' $mode)\n";			
        ($sp, $t) = commons::editPhase($sp, $t, 
            'oldphase' => $oldphase, 'newphase' => $newphase, 'mode' => $mode);	
    }
    
    my $fileout = "$zpdir/$file";
    #my $fileout = "$zpdir/100st135.h.rloc.pick" ;
    commons::writePicks($sp, $t, $fileout) if ($sp);
    
    # Clean edited picks
    print "model::edit() delete 'Ph$oldphase'\n";
    $model->getStation($st)->emptyTimes( keys => ['txin'], phases => [$oldphase]);
    $model->{time}->delete("Ph$oldphase");
    
		return;

}

sub getLayer {

=PROGhead2 C<model::getLayer($number)>

Returns the layer-object for $layernumber. Number is the same as in v.in

=cut

printf "(T) %s(@_)\n", commons::whoami() if $tree;

my ($model, $number) = @_;
    
    # Remove leading letters (B 1 -> 1)
    $number =~ s/[a-zA-Z]//g;
    
    # Floats do not have any number left
    #print "layer::getLayer -> number = $number\n" if $debug;
    
    #if ($number) {
        # Get layer-object
        my $l = $model->{layers}[$number-1];
        return $l;
    #} else {
        #return 0;
    #}
}

sub getStation {

=PROGhead2 C<model::getStation($number)>

Returns station object which matches C<$number> either as stationname or 
stationposition. Currently it's a bit of luck, which one is used. Look at 
position first but if a stationnumber accidently is the same as the position
of another, it gets confused. (This is a potential BUG, TODO)

=cut

printf "(T) %s(@_)\n", commons::whoami() if $tree;

# returns station-object
	# give either stationname or stationpositions
	my ($model, $number) = @_;
	#print "model::getStation($number)\n" if $debug;
	
	# Return station reference
	my $st;
	
	# Is number a station name? Then get profilekm for this station
	if (exists $model->{stationpositions}->{$number}) {
        #print "Input numer is a stationposition\n" if $debug;
		$number = $model->{stationpositions}->{$number};
	}	
	
	# Is number a stationkilometer? Then return the station object
	if (exists $model->{stations}->{$number}) {
        #print "Input number is a stationname\n" if $debug;
		 $st = $model->{stations}->{$number};
	}	
    
    # Check if station is still undefined
    unless ($st) {
        print "WARNING !!! Couldn't find a station for >$number<.\n";
    }
    
	# No check, if both might be right. It shouldn't anyway, because
	# stationkms have three decimals and stationnames don't
	#print "Return station $st->{name}\n";
	return $st;
}

sub getPhase {

=PROGhead2 C<model::getPhase(phase)>

Get alist of all station having a special phase picked

=cut

printf "(T) %s(@_)\n", commons::whoami() if $tree;

my ($model, @args) = @_;
    my $ph = $args[0];
    my @stations = sort { $a <=> $b } (keys(%{$model->{stations}}));
    
    my $s = ""; # Return string
    
    foreach my $st ( @stations ) {
        my $n = 0;
        
        if ( exists $model->getStation($st)->{txin}->{'1.000'}->{$ph} ) {
         $n += @{$model->getStation($st)->{txin}->{'1.000'}->{$ph}};
        }
        
        if ( exists $model->getStation($st)->{txin}->{'-1.000'}->{$ph} ) {
         $n += @{$model->getStation($st)->{txin}->{'-1.000'}->{$ph}};
        }
        
        # Only return string for stations with this phase
        if ( $n > 0 ) {
            my $name = $model->getStation($st)->{name};
            #print "Get station $name\n";
            $s .= sprintf("%3d: Ph %3d %4d picks\n", $name, $ph, $n);
        }
    }
    
    return $s;
}

sub getDensity {
#printf "(T) %s(@_)\n", commons::whoami() if $tree;

=PROGhead2 C<$model->getDensity (v)>

Returns density for given velocity C<v>. Conversion uses values from 
Barton(1986) "The relationship between seismic velocity and density in the continental crust - a useful constraint?"

=cut

my $model = shift;
    my $v = shift;
    
    $v = sprintf("%.1f",nearest(.1,$v));
    
    #my $conv = "ludwig"; 
    my $conv = "barton"; 
    #my $conv = $model->getConfig('densityconversion');

my $rho;
    
    if ( $conv eq "barton" ) {
        # Densities after Barton(1986) 
        # "The relationship between seismic velocity and density in the continental crust - a useful constraint?"
        my %densities = qw( 0.0 0.0 1.5	1.47 1.6	1.66 1.7	1.73 1.8	1.80 1.9	1.86 
            2.0	1.92 2.1	1.98 2.2	2.01 2.3	2.03 2.4	2.06 
            2.5	2.09 2.6	2.11 2.7	2.13 2.8	2.15 2.9	2.18 
            3.0	2.21 3.1	2.23 3.2	2.24 3.3	2.26 3.4	2.28 
            3.5	2.30 3.6	2.32 3.7	2.34 3.8	2.36 3.9	2.38 
            4.0	2.39 4.1	2.41 4.2	2.43 4.3	2.44 4.4	2.46 
            4.5	2.48 4.6	2.50 4.7	2.52 4.8	2.53 4.9	2.55 
            5.0	2.57 5.1	2.59 5.2	2.61 5.3	2.62 5.4	2.64 
            5.5	2.66 5.6	2.68 5.7	2.70 5.8	2.72 5.9	2.74 
            6.0	2.77 6.1	2.80 6.2	2.83 6.3	2.85 6.4	2.87 
            6.5	2.90 6.6	2.93 6.7	2.95 6.8	2.98 6.9	3.01 
            7.0	3.04 7.1	3.07 7.2	3.10 7.3	3.13 7.4	3.16 
            7.5	3.19 7.6	3.22 7.7	3.25 7.8	3.28 7.9	3.31 
            8.0	3.34 8.1	3.38 8.2	3.42 );

$rho = $densities{$v};    
    } 
    
    if ( $conv eq "funck" ) {
        # Densities after Ludwig(1970)
        # "Seismic Refraction" in THE SEA
        #a=-0.00283; b=0.0704; c=-0.598; d=2.23; e=-0.7;
        #printf("v= %f, rho= %8.2f\n", $1,a*$1^4+b*$1^3+c*$1^2+d*$1+e)
        
        # Funck 2004 has approximated ludwig/barton
        # But values are quite different!!
        $rho = -0.00283*$v**4 + 0.0704*$v**3 - 0.598*$v**2 + 2.23*$v -0.7;    
    }   
    
    unless ($rho) {
        #print "WARNING!! Cannot find density for v=$v\n";
        if ($v > 8.2 ) {
            $rho = 3.5
        }
        if ($v < 1.5 ) {
            $rho = 1.47;
        }
        #print "Use density rho = $rho\n";
    }
    
    return $rho;
}

sub read {

=PROGhead2 C<model::read('data')>

Routine starts individual model-subroutines for reading data.
Keys can be C<'vin', 'rays', 'times', 'str', 'vgrid'>

=cut

printf "(T) %s(@_)\n", commons::whoami() if $tree;

my ($model, @args) = @_; 
	my $r = ""; # Return message
	
	if (grep $_ eq "vin", @args) {
        # Delete contours, grid
        $model->reset();
		$model->_readVin;
        $model->_readRes();
	}
    
	if (grep $_ eq "fin", @args) {
        $model->reset();
		$model->_readFin;
	}
    
	if (grep $_ eq "rays", @args) {
		$r = $model->_readRays;
		
	}
	
	if (grep $_ eq "times", @args) {
		$model->_readTimes;
	}
	
	if (grep $_ eq "str", @args) {
		# Read Streamerpicks and convert them to depth
		$model->_readStr;
	}
	
	if (grep $_ eq "vgrid", @args){
		$model->_readVgrid;
	}
	return $r;
}

sub _readVin {

=PROGhead2 C<model::_readVin()>

Called by C<< $model->read(vin) >> and reads model file C<'v.in'>.

=cut

printf "(T) %s(@_)\n", commons::whoami() if $tree;

my ($model, %args) = @_;
    my $cns = $model->{space};
    my $xmax = $model->{rin}->{xmax}[0];
    my $xmin = $model->{rin}->{xmin}[0];
    my $file = "v.in";
    say "Reading file v.in" unless $quiet;
    
    open(FILE, $file) or croak "Can't open $file";
    my $l=0;
    my @line;
    my $currentlayer=0;
    my $pre; 	# prefix of line. holds layernumber or shows that line is contiuned
    my $layer; 	# current layer as layer-object
    my $v=0; 	# Flag for depth or velocity in second line of block
                # 0 for depth
                # 1 for last depth reading
                # 2 for velocity upper border
                # 3 for velocity lower border
                
    my @km;		# stores profilekm for more than one line
    my @d;		# stores depth corresponding to @km
    my @par;	# stores partial derivatives 
    
    #my (@uv, @lv);	# stores velocity-Nodes
    
    my $n = 0;		# Number of nodes for a layer. Used as flag for end of layer
    my $type = "depth"; # Type of data for second line in block: depth, upper or lower
                        # velocities
    
    my $missXmin = 0;   # Flag that gets set, is any xmin value for v is missing
    
    # Clean up old data
    $model->{layers} = (); # Remove any old layers (because of DESTROY function of layer-object)
    $model->{space}->delete('VNODE');
    $model->{space}->delete('ANNOTVNODE');
    
	while(<FILE>) {
		chomp;
        next unless length;     # Don't cry if line is empty
		$l += 1;
		($pre, my $tmp) = unpack "A2xA*";
        #print "(D) String $tmp contains of ".(length($tmp)/7)." packages \n" if $debug;
        if ( (length($tmp)/7) =~ m/^\d+.\d+$/ ) {
            print 
                 "\n\n\nWARNING!!! Something is wrong with your v.in file in line $.\n"
                ."I have problems to recognize the blocks. Please do not use\n"
                ."tabs to deliminate the numbers. Use spaces.\n\n\n";
		die;
        }
#		@line = split (' ', $tmp);
#		@line = unpack "A7" , $tmp;
		@line = unpack "(A7)".(length($tmp)/7) , $tmp;
		
        #print "(D) v.in $. type $type: ". join (',', @line)."\n" if $debug;
		
		# Blocks of three lines belong together
		# First line holds layernumber, profilekm
		if ( $l%3 == 1 ) {
			# Line holds kilometer. If one layer has more nodes than fit 
			# into one line, they still get pushed into the same array
			push @km, @line;
			
			# If layernumber has increased, create a new layer
			if ( $pre > $currentlayer ) {
				$currentlayer += 1;
                my $c = $model->{phasecolors}{sprintf('%d2',$pre-1)};
                # Layer color is determined by the reflected phasecolor from upper layer
                unless (defined $c) {
                    #print "(I) There is no specific color defined for the boundary of layer B$pre.\n I'm using black."
                    #."If you want the layer boundary in a different color \n";
                    $c = 'black';}
				$layer = new layer("model" => $model, "number" => $pre, "color" => $c);
				#print "Created a new Layer (B $pre).  My Profile is $xmax km long\n";
				#print "(B $pre) Color: $c\n";
			} 	
			
			# Check if end of profile is reached. km cannot be larger then
			# the whole profile
			if ( $line[-1] == $xmax ) {
				$v += 1;
				$n = @km;
				#print "End of layerkm for $pre, $n elements in km\n";
			} elsif ($line[-1] > $xmax) {
				print "Somethings wrong with your v.in-File. xmax = $xmax, but km = $line[-1]\n"; 
			}
		}
		
		# Second line of block holds depth/velocity for above km
		if ($l%3 == 2 ) {
			push @d, @line;
			#print "@line \n";
		}
	
		# Third line of block switches partial derivatives on/off
		if ($l%3 == 0 ) {
			push @par, @line;
			
            # Is depth array as large as km array?
			if ( $n == @d ) {
				# Amount of nodes in @d is the same as @km. And end of layer is reached
				# This is the last line, as km and d array have the same length
				#print "End of second line in block is reached\n";
				
                # is nodetype == depth?
				if ( $type =~ "depth" ) {
					#print "(D) This is last depth of layer @km, @d\n" if $debug;
                    #print "(D) addNodes vukm => \n@km, depth => \n@d, depthpar => \n@par\n" if $debug;
					$layer->addNodes("km" => \@km, "depth" => \@d, "depthpar" => \@par);
					@km = (); @d  = (); @par = (); $n = 0;	# Reset variables for next layer
					$type = "upVelo";   # Flag velocities for next line
					next;
				}
				## This shouldn't really happen. rayinvr wouldn't work. But just in case...
				## Test if number of km- and depth-nodes fit
				#if (@km != @d){
					#print "Length km: ".@km.", Length d: ".@d."\n";
					#die "Number of nodes for km and depth don't fit!!\n";
				#}
				#next;
			
				if ($type =~ "upVelo") {
					#print "velocity upper boundary\n";
					#print "k @km \n";
					#print "l @line \n";
					#print "d @d\n";
					#print "$currentlayer: readVin Partial derivatives upper layer >@par<\n";
                    #print "(D) addVNodes vukm => \n@km, vu => \n@d, vupar => \n@par\n" if $debug;
                    # Set missXmin-flag if no value for xmin found
                    #if ( $km[0] != $xmin){
                        #$missXmin = 1;
                        #unshift @km, $xmin;
                        #unshift @d,  $d[0];
                        #unshift @par, $par[0];
                        #print "(I) Added vl[km $xmin] nodes to $layer->{number}\n";
                    #}
                    $layer->addVNodes("vukm" => \@km, "vu" => \@d, "vupar" => \@par);
					@km = (); @d = (); @par = (); $n = 0;
					$type = "loVelo";
					next;
				}
				
				if ($type =~ "loVelo") {
					#print "velocity lower boundary \n";
					#print "@km \n";
					#print "@line \n";
					#print "$currentlayer: readVin Partial derivatives lower layer >@par<\n";
                    #print "(D) addVNodes vlkm => \n@km, vl => \n@d, vlpar => \n@par\n" if $debug;

# Set missXmin-flag if no value for xmin found
                    #if ( $km[0] != $xmin){
                        #$missXmin = 1;
                        #unshift @km, $xmin;
                        #unshift @d,  $d[0];
                        #unshift @par, $par[0];
                        #print "(I) Added vl[km $xmin] nodes to $layer->{number}\n";
                    #}
					$layer->addVNodes("vlkm" => \@km, "vl" => \@d, "vlpar" => \@par);
					push @{$model->{layers}}, $layer;
					$type = "depth";    # Next line is for depth again
					@km = (); @d = (); @par = (); $n = 0;
				}
			} else {
				#print "Number of second line elements hasn't reached km n $n d @d\n";
			}

}
		
	} #end of FILE-while
	close(FILE);
	# Adding last layer
	$layer->lastLayer("km" => \@km, "depth" => \@d);
	push @{$model->{layers}}, $layer;
    
    # Report missing xmin values!
    if ( $missXmin == 1 ) {
        my $t = "You are working with a v.in file that has no values for xmin. rayinvr works with this, but vin2xyz has problems. I'm adding nodes for xmin. Please write model to file before running dmplstsqr";
        my $m =   $model->{mainwindow}->Dialog(
                    -title => "Velocity nodes in v.in ",
                    -text => $t,
                    #-width => 75,
                    -wraplength => '6i',
                    -buttons => ['Ok']
                    );
     
        $m->Show;
    }
    
    #die;
    # Name layers:
    my %layernames = ( ' 1' => 'Water', ' 2' => 'Sediment', ' 3' => 'Sediment',
        ' 4' => 'Sediments ?', ' 5' => "Basement", ' 6' => "Upper crust",
        ' 7' => 'Middle crust', ' 8' => 'Lower crust', ' 9' => "Mantle");
    #print Dumper \%layernames;
    
    # Drawing layers
    foreach my $layer (@{$model->{layers}}){
        $layer->drawLayer;
        
        # TODO: change this naming behaviour
        if (exists $layernames{$layer->{number}}) {
            $layer->{name} = $layernames{$layer->{number}};
            #print "Name layer $layer->{number} $layer->{name}\n";
        } else {
            #print "Can't find a name for >$layer->{number}<\n";
        }
    }
    
    $model->writeV2xzyIn(); # Convert v.in to an vin2xyz-readable format
    $model->_readFin();
}

sub _readFin {

=PROGhead2 C<model::_readFin()>

Reads floating reflector file f.in

=cut

printf "(T) %s(@_)\n", commons::whoami() if $tree;

my ($model, %args) = @_;
    my $cns = $model->{space};
    my $xmax = $model->{xmax};
    my $file = "f.in";
    say "Reading file $file";

my $pre;    # prefix of line. holds layernumber or shows that line is contiuned
    my $layer;  # current layer as layer-object
    my $type='amount' ; # Flag linetype
                
    my @km;	    # stores profilekm for more than one line
    my @d;      # stores depth corresponding to @km
    my @par;    # stores partial derivatives 
    
    
    my $n = 0;		# Number of nodes for a layer. Used as flag for end of layer

open(FILE, $file) or (print "Can't open $file", return);
    while (<FILE>) {
        chomp;
        #s/^\s+//;               # no leading white
        #s/\s+$//;               # no trailing white
        
        next unless length;     # Don't cry if line is empty
        #my @line = split;
        
        # First line defines the number of nodes in the floating reflector
        if ( $type eq 'amount' ) {
            my @line = split;
            $n = $line[0];
            $type = 'km';
            $layer = new floatingReflector("model" => $model, 
                "number" => 'Float', "color" => "black");

next;
        }

# Input file: f.in contains the "floating" reflecting boundaries;
        # these are interfaces with no associated velocity discontinuity
        # (hence the word "floating") at which rays can reflect in
        # addition to the layer boundaries; the file format is similar to
        # that for v.in:
        # (1) the number of nodes defining the floating reflector
        # (format: I2)
        # (2) the x-coordinates (km) of the nodes defining the floating
        # reflector (format: 3X, <number of nodes>F7.2)
        # (3) the z-coordinates (km) of the nodes defining the floating
        # reflector (format: 3X, <number of nodes>F7.2)
        # (4) a 0 or 1 for each node listed above depending on whether (1)
        # or not (0) partial derivatives are to be calculated for
        # this particular node (format: 3X, <number of nodes>I7)
        # Lines (1) through (4) are repeated for each floating reflector
        
        my $line = $_;
        my @line = unpack "A3".("A7" x $n);
        #print "line: @line\n" if $debug;
        
        # Remove line prefix
        shift @line;
        if ( $type eq 'km' ) {

if (@line == $n) {
                # Got all nodes I expect for this reflector
                @km = @line;
                $type = 'depth';
                next;
            }
        }
        
        if ( $type eq 'depth' ) {
            @d = @line;
            $type = 'par';
            next;
        }
        
        if ( $type eq 'par' ) {
            @par = @line;
            
            #print "(D) addNodes \n"
                #."vukm     => @km\n"
                #."depth    => @d\n"
                #."depthpar => @par\n" if $debug;
            $layer->addNodes("km" => \@km, "depth" => \@d, "depthpar" => \@par);

@km = (); @d  = (); @par = (); $n = 0;	# Reset variables for next layer
            $type = 'amount';
            $layer->drawLayer;
            push @{$model->{layers}}, $layer;
            next;
        }
        
        
        
    }
    close(FILE);

}

sub _readRays{

=PROGhead2 C<model::_readRays()>

Large routine reading C<r1.out> and C<r2.out>. Before reading any data, old
ray-information is deleted from model-objects. Then phase information about rays
is read from C<r1.out> then data is read from C<r2.out> and connected with phases
from C<r1.out>.

Information from C<r1.out> is stored invariable C<@r1> for later use, when reading
C<r2.out>

=cut

printf "(T) %s(@_)\n", commons::whoami() if $tree;

my ($model, %args) = @_; 
    my $rayblock = 0; 	# Flag. If 1 than first part of file is read and rayinfo begins
    my $r1ray = -1; 	# raynumber in r1.out. should fit $raynumber from r2.out
    my @r1;				# Stores file r1.out, which keeps information about the number
                        # of points in a ray.
    my $phase;
    my $linenumber = 0;
    
    my $error = "";
    my $results = "";   # Store shot, phase and model information
    
    # TRY only deleting rays when the phase was recalculated. Otherwise keep informations
    ## Delete old rays
    # deleteRays
    #foreach (sort {$b <=> $a}(keys(%{$model->{stations}}))) {
        ## Safety maeasure for unblessed stations 
        #if  (ref $model->{stations}{$_} eq "station") {
            #$model->{stations}->{$_}->emptyRays;
        #} else {
            #print "WARNING !! There's an unblessed station in your model!! Call the priest !!\n";
        #}
    #}

####################################################################
    #                            READ R1.OUT                           #
    ####################################################################
	open (RAYS, "r1.out") or do {print "Can't open r1.out\n"; return 0};
	#shot  ray i.angle  f.angle   dist     depth red.time  npts code
	#  -1   1    1.838   -4.103    0.000    1.31  10.327     4   1.2
	#  -1   2    2.865   -5.802    0.000    0.74  10.359     5   1.2
	#
	# If npts is 500, ray tracing wasn't successfull
	# "dist" fits with last lines value for "x" for one ray in r2
	# "i.angle" = 90 - r2(ang2); fist line for one ray
	# "f.angle" = 90 - r2(ang1); last line for one ray
	
	my $status = 0;  # Has rayinvr finished? This status is set to one when 
		# rayinvrs output reached |----, which is followed by a tracing 
		# summary. If you have a bad model and rayinvr can't run through
		# this status will stay at 0.
	my ($sh, $raynr, $i, $f, $dist, $d, $redtime, $npts, $ph);
	$model->printStatusMessage("Reading r1.out..");
    print "Reading r1.out\n" unless $quiet;
	while (<RAYS>){
		$linenumber += 1;
		chomp;
		
        my $line = $_;
		#my @line = unpack("A4A4A9A9A9A8A8A6A6",$_);
		#   0      1      2   3   4     5     6        7      8
		($sh, $raynr, $i, $f, $dist, $d, $redtime, $npts, $ph) =  unpack("A4A4A9A9A9A8A8A6A6",$line);
		$ph =~ s/^\s+|\s+$//g; # remove trailing and leading spaces for phasecode
		
		# Read Phase numbers
		unless (defined $ph) {print "No phasecode defined >$ph<\n"; next;}
		if ($ph =~ /^\d\.?/){ # Only if there's a phasecode at the phasecodeposition

# Special treatment for floating reflectors
            # TODO 
            #if ($ph =~ /^\d\.4/ ) {
                #$ph =~ s/4$/2/;
                #print "A floating reflector has been traced. Change phase to $ph\n";
            
            #}

#print "Push phase $ph\n";
			#			raynr	npts	code/phase
			push @r1, [$raynr, $npts, $ph, $linenumber];
			#print "r1.out: $sh, $raynr, $i, $f, $dist, $d, $redtime, $npts, $ph\n";
		} else {
			print "r1.out: $line\n" unless $quiet;
            $results.="$line\n";
            
            # Find shot kilometer
            $line =~ s/^\s+//;               # no leading white
            my @line = split /\s+/, $line;
            
            # Find station for this km, if first element fits the required
            # format
            if ($line[0] && $line[0] =~ /\d\.\d/) {
                #print "Line array: $line[0] ";
                my $st = $model->getStation($line[0]);
                #print "Station = $st->{name}\n";
                chomp $results;
                $results .= "       st $st->{name}\n";
            }
			
		}
		
		if ($sh eq '|---') {
			#print "End of rays reached. Leaving loop .. \n";
			$status = 1;
			#last;
            #$model->printStatusMessage("\nThe model is talking to you!! ");
		}

if ( grep(/Number of data points used/, $line ) ) {
            $model->printStatusMessage("\n$line");
        }        
        
        if ( grep(/RMS/, $line ) ) {
            $model->printStatusMessage(",  $line");
        }
        
        if (  grep(/Normalized/, $line )) {
            $model->printStatusMessage(",  $line");
        }
        
        # Update Splash screen during start up
        if ( $model->_get("init") == 1 && $linenumber % 200 == 0 ) {
            #print "Update Splash\n";
            $model->_get("mainwindow")->update();
        }
	}
	close(RAYS);
    
    unless ($linenumber > 0 ) {
        print "ERROR: File r1.out is empty!!!\n";
        $model->printStatusMessage("\nERROR: File r1.out is empty!!! Check command line output for rayinvr error message");
        return 1;}
	print "r1.out read and closed\n" unless $quiet;
    
    $model->{results} = $results;   # Store modeled results
    
    # Check if there have been any errors during the rayinvr run
	if ($status == 0) {
        
        # Find last successfully traced station
        my $k = $model->{rin}->{xshot}[abs($sh)-1];
        my $s = $model->{stations}->{$k}->{name};
        print "------------------------------------------------------------------------------------\n";
		print "| WARNING!! RAYINVR HASN'T FINISHED                                                \n";
		print "| RAYS WILL BE MISSING IN DISPLAY                                                  \n";
		print "| LAST SHOT TRACED: SHOT$sh, RAY$raynr, DISTANCE$dist, DEPTH$d, PHASE $ph \n";
		print "|$sh, $raynr, $i, $f, $dist, $d, $redtime, $npts, $ph      \n";
                #4,   4, 79.256,-78.006,229.842,0.00,5.047, 46, 9.1

print "| SHOT$sh = station $s at km $k                                 \n";
		print "------------------------------------------------------------------------------------\n";
        
        $model->printStatusMessage(
            "\nERROR: rayinvr could trace all stations. Traveltimes are not updated."
            ." See command line for more information");
        return;
    }
    ####################################################################
    #                            READ R2.OUT                           #
    ####################################################################
	open (RAYS, "r2.out") or return "Can't open r2.out. No rays can be displayed. Run rayinvr.";
	print "Reading r2.out \n" unless $quiet;
    $model->printStatusMessage("                  .. reading r2.out ..");
	# To be written if you want to display velocity as gradient
	# Block definition in r2.out
	# number of blocks (nblk) is equal number of nodes in c.in
	#layer# 7  nblk=   2 (ivg,x1,x2,z11,z12,z21,z22,s1,b1,s2,b2,vp1,vs1,vp2,vs2,         vp3,vs3,vp4,vs4,c1,c2,...,c11)
	#ivg         1         1
	#x1     0.0000   37.9700
	#x2    37.9700  750.0000
	#z11   10.0000   10.0000
	#z12   10.0000   10.0000
	#z21   20.0000   20.0000
	#z22   20.0000   20.0000
	#s1     0.0000    0.0000
	#b1    10.0000   10.0000
	#s2	    0.0000    0.0000
	#b2    20.0000   20.0000
	#vp1    7.0000    7.0000
	#vs1    4.0415    4.0415
	#vp2    7.0000    7.0000
	#vs2    4.0415    4.0415
	#vp3    7.3000    7.3000
	#vs3    4.2147    4.2147
	#vp4    7.3000    7.3000
	#vs4    4.2147    4.2147
	#-0.477E-05 0.477E-05
	 #0.000E+00 0.000E+00
	 #0.114E+02 0.214E+03
	 #0.477E-06-0.477E-06
	 #0.254E+04 0.477E+05
	 #0.000E+00 0.000E+00
	 #0.380E+03 0.712E+04
	 #0.000E+00 0.000E+00
	 #0.000E+00 0.000E+00
	 #0.181E-03-0.340E-02
	#-0.181E-02 0.340E-01

$linenumber = 0;
	my $currentStation;	my $currentray = 0; my $currentgroup = 1; my $currentphase = 0;
	my $r2ray; my $pointnumber; my $group;
	my $km;
	$i=0;
	my $numLays; 	# number of layers, from line three in r2.out
	my $l = 0; 		# flags blocks of layer x
	my $nblk; 		# number of blocks within a layer
	my $skip = 0; 	# flag to skip one line
	my $unpack = "";
	my $ll = 29; 	# number of lines for layer blocks
					# if layer has more than hundred nodes, block spread
					# over two lines 
	my @blocks; 	# save data in matrix
	my $mod = -1; 	# Needed for layer block lines spreading over two lines in
					# file. It saves, if new line starts on even or odd linenumber
	while (<RAYS>){
		$linenumber += 1;
		chomp;
        
        # Update Splash screen
        if ( $model->_get("init") == 1 && $linenumber % 200 == 0 ) {
            #print "Update Splash\n";
            $model->_get("mainwindow")->update();
        }
        
		##################################################################
		# Reading block information
		##################################################################
        if ($rayblock == 0) {
			# Find the total number of layers in this model
			if ($linenumber == 3) {
				my @line = split /=/;
				$numLays = $line[-1];
                #print "Number of layers: $numLays\n";
				next;
			}
			# Find place with raypoints
			elsif ($_ =~ m/^gr*/) { 
			 	#print "Found start of rays at $linenumber: $_\n";
				$rayblock = 1;
				next;
			} 
			# Find description of modelblocks and set up variables for block 
			# reading
			elsif ($_ =~ m/^layer/ && $l >= 0 ) { 
                my @line = unpack("A6A2A7A4",$_);
				#my @line = split;
				$l = $line[1];
				$i = -1; # gives linenumber within layer-block
				$nblk   = $line[3];
				#print "Layer $l, number of blocks $nblk\n";
	
				#create unpack string from nblks
				$unpack = "";
				my $c = $nblk;
				if ($nblk > 100) {  # POSSIBLE BUG IF MORE THAN 200 HUNDRED
									# NODES ARE USED
					$c = 100;
					
					#$ll = 59;
					$ll = 29;
					$mod = $linenumber%2;
					#print "For this layer $l, nblk $nblk mod: $mod is too much line: $linenumber\n";
				} else {
					$ll = 29;
					$mod = -1;
				};
				for (my $i = 0; $i < $c; $i++){
					$unpack .= "A10";
				}
				#print "Unpack format: $unpack\n";
				
				next;
			}
			# Read Blockinformation for a layer
			elsif ($l > 0) {
				my @line = unpack($unpack,$_);
				#print "$linenumber: $i: @line\n";
				$i += 1; 
				if ($mod == $linenumber%2) {
					#print "Continue, paste line to previous\n >@line[-1]<\n";
					$i-=1; # This line is sequel of previus one
					# remove empty arrayplaces
					@line = grep {$_ ne ''} @line;
					#print "Sauber >@line[-1]<\n";
				} 
				
				######################################
				#if ($l == 2 ) {
					
					#if (defined @{$blocks[$i]}) {print "L2 in array @{$blocks[$i]}\n";}
					#print "L2 $i: @line\n";
				#}
				######################################
				push @{$blocks[$i]}, @line;
				
				if ($i == $ll ) { # if more than hundred nodes, it needs two lines
					
					#print "End of layer $l, blocks $nblk, lines ".@blocks." \n";
					$model->_addBlocks( "layer" => "$l", "blocks" => \@blocks ); 
					@blocks = (); # Clear array for next line
					$mod = -1; # Reset flag for long lines
					
                    #print "End of layer $l\n";
					if ($l == $numLays) { # All layers have been read once
						# more layer information is not neccessary at this 
						# stage of software
						#print "All layerblocks read\n";
						$l = -1;	
					}
				}
				next;
			}
		}
			
		next unless $rayblock == 1;

###################################################################################################
		# Reading raypathes
		#####################################
		# Ray is the same number as ray-column in r1. If ray equals zero, ray doesn't reach
		# the surface.
		# rays are numbered in column "ray". If that number increase, a new ray starts
		# if npt is 500, ray tracing wasn't successfull
		#   gr ray npt   x       z      ang1    ang2    v1     v2  lyr bk id iw
		#   1  0   1  20.285   2.000    0.00  -88.69   0.00   1.58  1  2 -1  1
		#   1  0   2  18.870   2.029  -88.93  -88.93   1.58   1.58  1  1 -1  1	
		
		# Problem occur if you use "split" because two columns can merge if numbers
		# become too large
		# "gr" is raygroup and means a group of rays with the same code, e.g. code 1.2 is group 1
		# code 2.2 is group 2
		my @line;

# rayinvr and xrayinvr have different output formats in r2.out. 
		# --> That seems to be a patch. Don't know where output-format is defined in rayinvr-source
		#if ( $model->{config}->{r2out} =~ m/[0-9]/ ) {
			# print "Config r2out <$model->{config}->{r2out}< matches a number!\n";
            my $pack;
            # Automatically check the r2out value!!
            if ( length $_ == 67 ) {
                $pack = "A2A3A4A8A8A8A8A7A7A3A3A3A3";
            } elsif ( length $_ == 69 ) {
                $pack = "A4A3A4A8A8A8A8A7A7A3A3A3A3";
            }else{
                # If none is correct, just skip this line
                print "Can't understand the format of r2.out. This is a programming error!";
                next;
                #die;
            }
			#@line = unpack("A".$model->{config}->{r2out}."A3A4A8A8A8A8A7A7A3A3A3A3",$_);
            @line = unpack($pack,$_);
		#} else {
			#print "Don't know your configfile-entry for r2.out! Got >$model->{config}->{r2out}<.\n".
			#"It should be the number of digits for first column of r2.out\n";
			#die;
		#}

$r2ray   	 = $line[1]; # equals $r1[0]	
		$pointnumber = $line[2]; # corresponds to (can't be bigger than) $r1[1] 
		$group 		 = $line[0]; # kind of correspondend to raycode in r1
		
		$km = sprintf "%.3f", $line[3];
		#$km =~ s/^\s+//;    # Remove leading whites
        
        # Check the format of group number. If too many rays are traced
        # it changes to '**'
		if ($group !~ m/\d+/){
			print 
                "\n\n\n________________________________________________________________\n"
                ."\nERROR\n"
                #."Gr >$group< is NO digit!\n"
                ."Line $linenumber: @line\n"
                ."Last station read: Station $model->{stations}->{$km}->{name} at "
                ."km $model->{stations}->{$km}->{position}, phase $currentphase\n";
            print "Your traced rays are too big to fit into r2.out. This is"
                ." a rayinvr problem and need to be solved there.\nYou need"
                ." to find a way, to reduce the number of rays (also failed)"
                ." tracings\n"
                ."________________________________________________________________\n\n\n"
                ;
            $model->printStatusMessage("\nYour r2.out is overflowed. Last station $model->{stations}->{$km}->{name} "
                ."km $model->{stations}->{$km}->{position}, phase $currentphase");
			return;
		}

# Second column in r2.out "ray" seems to be zero if raypath is not valid
		# Those rays are not listed in r1 and hence need to be ignored
		# You can have error messages in r2.out. They are sorted out with the $group-test
		# to be a digit.
		if ($rayblock == 1 && $group =~ m/\d+/ && $r2ray != 0) {
                       
            
			# A new ray starts
			if ($r2ray != $currentray || $group != $currentgroup){
				#print "New ray starts\n";
				#print "r2ray != currentray || group != currentgroup)";
				#print "$r2ray != $currentray || $group != $currentgroup)\n";
				#print "line:$. @line\n";

# If it's a new station, delete old rays in that station
                if ( ! $currentStation || $currentStation != $model->{stations}->{$km} ){
                    #print "Delete rays for station $model->{stations}->{$km}->{name} at km $km\n" if $debug;
                    #$model->{stations}->{$km}->emptyRays;
                    $model->getStation($km)->emptyRays;
                }
                #print "station $model->{stations}->{$km}->{name} at km $km\n";
            
				$currentray   = $r2ray;
				$currentgroup = $group;
				$currentStation = $model->{stations}->{$km};
				$r1ray+=1;
				
				$currentphase = $r1[$r1ray][2];
				#print "Phase is $currentphase. Previous phase is $phase\n" if ($phase != $currentphase);
				#$currentStation->emptyRays($currentphase) if ($phase != $currentphase);
				
				$phase = $r1[$r1ray][2];
				#print "r2.out:$linenumber: gr $group ray $r1ray, nbr: $r2ray, cur $currentray,".
				#" ph: $phase , km <$km> <@line\n";
				if (!defined $currentStation) {print "?? Which station? km $km in line $. of r2.out: @line\n"};
				
			}

if (!defined $r1[$r1ray][0] && $status == 0){ # there's no information in r1 for this ray
				# and rayinvr hasn't finished
				print "WARNING !! RAYINVR DIDN'T RUN THROUGH THE WHOLE TRACING\n";
				last;
			}
			
			# Check if r2rays fit:
			if ( $r1[$r1ray][0] != $r2ray){
				print "________________________________________________________________\n".
				      "\nERROR\n".
				  	  "r2rays of r1 and r2 don't fit!! index $r1ray $r1[$r1ray][0] != $line[1]\n".
					  "Station $currentStation->{name}, ray $currentray\nr2 line $linenumber: @line\n".
					  "r1 raynr, npts, ph, linenumber\n".
					  "r1 @{$r1[$r1ray]}\n".			
				      "________________________________________________________________\n\n".
				      "  What to do?\n".
					  "  - Check if rayinvr has finished without problems\n".
                      "      - PRay cannot recognize segmentation faults of rayinvr\n".
					  "      - if rayinvr has'nt finished, check v.in for crossing layers, large gradients, ..\n".
					  "        (use vmodel) and/or go to last working version and delete r2.out\n". 
					  "  - Check if format of r2.out matches your configuration:\n".
					  " 	 - get the number of digits of first column of raypart of r2.out\n".
					  " 	 - p.pl is expecting >$model->{config}->{r2out}< digits\n".
					  " 	 - Use 'r2out=NUMBER' in p.config to match r2.out with p.pl\n".
					  "_________________________________________________________________\n".
					  "\n\n";
                $model->printStatusMessage("\nERROR: rayinvr has problems to trace all data. See command line output for more information");
				close(RAYS);
                return 1;
				#die;
			}
            unless ($currentStation) {
                print "WARNING no station for ray: @line\n";
                $error = "\nERROR: Couldn't find station for some rays";
                next;
            }
            # Use real raycodes as 'phase'. e.g 1.2 instead of 12. In older
            # versions, this has been different and only the tx-code was used
            # to sort in rays
            # $phase = '5.1'
			$currentStation->addRays("phase" => $phase, "raynumber" => $r1ray,
									 "rays" =>  [$km, $line[4]] );
		} # if it is a ray
	} # end of while-loop for reading file
	
	close(RAYS);
	print "All rays read\n" unless $quiet;
    $model->printStatusMessage(".. Rays are read.");
	#die;
    $model->printStatusMessage($error);
}
sub _addBlocks {

=PROGhead2 C<model::_addBlocks()>

Add blocks read from C<v.in> to model

=cut

printf "(T) %s(@_)\n", commons::whoami() if $tree;

# Finds corresponding layer object and adds blocks there
	my ($model, %args) = @_;
	#print "model::addBlocks\n";
	my $layer = $args{layer};
	my $blocks = $args{blocks};
	#print "Layer $layer, blocks $blocks\n";
    if (defined $model->{layers}[$layer-1]) {
        $model->{layers}[$layer-1]->addBlocks($blocks);
    } else {
        warn "Undefined value for layer $layer-1\n";
    }
}

sub _readVgrid {

=PROGhead2 C<model::_readVgrid()>

Creates grid file from v.in, reads grid file and creates little colored 
squares to generate a smooth gradient.

vin2xyz -R$xmin/$xmax/$zmin/$zmax -x$dx -z$dz < v2xyz.in > v.grid

=cut

printf "(T) %s(@_)\n", commons::whoami() if $tree;

my ($model, %args) = @_;
    print "model::_readVgrid\n";
    $model->{space}->delete('VGRID');
    # Create gridfile with rayinvr-skript vin2xyz	
    my $xmin = $model->{xmin};
    my $xmax = $model->{xmax};
    my $zmin = $model->{zmin};
    my $zmax = $model->{zmax};
    my $dx = 1;
    my $dz = 0.5;
    
    my $command = "vin2xyz -R$xmin/$xmax/$zmin/$zmax -x$dx -z$dz < v2xyz.in > v.grid";
    print "$command\n";
    system ($command);
    
    # Read in vgrid
    print "Gridfile created\n";
    
    my $file = "v.grid";
    say "Reading file v.grid";
    my $cns = $model->{space};
    
    open(FILE, $file) or do {$model->{vgrid} = 0; die "Can't open $file"};
    while (<FILE>){
        chomp;
        my ($x, $z, $v) = split;
        my $color = _getColor($v);
        my @scaled = $model->model2screen([$x-$dx/2, $z-$dz/2, $x+$dx/2, $z+$dz/2], "space");
        
        $cns->createRectangle(@scaled, -fill => $color, -width=>0, -tags=>["VGRID", "km $x, d $z, v $v"]);
        
    }
    close (FILE);
    
    unless ($cns->find('withtag', 'VGRID')) {
        $model->{vgrid} = 0;
        $model->printStatusMessage("\nError in creating the velocity grid. Check the command manually:\n$command");
        }

$model->order;
    print "model::_readVgrid() is done\n";
}
sub tomoReadGrid{

=PROGhead2 C<model::tomoReadGrid()>

Reads and displays grid file from Tomo2D

=cut

printf "(T) %s(@_)\n", commons::whoami() if $tree;

my ($model, %args) = @_;
    print "model::tomoReadGrid\n";
    $model->{space}->delete('VGRID');
    my $file = $model->{config}->{tomoGrid};

my $cns = $model->{space};
    my $dx = 1;
    my $dz = 0.5;

# CLEANUP
    $cns->delete('TOMOGRID');
    $cns->delete('TOMOREFL');
    
    print "This is tomogrid trying to open $file\n";
    open(FILE, $file) or (print("Can't open $file\n"), return 0);
    print "This is tomogrid reading $file\n";
    my $line = 0;
    while (<FILE>){
        chomp;
        my ($x, $z, $v) = split;
        
        #print "Reading line $line: $x $z $v\n";
        my $color = _getColor($v);
        my @scaled = $model->model2screen([$x-$dx/2, $z-$dz/2, $x+$dx/2, $z+$dz/2], "space");
        $v = sprintf "%5.2f", $v;
        $cns->createRectangle(@scaled, -fill => $color, -width=>0, -tags=>["TOMOGRID", "km $x, d $z, v $v"]);
        $line++;
    }
    close (FILE);
    $model->{tomoGrid}=1;
    # Read Reflector file
    $model->_readxzfile($model->{config}->{tomoRefl});
    $model->_readxzfile('bathyxz.dat');
    $model->_readxzfile('basement.dat');
    
    return 1;
}

sub _readxzfile {

=PROGhead2 C<model::readxzfile($file [, diagram])>

Reads given file and creates a line in model space with tags C<'XZ'>.
Calls C<< $model->order >>.

=cut

printf "(T) %s(@_)\n", commons::whoami() if $tree;

# Read and plot given file
    my ($model, $file, $diagram) = @_;    
    $file =~ s/^\s+//;               # no leading white
    $file =~ s/\s+$//;               # no trailing white   
    #print "read xz