model.pm

###################################
# 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.pm;
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 file >$file<\n";

    my $color;

    ($file, $color) = split(/\s+/, $file,);
    $color = 'blue' unless ($color);

    $color = $model->_checkColor( $color );
    print "Draw file >$file<, color >$color<\n" if $verbose;



    my $cns = $model->{space};
    my $tag = "XZ";
    my $scalediagram = "space";

    if ($diagram && $diagram eq 'xt') {
        $cns = $model->{time};
        $tag = "XT";
        $scalediagram = "time";
    }

    my @refl;

    open(FILE, $file) or print("Can't open $file") , return 0;
        while (<FILE>){
            chomp;
            s/^\s*#.*//;            # whole line is commented
            s/^#.*//;               # whole line is commented
            next unless length;

            my ($x, $z) = split;
            push @refl, $x, $z;
        }
    close (FILE);

    my @scaled = $model->model2screen(\@refl, $scalediagram);
    # Line
    #$cns -> create('line',@scaled, -joinstyle=>"bevel", -fill=> $color, -width=>1,
            #tags => [$tag, "$file"]);
    my $s = .3; # size of circles
    for (my $i = 0; $i < $#scaled; $i += 2) {
        $cns -> create('oval',
            $scaled[$i]-$s, $scaled[$i+1]-$s ,$scaled[$i]+$s, $scaled[$i+1]+$s ,
            -fill=>$color, -outline=>$color, -width=>1,
            tags => [$tag, "$file"]);
    }


    if ($diagram && $diagram eq 'xt') {
        $model->{xt} = 1;
    } else {
        $model->{xz} = 1;
    }
}

sub tomoReadRays{
    # Reading in rays written by Korenaga's Tomo2D
    # Fileformat (ready for gmt plotting):
    # >
    # km d
    # .. ..
    # stkm std
    # >
    # .. next station

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

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

    my $file = $args{'file'};
    my $ph = $args{'ph'};

    print("tomoReadRays needs file >$file< and >$ph< \n"), return 0  unless (defined $file && defined $ph);
    print "model::readTomoRays ($file, $ph)\n";

    ### TODO TODO, only delete phases which shall get overwritten
    foreach (sort {$b <=> $a}(keys(%{$model->{stations}}))) {
        $model->{stations}->{$_}->emptyRays($ph);
    }


    my @data;   # array to store input file
    my $cns = $model->{space};
    my $st;
    my $ray=0;
    my $new = 1;    # Flag new ray
    my $oldkm = -1; # Check if new ray is for same station as last one
    open(FILE, $file) or  print("Can't open $file\n"), return 0;
    while (<FILE>){
        chomp;
        my @line = split;
        if ($new == 1) {
            my $km = $line[0];
            my $stkm = $model->getClosestStation($km);
            $st = $model->getStation($stkm);

            print "Draw tomo ray for station $st->{name} at km $km\n" if $debug;
            $new = 0;
            if ($oldkm != $km){
                $ray = 0;
                $oldkm = $km;
            }
        }

        if ($line[0] eq ">") {
            # This is when a new ray starts or ends
            $new = 1;
            @data = ();
            $ray++;
        } else {
            $st->addRays("phase" => $ph, "raynumber" => $ray,
                         "rays" =>  [$line[0], $line[1]] );
        }

    }
    close (FILE);
    print "Finished reading $file for Tomo2D-rays\n";


}
sub _makeContours {

=PROGhead2 C<< model::_makeContours( "source" => tomo/rayinvr) >>

Program creates contour file for rayinvr or tomo2D. The file name used
for tomo2D it's still hardwired to tomo.grd!!
For rayinvr the v.grd is created by vin2xyz and xyz2grd.

It's called by model->tomo for tomo and model->ordner for rayinvr

=cut

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

    my ($model, %args) = @_;
    print "model::_makeContours\n";

    die "Give datasource for Contours. Either rayinvr or tomo2d" unless (exists $args{source});

    # 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.1;

    my $grid = "v.grd";
    my $xyz = "v.xyz";

    if ($args{source} eq "rayinvr") {
        my $command =
        "vin2xyz -R$xmin/$xmax/$zmin/$zmax -x$dx -z$dz < v2xyz.in > $xyz";
        print "Make xyz: $command\n";
        system ($command);
        $command = $model->getConfig('gmt')." xyz2grd $xyz -R$xmin/$xmax/$zmin/$zmax -I1/0.1 -Gv.grd";
        print "Make grid:  $command\n";
        system ($command);

    }

    if ($args{source} eq "tomo") {
        $xyz = $model->getConfig("tomoGrid");
        $grid = "tomo.grd";
    }

    # -MC: means: M for use one file with multiple segments
    # C is a flag with marks a new segment
    #my $command = "grdcontour $grid -C0.2 -D -MC -JX1/1 > contours.ps";
    my $command =
       # Old gmt command (for gmt4)
       # "grdcontour v.grd -C0.2 -D -MC -JX25/-16 -B50/10:.'Version $model->{version}': ".
       # New configurable gmt-switch
       $model->getConfig('gmt')." grdcontour v.grd -C0.2 -D -MC -JX25/-16 -B50/10:.'Version $model->{version}': ".
       " -R$xmin/$xmax/$zmin/$zmax --HEADER_FONT_SIZE=9 > contours.ps";
    print "Make contours \n$command\n";
    system ($command);

    print "To create a postscript of your contours run:\n";
    print $model->getConfig('gmt')." grdimage v.grd -JX25/10 -R$xmin/$xmax/$zmin/$zmax -V -Ccolors_velocities.cpt > t.ps\n";
    print "model::_makeContours() is done\n";
}

sub _readContours {

=PROGhead2 C<model::_readContours>

Reads and draws contours with tag 'C<CLINE>'. No flag is set, no ordering
done. That's organized by C<order>.

=cut

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

    # Read file "contour" made by gmt programm grdcontour in sub _makeContours
    # and draws them.
    my ($model, %args) = @_;
    my $cns = $model->{space};
    #$model->_makeContours();

    print "_readContours\n";

    my $file = "contour";
    my @cline;  # Array for collecting point of a contour line
    my ($x, $z, $v);
    open(FILE, $file) or do {$model->printStatusMessage("\nCan't open file >$file<. Check the grdcontour command (see command lineoutput)"); return};
    while (<FILE>){
        chomp;

        # Skip first line
        if ($. == 1) {next};

        # If line matches "C.." then it's the start or end of a contourline
        unless ( m/C/) {
            ($x, $z, $v) = split;
            #print "Reading line $.: $x $z $v\n";
            push @cline, $x, $z;
        } else {
            #print "Draw Contour line $_, previous $v\n";
            my @line = split;
            #my $v = $line[1];
            my $color = _getColor($v);
            my $rho = $model->getDensity($v);

            my $w = 4;      # Contourwidth

            if ( $model->{contourcolor} == 0 ) {
                #$color = 'gray20';
                my $g = 110 ;
                $color = sprintf('#%02x%02x%02x', $g,$g,$g); # RGB
                $w = 1;
            }
            #print "Contour color for $v km/s is  $color\n";
            my @scaled = $model->model2screen(\@cline, "space");
            $cns -> create('line',@scaled, -joinstyle=>"bevel", -fill=>$color, -width=>$w, tags => ['CLINE', "$v rho $rho"]);
            @cline = ();
        }
}
    close (FILE);
    print "model::_makeContours() is done\n";
}

sub tomo {

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

C<model::tomo> handles buttons for drawing
C< $key =  tomorays, tomogrid>

Ends with calling C<< model->order() >>.

=cut

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

    my ($model, %args) = @_;
    print "model::tomo (%args)\n";
    print Dumper(\%args) ;

    my $cns = $model->{space};
    my $lzd = $model->{time};

    # ONLY TOMO2D FUNCTIONS LEFT
    # AT THE END: $model->order();
    if (exists $args{"tomorays"}) {


=PROGbegin html

<ol>
<li> Delete old rays <br>
    <code>$model->{stations}->{$_}->emptyRays;</code><br></li>

<li> Get version from file 'tomoversion' and update window title </li>

<li> Read rays and traveltimes <br>
    If 'tomoRaysPg' is defined, then Pg and PmP are read from different
    ray and tx-files <br>

    <code>
        $model->tomoReadRays('file' => $model->getConfig('tomoRaysPmP'),
        'ph' => $model->getConfig('tomoPhasePmP')); <br>
    $model->_readTx('file' => $model->getConfig('tomoTimesPg'),  key => 'txTomo');
    </code></li>

<li> Read tomo grid<br>
    <code>
    $model->tomoReadGrid;
    </code>

</ol>

=end html

=cut

        ## Just set value. Drawing is done by model::order

        #$model->{nodes} = 0; # Flaggin nodes not to show

        ####################
        # 1. Delete old rays
        ### TODO TODO, only delete phases which shall get overwritten
        foreach (sort {$b <=> $a}(keys(%{$model->{stations}}))) {
            $model->{stations}->{$_}->emptyRays;
        }

        # Get tomoversion
        my $version = "Tomo2D";
        my $file = "tomoversion";
        open(FILE, $file) or die "Can't open $file";
        while (<FILE>){
            chomp;
            #my ($x, $z, $v) = split;
            $version = $_;
            print "Version: $version\n";
        }
        close(FILE);
        $model->_set("version", $version);

        # Update windowtitle
        $model->{"mainwindow"}->title("PRay: Model: $version");
        $model->{"mainwindow"}->update;

        # There are two possibilities for input files. Either one file with
        # all rays for phases, or two files with Pg and PmP seperated
        if ( $model->getConfig('tomoPhasePg') ) {
            print "model->getConfig(tomoRaysPg) $model->getConfig('tomoRaysPg')\n";
            $model->tomoReadRays('file' => $model->getConfig('tomoRaysPg'),
                'ph' => $model->{config}->{'tomoPhasePg'});

            $model->tomoReadRays('file' => $model->getConfig('tomoRaysPmP'),
                'ph' => $model->getConfig('tomoPhasePmP'));

            $model->_readTx('file' => $model->getConfig('tomoTimesPg'),  key => 'txTomo');
            #$model->_readTx('file' => $model->getConfig('tomoTimesPmP'), key => 'txTomo');
            $model->{tomotimes} = 1
        } else {

            # One file with both phases. Cannot distinguish raypathes for
            # reflections and refractions. It doesn't matter for picked
            # times
            print "WARNING !! RAYS FOR PG AND PMP ARE UNDISTINGUISHABLE\n";
            if ($model->{config}->{'tomoRays'}) {
                $model->tomoReadRays('file' => $model->getConfig('tomoRays'),
                    'ph' => $model->{config}->{'tomoPhase'});
            }

            if ($model->{config}->{'tomoTimes'}){
                    $model->_readTx('file' => $model->getConfig('tomoTimes'), key => 'txTomo');
            }
            $model->{tomotimes} = 1
        }
        print "Read files\n";

=PROGpod

 model->{tomotimes} = 1, if traveltimes for tomo2D have been read
 No distinction between Pg, PmP or both is made in this flag

 otherwise model->{tomotimes} = 0

=cut

        my $e = $model->tomoReadGrid;
        if ( $e == 1 ) {
            print "Grid read successfully. Deactivate nodes\n";
            $model->{nodes} = 0;
        }
        #$model->{tomorays} = ${$args{'tomorays'}};
    } # end tomorays

    if (exists $args{'tomoGrid'}){
        # Only draw grid, if none is already there
        my $e = 1;
        unless ($cns->find('withtag', 'TOMOGRID')) {
            $e = $model->tomoReadGrid;}

        print "Error code is $e\n";
        # $e is set to zero, if reading the grid fails.
        # Then nodes should still be drawn
        if ( $e == 1 ) {
            print "Grid read successfully. Deactivate nodes\n";
            $model->{tomoGrid} = ${$args{'tomoGrid'}};
            $model->{nodes} = 0;
        }
    }

    if (exists $args{tomoContours}) {
        $model->_makeContours( "source" => "tomo");
        $model->_readContours;
        $model->{tomoContours} = ${$args{'tomoContours'}};
    }

    print "Order yourself!\n" if $debug;
    $model->order;
}

