HOME
ABOUT
- RESULTS
- differences
- BENEFITS
- HISTORY
- TEAM
- LOCATION
- FACILITIES
- BANKING
- MEMBERSHIPS
- APPROVALS
- LICENCES
- SUPPLIERS
- SPONSORSHIPS
- MEDIA
- PRIVACY
AUCTIONS
SHIPPING
FEES
- TS REWARDS
TOOLS
guides
FAQ
CONTACT
- CONNECT

VEHICLES
BRAND
- JAPANESE CARS
  - DAIHATSU
  - EUNOS
  - FORD
  - HONDA
  - ISUZU
  - LEXUS
  - MAZDA
  - MITSUBISHI
  - MITSUOKA
  - NISSAN
  - SUBARU
  - SUZUKI
  - TOYOTA
- GERMAN CARS
- AMERICAN CARS
- BRITISH CARS
- ITALIAN CARS
- FRENCH CARS
- SWEDISH CARS
- KOREAN CARS
TYPE
- mobility
- VENDING
- instruction
- TAXIS
- AMBULANCES
- FIRE ENGINES
- HEARSES
- LIMOUSINES
- COMMERCIAL
CLASS
FUEL
TRUCKS
minitrucks
- DAIHATSU
- HONDA
- MAZDA
- MITSUBISHI
- NISSAN
- SUBARU
- SUZUKI
- DUMP
- CRANE
- CAMPER
- REFRIGERATED
- 4WD
- NEW
BUSES
MOTORHOMES
- YAHOO!
- RAKUTEN
- DEALER

PARTS
- FREE REPORT
- PARTS CONTAINERS
- PARTS SYSTEMS
- PARTS PROTECTION
- BODY SHELLS
- DISMANTLING
- ONLINE PARTS
- NEW PARTS
- INTERIOR PARTS
- EXTERIOR PARTS
  - BONNETS
  - BUMPERS
  - GRILLES
  - FENDERS
  - DOORS
  - TRUNKS
  - SPOILERS
  - LIGHTS
  - EMBLEMS
  - CAMERAS
- ENGINES
- TRANSMISSIONS
- WHEELS & TYRES
  - WHEELS
  - TYRES
CUTS
PERFORMANCE PARTS
TRUCK PARTS
MOTORBIKE PARTS
- MOTORBIKE ENGINES
- MOTORBIKE ACCESSORIES

MOTORBIKES
MARINE
FORKLIFTS
MACHINERY
AGRICULTURAL
OTHER
COUNTRY
- AUSTRALIA
- CANADA
- KENYA
- MYANMAR
- NEW ZEALAND
- PAKISTAN
- TANZANIA
- UNITED STATES

CARVIEW

MOTORHOMES

Select Language

HTTP/2 200 date: Sat, 11 Oct 2025 22:08:04 GMT content-type: text/html; charset=UTF-8 server: cloudflare x-frame-options: DENY x-content-type-options: nosniff x-xss-protection: 1;mode=block vary: accept-encoding cf-cache-status: DYNAMIC content-encoding: gzip set-cookie: _csrf-frontend=c6d3e35fc2ec27800dcbbf021c40275183025a725c1cd6a98838a0111ed0334ea%3A2%3A%7Bi%3A0%3Bs%3A14%3A%22_csrf-frontend%22%3Bi%3A1%3Bs%3A32%3A%22BSxzmrFCzze-s_hm5jYUAFJsmN-BiM-g%22%3B%7D; HttpOnly; Path=/ cf-ray: 98d1ac87ba84a9c3-BLR #################################### THIS PROGRAMM IS DOCUMENTED WITH INLINE - Pastebin.com

Pastebin

API tools faq

paste

Login Sign up

SHARE

Untitled

johnkris

Oct 10th, 2025

88

0

Never

Add comment

Not a member of Pastebin yet? Sign Up, it unlocks many cool features!

C 143.56 KB | None | 0 0

raw download clone embed print report

