source: SH_SHM/trunk/util/locate_by_ttfit.csh @ 263

Revision 263, 8.4 KB checked in by marcus, 12 years ago (diff)

r136 | walther | 2010-08-20 09:29:00 +0200 (Fr, 20 Aug 2010) | 1 line

Adopted script to use GMT prefix if available.

  • Property svn:executable set to *
Line 
1#! /bin/csh
2#
3# file locate_by_ttfit.csh
4#      ===================
5#
6# $Revision: 136 $, $Date: 2010-08-20 09:29:00 +0200 (Fr, 20 Aug 2010) $
7#
8# Locate event via travel time fit.
9# K. Stammler, 5-Feb-2001
10
11if  ("$1" == "")  then
12        echo "Usage: $0 <evtfile> [<lat>] [<lon>] [<width>] [<stepsize>] [<dep1>] [<dep2>] [<ddep>] [<weight>] [<method>]"
13        exit
14endif
15
16# check for GMT command access. If GMT is defined, use it as prefix for
17# all commands.
18set testGMT=`which GMT`
19set GMT=""
20[ -x "$testGMT" ] && set GMT=$testGMT
21
22#set echo
23
24# get parameters
25set evtfile=$1
26set startlat=$2
27set startlon=$3
28set width=$4
29set stepsize=$5
30set dep1=$6
31set dep2=$7
32set ddep=$8
33set weight=$9
34set method=$10
35
36if  ("$width" == "")  set width=5
37if  ("$stepsize" == "")  set stepsize=0.5
38if  ("$weight" == "")  set weight=1.0
39if  ("$method" == "")  set method=simplex
40
41# set constants
42set wlat=$width
43set wlon=$width
44set dlat=$stepsize
45set dlon=$stepsize
46
47set resfile=locate_out.evt
48set relfile=ttfit_rel.000
49set fitxyz=ttfit_xyz.000
50set psfile=ttfit.ps
51set textresult=ttfit_edit.txt
52set cpt=ttfit.cpt
53set citycol=200/0/150
54#set fmt=-Jx1.0
55set imagsize=5
56set fmt=-JX$imagsize/$imagsize
57#set fmt=-JM$imagsize
58
59chdir $SH_SCRATCH
60# set units for plot, want to change this to cm
61$GMT gmtdefaults -Du >.gmtdefaults
62
63# delete previously existing result file
64if  (-e $resfile)  \rm $resfile
65
66set lab=1
67set res=`$SH_UTIL/floatop $width gt 5.0`
68if  ($res == 1)  set lab=2
69set res=`$SH_UTIL/floatop $width gt 10.0`
70if  ($res == 1)  set lab=5
71
72if  (! -e $evtfile)  then
73        echo "evt file $evtfile not found.  Abort."
74        exit
75endif
76
77# check number of phases
78set res=`grep -c '^Phase name' $evtfile`
79if  ($res < 4)  then
80        echo "$0 : need at least 4 phases.  Abort."
81        exit
82endif
83
84# find main phase
85set mainphase=""
86set tmp=`grep -c '^Phase name             : Pn$' $evtfile`
87if  ($tmp > 3)  set mainphase=Pn
88set tmp=`grep -c '^Phase name             : Pg$' $evtfile`
89if  ($tmp > 3)  set mainphase=Pg
90set tmp=`grep -c '^Phase name             : PP$' $evtfile`
91if  ($tmp > 3)  set mainphase=PP
92set tmp=`grep -c '^Phase name             : P$' $evtfile`
93if  ($tmp > 3)  set mainphase=P
94# check for PKP
95if  ("$mainphase" == "")  then
96        set tmp=`grep -c '^Phase name             : PKPdf$' $evtfile`
97        if  ($tmp > 4)  set mainphase=PKPdf
98        set tmp=`grep -c '^Phase name             : PKPab$' $evtfile`
99        if  ($tmp > 4)  set mainphase=PKPab
100        set tmp=`grep -c '^Phase name             : PKPbc$' $evtfile`
101        if  ($tmp > 4)  set mainphase=PKPbc
102endif
103
104echo mainphase $mainphase
105
106# create relative onset file
107if  (-e $relfile)  \rm $relfile
108$SH_UTIL/evt2reltimes $evtfile $relfile $mainphase
109
110# read preliminary location
111if  ("$startlat" == "")  then
112        set lat=`grep latitude $relfile`
113        if  ($#lat < 3)  then
114                echo "illegal output of evt2reltimes.  Abort."
115                #\rm $relfile
116                set lat=0
117                exit
118        else
119                set lat=$lat[3]
120        endif
121else
122        set lat=$startlat
123endif
124if  ("$startlon" == "")  then
125        set lon=`grep longitude $relfile`
126        if  ($#lon < 3)  then
127                echo "illegal output of evt2reltimes.  Abort."
128                #\rm $relfile
129                set lon=0
130                exit
131        else
132                set lon=$lon[3]
133        endif
134else
135        set lon=$startlon
136endif
137
138# make nice numbers
139set lat=`echo $lat | sed 's/\./ /'`
140set lat=$lat[1].5
141set lon=`echo $lon | sed 's/\./ /'`
142set lon=$lon[1].5
143
144#set echo
145
146# call travel time fit
147if  (-e $fitxyz)  \rm $fitxyz
148if  ("$method" == "grid")  then
149        $SH_UTIL/fit_travel -infoline -xyz -clat=$lat -clon=$lon -wlat=$wlat \
150                -wlon=$wlon -dlat=$dlat -dlon=$dlon -dep1=$dep1 -dep2=$dep2 \
151                -ddep=$ddep $relfile $weight >$fitxyz
152else
153        if  ("$dep1" == "$dep2")  then
154                set fixdepth="-fixdepth"
155        else
156                set fixdepth=""
157        endif
158        set xlat=`echo $lat | sed 's/-/s/'`
159        set xlon=`echo $lon | sed 's/-/w/'`
160        $SH_UTIL/fit_travel2 -trace $fixdepth -ftol=0.001 $relfile $xlat $xlon \
161                $dep1 $weight >$fitxyz
162endif
163#\rm $relfile
164
165set infoline=`tail -1 $fitxyz`
166if  ($#infoline < 11)  then
167        echo "infoline of $fitxyz has unknown syntax:"
168        echo "$infoline"
169        exit
170endif
171set minz=$infoline[3]
172set minlat=$infoline[5]
173set minlon=$infoline[7]
174set mindep=$infoline[9]
175set origtime=$infoline[11]
176if  ("$method" == "grid")  then
177        set maxz=`grep -v '^\!' $fitxyz | grep " $mindep depth" | awk '{print $3}' | sort -n | tail -1`
178else
179        set maxz=`grep -v '^\!' $fitxyz | awk '{print $3}' | sort -n | tail -1`
180endif
181
182# check for longitude larger than 180
183set xtmp=`echo $minlon | sed 's/\./ /'`
184set xtmp=$xtmp[1]
185if  ($xtmp >= 180)  then
186        set minlon=`echo "scale=3; $minlon - 360.0" | bc`
187endif
188
189# compute bounds of display
190set wlat2=`echo "scale=3; $wlat / 2" | bc`
191set wlon2=`echo "scale=3; $wlon / 2" | bc`
192set dlat2=`echo "scale=3; $dlat / 2" | bc`
193set dlon2=`echo "scale=3; $dlon / 2" | bc`
194set lolat=`echo "scale=3; $minlat - $wlat2" | bc`
195set hilat=`echo "scale=3; $minlat + $wlat2" | bc`
196set lolon=`echo "scale=3; $minlon - $wlon2" | bc`
197set hilon=`echo "scale=3; $minlon + $wlon2" | bc`
198
199set lolat=`echo "scale=3; $lolat - $dlat2" | bc`
200set hilat=`echo "scale=3; $hilat + $dlat2" | bc`
201set lolon=`echo "scale=3; $lolon - $dlon2" | bc`
202set hilon=`echo "scale=3; $hilon + $dlon2" | bc`
203
204# create text output
205if  (-e $textresult)  \rm $textresult
206touch $textresult
207echo "Best fit for"                   >>$textresult
208echo "$origtime"                      >>$textresult
209echo "lat=$minlat lon=$minlon"        >>$textresult
210if  ("$dep1" != "$dep2")  then
211        echo "Found depth $mindep km"      >>$textresult
212else
213        echo "Depth given as $mindep km"   >>$textresult
214endif
215echo "Min-rms=$minz (orange)"         >>$textresult
216echo "Max-rms=$maxz (blue)"           >>$textresult
217echo ""                               >>$textresult
218echo "Found residuals:"               >>$textresult
219echo ""                               >>$textresult
220grep '^\!' $fitxyz | grep -v 'minsq:' >>$textresult
221# $SH_TEXTEDIT $textresult &
222
223# create postscript output
224if  (-e $psfile)  \rm $psfile
225
226# create cpt file, make the best 3% turn orange
227set levz=`echo "scale=5; $maxz - $minz" | bc `
228set levz=`echo "scale=5; $levz * 0.03" | bc `
229set levz=`echo "scale=5; $levz + $minz" | bc `
230if  (-e $cpt)  \rm $cpt
231touch $cpt
232echo "$minz 255 150 0    $levz 255 255 0"  >>$cpt
233echo "$levz 255 255 0  $maxz 0 0 255"  >>$cpt
234
235# compute symbol size
236set symsiz=`echo "scale=5; $stepsize / $width * $imagsize" | bc`
237
238set xlat=`echo $minlat | sed 's/-/s/'`
239set xlon=`echo $minlon | sed 's/-/w/'`
240set areaname="`$SH_UTIL/fereg $xlat $xlon`"
241
242touch $psfile
243
244if  ("$method" == "grid")  then
245        grep -v '^\!' $fitxyz | grep " $mindep depth" | \
246                $GMT psxy -: -R$lolon/$hilon/$lolat/$hilat \
247                $fmt -Ss$symsiz -C$cpt -U \
248                -B${lab}:longitude:/${lab}:latitude::."$areaname":8 \
249                -K >>$psfile
250else
251        grep -v '^\!' $fitxyz | \
252                $GMT psxy -: -R$lolon/$hilon/$lolat/$hilat \
253                $fmt -Ss$symsiz -C$cpt -U \
254                -B${lab}:longitude:/${lab}:latitude::."$areaname":8 \
255                -K >>$psfile
256endif
257
258# recompute longitude from negative values to 180..360
259set xhilon=`echo $hilon | sed 's/-/m/'`
260set res=`$SH_UTIL/floatop $xhilon lt 0.0`
261if ($res == 1)  set xhilon=`echo "scale=3; $hilon + 360.0" | bc`
262set xlolon=`echo $lolon | sed 's/-/m/'`
263set res=`$SH_UTIL/floatop $xlolon lt 0.0`
264if ($res == 1)  set xlolon=`echo "scale=3; $lolon + 360.0" | bc`
265
266# plot coastlines and borders
267$GMT pscoast -O -R$xlolon/$xhilon/$lolat/$hilat $fmt -A10 \
268        -Df -N1/4/0/0/0/dashed -W1/0/0/0/solid -Ia/0.1/100 -K >>$psfile
269
270# plot cities
271awk '{print $1,$2}' $SH_UTIL/cities.dat \
272        | $GMT psxy -: -O -K -R$lolon/$hilon/$lolat/$hilat $fmt \
273        -W5/$citycol -G$citycol -Sc0.1 >>$psfile
274awk '{print $1,$2,"14","0","4","1","\ ",$3}' $SH_UTIL/cities.dat \
275        | $GMT pstext  -: -O -K -R$lolon/$hilon/$lolat/$hilat $fmt -G$citycol >>$psfile
276
277# integrate text into psfile
278set lcnt=`wc -l $textresult`
279set lcnt=$lcnt[1]
280set ypos=0
281set cont=-K
282while  ($lcnt > 0)
283        set line="`sed -n $lcnt'p' $textresult`"
284        if  ("$line" == "")  set line="\ "
285        if  ($lcnt == 1)  set cont=""
286        echo "6.0 $ypos.0 12 0.0 8 1  $line" | \
287                $GMT pstext -O $cont -JX10/8 -R0/10/0/100 >>$psfile
288        @ lcnt = $lcnt - 1
289        @ ypos = $ypos + 2
290end
291
292touch $resfile
293
294echo "Latitude               : $minlat"               >>$resfile
295echo "Longitude              : $minlon"               >>$resfile
296echo "Depth (km)             : $mindep"               >>$resfile
297echo "Origin time            : $origtime"             >>$resfile
298echo "Location method        : relative travel times" >>$resfile
299echo "--- End of Phase ---"                           >>$resfile
300
301if  ($?PSVIEW == 0)  setenv PSVIEW pageview
302
303# pop up window
304#$PSVIEW -right $PWD/$psfile & # this dies after the shell terminates (on Sun)
305
306$PSVIEW $PWD/$psfile &
307
308#\rm $fitxyz
309
310
Note: See TracBrowser for help on using the repository browser.