sub order {

=PROGhead2 C<model::order()>

Change the sequence of object on canvases. Does not draw new objects
(apart from contours (should that be changed in future?)). Reads config
of objects from model and raises or lowers the object accordingly. E.g.
if nodes should not be displayed the are lowered, so they lie behind the
model layers and cannot be seen. If they should be shown later, they
don't need to be redrawn

=cut

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

    # Make order of objects according to settings
    # Like raise or lower bounds, nodes, ..
    my ($model) = shift;
    my $cns = $model->{space};
    my $lzd = $model->{time};

    print "model::order\n" if $tree;


    if ($model->{tomoGrid} == 1) {
        $cns->raise('TOMOGRID');
        $cns->raise('TOMOREFL');
    } else {
        $cns->lower('TOMOGRID');
        $cns->lower('TOMOREFL');
    }

    if ($model->{vgrid} == 1) {
        #print "Raise VGRID\n";
        unless ($cns->find('withtag', 'VGRID')) {
            $model->_readVgrid;
        }
        $cns->raise('VGRID');
    } else {
        #print "Lower VGRID\n";
        $cns->lower('VGRID');
    }

    if ($model->{contours} == 1){
        unless ($cns->find('withtag', 'CLINE')) {
            print "There should be contours, but there aren't. So I make them\n";
            $model->_makeContours("source" => 'rayinvr');
            $model->_readContours();
        }
        $cns->raise('CLINE');
    } else {
        $cns->lower('CLINE')
    }

    #print "model->{xz} = $model->{xz}\n";
    if ($model->{xz} == 1){
        $cns->raise('XZ');
    } else {
        $cns->lower('XZ')
    }

    #print "model->{xt} = $model->{xt}\n";
    if ($model->{xt} == 1){
        $lzd->raise('XT');
    } else {
        $lzd->lower('XT')
    }

    if ($model->{blocks} == 1 ){
        $cns->raise('BLOCK');
    } else {
        $cns->lower('BLOCK');
    }

    #if ($model->{tomorays} == 1){
        #unless ($cns->find('withtag', 'TOMORAYS')) {
            #print "There should be tomorays, but there aren't. So I make them\n";
            #$model->tomoReadRays();
        #}
        ##$cns->raise('TOMORAYS');
    #} else {
        ##$cns->lower('TOMORAYS')
    #}

    $cns->raise('AXES');
    #$cns->raise('RAYS');
    if ($model->{ShowRays} == 1) {
        print "order::Raise RAYS\n" if $debug;
        $cns->raise('RAYS');
    } else {
        print "order::Lower RAYS\n" if $debug;
        $cns->lower('RAYS');
    }


    #if ($model->{nodes} == 1 ){
        #$cns->raise('BOUND'); $cns->raise('NODE');
    #} else {
        #$cns->lower('BOUND'); $cns->lower('NODE');
    #}
    $cns->raise('BOUND');
    if ($model->{nodes} == 1 ){
        $cns->raise('NODE');
    } else {
        $cns->lower('NODE');
    }

    if ($model->{vnodes} == 1 ){
        unless ($cns->find('withtag', 'VNODE')) {
            foreach my $layer (@{$model->{layers}}){
                $layer->drawVNodes;
            }
        }
        $cns->raise('VNODE');
    } else {
        $cns->lower('VNODE');
    }

    if ($model->{annotvnodes} == 1 && $model->{vnodes} == 1 ){
        $cns->raise('ANNOTVNODE');
    } else {
        $cns->lower('ANNOTVNODE');
    }

    if ($model->{densities} == 1 && $model->{vnodes} == 1 ){
        $cns->raise('ANNOTRHO');
        $cns->lower('ANNOTVNODE');
    } else {
        $cns->lower('ANNOTRHO');
    }

    if ($model->{resolution} == 1 && $model->{vnodes} == 1 ){
        $cns->raise('ANNOTRES');
        $cns->lower('ANNOTVNODE');
    } else {
        $cns->lower('ANNOTRES');
    }

    $cns->raise('STATION');
    $lzd->raise('PICK');

    if ($model->{PicksMan} == 0) {
        $lzd->lower('txin');
    }
    if ($model->{PicksCal} == 0) {
        $lzd->lower('txout');
    }
}

sub drawSingleStationData {

=PROGhead2 C<model::drawSingleStationData( $station, $value)>

Called by p::bc_drawstation and sets ONE station. It's either added or
removed from array C<< $model->{drawnStations} >>. Calls C<<station::drawSingleStationData(value)>
for drawing.

 - $station->draw(value);
 - Update $model->{drawnStations}
 - $model->order;

=cut

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

    my $model = shift;
    my $arg = shift;
    my $value = shift;
    #print "model::drawSingleStationData($arg) <$value> on/off\n" if ($model->tree());

    my $pos = $model->{stationpositions}->{$arg};
    $model->{stations}{$pos}->draw($value);     # <----- DRAWING

    # Update array $model->{drawnStations}
    if ($value == 1 ){
        push @{$model->{drawnStations}}, $pos;
    } else {
        # Delete station picks, if no rays are there
        my @stations = grep {$_ ne $pos} @{$model->{drawnStations}};
        $model->{drawnStations} = \@stations;
    }
    #print "drawn Stations: @{$model->{drawnStations}}\n";

    $model->order;
}

sub drawPhaseStationList{

=PROGhead2 C<< model::drawPhaseStationList("phases" => [@phaselist], "stations" => [@stationlist]) >>

This sub is ONLY called by drawOOP from p.pl. It delete all drawings and
draws the given lists of phases and rays.

 $model->{drawnStations}=$args{stations}
 $st->draw;
 $model->order;

Also takes the 'redraw' subroutine. If no arguments are given, drawn Phases and stations
are redrawn.

=cut

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

    #say "model::drawPhaseStationList";
    my ($model, %args) = @_;

    # Set new stations and phases
    $model->{drawnPhases}   = $args{phases}   if (exists $args{phases});
    $model->{drawnRays}     = $args{rays}     if (exists $args{rays});
    $model->{drawnStations} = $args{stations} if (exists $args{stations});

    my @positions = @{$model->{drawnStations}}; # Stationpositions to be drawn
    my @phases    = @{$model->{drawnPhases}};   # Phases to be drawn

    # Draw rays and times for given stationlist
    for (my $i = 0; $i <= $#positions; $i++){
        my $st = $model->getStation($positions[$i]);
        #say "Draw rays and times for Station at position >$positions[$i]<, $st->{name}";
        $model->{stations}{$positions[$i]}->draw(1);
    }
    $model->order;
}

