###################################
# THIS PROGRAMM IS DOCUMENTED WITH INLINE PODs
# There are different comments for users and programmers
# Extract pod with
# USERs: sed '/=PROG/,/=cut/d' p.pl | sed 's/=USER/=/g' | pod2html
# PROGRAMMERs: sed '/=USER/,/=cut/d' p.pl | sed 's/=PROG/=/g' | pod2html
###################################
=PROGhead1 NAME model.pm
model - containing objects to represent a rayinvr model used by p.pl
=cut
{
########################################################################
package model;
use strict;
use warnings;
use Carp;
use lib '/opt/pray/lib';
use Carp qw(cluck);
use feature qw(switch say);
use POSIX ('ceil', 'floor');
use Data::Dumper;
use Graphics::ColorUtils;
use Math::Round;
use File::Copy;
use layer;
use floatingReflector;
use station;
use codes;
#use Tk::Animation;
#use Tk::Splashscreen;
=PROGhead1 C
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
model->tomo [\&model::tomo, $model, "tomorays", \"1"]
=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
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
Creates startup-splash (if configured) and reads model data. Currently
started from C. 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
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
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
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
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
Returns position of closest station for given km. Currently only used when
C.
=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
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
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
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
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 = ) {
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
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
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
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. 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
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
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() {
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
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 () {
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, F7.2)
# (3) the z-coordinates (km) of the nodes defining the floating
# reflector (format: 3X, 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, 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
Large routine reading C and C. Before reading any data, old
ray-information is deleted from model-objects. Then phase information about rays
is read from C then data is read from C and connected with phases
from C.
Information from C is stored invariable C<@r1> for later use, when reading
C
=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 (){
$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 (){
$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
Add blocks read from C 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
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 (){
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
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 (){
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
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