###################################
# THIS PROGRAMM IS DOCUMENTED WITH INLINE PODs
# There are different comments for users and programmers
# Extract pod with
# USERs: sed '/=PROG/,/=cut/d' p.pl | sed 's/=USER/=/g' | pod2html
# PROGRAMMERs: sed '/=USER/,/=cut/d' p.pl | sed 's/=PROG/=/g' | pod2html
###################################
=PROGhead1 NAME model.pm
model - containing objects to represent a rayinvr model used by p.pl
=cut
{
########################################################################
package model;
use strict;
use warnings;
use Carp;
use lib '/opt/pray/lib';
use Carp qw(cluck);
use feature qw(switch say);
use POSIX ('ceil', 'floor');
use Data::Dumper;
use Graphics::ColorUtils;
use Math::Round;
use File::Copy;
use layer;
use floatingReflector;
use station;
use codes;
#use Tk::Animation;
#use Tk::Splashscreen;
=PROGhead1 C<model>
Draws rayinvr-model for seismic refraction data, keeps data storage, do
administration
=PROGhead2 Tomo2D
=begin html
'model' can also be used for viewing data from tomo2D. Reading tomo-data
is started with button 'Tomo2D' from main PRay GUI which invokes
<code>
model->tomo [\&model::tomo, $model, "tomorays", \"1"]
</code>
=end html
=PROGhead2 Model keys
contours
d2
km1
debug
profilelength
stationpositions
raycode
balloon
layers
d1
tmin
km2
tmax
stations
vredstate
layer
vgrid
zmax
tomotimes
ytscale
vred
foo
nodes
box
phasecodes r
ays
PicksCal
readStr
time
t1
t2
PicksMan
zmin
drawnStations
canvasheigth
xscale
canvaswidth
yscale
ShowRays
version
vnodes
xmin
space
drawnPhases
xmax
config
tomoGrid
tomorays
phasecolors
DEBUGING
All station objects are stored in the hash $model->{stations}
x keys(%{$model->{stations}}) with position-km as hashkey.
e.g.:
DB<17> x keys(%{$model->{stations}})
0 250.000
1 89.300
2 175.000
3 25.000
4 225.000
5 50.000
6 275.000
7 200.000
8 100.000
9 125.000
10 150.000
DB<6> x $model->getStation(127) Gives you the stationreference
=cut
our ($tree, $debug, $dev, $verbose, $quiet); # Debugging and verbosity variables
sub new {
=PROGhead2 C<model::new(%args)>
A new model is born, blessed and issued with configuration values. Only argements,
which are also defined in %defaults are allowed.
my $model = new model('space'=> $cns, "time" => $lzd,
"balloon" => $balloon,
"zmin" => $config{zmin}, "zmax" => $config{zmax},
"yscale" => $yscale, "ytscale" => $ytscale, "xscale" => $xscale,
);
In Future it shall be:
my $model = new model(
["reverseStation" => (1/0)], stationfile => statxz,
"working dir" => ... , tmin, tmax
"zmin" => $config{zmin}, "zmax" => $config{zmax},
);
=cut
#print "model::new()\n";
my ($class, %args) = @_;
my $model = bless ({}, $class);
$model->{rays}=undef;
#$model->init(%args);
# old $model->init function:
# Define default parameters for a model object
my %defaults = ( "space" => undef, # canvas for model
"time" => undef, # canvas for ttdiagramm
"statusline" => undef, # GUI status line
"mainwindow" => undef, # MAIN WINDOW
"splash" => undef,
"icons" => undef,
"balloon" => undef, # window, needed for balloons
"init" => -1, # Is this model initialized? -1 no, 0 yes, 1 currently initializing
"layer" => 0, # ??? obsolete??
"layers" => undef, # array for storing layer objects
"profilelength"=> undef,
"box" => undef,
"foo" => "bar",
"zmin" => undef,
"zmax" => undef,
"xscale" => undef,
"yscale" => undef,
"ytscale" => undef,
"vred" => 3,
"vredstate" => 1,
"canvaswidth" => 1600, # Width of cavas for time and model
"canvasheigth" => 400, # Heigth of cavas for time and model
"tmin" => 0,
"tmax" => 12,
"stations" => undef, # refernce hash with links to stations
"stationpositions" => undef, # hash with names and positions of stations
"stationkm" => undef, # hash with positions and names of stations
"drawnStations" => [],
"phasecolors" => { "00" => "darkgrey",
"12" => "red",
"21" => "green",
"22" => "blue",
"31" => "orange",
"32" => "black",
"41" => "yellow",
"42" => "springgreen",
"51" => "green",
"52" => "blue",
"61" => "orange",
"62" => "black",
"71" => "yellow",
"72" => "springgreen",
},
# colors for phases and layer boundarys.
"nodes" => 1, # Show and edit nodes
"vnodes" => 0, # Show and edit vnodes
"annotvnodes" => 1,
"drawnPhases" => undef, # which phases should be drawn
"drawnRays" => undef, # which rays should be drawn
"PicksMan" => 1, # Show manual picks?
"PicksCal" => 1, # Show calculated picks?
"vgrid" => 0,
"config" => undef, # Reference to configurationhash
"rin" => undef, # Reference to rin
"contours" => 0,
"contourcolor" => 0,
"blocks" => 0,
"debug" => undef,
"ShowRays" => 1, # Show rays? Usefull when too many raypathes block view on model
"phasecodes" => { 21 => 2.1}, #conversion of phasecodes between ZP and RI
"raycode" => { 2.1 => 21}, #conversion of phasecodes between ZP and RI
"tomorays" => 0, # TEMPORARY
"tomoGrid" => 0, # Inital value for displaying tomogrid
"tomotimes" => 0, # Flag for read tomo2D traveltime data
"version" => '', # Save model version
"glueNodes" => 1, # Moving pinched nodes together?
"xz" => 0, # Show xz overlay files?
"xt" => 0, # Show xt overlay files?
"igmasCoords" => {}, # Hash to store igmas coordinates (see layer->getIgmasCoordintes()
"densities" => 0, # Annotate vnodes with densities after Barton
"resolution" => 0,
"codes" => undef, # Store codes object
);
foreach (keys(%defaults)){
#print "Setting default for $_ -> $defaults{$_}\n";
$model->{$_}=$defaults{$_};
}
foreach (keys(%args)){
say "initModel: $_ $args{$_}" if $debug;
croak "Unknown attribute <$_> for station initialization" unless exists $defaults{$_};
$model->{$_} = $args{$_};
}
#foreach (keys(%{$model->{config}})) {
#say "initModel: apply config $_ $model->{config}->{$_}";
#$model->{$_} = $model->{config}->{$_} if exists $defaults{$_};
#}
$model->{vred} = $model->{config}->{vred};
$model->{box} = [0, 0, $model->{canvaswidth}, $model->{canvasheigth}]; # Size of model and runtimediagram in px
# Setup scale parameters
$model->{km1} = $model->{config}->{xmin};
$model->{km2} = $model->{config}->{xmax} ;
$model->{d1} = $model->{config}->{zmin};
$model->{d2} = $model->{config}->{zmax};
$model->{t1} = $model->{config}->{tmin};
$model->{t2} = $model->{config}->{tmax};
$model->{xmin} = $model->{config}->{xmin};
$model->{xmax} = $model->{config}->{xmax};
$model->{profilelength} = $model->{xmax} - $model->{xmin};
my $cns = $model->{space};
#say "Exit initalization of your model";
#print Dumper \$model->{color};
# Setup debugging variables
$tree = $model->{debug}->{tree};
$debug = $model->{debug}->{debug};
$dev = $model->{debug}->{dev};
$verbose = $model->{debug}->{verbose};
$quiet = $model->{debug}->{quiet};
return $model;
}
sub init {
=PROGhead2 C<model::init(%args)>
Creates startup-splash (if configured) and reads model data. Currently
started from C<p.pl>. Cannot be run by model->new because stations need
to be defined in the model.
=cut
printf "(T) %s(@_)\n", commons::whoami() if $tree;
my ($model, %args) = @_;
my $splash;
my $mw;
# Display splash screen if switched on in configuration file
if ( $model->{config}->{splash} == 1) {
print "Display splash screen\n";
$model->_set("init", 1);
$mw = $model->_get("mainwindow");
# Make splash screen
$splash = $mw->Splashscreen;
$mw->withdraw;
$splash->Splash; # show Splashscreen
$splash->Label(-text => 'Building your Model ...', -bg => 'blue')->
pack(qw/-fill both -expand 1/);
my $animate;
if (@ARGV) {
$animate = $mw->Animation;
foreach (@ARGV) {
$animate->add_frame($mw->Photo(-file => $_));
}
} else {
#my $gif89 = Tk->findINC('anim.gif');
my $icons = $model->_get("icons");
my $gif89 = "$icons/Praysplash.gif";
#my $gif89 = '/projects/nam2011/bin/pray/icons/Praysplash.gif';
print "Take gif $gif89\n";
$animate = $mw->Animation(-format => 'gif', -file => $gif89);
}
$animate->set_image(0);
$animate->start_animation;
$splash->Label(-text => "Button me",
#-image => $mw->Photo(-file => "$icons/reload.gif")
-image => $animate
)-> pack(qw/-fill both -expand 1/);
$model->_set("splash" => $splash);
$mw->update;
}
## READ DATA
$model->read("vin", "rays", "times");
$model->_set("init", 0);
$model->drawAxes();
# Draw stations
foreach (keys(%{$model->{stations}})){
if (ref $model->{stations}{$_} eq "station") {
$model->{stations}{$_}->drawStation;
} else {
print "WARNING !! There's an unblessed station in your model!! Call the priest !!\n";
}
}
# Draw user defined xz if configured
if ($model->{config}->{xz}) {
foreach my $file (@{$model->{config}->{xz}}) {
$model->_readxzfile("$file");
}
}
# Draw user defined xz if configured
if ($model->{config}->{xt}) {
foreach my $file (@{$model->{config}->{xt}}) {
$model->_readxzfile("$file", 'xt');
}
}
if ( $args{splash} == 1) {
$splash->Destroy; # tear down Splashscreen
$mw->deiconify; # show main window
}
}
sub _get{
=PROGhead2 C<model::_get($key, ..)>
function returns any value stored in model.
if a key is not found, undef is returned
THIS FUNCTION IS INTENDED FOR INTERNAL USE
external is C<model::get>
arguments are a list of keys to retrieve
=cut
#printf "(T) %s(@_)\n", commons::whoami() if $tree;
my ($model, @args) = @_;
return undef unless @args;
my ($key, $value);
my @return;
foreach $key (@args){
if (exists $model->{$key}){
$value = $model->{$key};
#print "Get ".$key." = ".$value."\n";
} else {
print "ERROR: no such key as $key in model!!\n";
$value = undef;
}
push @return, $value;
}
#D Return scalar, if only one argument
#D and referece to array if more
if ( @return == 1) {
return $value;
} else {
return \@return;
}
} #D model->_get END
sub _set{
=PROGhead2 C<< model::_set( $key => $value, .. ) >>
model->_set can set a reference to data for any existing key in
a model
THIS FUNCTION IS INTENDED FOR INTERNAL USE
external is set
=cut
printf "(T) %s(@_)\n", commons::whoami() if $tree;
my ($model, %args) = @_;
my ($key, $value);
while (($key, $value) = each(%args)){
if (exists $model->{$key}){
#print "Set ".$key.", ".$value."\n";
$model->{$key} = $value;
} else {
print "ERROR: no such key as $key in model!!\n";
die;
}
}
}
sub status {
=PROGhead2 C<model::status()>
Prints number of layers, average, minimum and maximum depth and velocities, number of stations with phases.
$model->status;
=cut
printf "(T) %s(@_)\n", commons::whoami() if $tree;
say "--------------------------";
say "This is your status-report";
my ($model, %args) = @_;
my $s =" B min max mean variance stddev\n"; # Return string
my $n = keys(%{$model->{stations}});
say "Your model has <$n> stations\n";
#foreach (sort {$b <=> $a}(keys(%{$model->{stations}}))) {
#$model->{stations}->{$_}->statusstation;
#}
my @layers = @{$model->{layers}};
say "Your model has <".@layers."> layers";
for (my $i=0; $i<=$#layers; $i++){
#say "Status layer $i";
$s .= $model->{layers}[$i]->status(%args);
}
say "--------------------------";
print "Modelparameter : ";
my @par = keys (%$model);
print "@par\n";
#print "ISA: \n";
#$model->DESTROY;\
return $s;
}
sub printStatusMessage{
=PROGhead2 C<model::printStatusMessage("\nmsg")>
Prints msg to GUI. Same like main prog in p.pl
=cut
printf "(T) %s(@_)\n", commons::whoami() if $tree;
my ($model, $msg) = @_;
my $stline = $model->{"statusline"};
$stline->insert ('end', "$msg");
$stline->see('end');
$model->_get("mainwindow")->update();
}
sub getClosestStation{
=PROGhead2 C<model::getClosestStation($km)>
Returns position of closest station for given km. Currently only used when
C<model::tomoReadRays>.
=cut
printf "(T) %s(@_)\n", commons::whoami() if $tree;
my ($model, @args) = @_;
# Used if you want to calculate moveout velocity.
# coding not complete
my $x = $args[0];
my @positions = sort{$a <=> $b}(keys(%{$model->{stations}}));
my $dx;
my $min = 900;
my $st;
foreach (@positions) {
$dx = abs($_ - $x);
print "pos: $_, dx $dx, min $min\n" if ( $model->debug());
if ($dx > $min) {
print "Closest Position:dx $dx min $min pos $st\n" if ( $model->debug());
last;
} else {
#print "Closest Position: $_\n";
$min = $dx;
$st = $_;
print "Closest Position: $min, $dx, $st\n" if ( $model->debug());
}
}
print "Returning $st\n" if ( $model->debug());
return $st;
}
sub addStation{
=PROGhead2 C<model::addStation()>
Initializes new station with C<< new station(%args, "model" => $model) >>
and fills information in model-object:
# Storage place for stations: $model->stations{position}
$model->{stations}{$st->{position}} = $st;
# Information to find station
# name <-> position
$model->{stationpositions}->{$st->{name}} = $st->{position};
$model->{stationkm}->{$st->{position}} = $st->{name};
Stations are stored at hash key C<< $model->{stations}{$st->{position}} = $st >> !!
Station position is an integer with three numbers behind comma.
Check stations in model in perls debuging console:
x keys(%{$model->{stations}})
=cut
printf "(T) %s(@_)\n", commons::whoami() if $tree;
my ($model, %args) = @_;
my $st = new station(%args, "model" => $model);
#say "at ".$st->{position};
# Storing station!!
$model->{stations}{$st->{position}} = $st;
# Add information about station. Enables finding the station with
# stationname or stationposition
$model->{stationpositions}->{$st->{name}} = $st->{position};
$model->{stationkm}->{$st->{position}} = $st->{name};
}
sub edit{
=PROGhead2 C<model::edit()>
used for all kind of model editing like adding, deleting, moving nodes
C<< $model->edit("op" => "edit", "tags" => \@tags, "value" => [$km, $v]); >>
=over
=item 'editPhase'
# For SPrange selection: tags = [phase, ["tagsofPick1"], ["tagsofPick2"]
# For changing all picks: tags = [phase, st]
=item 'allV'
=item other
calls
layer->edit(%args);
=item return message
=back
=cut
printf "(T) %s(@_)\n", commons::whoami() if $tree;
my ($model, %args) = @_;
my @tags = @{$args{"tags"}};
my $op = $args{"op"};
my $msg = ""; # Return message to model
#print "model::edit tags: @tags\n";
#print Dumper(\%args);
my $type = $tags[0];
print "My type is $type\n" if $debug;
if ($op eq "editPhase") {
$model->editPhase(%args);
}
if ($op =~ m/.*Shot/) {
#Find out, which station was clicked
my $station = $model->getStation($tags[2]);
$station->set(%args);
}
#####
# THIS HANDLES ALL OTHER OPERATIONS !!!
if ($type ne "allV" && $op ne "editPhase" && $op !~ m/.*Shot/){
#print "model::edit($type) args = \n";
#print Dumper \%args;
my $layer = $tags[1];
my $node = $tags[2];
# TODO: Bug: Throws error for floating refl. Seems to work with
# just one refl, but probably causes error in the future
$model->getLayer($layer)->edit(%args);
}
if ($type eq "allV"){
for(my $i = 0; $i < $#{$model->{layers}}; $i++){
#print "Save v nodes into layer $i\n";
my $km = $args{value}->[0]->[$i];
my $vu = $args{value}->[1]->[$i];
my $vl = $args{value}->[2]->[$i];
my $vupar = $args{value}->[3]->[$i];
my $vlpar = $args{value}->[4]->[$i];
#print "Editging km @$vl\n";
# Tags are not necessary, but $layer->edit needs some values there
$model->{layers}->[$i]
->edit( "op" => "editAllV", "tags" => ["TYPE", "LAYER", "NODE"], "value" => [$km, $vu, $vl, $vupar, $vlpar]);
}
}
$model->order;
}
sub editAllParDerivs {
=PROGhead2 C<model::editAllParDerivs()>
Set all partial derivatives to given value. This sub is used, when calculating
the resolution for a model. So far no other use.
=cut
printf "(T) %s(@_)\n", commons::whoami() if $tree;
my ($model, $switch) = @_;
my $op;
$op = 'setPar' if ($switch == 1) ;
$op = 'unsetPar' if ($switch == 0);
print "model:: editAllParDerivs: $op\n" if $debug;
my @coords = (0,0); # Fake, not needed
foreach my $layer (@{$model->{layers}}) {
my @tags = ('LAYER', "B$layer->{number}");
$layer->edit("tags", \@tags, "value", \@coords, "op" , $op);
}
$model->order();
}
sub editPhase {
=PROGhead2 C<model::editPhase()>
Edit phasecodes in zp-Files. Called with right click on a phase in the
traveltime canvas.
=cut
printf "(T) %s(@_)\n", commons::whoami() if $tree;
my ($model, %args) = @_;
my @tags = @{$args{"tags"}};
my $mode = $args{mode};
(my $oldphase = $args{value}[0]) =~ s/\.//;
(my $newphase = $args{value}[1]) =~ s/\.//;
# For SPrange selection: tags = [phase, ["tagsofPick1"], ["tagsofPick2"]
# For changing all picks: tags = [phase, st]
print Dumper(\@tags) if $debug;
# Environment
my $config = $model->{config};
my $zpdir = $config->{zpdir};
my ($file, $st, $sp1, $sp2);
# If only selected picks should change, SPs have to be found for them
if (@tags == 3){ # If there are not 3 elements in tags array, then there's no selection of SPrange
# Get variables for first pick
(my $st1 = $tags[1][1]) =~ s/RAYS//;
(my $ph1 = $tags[1][2]) =~ s/Ph//;
my ($pkm1, $t1, $unc) = split (/ /, $tags[1][3]);
$pkm1 =~ s/km//;
$t1 =~ s/t//;
$unc =~ s/unc//;
$ph1 =~ s/\.//; # Phases in Pickfile are without decimal ### TODO CHANGE THIS TO USE RIN MATCHES
# Calculate offset of this pick
my $off1 = sprintf ("%.3f", $pkm1 - $model->{stationpositions}->{$st1});
print "Pickoffset is $off1\n" if $debug;
# Get variables for second pick
(my $st2 = $tags[2][1]) =~ s/RAYS//;
(my $ph2 = $tags[2][2]) =~ s/Ph//;
my ($pkm2, $t2, $unc2) = split (/ /, $tags[2][3]);
$pkm2 =~ s/km//;
$t2 =~ s/t//;
$unc2 =~ s/unc//;
$ph2 =~ s/\.//; # Phases in Pickfile are without decimal ### TODO CHANGE THIS TO USE RIN MATCHES
# Calculate offset of this pick
my $off2 = sprintf ("%.3f", $pkm2 - $model->{stationpositions}->{$st2});
print "Pickoffset is $off2\n" if $debug;
$file = $model->getStation($st1)->{zpfile};
$file =~ s/head/unc/;
# Find picks in unc file to get shotpoints
my $line;
open (FILE ,"<$zpdir/$file") or die "Can't open myfile: $_!\n";
while ($line = <FILE>) {
if ($line =~ /$off1\s+$ph1/i) {
my @line = split ( /\s+/ , $line);
$sp1 = $line[1];
print $line };
if ($line =~ /$off2\s+$ph2/i) {
my @line = split ( /\s+/ , $line);
$sp2 = $line[1];
print $line };
}
print "SP=$sp1, SP = $sp2\n" if $debug;
if (!defined $sp1 || !defined $sp2) {
print "COULDN'T FIND BOTH SHOTPOINTS\n";
return undef;}
# Sort shotpoints so sp2 is larger than sp1
if ($sp2 < $sp1 ) {
print "Swap sp1 $sp1 with sp2 $sp2\n" if $debug;
my $t = $sp2;
$sp2 = $sp1;
$sp1 = $t;
print "Now: sp1 $sp1, sp2 $sp2\n" if $debug;}
#print "My stationnumber is $st1, ph $ph1 km $pkm1, t $t1, unc $unc\nOpen file $file\n";
$st = $st1;
} # Now SPs for picks are found
else {
$st = $tags[1];
print "CHange all picks for station $st\n" if $debug;
}
$file = $model->getStation($st)->{zpfile};
$file =~ s/head/pick/;
#$file = "100st135.h.rloc.pickbak";
my ($sp, $t) = commons::readPicks("$zpdir/$file");
die "No file $file found\n" if (!defined $sp );
if ($sp1) {
print "CALL editPhase($sp, $t,$sp, $t, 'sp1' => $sp1 'sp2' => $sp2,".
"'oldphase' => $oldphase, 'newphase' => $newphase, 'mode' => $mode \n";
($sp, $t) = commons::editPhase($sp, $t, 'sp1' => $sp1, 'sp2' => $sp2,
"oldphase" => $oldphase, "newphase" => $newphase, "mode" => $mode);
} else {
print "CALL editPhase($sp, $t, 'oldphase' $oldphase 'newphase' $newphase 'mode' $mode)\n";
($sp, $t) = commons::editPhase($sp, $t,
'oldphase' => $oldphase, 'newphase' => $newphase, 'mode' => $mode);
}
my $fileout = "$zpdir/$file";
#my $fileout = "$zpdir/100st135.h.rloc.pick" ;
commons::writePicks($sp, $t, $fileout) if ($sp);
# Clean edited picks
print "model::edit() delete 'Ph$oldphase'\n";
$model->getStation($st)->emptyTimes( keys => ['txin'], phases => [$oldphase]);
$model->{time}->delete("Ph$oldphase");
return;
}
sub getLayer {
=PROGhead2 C<model::getLayer($number)>
Returns the layer-object for $layernumber. Number is the same as in v.in
=cut
printf "(T) %s(@_)\n", commons::whoami() if $tree;
my ($model, $number) = @_;
# Remove leading letters (B 1 -> 1)
$number =~ s/[a-zA-Z]//g;
# Floats do not have any number left
#print "layer::getLayer -> number = $number\n" if $debug;
#if ($number) {
# Get layer-object
my $l = $model->{layers}[$number-1];
return $l;
#} else {
#return 0;
#}
}
sub getStation {
=PROGhead2 C<model::getStation($number)>
Returns station object which matches C<$number> either as stationname or
stationposition. Currently it's a bit of luck, which one is used. Look at
position first but if a stationnumber accidently is the same as the position
of another, it gets confused. (This is a potential BUG, TODO)
=cut
printf "(T) %s(@_)\n", commons::whoami() if $tree;
# returns station-object
# give either stationname or stationpositions
my ($model, $number) = @_;
#print "model::getStation($number)\n" if $debug;
# Return station reference
my $st;
# Is number a station name? Then get profilekm for this station
if (exists $model->{stationpositions}->{$number}) {
#print "Input numer is a stationposition\n" if $debug;
$number = $model->{stationpositions}->{$number};
}
# Is number a stationkilometer? Then return the station object
if (exists $model->{stations}->{$number}) {
#print "Input number is a stationname\n" if $debug;
$st = $model->{stations}->{$number};
}
# Check if station is still undefined
unless ($st) {
print "WARNING !!! Couldn't find a station for >$number<.\n";
}
# No check, if both might be right. It shouldn't anyway, because
# stationkms have three decimals and stationnames don't
#print "Return station $st->{name}\n";
return $st;
}
sub getPhase {
=PROGhead2 C<model::getPhase(phase)>
Get alist of all station having a special phase picked
=cut
printf "(T) %s(@_)\n", commons::whoami() if $tree;
my ($model, @args) = @_;
my $ph = $args[0];
my @stations = sort { $a <=> $b } (keys(%{$model->{stations}}));
my $s = ""; # Return string
foreach my $st ( @stations ) {
my $n = 0;
if ( exists $model->getStation($st)->{txin}->{'1.000'}->{$ph} ) {
$n += @{$model->getStation($st)->{txin}->{'1.000'}->{$ph}};
}
if ( exists $model->getStation($st)->{txin}->{'-1.000'}->{$ph} ) {
$n += @{$model->getStation($st)->{txin}->{'-1.000'}->{$ph}};
}
# Only return string for stations with this phase
if ( $n > 0 ) {
my $name = $model->getStation($st)->{name};
#print "Get station $name\n";
$s .= sprintf("%3d: Ph %3d %4d picks\n", $name, $ph, $n);
}
}
return $s;
}
sub getDensity {
#printf "(T) %s(@_)\n", commons::whoami() if $tree;
=PROGhead2 C<$model->getDensity (v)>
Returns density for given velocity C<v>. Conversion uses values from
Barton(1986) "The relationship between seismic velocity and density in the continental crust - a useful constraint?"
=cut
my $model = shift;
my $v = shift;
$v = sprintf("%.1f",nearest(.1,$v));
#my $conv = "ludwig";
my $conv = "barton";
#my $conv = $model->getConfig('densityconversion');
my $rho;
if ( $conv eq "barton" ) {
# Densities after Barton(1986)
# "The relationship between seismic velocity and density in the continental crust - a useful constraint?"
my %densities = qw( 0.0 0.0 1.5 1.47 1.6 1.66 1.7 1.73 1.8 1.80 1.9 1.86
2.0 1.92 2.1 1.98 2.2 2.01 2.3 2.03 2.4 2.06
2.5 2.09 2.6 2.11 2.7 2.13 2.8 2.15 2.9 2.18
3.0 2.21 3.1 2.23 3.2 2.24 3.3 2.26 3.4 2.28
3.5 2.30 3.6 2.32 3.7 2.34 3.8 2.36 3.9 2.38
4.0 2.39 4.1 2.41 4.2 2.43 4.3 2.44 4.4 2.46
4.5 2.48 4.6 2.50 4.7 2.52 4.8 2.53 4.9 2.55
5.0 2.57 5.1 2.59 5.2 2.61 5.3 2.62 5.4 2.64
5.5 2.66 5.6 2.68 5.7 2.70 5.8 2.72 5.9 2.74
6.0 2.77 6.1 2.80 6.2 2.83 6.3 2.85 6.4 2.87
6.5 2.90 6.6 2.93 6.7 2.95 6.8 2.98 6.9 3.01
7.0 3.04 7.1 3.07 7.2 3.10 7.3 3.13 7.4 3.16
7.5 3.19 7.6 3.22 7.7 3.25 7.8 3.28 7.9 3.31
8.0 3.34 8.1 3.38 8.2 3.42 );
$rho = $densities{$v};
}
if ( $conv eq "funck" ) {
# Densities after Ludwig(1970)
# "Seismic Refraction" in THE SEA
#a=-0.00283; b=0.0704; c=-0.598; d=2.23; e=-0.7;
#printf("v= %f, rho= %8.2f\n", $1,a*$1^4+b*$1^3+c*$1^2+d*$1+e)
# Funck 2004 has approximated ludwig/barton
# But values are quite different!!
$rho = -0.00283*$v**4 + 0.0704*$v**3 - 0.598*$v**2 + 2.23*$v -0.7;
}
unless ($rho) {
#print "WARNING!! Cannot find density for v=$v\n";
if ($v > 8.2 ) {
$rho = 3.5
}
if ($v < 1.5 ) {
$rho = 1.47;
}
#print "Use density rho = $rho\n";
}
return $rho;
}
sub read {
=PROGhead2 C<model::read('data')>
Routine starts individual model-subroutines for reading data.
Keys can be C<'vin', 'rays', 'times', 'str', 'vgrid'>
=cut
printf "(T) %s(@_)\n", commons::whoami() if $tree;
my ($model, @args) = @_;
my $r = ""; # Return message
if (grep $_ eq "vin", @args) {
# Delete contours, grid
$model->reset();
$model->_readVin;
$model->_readRes();
}
if (grep $_ eq "fin", @args) {
$model->reset();
$model->_readFin;
}
if (grep $_ eq "rays", @args) {
$r = $model->_readRays;
}
if (grep $_ eq "times", @args) {
$model->_readTimes;
}
if (grep $_ eq "str", @args) {
# Read Streamerpicks and convert them to depth
$model->_readStr;
}
if (grep $_ eq "vgrid", @args){
$model->_readVgrid;
}
return $r;
}
sub _readVin {
=PROGhead2 C<model::_readVin()>
Called by C<< $model->read(vin) >> and reads model file C<'v.in'>.
=cut
printf "(T) %s(@_)\n", commons::whoami() if $tree;
my ($model, %args) = @_;
my $cns = $model->{space};
my $xmax = $model->{rin}->{xmax}[0];
my $xmin = $model->{rin}->{xmin}[0];
my $file = "v.in";
say "Reading file v.in" unless $quiet;
open(FILE, $file) or croak "Can't open $file";
my $l=0;
my @line;
my $currentlayer=0;
my $pre; # prefix of line. holds layernumber or shows that line is contiuned
my $layer; # current layer as layer-object
my $v=0; # Flag for depth or velocity in second line of block
# 0 for depth
# 1 for last depth reading
# 2 for velocity upper border
# 3 for velocity lower border
my @km; # stores profilekm for more than one line
my @d; # stores depth corresponding to @km
my @par; # stores partial derivatives
#my (@uv, @lv); # stores velocity-Nodes
my $n = 0; # Number of nodes for a layer. Used as flag for end of layer
my $type = "depth"; # Type of data for second line in block: depth, upper or lower
# velocities
my $missXmin = 0; # Flag that gets set, is any xmin value for v is missing
# Clean up old data
$model->{layers} = (); # Remove any old layers (because of DESTROY function of layer-object)
$model->{space}->delete('VNODE');
$model->{space}->delete('ANNOTVNODE');
while(<FILE>) {
chomp;
next unless length; # Don't cry if line is empty
$l += 1;
($pre, my $tmp) = unpack "A2xA*";
#print "(D) String $tmp contains of ".(length($tmp)/7)." packages \n" if $debug;
if ( (length($tmp)/7) =~ m/^\d+.\d+$/ ) {
print
"\n\n\nWARNING!!! Something is wrong with your v.in file in line $.\n"
."I have problems to recognize the blocks. Please do not use\n"
."tabs to deliminate the numbers. Use spaces.\n\n\n";
die;
}
# @line = split (' ', $tmp);
# @line = unpack "A7" , $tmp;
@line = unpack "(A7)".(length($tmp)/7) , $tmp;
#print "(D) v.in $. type $type: ". join (',', @line)."\n" if $debug;
# Blocks of three lines belong together
# First line holds layernumber, profilekm
if ( $l%3 == 1 ) {
# Line holds kilometer. If one layer has more nodes than fit
# into one line, they still get pushed into the same array
push @km, @line;
# If layernumber has increased, create a new layer
if ( $pre > $currentlayer ) {
$currentlayer += 1;
my $c = $model->{phasecolors}{sprintf('%d2',$pre-1)};
# Layer color is determined by the reflected phasecolor from upper layer
unless (defined $c) {
#print "(I) There is no specific color defined for the boundary of layer B$pre.\n I'm using black."
#."If you want the layer boundary in a different color \n";
$c = 'black';}
$layer = new layer("model" => $model, "number" => $pre, "color" => $c);
#print "Created a new Layer (B $pre). My Profile is $xmax km long\n";
#print "(B $pre) Color: $c\n";
}
# Check if end of profile is reached. km cannot be larger then
# the whole profile
if ( $line[-1] == $xmax ) {
$v += 1;
$n = @km;
#print "End of layerkm for $pre, $n elements in km\n";
} elsif ($line[-1] > $xmax) {
print "Somethings wrong with your v.in-File. xmax = $xmax, but km = $line[-1]\n";
}
}
# Second line of block holds depth/velocity for above km
if ($l%3 == 2 ) {
push @d, @line;
#print "@line \n";
}
# Third line of block switches partial derivatives on/off
if ($l%3 == 0 ) {
push @par, @line;
# Is depth array as large as km array?
if ( $n == @d ) {
# Amount of nodes in @d is the same as @km. And end of layer is reached
# This is the last line, as km and d array have the same length
#print "End of second line in block is reached\n";
# is nodetype == depth?
if ( $type =~ "depth" ) {
#print "(D) This is last depth of layer @km, @d\n" if $debug;
#print "(D) addNodes vukm => \n@km, depth => \n@d, depthpar => \n@par\n" if $debug;
$layer->addNodes("km" => \@km, "depth" => \@d, "depthpar" => \@par);
@km = (); @d = (); @par = (); $n = 0; # Reset variables for next layer
$type = "upVelo"; # Flag velocities for next line
next;
}
## This shouldn't really happen. rayinvr wouldn't work. But just in case...
## Test if number of km- and depth-nodes fit
#if (@km != @d){
#print "Length km: ".@km.", Length d: ".@d."\n";
#die "Number of nodes for km and depth don't fit!!\n";
#}
#next;
if ($type =~ "upVelo") {
#print "velocity upper boundary\n";
#print "k @km \n";
#print "l @line \n";
#print "d @d\n";
#print "$currentlayer: readVin Partial derivatives upper layer >@par<\n";
#print "(D) addVNodes vukm => \n@km, vu => \n@d, vupar => \n@par\n" if $debug;
# Set missXmin-flag if no value for xmin found
#if ( $km[0] != $xmin){
#$missXmin = 1;
#unshift @km, $xmin;
#unshift @d, $d[0];
#unshift @par, $par[0];
#print "(I) Added vl[km $xmin] nodes to $layer->{number}\n";
#}
$layer->addVNodes("vukm" => \@km, "vu" => \@d, "vupar" => \@par);
@km = (); @d = (); @par = (); $n = 0;
$type = "loVelo";
next;
}
if ($type =~ "loVelo") {
#print "velocity lower boundary \n";
#print "@km \n";
#print "@line \n";
#print "$currentlayer: readVin Partial derivatives lower layer >@par<\n";
#print "(D) addVNodes vlkm => \n@km, vl => \n@d, vlpar => \n@par\n" if $debug;
# Set missXmin-flag if no value for xmin found
#if ( $km[0] != $xmin){
#$missXmin = 1;
#unshift @km, $xmin;
#unshift @d, $d[0];
#unshift @par, $par[0];
#print "(I) Added vl[km $xmin] nodes to $layer->{number}\n";
#}
$layer->addVNodes("vlkm" => \@km, "vl" => \@d, "vlpar" => \@par);
push @{$model->{layers}}, $layer;
$type = "depth"; # Next line is for depth again
@km = (); @d = (); @par = (); $n = 0;
}
} else {
#print "Number of second line elements hasn't reached km n $n d @d\n";
}
}
} #end of FILE-while
close(FILE);
# Adding last layer
$layer->lastLayer("km" => \@km, "depth" => \@d);
push @{$model->{layers}}, $layer;
# Report missing xmin values!
if ( $missXmin == 1 ) {
my $t = "You are working with a v.in file that has no values for xmin. rayinvr works with this, but vin2xyz has problems. I'm adding nodes for xmin. Please write model to file before running dmplstsqr";
my $m = $model->{mainwindow}->Dialog(
-title => "Velocity nodes in v.in ",
-text => $t,
#-width => 75,
-wraplength => '6i',
-buttons => ['Ok']
);
$m->Show;
}
#die;
# Name layers:
my %layernames = ( ' 1' => 'Water', ' 2' => 'Sediment', ' 3' => 'Sediment',
' 4' => 'Sediments ?', ' 5' => "Basement", ' 6' => "Upper crust",
' 7' => 'Middle crust', ' 8' => 'Lower crust', ' 9' => "Mantle");
#print Dumper \%layernames;
# Drawing layers
foreach my $layer (@{$model->{layers}}){
$layer->drawLayer;
# TODO: change this naming behaviour
if (exists $layernames{$layer->{number}}) {
$layer->{name} = $layernames{$layer->{number}};
#print "Name layer $layer->{number} $layer->{name}\n";
} else {
#print "Can't find a name for >$layer->{number}<\n";
}
}
$model->writeV2xzyIn(); # Convert v.in to an vin2xyz-readable format
$model->_readFin();
}
sub _readFin {
=PROGhead2 C<model::_readFin()>
Reads floating reflector file f.in
=cut
printf "(T) %s(@_)\n", commons::whoami() if $tree;
my ($model, %args) = @_;
my $cns = $model->{space};
my $xmax = $model->{xmax};
my $file = "f.in";
say "Reading file $file";
my $pre; # prefix of line. holds layernumber or shows that line is contiuned
my $layer; # current layer as layer-object
my $type='amount' ; # Flag linetype
my @km; # stores profilekm for more than one line
my @d; # stores depth corresponding to @km
my @par; # stores partial derivatives
my $n = 0; # Number of nodes for a layer. Used as flag for end of layer
open(FILE, $file) or (print "Can't open $file", return);
while (<FILE>) {
chomp;
#s/^\s+//; # no leading white
#s/\s+$//; # no trailing white
next unless length; # Don't cry if line is empty
#my @line = split;
# First line defines the number of nodes in the floating reflector
if ( $type eq 'amount' ) {
my @line = split;
$n = $line[0];
$type = 'km';
$layer = new floatingReflector("model" => $model,
"number" => 'Float', "color" => "black");
next;
}
# Input file: f.in contains the "floating" reflecting boundaries;
# these are interfaces with no associated velocity discontinuity
# (hence the word "floating") at which rays can reflect in
# addition to the layer boundaries; the file format is similar to
# that for v.in:
# (1) the number of nodes defining the floating reflector
# (format: I2)
# (2) the x-coordinates (km) of the nodes defining the floating
# reflector (format: 3X, <number of nodes>F7.2)
# (3) the z-coordinates (km) of the nodes defining the floating
# reflector (format: 3X, <number of nodes>F7.2)
# (4) a 0 or 1 for each node listed above depending on whether (1)
# or not (0) partial derivatives are to be calculated for
# this particular node (format: 3X, <number of nodes>I7)
# Lines (1) through (4) are repeated for each floating reflector
my $line = $_;
my @line = unpack "A3".("A7" x $n);
#print "line: @line\n" if $debug;
# Remove line prefix
shift @line;
if ( $type eq 'km' ) {
if (@line == $n) {
# Got all nodes I expect for this reflector
@km = @line;
$type = 'depth';
next;
}
}
if ( $type eq 'depth' ) {
@d = @line;
$type = 'par';
next;
}
if ( $type eq 'par' ) {
@par = @line;
#print "(D) addNodes \n"
#."vukm => @km\n"
#."depth => @d\n"
#."depthpar => @par\n" if $debug;
$layer->addNodes("km" => \@km, "depth" => \@d, "depthpar" => \@par);
@km = (); @d = (); @par = (); $n = 0; # Reset variables for next layer
$type = 'amount';
$layer->drawLayer;
push @{$model->{layers}}, $layer;
next;
}
}
close(FILE);
}
sub _readRays{
=PROGhead2 C<model::_readRays()>
Large routine reading C<r1.out> and C<r2.out>. Before reading any data, old
ray-information is deleted from model-objects. Then phase information about rays
is read from C<r1.out> then data is read from C<r2.out> and connected with phases
from C<r1.out>.
Information from C<r1.out> is stored invariable C<@r1> for later use, when reading
C<r2.out>
=cut
printf "(T) %s(@_)\n", commons::whoami() if $tree;
my ($model, %args) = @_;
my $rayblock = 0; # Flag. If 1 than first part of file is read and rayinfo begins
my $r1ray = -1; # raynumber in r1.out. should fit $raynumber from r2.out
my @r1; # Stores file r1.out, which keeps information about the number
# of points in a ray.
my $phase;
my $linenumber = 0;
my $error = "";
my $results = ""; # Store shot, phase and model information
# TRY only deleting rays when the phase was recalculated. Otherwise keep informations
## Delete old rays
# deleteRays
#foreach (sort {$b <=> $a}(keys(%{$model->{stations}}))) {
## Safety maeasure for unblessed stations
#if (ref $model->{stations}{$_} eq "station") {
#$model->{stations}->{$_}->emptyRays;
#} else {
#print "WARNING !! There's an unblessed station in your model!! Call the priest !!\n";
#}
#}
####################################################################
# READ R1.OUT #
####################################################################
open (RAYS, "r1.out") or do {print "Can't open r1.out\n"; return 0};
#shot ray i.angle f.angle dist depth red.time npts code
# -1 1 1.838 -4.103 0.000 1.31 10.327 4 1.2
# -1 2 2.865 -5.802 0.000 0.74 10.359 5 1.2
#
# If npts is 500, ray tracing wasn't successfull
# "dist" fits with last lines value for "x" for one ray in r2
# "i.angle" = 90 - r2(ang2); fist line for one ray
# "f.angle" = 90 - r2(ang1); last line for one ray
my $status = 0; # Has rayinvr finished? This status is set to one when
# rayinvrs output reached |----, which is followed by a tracing
# summary. If you have a bad model and rayinvr can't run through
# this status will stay at 0.
my ($sh, $raynr, $i, $f, $dist, $d, $redtime, $npts, $ph);
$model->printStatusMessage("Reading r1.out..");
print "Reading r1.out\n" unless $quiet;
while (<RAYS>){
$linenumber += 1;
chomp;
my $line = $_;
#my @line = unpack("A4A4A9A9A9A8A8A6A6",$_);
# 0 1 2 3 4 5 6 7 8
($sh, $raynr, $i, $f, $dist, $d, $redtime, $npts, $ph) = unpack("A4A4A9A9A9A8A8A6A6",$line);
$ph =~ s/^\s+|\s+$//g; # remove trailing and leading spaces for phasecode
# Read Phase numbers
unless (defined $ph) {print "No phasecode defined >$ph<\n"; next;}
if ($ph =~ /^\d\.?/){ # Only if there's a phasecode at the phasecodeposition
# Special treatment for floating reflectors
# TODO
#if ($ph =~ /^\d\.4/ ) {
#$ph =~ s/4$/2/;
#print "A floating reflector has been traced. Change phase to $ph\n";
#}
#print "Push phase $ph\n";
# raynr npts code/phase
push @r1, [$raynr, $npts, $ph, $linenumber];
#print "r1.out: $sh, $raynr, $i, $f, $dist, $d, $redtime, $npts, $ph\n";
} else {
print "r1.out: $line\n" unless $quiet;
$results.="$line\n";
# Find shot kilometer
$line =~ s/^\s+//; # no leading white
my @line = split /\s+/, $line;
# Find station for this km, if first element fits the required
# format
if ($line[0] && $line[0] =~ /\d\.\d/) {
#print "Line array: $line[0] ";
my $st = $model->getStation($line[0]);
#print "Station = $st->{name}\n";
chomp $results;
$results .= " st $st->{name}\n";
}
}
if ($sh eq '|---') {
#print "End of rays reached. Leaving loop .. \n";
$status = 1;
#last;
#$model->printStatusMessage("\nThe model is talking to you!! ");
}
if ( grep(/Number of data points used/, $line ) ) {
$model->printStatusMessage("\n$line");
}
if ( grep(/RMS/, $line ) ) {
$model->printStatusMessage(", $line");
}
if ( grep(/Normalized/, $line )) {
$model->printStatusMessage(", $line");
}
# Update Splash screen during start up
if ( $model->_get("init") == 1 && $linenumber % 200 == 0 ) {
#print "Update Splash\n";
$model->_get("mainwindow")->update();
}
}
close(RAYS);
unless ($linenumber > 0 ) {
print "ERROR: File r1.out is empty!!!\n";
$model->printStatusMessage("\nERROR: File r1.out is empty!!! Check command line output for rayinvr error message");
return 1;}
print "r1.out read and closed\n" unless $quiet;
$model->{results} = $results; # Store modeled results
# Check if there have been any errors during the rayinvr run
if ($status == 0) {
# Find last successfully traced station
my $k = $model->{rin}->{xshot}[abs($sh)-1];
my $s = $model->{stations}->{$k}->{name};
print "------------------------------------------------------------------------------------\n";
print "| WARNING!! RAYINVR HASN'T FINISHED \n";
print "| RAYS WILL BE MISSING IN DISPLAY \n";
print "| LAST SHOT TRACED: SHOT$sh, RAY$raynr, DISTANCE$dist, DEPTH$d, PHASE $ph \n";
print "|$sh, $raynr, $i, $f, $dist, $d, $redtime, $npts, $ph \n";
#4, 4, 79.256,-78.006,229.842,0.00,5.047, 46, 9.1
print "| SHOT$sh = station $s at km $k \n";
print "------------------------------------------------------------------------------------\n";
$model->printStatusMessage(
"\nERROR: rayinvr could trace all stations. Traveltimes are not updated."
." See command line for more information");
return;
}
####################################################################
# READ R2.OUT #
####################################################################
open (RAYS, "r2.out") or return "Can't open r2.out. No rays can be displayed. Run rayinvr.";
print "Reading r2.out \n" unless $quiet;
$model->printStatusMessage(" .. reading r2.out ..");
# To be written if you want to display velocity as gradient
# Block definition in r2.out
# number of blocks (nblk) is equal number of nodes in c.in
#layer# 7 nblk= 2 (ivg,x1,x2,z11,z12,z21,z22,s1,b1,s2,b2,vp1,vs1,vp2,vs2, vp3,vs3,vp4,vs4,c1,c2,...,c11)
#ivg 1 1
#x1 0.0000 37.9700
#x2 37.9700 750.0000
#z11 10.0000 10.0000
#z12 10.0000 10.0000
#z21 20.0000 20.0000
#z22 20.0000 20.0000
#s1 0.0000 0.0000
#b1 10.0000 10.0000
#s2 0.0000 0.0000
#b2 20.0000 20.0000
#vp1 7.0000 7.0000
#vs1 4.0415 4.0415
#vp2 7.0000 7.0000
#vs2 4.0415 4.0415
#vp3 7.3000 7.3000
#vs3 4.2147 4.2147
#vp4 7.3000 7.3000
#vs4 4.2147 4.2147
#-0.477E-05 0.477E-05
#0.000E+00 0.000E+00
#0.114E+02 0.214E+03
#0.477E-06-0.477E-06
#0.254E+04 0.477E+05
#0.000E+00 0.000E+00
#0.380E+03 0.712E+04
#0.000E+00 0.000E+00
#0.000E+00 0.000E+00
#0.181E-03-0.340E-02
#-0.181E-02 0.340E-01
$linenumber = 0;
my $currentStation; my $currentray = 0; my $currentgroup = 1; my $currentphase = 0;
my $r2ray; my $pointnumber; my $group;
my $km;
$i=0;
my $numLays; # number of layers, from line three in r2.out
my $l = 0; # flags blocks of layer x
my $nblk; # number of blocks within a layer
my $skip = 0; # flag to skip one line
my $unpack = "";
my $ll = 29; # number of lines for layer blocks
# if layer has more than hundred nodes, block spread
# over two lines
my @blocks; # save data in matrix
my $mod = -1; # Needed for layer block lines spreading over two lines in
# file. It saves, if new line starts on even or odd linenumber
while (<RAYS>){
$linenumber += 1;
chomp;
# Update Splash screen
if ( $model->_get("init") == 1 && $linenumber % 200 == 0 ) {
#print "Update Splash\n";
$model->_get("mainwindow")->update();
}
##################################################################
# Reading block information
##################################################################
if ($rayblock == 0) {
# Find the total number of layers in this model
if ($linenumber == 3) {
my @line = split /=/;
$numLays = $line[-1];
#print "Number of layers: $numLays\n";
next;
}
# Find place with raypoints
elsif ($_ =~ m/^gr*/) {
#print "Found start of rays at $linenumber: $_\n";
$rayblock = 1;
next;
}
# Find description of modelblocks and set up variables for block
# reading
elsif ($_ =~ m/^layer/ && $l >= 0 ) {
my @line = unpack("A6A2A7A4",$_);
#my @line = split;
$l = $line[1];
$i = -1; # gives linenumber within layer-block
$nblk = $line[3];
#print "Layer $l, number of blocks $nblk\n";
#create unpack string from nblks
$unpack = "";
my $c = $nblk;
if ($nblk > 100) { # POSSIBLE BUG IF MORE THAN 200 HUNDRED
# NODES ARE USED
$c = 100;
#$ll = 59;
$ll = 29;
$mod = $linenumber%2;
#print "For this layer $l, nblk $nblk mod: $mod is too much line: $linenumber\n";
} else {
$ll = 29;
$mod = -1;
};
for (my $i = 0; $i < $c; $i++){
$unpack .= "A10";
}
#print "Unpack format: $unpack\n";
next;
}
# Read Blockinformation for a layer
elsif ($l > 0) {
my @line = unpack($unpack,$_);
#print "$linenumber: $i: @line\n";
$i += 1;
if ($mod == $linenumber%2) {
#print "Continue, paste line to previous\n >@line[-1]<\n";
$i-=1; # This line is sequel of previus one
# remove empty arrayplaces
@line = grep {$_ ne ''} @line;
#print "Sauber >@line[-1]<\n";
}
######################################
#if ($l == 2 ) {
#if (defined @{$blocks[$i]}) {print "L2 in array @{$blocks[$i]}\n";}
#print "L2 $i: @line\n";
#}
######################################
push @{$blocks[$i]}, @line;
if ($i == $ll ) { # if more than hundred nodes, it needs two lines
#print "End of layer $l, blocks $nblk, lines ".@blocks." \n";
$model->_addBlocks( "layer" => "$l", "blocks" => \@blocks );
@blocks = (); # Clear array for next line
$mod = -1; # Reset flag for long lines
#print "End of layer $l\n";
if ($l == $numLays) { # All layers have been read once
# more layer information is not neccessary at this
# stage of software
#print "All layerblocks read\n";
$l = -1;
}
}
next;
}
}
next unless $rayblock == 1;
###################################################################################################
# Reading raypathes
#####################################
# Ray is the same number as ray-column in r1. If ray equals zero, ray doesn't reach
# the surface.
# rays are numbered in column "ray". If that number increase, a new ray starts
# if npt is 500, ray tracing wasn't successfull
# gr ray npt x z ang1 ang2 v1 v2 lyr bk id iw
# 1 0 1 20.285 2.000 0.00 -88.69 0.00 1.58 1 2 -1 1
# 1 0 2 18.870 2.029 -88.93 -88.93 1.58 1.58 1 1 -1 1
# Problem occur if you use "split" because two columns can merge if numbers
# become too large
# "gr" is raygroup and means a group of rays with the same code, e.g. code 1.2 is group 1
# code 2.2 is group 2
my @line;
# rayinvr and xrayinvr have different output formats in r2.out.
# --> That seems to be a patch. Don't know where output-format is defined in rayinvr-source
#if ( $model->{config}->{r2out} =~ m/[0-9]/ ) {
# print "Config r2out <$model->{config}->{r2out}< matches a number!\n";
my $pack;
# Automatically check the r2out value!!
if ( length $_ == 67 ) {
$pack = "A2A3A4A8A8A8A8A7A7A3A3A3A3";
} elsif ( length $_ == 69 ) {
$pack = "A4A3A4A8A8A8A8A7A7A3A3A3A3";
}else{
# If none is correct, just skip this line
print "Can't understand the format of r2.out. This is a programming error!";
next;
#die;
}
#@line = unpack("A".$model->{config}->{r2out}."A3A4A8A8A8A8A7A7A3A3A3A3",$_);
@line = unpack($pack,$_);
#} else {
#print "Don't know your configfile-entry for r2.out! Got >$model->{config}->{r2out}<.\n".
#"It should be the number of digits for first column of r2.out\n";
#die;
#}
$r2ray = $line[1]; # equals $r1[0]
$pointnumber = $line[2]; # corresponds to (can't be bigger than) $r1[1]
$group = $line[0]; # kind of correspondend to raycode in r1
$km = sprintf "%.3f", $line[3];
#$km =~ s/^\s+//; # Remove leading whites
# Check the format of group number. If too many rays are traced
# it changes to '**'
if ($group !~ m/\d+/){
print
"\n\n\n________________________________________________________________\n"
."\nERROR\n"
#."Gr >$group< is NO digit!\n"
."Line $linenumber: @line\n"
."Last station read: Station $model->{stations}->{$km}->{name} at "
."km $model->{stations}->{$km}->{position}, phase $currentphase\n";
print "Your traced rays are too big to fit into r2.out. This is"
." a rayinvr problem and need to be solved there.\nYou need"
." to find a way, to reduce the number of rays (also failed)"
." tracings\n"
."________________________________________________________________\n\n\n"
;
$model->printStatusMessage("\nYour r2.out is overflowed. Last station $model->{stations}->{$km}->{name} "
."km $model->{stations}->{$km}->{position}, phase $currentphase");
return;
}
# Second column in r2.out "ray" seems to be zero if raypath is not valid
# Those rays are not listed in r1 and hence need to be ignored
# You can have error messages in r2.out. They are sorted out with the $group-test
# to be a digit.
if ($rayblock == 1 && $group =~ m/\d+/ && $r2ray != 0) {
# A new ray starts
if ($r2ray != $currentray || $group != $currentgroup){
#print "New ray starts\n";
#print "r2ray != currentray || group != currentgroup)";
#print "$r2ray != $currentray || $group != $currentgroup)\n";
#print "line:$. @line\n";
# If it's a new station, delete old rays in that station
if ( ! $currentStation || $currentStation != $model->{stations}->{$km} ){
#print "Delete rays for station $model->{stations}->{$km}->{name} at km $km\n" if $debug;
#$model->{stations}->{$km}->emptyRays;
$model->getStation($km)->emptyRays;
}
#print "station $model->{stations}->{$km}->{name} at km $km\n";
$currentray = $r2ray;
$currentgroup = $group;
$currentStation = $model->{stations}->{$km};
$r1ray+=1;
$currentphase = $r1[$r1ray][2];
#print "Phase is $currentphase. Previous phase is $phase\n" if ($phase != $currentphase);
#$currentStation->emptyRays($currentphase) if ($phase != $currentphase);
$phase = $r1[$r1ray][2];
#print "r2.out:$linenumber: gr $group ray $r1ray, nbr: $r2ray, cur $currentray,".
#" ph: $phase , km <$km> <@line\n";
if (!defined $currentStation) {print "?? Which station? km $km in line $. of r2.out: @line\n"};
}
if (!defined $r1[$r1ray][0] && $status == 0){ # there's no information in r1 for this ray
# and rayinvr hasn't finished
print "WARNING !! RAYINVR DIDN'T RUN THROUGH THE WHOLE TRACING\n";
last;
}
# Check if r2rays fit:
if ( $r1[$r1ray][0] != $r2ray){
print "________________________________________________________________\n".
"\nERROR\n".
"r2rays of r1 and r2 don't fit!! index $r1ray $r1[$r1ray][0] != $line[1]\n".
"Station $currentStation->{name}, ray $currentray\nr2 line $linenumber: @line\n".
"r1 raynr, npts, ph, linenumber\n".
"r1 @{$r1[$r1ray]}\n".
"________________________________________________________________\n\n".
" What to do?\n".
" - Check if rayinvr has finished without problems\n".
" - PRay cannot recognize segmentation faults of rayinvr\n".
" - if rayinvr has'nt finished, check v.in for crossing layers, large gradients, ..\n".
" (use vmodel) and/or go to last working version and delete r2.out\n".
" - Check if format of r2.out matches your configuration:\n".
" - get the number of digits of first column of raypart of r2.out\n".
" - p.pl is expecting >$model->{config}->{r2out}< digits\n".
" - Use 'r2out=NUMBER' in p.config to match r2.out with p.pl\n".
"_________________________________________________________________\n".
"\n\n";
$model->printStatusMessage("\nERROR: rayinvr has problems to trace all data. See command line output for more information");
close(RAYS);
return 1;
#die;
}
unless ($currentStation) {
print "WARNING no station for ray: @line\n";
$error = "\nERROR: Couldn't find station for some rays";
next;
}
# Use real raycodes as 'phase'. e.g 1.2 instead of 12. In older
# versions, this has been different and only the tx-code was used
# to sort in rays
# $phase = '5.1'
$currentStation->addRays("phase" => $phase, "raynumber" => $r1ray,
"rays" => [$km, $line[4]] );
} # if it is a ray
} # end of while-loop for reading file
close(RAYS);
print "All rays read\n" unless $quiet;
$model->printStatusMessage(".. Rays are read.");
#die;
$model->printStatusMessage($error);
}
sub _addBlocks {
=PROGhead2 C<model::_addBlocks()>
Add blocks read from C<v.in> to model
=cut
printf "(T) %s(@_)\n", commons::whoami() if $tree;
# Finds corresponding layer object and adds blocks there
my ($model, %args) = @_;
#print "model::addBlocks\n";
my $layer = $args{layer};
my $blocks = $args{blocks};
#print "Layer $layer, blocks $blocks\n";
if (defined $model->{layers}[$layer-1]) {
$model->{layers}[$layer-1]->addBlocks($blocks);
} else {
warn "Undefined value for layer $layer-1\n";
}
}
sub _readVgrid {
=PROGhead2 C<model::_readVgrid()>
Creates grid file from v.in, reads grid file and creates little colored
squares to generate a smooth gradient.
vin2xyz -R$xmin/$xmax/$zmin/$zmax -x$dx -z$dz < v2xyz.in > v.grid
=cut
printf "(T) %s(@_)\n", commons::whoami() if $tree;
my ($model, %args) = @_;
print "model::_readVgrid\n";
$model->{space}->delete('VGRID');
# Create gridfile with rayinvr-skript vin2xyz
my $xmin = $model->{xmin};
my $xmax = $model->{xmax};
my $zmin = $model->{zmin};
my $zmax = $model->{zmax};
my $dx = 1;
my $dz = 0.5;
my $command = "vin2xyz -R$xmin/$xmax/$zmin/$zmax -x$dx -z$dz < v2xyz.in > v.grid";
print "$command\n";
system ($command);
# Read in vgrid
print "Gridfile created\n";
my $file = "v.grid";
say "Reading file v.grid";
my $cns = $model->{space};
open(FILE, $file) or do {$model->{vgrid} = 0; die "Can't open $file"};
while (<FILE>){
chomp;
my ($x, $z, $v) = split;
my $color = _getColor($v);
my @scaled = $model->model2screen([$x-$dx/2, $z-$dz/2, $x+$dx/2, $z+$dz/2], "space");
$cns->createRectangle(@scaled, -fill => $color, -width=>0, -tags=>["VGRID", "km $x, d $z, v $v"]);
}
close (FILE);
unless ($cns->find('withtag', 'VGRID')) {
$model->{vgrid} = 0;
$model->printStatusMessage("\nError in creating the velocity grid. Check the command manually:\n$command");
}
$model->order;
print "model::_readVgrid() is done\n";
}
sub tomoReadGrid{
=PROGhead2 C<model::tomoReadGrid()>
Reads and displays grid file from Tomo2D
=cut
printf "(T) %s(@_)\n", commons::whoami() if $tree;
my ($model, %args) = @_;
print "model::tomoReadGrid\n";
$model->{space}->delete('VGRID');
my $file = $model->{config}->{tomoGrid};
my $cns = $model->{space};
my $dx = 1;
my $dz = 0.5;
# CLEANUP
$cns->delete('TOMOGRID');
$cns->delete('TOMOREFL');
print "This is tomogrid trying to open $file\n";
open(FILE, $file) or (print("Can't open $file\n"), return 0);
print "This is tomogrid reading $file\n";
my $line = 0;
while (<FILE>){
chomp;
my ($x, $z, $v) = split;
#print "Reading line $line: $x $z $v\n";
my $color = _getColor($v);
my @scaled = $model->model2screen([$x-$dx/2, $z-$dz/2, $x+$dx/2, $z+$dz/2], "space");
$v = sprintf "%5.2f", $v;
$cns->createRectangle(@scaled, -fill => $color, -width=>0, -tags=>["TOMOGRID", "km $x, d $z, v $v"]);
$line++;
}
close (FILE);
$model->{tomoGrid}=1;
# Read Reflector file
$model->_readxzfile($model->{config}->{tomoRefl});
$model->_readxzfile('bathyxz.dat');
$model->_readxzfile('basement.dat');
return 1;
}
sub _readxzfile {
=PROGhead2 C<model::readxzfile($file [, diagram])>
Reads given file and creates a line in model space with tags C<'XZ'>.
Calls C<< $model->order >>.
=cut
printf "(T) %s(@_)\n", commons::whoami() if $tree;
# Read and plot given file
my ($model, $file, $diagram) = @_;
$file =~ s/^\s+//; # no leading white
$file =~ s/\s+$//; # no trailing white
#print "read xz 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
<code>$model->{stations}->{$_}->emptyRays;</code> </li>
<li> Get version from file 'tomoversion' and update window title </li>
<li> Read rays and traveltimes
If 'tomoRaysPg' is defined, then Pg and PmP are read from different
ray and tx-files
<code>
$model->tomoReadRays('file' => $model->getConfig('tomoRaysPmP'),
'ph' => $model->getConfig('tomoPhasePmP'));
$model->_readTx('file' => $model->getConfig('tomoTimesPg'), key => 'txTomo');
</code></li>
<li> Read tomo grid
<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
<code>$model->_readTx( file => $file, key => 'txin');</code>
=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
<code>$model->_readTx( file => $file, key => 'txout');</code>
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
tx-times for tomo are stored in <code>$station->{$key}</code>
format: <code>push @{$data{$pos}{$richtung}{$ph}}, [ $km, $t, $unc ];
$station->{$key}->{$richtung}->{$ph} = @{$data{$pos}{$richtung}{$ph}}
</code>
=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

Add Comment

Please, Sign In to add comment

Public Pastes

📌 Swapzone +37% glitch ⭐ B
JavaScript | 5 sec ago | 0.25 KB
📝 Crypto Swap Glitch ✅ Working
JavaScript | 9 sec ago | 0.24 KB
✅⭐ Make huge profits on trading ⭐⭐ 7
JavaScript | 15 sec ago | 0.25 KB
📝 EASY MONEY GUIDE ✅ Working
JavaScript | 18 sec ago | 0.24 KB
⭐✅ Marketplace Glitch ✅ Working ✅ NEVER SEEN...
JavaScript | 26 sec ago | 0.25 KB
📌 ChangeNOW Exploit
JavaScript | 27 sec ago | 0.24 KB
⭐✅ Exploit 2500$ in 15 Minutes⭐⭐⭐ 1
JavaScript | 34 sec ago | 0.25 KB
📌 Instant BTC Profit Method ✅ Working
JavaScript | 38 sec ago | 0.24 KB

create new paste / syntax languages / archive / faq / tools / night mode / api / scraping api / news / pro
privacy statement / cookies policy / terms of service / security disclosure / dmca / report abuse / contact

By using Pastebin.com you agree to our cookies policy to enhance your experience.
Site design & logo © 2025 Pastebin

We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand

Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!

HOME
ABOUT
AUCTIONS
SHIPPING
FEES
TOOLS
HOW
FAQ
CONTACT

Original Source | Taken Source