sub drawPhase {

=PROGhead2 C<model::drawPhase($ph, $value)>

Called by p::bc_drawPhase when phase/raybutton has changed.


 - Update $model->{drawnRays}
 - Update $model->{drawnPhases}
 - Draws rays and times for all stations in $model->{drawnStations}
 - $model->order

=cut

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

   #print Dumper \@_;
    #print @_;

    my ($model, %args) = @_;
    my $cns = $model->{space};
    my $lzd = $model->{time};

    if (exists $args{ray}) {
        my $rc = $args{ray}[0];
        my $value = $args{ray}[1];
        print "model::drawPhase() set ray $rc to $value\n" if $debug;

        if ($value == 0) {
            print "model::drawPhase() Delete phase $rc from @{$model->{drawnRays}}\n" if $debug;
            $cns->delete("Ph$rc");

            my @raycodes = grep {$_ ne $rc} @{$model->{drawnRays}};
            $model->{drawnRays} = \@raycodes;
        } else {
            push @{$model->{drawnRays}}, $rc;
            say "model::drawPhase() Added $rc to >@{$model->{drawnRays}}<" if $debug;
        }
    } # if exists $args{ray}

    if (exists $args{phase}){
        my $ph = $args{phase}[0];
        my $value = $args{phase}[1];
        print "model::drawPhase() set phase $ph, to $value\n" if $debug;
        if ($value == 0) {
            print "model::drawPhase() Delete phase $ph from @{$model->{drawnPhases}}\n" if $debug;
            $lzd->delete("Ph$ph");

            my @phases = grep {$_ ne $ph} @{$model->{drawnPhases}};
            $model->{drawnPhases} = \@phases;
        } else {
            push @{$model->{drawnPhases}}, $ph;
            say "model::drawPhase() Added phase $ph to >@{$model->{drawnPhases}}<" if $debug;
        }
    } # if exists $args{phase}

    if (defined $model->{drawnStations}){
    # Is phase button pushed before stationbuttons?
        my @positions = @{$model->{drawnStations}};
        for (my $i = 0; $i <= $#positions; $i++){
            say "model::drawPhase() Draw for Station $positions[$i]" if $debug;

            $model->{stations}{$positions[$i]}->draw;
        }
    }

    $model->order;
}

sub drawAxes{

=PROGhead2 C<model::drawAxes()>

Deletes existing axes and draws new axes for model and traveltime diagram

=cut

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

    my $model = shift;
    my $cns = $model->{space};
    my $lzd = $model->{time};
    my $box = $model->{box};

    #say "model::drawAxes()";

    ## Delete existing axes
    $cns->delete('AXES');
    $lzd->delete('AXES');
    $lzd->delete('GRID');

    my $km1 = $model->{km1};
    my $km2 = $model->{km2};
    my $xmin = $model->{xmin};
    my $xmax = $model->{xmax};

    my $d1 = $model->{d1};
    my $d2 = $model->{d2};

    ######
    # PROFILE LENGTH
    ## Draw x-axes
    #print "Draw xaxes\n";
    my @tick = $model->model2screen([$km1, $d2, $km2, $d2], "space");
    $cns -> create('line',
            $tick[0], $tick[1]-2, $tick[2], $tick[3]-2,
            -joinstyle=>"bevel", -fill=>"black", -width=>2, -tags => ["AXES"]);

    #print "Draw ticks\n";
    my $stepx = 50; # default
    my $xdist = $km2-$km1;
    if ($xdist < 10) {
        $stepx = 1;}
    elsif ($xdist < 40) {
        $stepx = 5;}
    elsif ($xdist < 210) {
        $stepx = 10;}
    elsif ($xdist < 500) {
        $stepx = 25;}

    # Start at a nice number like 0,10, or so
    my $a = sprintf "%d",$xmin/$stepx;
    my $x0 = $a*$stepx;

    #print "step $stepx, xdist $xdist, km= $km2 - $km1 start from >$x0<\n";

    ###############################
    ## Labels and ticks for x-axes

    for (my $i=$x0; $i <= $xmax; $i+=$stepx) {
        #print "Mark tick $i, $d2 , from ".($xmin+$xmin%$stepx)."\n";
        @tick = $model->model2screen([$i, $d2], "space");
        $cns -> create('line',
            $tick[0], $tick[1]-5, $tick[0], $tick[1],
            -joinstyle=>"bevel", -fill=>"black", -width=>1, -tags => ["AXES"]);

        $cns -> createText(
            $tick[0], $tick[1]-10,
            -fill=>"black", -text => "$i", justify => "right", -tags => ["AXES"]);

        # Same for TT-diagramm
        $lzd -> create('line',
            $tick[0], $tick[1]-5, $tick[0], $tick[1],
            -joinstyle=>"bevel", -fill=>"black", -width=>1, -tags => ["AXES"]);

        $lzd -> createText(
            $tick[0], $tick[1]-10,
            -fill=>"black", -text => "$i", justify => "right", -tags => ["AXES"]);

    }

    #####################
    # DEPTH
    # Draw y-axes
    my $dy = ($d2-$d1);
    my $stepy = 10; #default
    my $y1 = $d1;
    my $y2 = $d2;

    if ($dy < 10) {
        $stepy = 1;}
    elsif ($dy < 26) {
        $stepy = 5;}

    # Draw Y-axes
    $cns -> create('line',
            $model->model2screen([$km1, $d1, $km1, $d2], "space"),
            -joinstyle=>"bevel", -fill=>"black", -width=>4, -tags => ["AXES"]);

    # Tick and Label Y-AXES
    # Start at a nice point
    $a = sprintf "%d",$model->{zmin}/$stepy;
    my $y0 = $a*$stepy;

    for (my $i=$y0; $i <= $model->{zmax}; $i+=$stepy) {
        @tick = $model->model2screen([$km1, $i, $km1, $i], "space");
        $cns ->
            create('line',
             $tick[0], $tick[1], $tick[2]+5, $tick[3],
            -joinstyle=>"bevel", -fill=>"black", -width=>1, -tags => ["AXES"]);
        $cns ->
            createText(
            $tick[0]+15, $tick[1],
            -text => $i,
            -fill=>"black", -tags => ["AXES"]);
    }


    #############
    # Working on TTDiagramm
    # Draw x-axes

    my $tmin = $model->{tmin};
    my $tmax = $model->{tmax};
    my $t1 = $model->{t1};
    my $t2 = $model->{t2};
    my $dt = $t2-$t1;
    $stepy = 1;
    if ($dt <= 2) {
        $stepy = 0.1;
    } elsif ($dt <= 4) {
        $stepy = 0.5;
    }
    #print "drawAxes:: tmin = $tmin, tmax = $tmax, t1 = $t1, t2 = $t2\n";

    # Draw axes
    @tick = $model->model2screen([$km1, $tmin, $km2, $tmax], "time", $km1);
    $lzd -> create('line',
            $tick[0], $tick[1], $tick[0], $tick[3],
            -joinstyle=>"bevel", -fill=>"black", -width=>2, -tags => ["AXES"]);

    # Draw ticks
    for (my $i=$tmin; $i <= $tmax; $i+=$stepy) {
        # $km2 is used for lenght of gridline
        # $i+stepy/2 is used for minor ticks and gridlines
        my $minortick = $i+$stepy/2;
        #print "Minortick = $minortick. major = $i\n";
        # need $km1 again for correct calculation of minortick
        #                       0         2                  4                  6         8
        @tick = $model->model2screen([$km1, $i, $km2, $minortick, $xmin, $minortick, $xmax, $i, $km1, $minortick], "time", $km1);


         #Draw gridlines (major)
        $lzd -> create('line',
                $tick[4], $tick[1], $tick[6], $tick[1],
                -joinstyle=>"bevel", -fill=>"grey", -width=>1, -tags => ["GRID"]);

        # Draw minor gridline
        $lzd -> create('line',
                $tick[4], $tick[9], $tick[6], $tick[9],
                -joinstyle=>"bevel", -fill=>"lightgrey", -width=>1, -tags => ["GRID"]);
        #print "Major grid: $tick[4], $tick[1], $tick[6], $tick[1]\n";
        #print "Minor grid: $tick[4], $tick[9], $tick[6], $tick[5]\n";

        # Major ticks
        $lzd -> create('line',
                $tick[0], $tick[1], $tick[0]+5, $tick[1],
                -joinstyle=>"bevel", -fill=>"black", -width=>1, -tags => ["AXES"]);
        # Minor ticks
        $lzd -> create('line',
                $tick[0], $tick[5], $tick[0]+2, $tick[5],
                -joinstyle=>"bevel", -fill=>"black", -width=>1, -tags => ["AXES"]);


        $lzd -> createText(
                $tick[0]+15, $tick[1],
                -fill=>"black", -text => "$i", -tags => ["AXES"]);
    }
}


sub _readTimes {

=PROGhead2 C<model::_readTimes()>

Function reads traveltimes from rayinvr. C<tx.in> und C<tx.out>.

=cut

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

    my $model = shift;
    #print "model::readTimes()\n" if ($model->tree());
    my $file = "tx.in";

    #say "Read picked times from $file";
    my $r = $model->_readTx( file => $file, key => 'txin');

=PROGfor html
stores txin-data using<br>
<code>$model->_readTx( file => $file, key => 'txin');</code>
<br>

=cut

    #print "Finished reading tx.in\n";

    #say "Reading calculated Traveltimes";
    $file = "tx.out";
    $r = $model->_readTx( file => $file, key => 'txout');

=PROGfor html
stores txout-data using<br>
<code>$model->_readTx( file => $file, key => 'txout');</code>
<br>
Times for tomo2D are currently read in function <code>tomo</code> and stored in key
<code>tomoTimesPg</code>, <code>tomoTimesPmP</code> or <code>tomoTimes</code>

=cut

    #print "Finished reading tx.out\n";
}

