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

Revision 46, 8.2 KB checked in by marcus, 14 years ago (diff)

r31 | svn | 2008-03-03 20:47:24 +0100 (Mo, 03 Mär 2008) | 1 line

first working version for reading ingres sfdb; new check for origin in 1970; fixed ps output of locate_by_ttfit

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