sub _readTx {

=PROGhead2 C<< model::_readTx( file => 'txfile', key => 'key') >>

Reads given C<txfile> and stores it in C<< $model->{key} >>.

=cut

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

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

    unless (exists $args{file} && exists $args{key}) {
        warn "File AND key are needed for _readTX!!\n".
             "Got %args\n";
        return 1;
    }

    my $file = $args{'file'};
    my $key = $args{'key'}; # Where to store data, model-hash key

    print "Reading tx-file >$file<\n" unless $quiet;

    open(T, $file) or print "Can't open $file\n", return 0;
    my $pos;
    my $station;
    my $richtung; # Welche Seite der Station ist der Pick? -1 fuer links, 1 fuer rechts

    my @readStations = (); # This arrray is a list of stations in your tx.in
    # If a station is twice in your tx.in it won't get cleared in the
    # second run

    #############################
    # Reading picks
    my %data;   # Hash of hashes with stations and phases as key and datapicks as values
    while(<T>){
        chomp;
        #print "$.: $_\n";
        s/^\s*#.*//;            # whole line is commented
        s/^#.*//;               # whole line is commented

        next unless length;

        my ($km, $t, $unc, $ph) = split;
        ## Lines like this mark new branches
        ##   188.959    -1.000     0.000     0
        ##  station-km  shot dir    mark    mark

        next unless ($t);

        $ph = 0 unless ($ph);
        #print "Read >$km, $t, $unc, $ph<\n";
        if (($t == -1 || $t == 1) && $ph == 0) {
        # if (($t == -1 || $t == 1) && $unc == 0 && $ph == 0) { # old conditions.
            # But $unc does not need to be zero. rayinvr works with non-zero values
            # New branch starts. First picks are usually left
            $richtung = $t;
            $pos = sprintf("%.3f",$km);
            #print "Reading arrivals for station at $pos, shot direction $t, $richtung\n";
            next;
        }

        unless ($pos ) {
            print "$.: $_\n";
            print "Station position undefined!!\n";
            next;
        }

        # tx-files all end with line
        #     0.000     0.000     0.000        -1
        if ($t == 0 && $unc == 0 && $ph == -1) {
            #print "Reached last line\n";
            next;
        }

        # This if-clause catches files with all picks in one branch, even
        # if they are on different sides of the station
        if ($km < $pos) {
            $richtung = '-1.000';
        } elsif ($km > $pos) {
            $richtung = '1.000';
        }

        if ($ph > 0) {
            # This is a pick. Save it to station
            # phase = 0     New station or branch
            # phase = -1    End of File

            # Store values
            # Pushing in unscaled values
            # print "Add pick to station at $pos, $richtung, $ph Data: $km, $t, $unc \n" if $dev;
            push @{$data{$pos}{$richtung}{$ph}}, [ $km, $t, $unc ];

=PROGfor html
<p> tx-times for tomo are stored in <code>$station->{$key}</code><br>
format: <code>push @{$data{$pos}{$richtung}{$ph}}, [ $km, $t, $unc ];
$station->{$key}->{$richtung}->{$ph} = @{$data{$pos}{$richtung}{$ph}}
</code><br>

=cut

        }
    } # END OF READING tx-formated file
    #print "Finished reading tx-file $file\n";
    close(T);

    #####
    # Storing data in station-object
    # Loop through stationpositions
    #print "Add pick data to stations\n";
    foreach my $pos (keys %data){
        my $station = $model->{stations}{$pos};
        #my $station = $model->{stations}{int($pos)};
        # Use int($pos) to remove unneeded numbers behind comma

        #print "Station >$station->{name}< at >$pos<\n";
        unless ( $station ) {
            print "\nError, no station defined for $pos!! This is likely a bug\n";
            $model->printStatusMessage("\nError, no station defined for $pos. This is likely a bug!!");
            next;
        }
        #print "Station $station->{name} at $pos is defined\n";
        $station->{$key} = undef; # Deleting ALL old traveltimes
        $model->{time}->delete('withtag',"RAYS$station->{name} "); # Clean up old data

        # Loop through directions
        foreach my $richtung (keys %{$data{$pos}}){
            # Loop through phases
            foreach my $ph (keys %{$data{$pos}{$richtung}}){
                #print "Station $station->{name}: Use phase $ph";
                my @picks = @{$data{$pos}{$richtung}{$ph}};
                #print " with >".@picks."< picks\n";

                $model->{time}->delete('withtag',"Ph$ph"); # Clean up old data
                #print "model::_readTx: Add >".@picks."< picks to station $station->{name}, key = $key, richtung = $richtung, ph = $ph\n";
                #$station->{$key}->{$richtung}->{$ph} = \@picks;    # This removes old existing data for that phase
                $station->addPhase(key => $key, direction => $richtung, phase => $ph, data => \@picks );
            }
        }
    }
    return \%data;

}

#sub getCode {

#=PROGhead2 C<model::getCode(ph)> depricated!!


#Connects raycodes (ray) with phase numbering in tx-file (ivray) and returns the opposite
#code to the input code, e.g.: gets 1.2 (a ray-code) and returns 12
#( the tx code)

#NEW: new object CODES hold this information.

#=cut

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

    #my ($model, %args) = @_;
    #my $phase = "0";

    ## phasecodes has ZP phases as keys (e.g. 12)
    ## raycode has RI phases as keys (e.g. 1.2)

    #return $model->{codes}->get(%args);


    ## Got tx-code, get ray-code
    #if (exists $args{phase} ) {
        #my $ph = $args{phase};
        #if (exists $model->{phasecodes}->{$ph}) {
            #$phase = $model->{phasecodes}->{$ph};
        #} else {
            #print "Cannot find raycode for phase >$ph< in >".join(" ",sort(keys(%{$model->{phasecodes}})))."<\n";
        #}
    #}

    ## Got raycode, Get tx-phasecode
    #elsif (exists $args{ray} ) {
        #my $ray = $args{ray};

        ## Phase is reflected off floating reflector and has a an ambigous
        ## phase naming. (ray = X.2 in r.in, but X.4 in r1.out)
        #if ($ray =~ /.4$/ ) {
            #$ray =~ s/.4$/.2/;
        #}
        ##$model->{codes}->getCode(
        ##my @phases =  grep { $model->{raycode}[$_]}  0 .. $#{$model->{raycode}};
        ##if () {
            ##$phase = $model->{raycode}->{$ph};
        ##} else {
            ##print "Cannot find tx phasecode for phase >$ph< in >".join(" ",sort(keys(%{$model->{phasecodes}})))."<\n";
        ##}
    #}

    ## Didn't get specification for code
    #else {
        #print "I don't know, what you gave me\n";
        #print Dumper \%args;
    #}

    #return $phase;
#}

sub getConfig {

=PROGhead2 C<model::getConfig(key)>

Returns corresponding value for given key or undef if not defined

=cut

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

    my $model = shift;
    my $key = shift;
    my $value = undef;

    if ( exists $model->{config}->{$key} && defined  $model->{config}->{$key}) {
        $value = $model->{config}->{$key};
    } else {
        print "model->getConfig($key) - No entry for $key\n";
        return undef;
    }

    return $value;
}

sub screen2model {

=PROGhead2 C<model::screen2model([$x, $y], "space/time", $pos)>

Converts screen coordinates in km. C<$pos> is optional, if given,
times will be unreduced. Otherwise they are reduced.

=cut

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

    my $model = shift;
    my $coords = shift;
    my @coords = @$coords;
    my $diagramm = shift;   # Convert in timespace or modelspace?
    my $pos = shift if @_;

    # Error check for typos
    unless ( $diagramm =~ "time" || $diagramm =~ "space") {
        confess "model::model2screen: Don't know which diagram to use for >$diagramm<";
    }

    my @scaled = ();

    # Get current scaling values
    my $yscale  = $model->{yscale};
    my $ytscale  = $model->{ytscale};
    my $xscale  = $model->{xscale};

    # Get model dimensions
    my $zmax = $model->{zmax};
    my $zmin = $model->{zmin};
    my $xmin = $model->{xmin};
    my $tmin = $model->{tmin};
    my $tmax = $model->{tmax};

    my $box = @{$model->{time}->configure(-scrollregion)}[-1];

    # Loop through given data and convert
    for (my $j=0; $j<= $#coords; $j++){
        # data comes in [x1,y1,x2,y2,..]
        #print "Convert $j: $coords[$j]\n";
        if ($j%2 != 0) {

            # y-value. Conversion depends on type of space: model or time
            if ($diagramm eq "space"){
                # y: model depth (z) in km
                push @scaled, $coords[$j]/$yscale+$zmin;
                #print "Depth $scaled[-1]\n";

            } elsif ($diagramm eq "time"){
                # y: traveltime value in s
                            # Reverse time?

                if ( $model->{config}->{reverseTime} ) {
                    #print "(DEV) reverse coord $coords[$j] to time with $box->[3] to " if $dev;
                    $coords[$j] = $box->[3]-$coords[$j];
                    #print "(DEV) Reversed coord: $coords[$j] \n" if $dev;
                }

                my $t = $coords[$j]/$ytscale+$tmin;

                if ($pos && $model->{vredstate} == 1) {
                    # A station position is defined. Scale coordinates
                    # without reduction
                    my $vred = $model->{vred};
                    $t = $t + abs($pos - abs($scaled[-1]))/$vred + $tmin;
                }
                push @scaled, $t;
            }
        } else {
            # x: Profilkm
            push @scaled, $coords[$j]/$xscale+$xmin;

        }

    }

    return @scaled;
}

sub model2screen{

=PROGhead2 C<model::model2screen([$x, $y], "space/time", $pos)>

Converts model coordinates in screen coordinates and reduces traveltimes
if needed (if C<$pos> is given and C<< $model->{vredstate} >> is set)

=cut

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

    # call: model2screen( model, coords, diagramm, stationkm)
    #print "Scale Got @_\n";
    my $model = shift; # Model

    my $zmax = $model->{zmax};
    my $zmin = $model->{zmin};
    my $xmin = $model->{xmin};
    #my $box = $model->{box};
    my $box = @{$model->{time}->configure(-scrollregion)}[-1];
    my $profilelength = $model->{profilelength};
    my $tmin = $model->{tmin};
    my $tmax = $model->{tmax};

    # TRY OOP SCALE
    my $yscale  = $model->{yscale};     # for model
    my $ytscale = $model->{ytscale};    # for time
    my $xscale  = $model->{xscale};

    my $vredbutton = $model->{vredstate};
    my $vred = $model->{vred};

    my @scaled; # = (1,2,3);
    my @in=$_[0];
    my $diagramm=$_[1];

    # Error check for typos
    unless ( $diagramm =~ "time" || $diagramm =~ "space") {
        confess "model::model2screen: Don't know which diagram to use for >$diagramm<";
    }

    my $pos = $_[2];

    for (my $j=0; $j<= $#{$in[0]}; $j++){

        # Every second value is a time value
        # data: km1, t1, km2, t2, ..
        if( $j%2 != 0 && $diagramm eq "time"){

            my $t;  # time in screen coordinates

            # reduce time
            if ($vredbutton == 1 ){

                # no offset for reduction known
                if (!defined $pos) {
                    #print "No position for reducing given, push in ".(($in[0][$j])*$ytscale)."\n" ;
                    # set position to current km, which basically undo the reduction
                    #$pos = $in[0][$j-1] if (!defined $pos);
                    #print "Set position to $pos\n";
                    $t = ($in[0][$j] - $tmin)*$ytscale;
                    #push @scaled, ($in[0][$j] - $tmin)*$ytscale;

                } else {
                    # Reduce time
                    #say " My position for reducing traveltime is $pos. data: >$in[0][$j-1]<";
                    #               Zeit     - |     x0-x            / v   |
                    #push @scaled, ($in[0][$j]-(abs($pos-$in[0][$j-1])/$vred) - $tmin)*$ytscale;
                    $t = ($in[0][$j]-(abs($pos-$in[0][$j-1])/$vred) - $tmin)*$ytscale;

                    #push @scaled, ($tmax - $in[0][$j]-(abs($pos-$in[0][$j-1])/$vred))*$ytscale;
                    #print "Scaled; @scaled\n";
                }
            } else {
                # Time value is not reduced

                #print "Button $vredbutton \n";
                $t = ($in[0][$j] - $tmin)*$ytscale;
                #push @scaled, ($in[0][$j] - $tmin)*$ytscale;

                #push @scaled, ($tmax - $in[0][$j])*$ytscale;
                #print "Zeit skaliert mit $in[0][$j]*$ytscale=". ($in[0][$j])*$ytscale."\n";
            }

            # Reverse time?
            if ( $model->{config}->{reverseTime} ) {
                #print "(DEV) reverse time $in[0][$j] to coordinate $t with $box->[3] to " if $dev;
                $t = $box->[3]-$t;
                #print "$t\n" if $dev;
            }
            push @scaled, $t;

        } elsif ($j%2 != 0 && $diagramm eq "space"){
            # DEPTH IN MODEL (y-value)
            push @scaled, ($in[0][$j]-$zmin)*$yscale;
            #print "Scaled tiefe  $xscale, $yscale or $ytscale is @scaled \n";

        } else {
            # Profilkilometer
            push @scaled,
            ($in[0][$j]-$xmin)*
            $xscale;

            #print "(DEV) m2s ( $in[0][$j] - $xmin) * $xscale\n" if $dev;
            #print "Scale with x>$xscale<, y>$yscale< or yt>$ytscale< is @scaled \n";
        }
    }
    #print "Returning array with ".@scaled." elements\n";
    return @scaled;
    }

sub reduce {

=PROGhead2 C<model::reduce(km,time,pos,v)>

Reduce given C<<time>> for C<<pos-km>> with current model reduction velocity.

=cut

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

    my ($model) = shift;
    my $km = shift;
    my $t = shift;
    my $pos = shift;
    my $v = $model->{vred};

    my $tred = $t - ( (abs($pos-$km) )/ $v) if ($v > 0);

    # Don't reduce if reduction is switched off.
    $tred = $t if ($model->{vredstate} == 0 );
    #print "Don't reduce time. It's switched off\n" if ($model->{vredstate} == 0   && $debug);

    #print "Got $t s, $pos - $km = ".abs($pos-$km)."km , $v , return reduced time t = $tred\n";
    return $tred;
}

sub unreduce {

=PROGhead2 C<model::unreduce(km,time,pos)>

Removes reduction for given C<<time>> for C<<pos-km>> with current model reduction velocity.

=cut

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

    my ($model) = shift;
    my $km = shift;
    my $tred = shift;
    my $pos = shift;
    my $v = $model->{vred};

    my $t = $tred + ( (abs($pos-$km) )/ $v);
    #print "Got reduced time t = $tred s, $pos - $km = ".abs($pos-$km)."km , $v , return $t\n";
    return $t;
}

sub updateScale {

=PROGhead2 C<model::updateScale()>

Get's new values for y-,yt- and xscale. Deletes and redraws velocity and
depth nodes to keep a circular shape

=cut

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

    # Call: $model->updateScale("yscale" => $yscale, "ytscale" => $ytscale, "xscale" => $xscale);
    my ($model, %args) = @_;
    my $cns = $model->{space};
    my $lzd = $model->{time};

    #print "Updating to: ".%args."\n";
    while (my ($key, $value) = each(%args)){
        #say "$key => $value ";
        $model->{$key}= $value;

    }
    $model->drawAxes;

    $model->{space}->delete('VNODE');
    $model->{space}->delete( "NODE");
    foreach (@{$model->{layers}}){
        $_->drawNodes;
        if ($model->{'vnodes'} == 1){
            $_->drawVNodes;}
    }

    # Delete traveltimes
    $lzd->delete('PICK');
    if (defined $model->{drawnStations}){
    # Is phase button pushed before stationbuttons?
        my @positions = @{$model->{drawnStations}};
        for (my $i = 0; $i <= $#positions; $i++){
            say "model::drawPhase() Draw for Station $positions[$i]" if $debug;

            $model->{stations}{$positions[$i]}->draw;
        }
    }
    $model->order;
}

sub set {

=PROGhead2 C<< model::set( $key => $value ) >>

Interface for external setting of model configuration. C<$key> is checked
before applied. Only C<vredstate vred PicksMan PicksCal ShowRays glueNodes
version xz nodes annotvnodes contours>
are allowed. C<< $model->order >> does the changing.

=cut

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

    my ($model, %args) = @_;
    #say "model::set(%args)"; print Dumper(\%args);

    # Only options in this array are allowed to be set.
    my @options = (qw/vredstate vred PicksMan PicksCal ShowRays glueNodes
        version xz xt nodes annotvnodes contours contourcolor blocks vgrid vnodes densities
        resolution/);

    # Safe all valid keys and values in model
    foreach my $key (keys(%args)) {
        #print "model::set $key = $args{$key}\n";
        if (grep {$_ eq $key} @options){
        # Is key allowed to be set?

            if (ref($args{$key})) {
            # is argument a reference?

                $model->{$key} = ${$args{$key}};
                print "model::set $key = ${$args{$key}}\n" if $debug;
            } else {
            # argument is no reference

                $model->{$key} = $args{$key};
                print "model::set $key = $args{$key}\n" if $debug;
            }
        } else {
        # key is not allowed to be set
            print "Unknown option >$key< in model::set\nOptions are @options";
        }
    }

    ###################################################################
    # Some keys require action like drawing objects
    # - blocks
    # - vred

    # If traveltime reduction has been changed, picks have to be redrawn
    if ( exists $args{"vredstate"} || exists $args{"vred"} || exists $args{"PicksMan"}
        || exists $args{"PicksCal"} ){
            #TODO: find a more elegant way to remove observed and calculated picks
            # from station. fix this with drawPhaseStationList
            $model->{time}->delete('withtag', 'PICK');
        $model->drawPhaseStationList();
    }

=PROGpod

 if (exists $args{"blocks"} && == 1 )
 delete BLOCK
 foreach layer->_drawBlock

=cut

    if (exists $args{"blocks"} && ${$args{"blocks"}} == 1 ) {
        $model->{space}->delete('BLOCK');
        foreach my $layer (@{$model->{layers}}){
            $layer->_drawBlocks;
        }
    }

    $model->order;
}

sub reset {

=PROGhead2 C<model::reset()>

This routine resets the model and is run, when the model is changed.
So all things, that change as consequence are deleted.
That are contours and grid. It deletes them, so they have to be
regenerated by model:: order.
All rays are deleted from stations.
All calculated traveltimes.

=cut

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

    my ($model, %args) = @_;
    my $lzd = $model->{time};
    my $cns = $model->{space};

    $cns->delete('CLINE');
    $cns->delete('VGRID');

    ## Delete old data
    foreach (sort {$b <=> $a}(keys(%{$model->{stations}}))) {
        # Safety maeasure for unblessed stations
        if  (ref $model->{stations}{$_} eq "station") {
            $model->{stations}->{$_}->emptyRays;
            $model->{stations}->{$_}->emptyTimes(keys => ['txout']);
        } else {
            print "WARNING !! There's an unblessed station in your model!! Call the priest !!\n";
        }
    }
}

sub writeVin {

=PROGhead2 C<model::writeVin()>

Called by p and controls writing C<v.in> formatted file by each layer calling
C<< layer->writeVin >> and resets model.

=cut

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

    my ($model, %args) = @_;
    my @layers = @{$model->{layers}};

    my $file = "v.in";
    my $output = "v";   # Default output are velocities
    if (exists $args{file}) {
        if ( $args{file} eq "rho.in") {
        # Write densities?
            $output = 'rho';
            $file = "rho.in";
        } else {
        # Write resolutions
            $file = $args{file};
            $output = "res";
        }
        print "Write file $file\n";
    }

    #print "model:: Writing v.in\n";

    # Start writing for each layer
    # Desired output is
    # v     velocities
    # rho   densities
    # res   resolution values

    # Delete old v.in file
    #print "(DEV) remove v.in ($file)\n";
    #unlink $file;

    open(FILE, ">$file") or die "Can't open $file";
    foreach (@layers){
        printf FILE $_->writeVin( 'output' => $output );
    }
    print "(D) Wrote file >$file< with $output-data\n" if $debug;
    close(FILE);

    # If v.in gets written, write an additional file for conversion
    # with vin2xyz (that v.in file needs values at both boundaries
    if ( $file eq "v.in" ) {
        #$model->writeV2xyzIn();
        $model->writeV2xzyIn;
    }

    # Copy file to export path (do keep a copy in rayinvr folder for
    # running plots and stuff
    print "copy($file, $model->{config}->{exportpath}$file\n";
    copy("$file", "$model->{config}->{exportpath}$file");

    # Why resetting?
    $model->reset;
}

sub writeV2xzyIn {

    my ($model, %args) = @_;
    my @layers = @{$model->{layers}};

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

    my $file = "v2xyz.in";
    my $output = "v2xyz";
    open(FILE, ">$file") or die "Can't open $file";
    foreach (@layers){
        printf FILE $_->writeVin( 'output' => $output );
    }
    close(FILE);
    print "(D) Wrote file >$file< with $output-data\n" if $debug;

}

sub get {

=PROGhead2 C<model::get(key)>

 model->get()
 interface function to p.pl to get information about model

 Extract information out of tags for vnodes

 key can be: vnode, vnodes, 1d or any key from model hash

=cut

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

    my $model = shift;
    my $arg = shift;
    my $cns = $model->{space};

    if ($arg eq "vnode"){
        # Get information for a single node

        my $id  = $cns->find(qw/withtag current/);
        my @tags = $cns->gettags($id);
        my $layer = $tags[1];
        #my ($km, $v, $vu, $vl, $vupar, $vlpar);
        my ($km, $v, $pos) = split (/ /, $tags[4]);
        # $pos is either vu or vl
        print "I'm .$arg. >$tags[1]<, id @$id, @tags, $v = v, $km = km\n";
        # get(vu/vl, index)
        #($km, $v, $vu, $vl, $vupar, $vlpar) = $model->getLayer($layer)->get("vnode", $tags[5], $tags[2]);
        #return ($v, $km, $vu, $vl, $vupar, $vlpar);
        my @return  = $model->getLayer($layer)->get("vnode", $tags[5], $tags[2]);
        return @return;
    }

    if ($arg eq "vnodes"){
        # Get information for ALL vnodes

        my (@km, @vu, @vl, @vupar, @vlpar);
        # Get velocity nodes for each layer
        foreach (@{$model->{layers}}) {
            print "(DEV) Get nodes for $_->{number}\n" if $dev;
            #last layer doesn't have any velocity nodes
            last if ($_->{last} == 1);

            # lkm: layerkm
            my ($lkm, $lvu, $lvl, $lvupar, $lvlpar) = $_->get("vnodes");
            #print "Got back @$lkm\n@$lvu\n@$lvl\n";
            push @km, $lkm;
            push @vu, $lvu;
            push @vl, $lvl;
            push @vupar, $lvupar;
            push @vlpar, $lvlpar;

        }
        # Return references to arrays containing references
        return \@km, \@vu, \@vl, \@vupar, \@vlpar;
    }

    if ($arg eq "1d") {
        # Print out 1D velocity profile

        my $kmstring = shift;
        my @km = split(/,/, $kmstring);

        # Clean directory from old files
        #opendir(DIR, "") or die $!;
        #my @files = grep { /^vd_\./ && -f "$_"}  readdir(DIR);
        #closedir(DIR);
        my @files = glob("km*_v*.vz");
        print "Delete files @files\n";

        foreach (@files){
            unlink("$_");   # Potential Problem if you run it on Windows!!
        }

        foreach my $km (@km) {
            $km =~ s/\s//g;
            #$km =~ sprintf ("%03f", $km);
            print "Extract for km >$km<\n";
            #my $file = "km${km}_v$model->{version}.vz";
            my $file = sprintf ("km%03d_v%d.vz", $km, $model->{version});
            open(FILE, ">$file") or die "Cannot open $file";
            foreach my $layer (@{$model->{layers}}){
                my $vd = $layer->get($arg, $km);
                print $vd;
                printf FILE $vd;
            }
        close(FILE);
        }

    }

    if (exists $model->{$arg}) {
        return $model->_get($arg);
    }

}
sub _getColor {

=PROGhead2 C<model::_getColor(v)>

Returns color code depending on given velocity.Used for
gradient grid, contour lines and layer colors.

=cut

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

    my $v = shift;

    ######
    # Create color-palette
    my $vmin = 1.2;     # Lowest velocity to color
    my $vmax = 8;       # Highest velocity

    # Make sure, you don't overstep your limits
    if ( $v > $vmax ) { $v = $vmax }

    # Change color in steps of 0.1 km/s
    my $hmin = 240;     # "blue" in HSV-colorsystem
    #my $hmin = 340;        # "blue" in HSV-colorsystem
    my $hmax = 0;       # "red"

    my $i = sprintf "%3i", ($vmax - $v)/(($vmax - $vmin))*$hmin;
                               #h   s    v

    my $color = hsv2rgb( $i, 0.50, 0.80 ); # s und v < 1 !!
    #print "Color $i, $color, \n";

    return $color;
}

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

    my $model = shift;
    my $color = shift;
    my $cns = $model->{space};

    my @coords = ( qw/0 0 10 10/ );
    eval {$cns -> create('line', @coords, -joinstyle=>"bevel",
        -fill=> $color, -width=>1, tags => ['COLOR'])};

    if ( $@ ) {
        #print "ERROR $@\n";
        $color = "#$color";
    }

    $cns->delete('COLOR');
    return $color;
}

sub exportRays {

=PROGhead2 C<model::exportRays()>

 1.  Delete old exported data in output directory (this can be switched
     off in p.config)
 2a. Loop through all stations
 2b. Loop through all phases
     Export rays
     Export picked phases
     Export calculated phases

=cut

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

    my ($model, %args) = @_;
    my $dir = $model->{config}->{exportpath};

    print "Export Rays and picks to $dir\n";

    ############################################
    # 1. Remove old rays and picks in export directory
    unless (exists $model->{config}->{deleteExported} && $model->{config}->{deleteExported} == 0) {
        print "Deleting old exported data\n";
        if (-d $dir) {
            opendir(DIR, $dir) or die $!;
            my @files = grep { -f "$dir/$_"}  readdir(DIR);
            closedir(DIR);

            # Delete only pick, rays and refl-files
            foreach (@files){
                if ($_ =~ m/(pick_st*)/ || $_ =~ m/(ray_st*)/ || $_ =~ m/(refl_st*)/ ){
                    # $1 represents the digit selected by regex (\d+)
                    #print "Delete $dir$_\n";
                    unlink("$dir/$_");   # Potential Problem if you run it on Windows!!
                }
            }
        }
    }

    # Export absolute AND relative output
    for ( my $absolute = 0 ; $absolute <=1; $absolute++) {
    ############################################
    # 2a. Loop through all stations
    foreach (keys(%{$model->{stations}})){
        my $pos = $_;
        my $station = $model->{stations}{$pos};
        say "----- Get data for station $station->{name} at $pos";

        # Relative position used for output x-position
        my $relpos = $pos;  # x position is relavitve to obs

        my $fileRays = $dir."/ray_st$_.dat";
        my $fileReflected = "$dir/refl_st$_.dat";
        my $fileOutPicks = $dir."/pick_st$_"."out.dat";
        my $fileInPicks = "$dir/pick_st$_"."in.dat";

        #if ($model->getConfig("outpos") eq "absolute" ) {
        if ($absolute ) {
            $relpos = 0;    # x position is absolute model km
            $fileRays        = "$dir/ray_st${_}_a.dat";
            $fileReflected   = "$dir/refl_st${_}_a.dat";
            $fileOutPicks    = "$dir/pick_st${_}out_a.dat";
            $fileInPicks     = "$dir/pick_st${_}in_a.dat";
        } else {
            $fileRays        = "$dir/ray_st${_}_r.dat";
            $fileReflected   = "$dir/refl_st${_}_r.dat";
            $fileOutPicks    = "$dir/pick_st${_}out_r.dat";
            $fileInPicks     = "$dir/pick_st${_}in_r.dat";
        }

        # Get rays
        my $rays = $model->{stations}{$pos}->getRays;
        #my $picks = $model->{stations}{$pos}->getPicks("direction" => '-1.000');
        my $manpicks = $model->{stations}{$pos}->{txin};
        next unless (defined $rays);

        my @rcodes = sort(keys(%$rays));
        print "(D) Got rays: $rays, ph: @rcodes\n" if $debug;

        next unless ( @rcodes);


        open (FILERAYS,     ">$fileRays")     or die "Can't open $fileRays\n";
        open (FILEREFL,     ">$fileReflected")or die "Can't open $fileReflected\n";
        open (FILEOUTPICKS, ">$fileOutPicks") or die "Can't open $fileOutPicks\n";
        open (FILEINPICKS,  ">$fileInPicks")  or die "Can't open $fileInPicks\n";
        #print "Write files $fileRays, $fileReflected, $fileOutPicks, $fileInPicks\n";

        ############################################
        # 2b. Loop through all raycodes
        foreach my $rc (@rcodes) {
            print "Print raycode $rc\n" if $debug;
            my $newDir = -1; # Start with rays left of station

            ###########################
            # RAYS
            foreach (@{$rays->{$rc}}) {
                next unless (defined $_ );
                #print "Ray: $_\n" ;

                # Export turning points of reflections
                if ($rc =~ m/2$/) {
                    #$turningPoints = 1;
                    my $i = 2;
                    my @a = @{$_};
                    my %h = @a;
                    my %i = reverse %h;
                    my $max = (sort(values(%h)))[-1];
                    #$i{$max};
                    #print "Print raycode $rc to file\n";
                    if ($newDir == -1 &&  $pos < $i{$max}) {
                        #print "Change direction at $i{$max} $max\n";
                        print FILEREFL ">\n";
                        $newDir = 1;
                    }
                    print FILEREFL "$i{$max} $max\n";
                }

                # Export raypath
                for (my $i = 0; $i < $#{$_}; $i += 2) {
                    # 0 is just printed for compatibility
                    # with make_picks_rays.csh. There's a -1 for left shots and 1 for
                    # right shots. But I don't make this distinction here
                    #print "$_->[$i] $_->[$i+1] 0 $rc\n";
                    print FILERAYS "$_->[$i] $_->[$i+1] 0 $rc\n";

                }
                #print ">\n";
                print FILERAYS ">\n";

            } # end of rays loop
            print FILEREFL ">\n";
        } # end of rcodes loop

            ###########################
            # PICKS
            my $key = $model->{tomotimes} ? 'txTomo' : 'txout' ;


        foreach (qw/ -1.000 1.000/){
            ###########################
            # CALCULATED TRAVELTIMES
            #my $ph = $model->{codes}->get('ray'=>$rc);
            # Get all phases for a station
            foreach my $ph (@{$station->getPhases()}) {

                print "(D) Try to find out if picks are defined for key $key : $_, phasecode $ph\n" if $debug;
                my $picks = $station->getPicks( 'direction' => $_, 'phase' => $ph, 'key' => $key );

                if (defined $picks) {
                    print "(D) Got picks $picks\n" if $debug;
                    #print "txout for station at $pos, ph $ph, dir $_ is defined \n";
                    #print "Got an array $_\n";
                    #my $picks = @{$picks->{$_}->{$ph}};
                    #print "There are <".@{$picks}."> calculated picks for phase $ph\n";

                    for (my $i = 0; $i <= $#{$picks}; $i += 1) {
                        # Copy values to variabls, otherwise data is changed
                        # in the object!!
                        my $x = $picks->[$i][0];
                        my $t = $picks->[$i][1];
                        #print "my x: $x ref($x), t=$t\n";
                        printf FILEOUTPICKS "%.3f %.3f 0 %d\n",
                            ($x-$relpos),
                            $model->reduce($x,$t, $pos),
                            $ph;
                    }
                    print FILEOUTPICKS ">\n";
                }

                ###########################
                # PICKED TRAVELTIMES
                $picks = $station->getPicks
                    ( 'direction' => $_, 'phase' => $ph, 'key' => 'txin' );
                print "(D) Get picks for direction $_\n" if $debug;
                next unless (defined $picks);

                #print "There are >$#{$picks}< for Phase $ph and $station->{name} and direction $_\n";
                for (my $i = 0; $i <= $#{$picks}; $i += 1) {
                    # Copy values to variabls, otherwise data is changed
                    # in the object!!
                    my $x = $picks->[$i][0];
                    my $t = $picks->[$i][1];
                    my $unc = $picks->[$i][2];

                    printf FILEINPICKS "%.3f %.3f 0 %d\n",
                        ($x-$relpos),
                        $model->reduce($x,$t-$unc, $pos),
                        $ph;
                    printf FILEINPICKS "%.3f %.3f 0 %d\n",
                        ($x-$relpos),
                        $model->reduce($x,$t+$unc, $pos),
                        $ph;
                    print FILEINPICKS ">\n";
                }

            } # foreach direction
        } # foreach phase
            #} else {
                ## EXPORT PICKS FROM TOMO2D
                #my $picks = $model->{stations}{$pos}->{txTomo};
                #next unless $picks;
                #print "Got $picks \n";
                #foreach (keys(%{$picks})) {
                    #my $ph = $_;

                    #for (my $i = 0; $i < $#{$picks->{$ph}}; $i++) {
                        ## Copy values to variabls, otherwise data is changed
                        ## in the object!!
                        #next unless (defined $picks->{$ph}->[$i]);
                        #my $x = $picks->{$ph}->[$i][0];
                        #my $t = $picks->{$ph}->[$i][1];
                        ##print "my x: $x ref($x)\n";

                        #print FILECALPICKS ($x-$pos)." "
                            #.$model->reduce($x,$t, $pos)
                            #." 0 $ph\n";

                        ## if this x and next x don't have same sign
                        ## insert ">\n" for gmt to start a new line
                        #unless ( $picks->{$ph}->[$i+1][0] && commons::eqsig ($x-$pos, $picks->{$ph}->[$i+1][0]-$pos) ) {
                            #print FILECALPICKS ">\n";
                        #}

                    #} # end of for pick loop
                    #print FILECALPICKS ">\n";
                #} # end of foreach keys(picks)
            #} # end if-else tomo2d picks


            # Get picks
            # OUTPUT, ready for plotting with gmt, SHOULD LOOK LIKE:
            #-18.077 3.731375 0 0.0
            #-18.077 3.831375 0 0.0
            #>
            #-17.254 3.66325 0 0.0
            #-17.254 3.76325 0 0.0
            #
            # Picks are seperated by ">" to mark new segment for gmts multiple
            # segments plotting
            # zeroes at the end are inserted for compartibility with make_rays
            # script by Tobi


        close(FILERAYS);
        close(FILEREFL);
        close(FILEOUTPICKS);
        close(FILEINPICKS);

    } # end of station loop
} # for absolute = 0/1

} # end sub export rays

sub resolution {

=PROGhead2 C<model::resolution>

Control subs for reading and writing resolution. Called from p.pl,
export resolution button

=cut

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

    my $model = shift;

    # Read resolution
    $model->_readRes();
    # Now print res.in to file
    $model->writeVin( 'file' => "res.in");
    $model->printStatusMessage(" .. wrote resolution to res.in ..");
    $model->writeXZV();
    $model->order;

} # sub resolution


sub _readRes {

=PROGhead2 C<model::_readRes>

Reads d.out and stores data in model->layers.
d.out contains a list with nodes which have been selected for adjustments
and an updated velocity model (same that is written to v.in).
The node list contains only the original value, adjustment and new value
but no position or index to identify which node it is.
You've got to look into the model to find the nodes position.

For resolution plots depth resolution makes no sense.

=cut

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

    my $model = shift;

    my @type;
    my @value;
    my @resolution;
    my @stddev;
    my @results;    # Save result-table (all information for nodes)
    my $paragraph = 'start'; # Flag

    # Keep current position infos
    my $currentLayer = 0;
    my $currentNode = 0;

    # Save data in an array of arrays and add them to each layer
    # when finished reading
    my @res;

    ####################################################################
    # Read resolution file
    my $file = "d.out";
    open(FILE, $file) or print "Can't open $file\n", return 0;
    while(<FILE>){
        chomp;
        s/^\s+//;               # no leading white
        s/\s+$//;               # no trailing white
        next unless length;

        my $line = $_;
        my @line = split(/\s+/,$line);
        #                  t or un ad nw re std
        #print "My line is >@line<\n";
        #print "Line $. ".join(';',@line)."\n";

        if ($line =~ m/overall damping factor/) {
            print "Damping factor = $line[-1]\n";
            next;
        }

        # Set flag for reading adjustments, uncertainties, resolution, .. values
        if ($line[0] eq 'type'){
            #print "Reading resolution values\n";
            $paragraph = 'resolution';
            next;
        }

        # set flag for reading velocity model
        if ($line[0] eq 'velocity'){
            #print "Reading updated model\n";
            $paragraph = 'velocity';
            last;
        }

        # Reading resolution values
        if ($paragraph eq 'resolution') {
            # Overwrite line. Split does not work proberly when values
            # become too large
            @line = unpack("A1A11A11A11A11A11A11",$line);
            # Remove whitespaces
            @line = split('\s+',join(' ', @line));
            #print "Line $. ".join(';',@line)."\n" if $debug;

            push @type, "$line[0]";
            push @value, "$line[4]";
            push @resolution, "$line[5]";
            push @stddev, "$line[6]";

            print "$. resolution $resolution[-1] $line[5]\n" if $debug; # Print resolution value

            # save linenumber
            unshift @line, $.;
            push @results, \@line;
            #print "READ: ".join('-',@line)."\n";
        } # endif reading resolution

    } # end while FILE
    close(FILE);

    #print "Resolution values: @resolution\n";
    #print "Standard deviation: @stddev\n"
        #."Value: @value\n\n\n";
    ####################################################################

    my $n = @{$model->{layers}}; # Get number of layers

    # Find nodes in layer-objects
    my $i = 0;  # Index for nodes
    my $j = 0;  # Index for layers

    # nodetype are the old values that can be compared with the information
    # given in the nodes list to see if the node fits
    my @nodetype = qw(depth vu vl);   # Do NOT change sequence!!!
    my %nodepos = ( depth => 'km',
                    vu => 'vukm',
                    vl => 'vlkm');


    #print "Start looking for nodes\n"
        #."Node type value resolution stddev\n";
    while ( $i < $#type ) { # Loop through all nodes
        my $type = $type[$i];

        print "$i $type[$i] $value[$i] $resolution[$i] $stddev[$i]\n" if $debug;

        while ($j < $n ) { # loop through n-layers
            last if ($i > $#type);

            my $layer = $model->{layers}[$j];
            my @dres;

            print "(D) Look in layer $layer->{number} for $i "
                ."$type[$i] $value[$i] $resolution[$i] $stddev[$i]\n" if $debug;

            # Get partial derivative arrays of current layer for depth,
            # upper and lower velocity in this sequence.
            foreach my $nodetype (@nodetype){ # loop through depth, vu, vl
                my $derivatives = $nodetype.'par';
                #print "Look for derivativekey $derivatives and values $nodetype\n";

                # Find indizes of elements with activated partial derivs
                my @activated = grep { $layer->{$derivatives}[$_] == 1 } 0 .. $#{$layer->{$derivatives}};

                # Create with same length for resolution
                @dres = (-1)x@{$layer->{$derivatives}};

                #if (@activated) {
                    print "(D) Found n = ".@activated." nodes with partial derivatives = 1 in layer $layer->{number} $nodetype: >@activated< \n" if $debug;
                    #print "Created resolution array >@dres<  with ".@dres." elements \n";
                #}

                # Check if values fit and save resolution at same index
                foreach my $active (@activated) {
                    # Check if values are equal
                    # Due to differences in rounding procedures it's necessary
                    # to allow small differences between values
                    #print "Get node $active from layer $layer->{number} in array @{$layer->{$nodetype}}\n";
                    unless (defined $value[$i]) {
                        my $msg =
                        "\n\n WARNING !!!\n\n"
                        ."There's no new value in d.in for $nodetype-node $layer->{$nodetype}[$active] at "
                        ."km $layer->{$nodepos{$nodetype}}[$active]\n"
                        ."\n\nThis is a rayinvr-problem!\n\n"
                        ;
                        print $msg;
                        $model->printStatusMessage($msg);
                        last;
                    }

                    # Ignore nodes with value == 0
                    # rayinvr interpolates them with value from other nodes and
                    # does not calculate a resolution for them
                    if ( $layer->{$nodetype}[$active] == 0 ) {
                        print "(D) ignore node index $nodetype [ $active ] = $layer->{$nodetype}[$active] "
                            ."for layer $layer->{number} at km $layer->{$nodepos{$nodetype}}[$active].\n";
                        next;
                    }

                    #print "Look for value $value[$i]\n";
                    my $diff = $value[$i] - $layer->{$nodetype}[$active];
                    my $eps = 0.0051;

                    if ( $resolution[$i] !~ /\d/ ) {
                        my $msg = "WARNING !! your resolution is no number!! res = $resolution[$i]. Set to zero\n";
                        print $msg;
                        $resolution[$i] = 0;
                    }

                    if ($resolution[$i] < 0 || $resolution[$i] > 1){

                        my $msg = "\n\n!! WARNING !!\n\n"
                            ."Problem with node in line @{$results[$i]}\n"
                            ."Resolution >$resolution[$i]< for node $value[$i] makes no sense!! Your rayinvr seems to have problems\n";

                        # Setting to -1 will remove this value, when writing vin.
                        # Setting to 0 will keep it.
                        $resolution[$i]=-1;
                        #$resolution[$i]=0;

                        $msg.="Set value to $resolution[$i] for output resolution file\n"
                            ."Node position: B$layer->{number}, block '$nodetype', node number $active "
                            ."km $layer->{$nodepos{$nodetype}}[$active]\n";
                        print $msg;
                        $model->printStatusMessage($msg);

                        #return;
                    }

                    # The node from the list is found in the model
                    if ( abs($diff) < $eps){
                        print "Found node $i (line $results[$i][0]) in B$layer->{number}, $nodetype, $active "
                         ."Values are equal $value[$i] == $layer->{$nodetype}[$active], resolution = $resolution[$i]\n" if $debug;
                        $dres[$active] = $resolution[$i];
                        $i++;   # Work on  next value
                        #print "(D) i=$i, total number of nodes $#type\n" if $debug;
                        last if ($i > $#type);
                    } else {
                        my $msg =
                            "\n\n!! WARNING !!\n\n"
                            ."Problem with node in line @{$results[$i]}\n"
                            ."Index $active, Values are not equal $value[$i] != $layer->{$nodetype}[$active], resolution = $resolution[$i]\n"
                            ."\nDiff is $diff, a difference of $eps is allowed for 'equal' values\n"

                            ."'new val.' in d.out for $i.th node is >$value[$i]<\n"
                            ."current value in v.in for $i.th activated node is $layer->{$nodetype}[$active]\n"
                            ."Node position: B$layer->{number}, block '$nodetype', node number $active "
                            ."km $layer->{$nodepos{$nodetype}}[$active]"
                            ."\nblock abbreviations:\n"
                            ."depth = depth nodes (first block)\n"
                            ."vu = upper velocity nodes (second block)\n"
                            ."vl = lower velocity nodes (third block)\n"
                            ."\n\nYour v.in is not the v.in written by dmplsqr."
                            ." Resolutions from d.out are not valid for this v.in!\n"
                            ."Run dmplsqr and then try again to export resolutions\n"
                            ;
                            print $msg;
                            $model->printStatusMessage($msg);

                            return;
                    }
                } # end foreach array with partial derivatives

                # Store resolutionarray in layer
                my $reskey = $nodetype.'res';
                print "Save resolutions >@dres< to >$layer->{number}< at >$reskey<\n" if $debug;
                $layer->{$reskey} = [@dres];
            } # foreach derivativekeys

            $j++;   # next layer
        } #end while $j < $n
        if ($j > $n) {
            my $msg = "Error, some values haven't been found\n";
            print $msg;
            $model->printStatusMessage($msg);
            return;
        }
        if ($i < $#type && $j == $n) {

            my $msg =  "\n\n WARNING !! \n\n"
                ."Could not find all nodes with partial derivative set in"
                ." your v.in !! d.out is not valid for this model\n";

            print $msg;
            $model->printStatusMessage($msg);
            return;
        }
        print "End while $i. $j\n";
        #$i++;

    } # while i < type

} # sub _readRes

sub writeXZV {

=PROGhead2 C<model::writeXZV()>

Writes a gmt-readable file with velocity nodes.
Function that calls writeXZV for each layer and writes the returned string
to file C<< $model->{config}->{exportpath}v.xzv >>

=cut

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

    my $model = shift;
    #my $file = "$model->{config}->{exportpath}/v.xzv";
    my $file = "v.xzv";

    open(FILE, ">$file") or die "Can't open $file";
    foreach my $layer (@{$model->{layers}}) {
        #print "Write layer $layer->{number}\n";
        printf FILE $layer->writeXZV();
    }
    close(FILE);

    # Copy file to export path (do keep a copy in rayinvr folder for
    # running plots and stuff
    print "copy($file, $model->{config}->{exportpath}$file\n";
    copy("$file", "$model->{config}->{exportpath}$file");

    $model->printStatusMessage("exported velocity nodes to $file ");
}

sub exportPolygons {

=PROGhead2 C<model::exportPolygons()>

Writes a gmt-readable file with layer polygons.
Function that calls C<writePoly()> for each layer and writes the returned string
to file C<< $model->{config}->{exportpath}v.poly >>

=cut


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

    my $model = shift;
    #my $file = "$model->{config}->{exportpath}/v.xzv";
    my $file = "v.poly";

    open(FILE, ">$file") or die "Can't open $file";
    foreach my $layer (@{$model->{layers}}) {
        print "(DEV) Write layer $layer->{number}\n" if $dev;
        printf FILE $layer->writePoly();
    }
    close(FILE);

    # Copy file to export path (do keep a copy in rayinvr folder for
    # running plots and stuff
    print "copy($file, $model->{config}->{exportpath}$file\n";
    copy("$file", "$model->{config}->{exportpath}$file");

    $model->printStatusMessage("\nExported layer polygons to $file ");

}


sub exportIgmas {

=PROGhead2 C<model::exportIgmas()>

Function that calls writeIgmas for each layer and writes the returned string
to file C<< $model->{config}->{exportpath}pVersion.mod >>

=cut

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

    my $model = shift;
    my $file = "$model->{config}->{exportpath}/igmas.mod";
    #my $file = "v.xzv";

    print "Export model to igmas format in file $file\n";

    open(FILE, ">$file") or die "Can't open $file";

    printf FILE "#gmdf_ascii\n\n\n";

    # Get densities
    foreach my $layer (@{$model->{layers}}) {
        #next if ( $layer->{number} ==  10 ); #????
        printf FILE $layer->igmasGetDensity();
    }

    printf FILE
         "\nreference_density    3.21\n\n";

    printf FILE $model->_igmasWriteSection('back');
    printf FILE $model->_igmasWriteSection('model');
    printf FILE $model->_igmasWriteSection('front');

    close(FILE);

    $model->printStatusMessage("\nExported igmas model to $file");
    $model->printStatusMessage("\nNodes are reduced for clarity. Coordinates"
        ." for igmas model are still hardwired in this code. "
        ."Please manually edit your model file.");

}

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

    my $model = shift;
    my $section = shift;

    # Define two origin points
    #               x1       y1        x2       y2
    my @point = ( [453799, 7745287],[846635, 7902328] );    # P100
    #my @point = ( [ 67348, 7763403],[509778, 7563694] );    # P150

    print "# (I) Coordinates for igmas model are still hardwired in this code. "
        ."Please manually edit your model file.\n";

    if ( $section eq "back" ) {
        $point[0][1] += 10000; # 10 km north
        $point[1][1] += 10000; # 10 km north
        $point[0][1] += 200000; # 200 km north
        $point[1][1] += 200000; # 200 km north
        #$point[0][1] += 1000000; # 1000 km north
        #$point[1][1] += 1000000; # 1000 km north

    }

    if ( $section eq "front" ) {
        $point[0][1] -= 10000; # 10 km south
        $point[1][1] -= 10000; # 10 km south
        $point[0][1] -= 200000; # 200 km south
        $point[1][1] -= 200000; # 200 km south
       #$point[0][1] -= 1000000; # 1000 km south
        #$point[1][1] -= 1000000; # 1000 km south
    }

    my $s =
         "\n\n"
        ."#####################################################\n"
        ."# SECTION $section\n"
        ."#####################################################\n"
        ."\nsection $section\n"
        ."point       $point[0][0]    $point[0][1]\n"
        ."point       $point[1][0]    $point[1][1]\n"
        ."\n\n"
        ."coordinates\n";

    # Get coordintes of nodes
    my $n = @{$model->{layers}};
    #$n = 3;
    for (my $i = 1; $i <= $n ; $i++ ) {
        #print "########################################\nGet coordinates for layer $i\n" if ($debug);
        $s .= $model->getLayer($i)->igmasGetCoordinates($section);
    }

    #print "###############\nGET VERTICES for section $section\n" if ($debug);
    # Get vertices for each layer
    for (my $i = 1; $i <= $n-1; $i++ ) {
        $s .=   $model->getLayer($i)->igmasGetVertices($section);
    }

    # Clean up
    $model->{igmasCoords} = {};

    return $s;
}


1;
} # end of model package