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 | |
---|
11 | if ("$1" == "") then |
---|
12 | echo "Usage: $0 <evtfile> [<lat>] [<lon>] [<width>] [<stepsize>] [<dep1>] [<dep2>] [<ddep>] [<weight>] [<method>]" |
---|
13 | exit |
---|
14 | endif |
---|
15 | |
---|
16 | # check for GMT command access. If GMT is defined, use it as prefix for |
---|
17 | # all commands. |
---|
18 | set testGMT=`which GMT` |
---|
19 | set GMT="" |
---|
20 | [ -x "$testGMT" ] && set GMT=$testGMT |
---|
21 | |
---|
22 | #set echo |
---|
23 | |
---|
24 | # get parameters |
---|
25 | set evtfile=$1 |
---|
26 | set startlat=$2 |
---|
27 | set startlon=$3 |
---|
28 | set width=$4 |
---|
29 | set stepsize=$5 |
---|
30 | set dep1=$6 |
---|
31 | set dep2=$7 |
---|
32 | set ddep=$8 |
---|
33 | set weight=$9 |
---|
34 | set method=$10 |
---|
35 | |
---|
36 | if ("$width" == "") set width=5 |
---|
37 | if ("$stepsize" == "") set stepsize=0.5 |
---|
38 | if ("$weight" == "") set weight=1.0 |
---|
39 | if ("$method" == "") set method=simplex |
---|
40 | |
---|
41 | # set constants |
---|
42 | set wlat=$width |
---|
43 | set wlon=$width |
---|
44 | set dlat=$stepsize |
---|
45 | set dlon=$stepsize |
---|
46 | |
---|
47 | set resfile=locate_out.evt |
---|
48 | set relfile=ttfit_rel.000 |
---|
49 | set fitxyz=ttfit_xyz.000 |
---|
50 | set psfile=ttfit.ps |
---|
51 | set textresult=ttfit_edit.txt |
---|
52 | set cpt=ttfit.cpt |
---|
53 | set citycol=200/0/150 |
---|
54 | #set fmt=-Jx1.0 |
---|
55 | set imagsize=5 |
---|
56 | set fmt=-JX$imagsize/$imagsize |
---|
57 | #set fmt=-JM$imagsize |
---|
58 | |
---|
59 | chdir $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 |
---|
64 | if (-e $resfile) \rm $resfile |
---|
65 | |
---|
66 | set lab=1 |
---|
67 | set res=`$SH_UTIL/floatop $width gt 5.0` |
---|
68 | if ($res == 1) set lab=2 |
---|
69 | set res=`$SH_UTIL/floatop $width gt 10.0` |
---|
70 | if ($res == 1) set lab=5 |
---|
71 | |
---|
72 | if (! -e $evtfile) then |
---|
73 | echo "evt file $evtfile not found. Abort." |
---|
74 | exit |
---|
75 | endif |
---|
76 | |
---|
77 | # check number of phases |
---|
78 | set res=`grep -c '^Phase name' $evtfile` |
---|
79 | if ($res < 4) then |
---|
80 | echo "$0 : need at least 4 phases. Abort." |
---|
81 | exit |
---|
82 | endif |
---|
83 | |
---|
84 | # find main phase |
---|
85 | set mainphase="" |
---|
86 | set tmp=`grep -c '^Phase name : Pn$' $evtfile` |
---|
87 | if ($tmp > 3) set mainphase=Pn |
---|
88 | set tmp=`grep -c '^Phase name : Pg$' $evtfile` |
---|
89 | if ($tmp > 3) set mainphase=Pg |
---|
90 | set tmp=`grep -c '^Phase name : PP$' $evtfile` |
---|
91 | if ($tmp > 3) set mainphase=PP |
---|
92 | set tmp=`grep -c '^Phase name : P$' $evtfile` |
---|
93 | if ($tmp > 3) set mainphase=P |
---|
94 | # check for PKP |
---|
95 | if ("$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 |
---|
102 | endif |
---|
103 | |
---|
104 | echo mainphase $mainphase |
---|
105 | |
---|
106 | # create relative onset file |
---|
107 | if (-e $relfile) \rm $relfile |
---|
108 | $SH_UTIL/evt2reltimes $evtfile $relfile $mainphase |
---|
109 | |
---|
110 | # read preliminary location |
---|
111 | if ("$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 |
---|
121 | else |
---|
122 | set lat=$startlat |
---|
123 | endif |
---|
124 | if ("$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 |
---|
134 | else |
---|
135 | set lon=$startlon |
---|
136 | endif |
---|
137 | |
---|
138 | # make nice numbers |
---|
139 | set lat=`echo $lat | sed 's/\./ /'` |
---|
140 | set lat=$lat[1].5 |
---|
141 | set lon=`echo $lon | sed 's/\./ /'` |
---|
142 | set lon=$lon[1].5 |
---|
143 | |
---|
144 | #set echo |
---|
145 | |
---|
146 | # call travel time fit |
---|
147 | if (-e $fitxyz) \rm $fitxyz |
---|
148 | if ("$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 |
---|
152 | else |
---|
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 |
---|
162 | endif |
---|
163 | #\rm $relfile |
---|
164 | |
---|
165 | set infoline=`tail -1 $fitxyz` |
---|
166 | if ($#infoline < 11) then |
---|
167 | echo "infoline of $fitxyz has unknown syntax:" |
---|
168 | echo "$infoline" |
---|
169 | exit |
---|
170 | endif |
---|
171 | set minz=$infoline[3] |
---|
172 | set minlat=$infoline[5] |
---|
173 | set minlon=$infoline[7] |
---|
174 | set mindep=$infoline[9] |
---|
175 | set origtime=$infoline[11] |
---|
176 | if ("$method" == "grid") then |
---|
177 | set maxz=`grep -v '^\!' $fitxyz | grep " $mindep depth" | awk '{print $3}' | sort -n | tail -1` |
---|
178 | else |
---|
179 | set maxz=`grep -v '^\!' $fitxyz | awk '{print $3}' | sort -n | tail -1` |
---|
180 | endif |
---|
181 | |
---|
182 | # check for longitude larger than 180 |
---|
183 | set xtmp=`echo $minlon | sed 's/\./ /'` |
---|
184 | set xtmp=$xtmp[1] |
---|
185 | if ($xtmp >= 180) then |
---|
186 | set minlon=`echo "scale=3; $minlon - 360.0" | bc` |
---|
187 | endif |
---|
188 | |
---|
189 | # compute bounds of display |
---|
190 | set wlat2=`echo "scale=3; $wlat / 2" | bc` |
---|
191 | set wlon2=`echo "scale=3; $wlon / 2" | bc` |
---|
192 | set dlat2=`echo "scale=3; $dlat / 2" | bc` |
---|
193 | set dlon2=`echo "scale=3; $dlon / 2" | bc` |
---|
194 | set lolat=`echo "scale=3; $minlat - $wlat2" | bc` |
---|
195 | set hilat=`echo "scale=3; $minlat + $wlat2" | bc` |
---|
196 | set lolon=`echo "scale=3; $minlon - $wlon2" | bc` |
---|
197 | set hilon=`echo "scale=3; $minlon + $wlon2" | bc` |
---|
198 | |
---|
199 | set lolat=`echo "scale=3; $lolat - $dlat2" | bc` |
---|
200 | set hilat=`echo "scale=3; $hilat + $dlat2" | bc` |
---|
201 | set lolon=`echo "scale=3; $lolon - $dlon2" | bc` |
---|
202 | set hilon=`echo "scale=3; $hilon + $dlon2" | bc` |
---|
203 | |
---|
204 | # create text output |
---|
205 | if (-e $textresult) \rm $textresult |
---|
206 | touch $textresult |
---|
207 | echo "Best fit for" >>$textresult |
---|
208 | echo "$origtime" >>$textresult |
---|
209 | echo "lat=$minlat lon=$minlon" >>$textresult |
---|
210 | if ("$dep1" != "$dep2") then |
---|
211 | echo "Found depth $mindep km" >>$textresult |
---|
212 | else |
---|
213 | echo "Depth given as $mindep km" >>$textresult |
---|
214 | endif |
---|
215 | echo "Min-rms=$minz (orange)" >>$textresult |
---|
216 | echo "Max-rms=$maxz (blue)" >>$textresult |
---|
217 | echo "" >>$textresult |
---|
218 | echo "Found residuals:" >>$textresult |
---|
219 | echo "" >>$textresult |
---|
220 | grep '^\!' $fitxyz | grep -v 'minsq:' >>$textresult |
---|
221 | # $SH_TEXTEDIT $textresult & |
---|
222 | |
---|
223 | # create postscript output |
---|
224 | if (-e $psfile) \rm $psfile |
---|
225 | |
---|
226 | # create cpt file, make the best 3% turn orange |
---|
227 | set levz=`echo "scale=5; $maxz - $minz" | bc ` |
---|
228 | set levz=`echo "scale=5; $levz * 0.03" | bc ` |
---|
229 | set levz=`echo "scale=5; $levz + $minz" | bc ` |
---|
230 | if (-e $cpt) \rm $cpt |
---|
231 | touch $cpt |
---|
232 | echo "$minz 255 150 0 $levz 255 255 0" >>$cpt |
---|
233 | echo "$levz 255 255 0 $maxz 0 0 255" >>$cpt |
---|
234 | |
---|
235 | # compute symbol size |
---|
236 | set symsiz=`echo "scale=5; $stepsize / $width * $imagsize" | bc` |
---|
237 | |
---|
238 | set xlat=`echo $minlat | sed 's/-/s/'` |
---|
239 | set xlon=`echo $minlon | sed 's/-/w/'` |
---|
240 | set areaname="`$SH_UTIL/fereg $xlat $xlon`" |
---|
241 | |
---|
242 | touch $psfile |
---|
243 | |
---|
244 | if ("$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 |
---|
250 | else |
---|
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 |
---|
256 | endif |
---|
257 | |
---|
258 | # recompute longitude from negative values to 180..360 |
---|
259 | set xhilon=`echo $hilon | sed 's/-/m/'` |
---|
260 | set res=`$SH_UTIL/floatop $xhilon lt 0.0` |
---|
261 | if ($res == 1) set xhilon=`echo "scale=3; $hilon + 360.0" | bc` |
---|
262 | set xlolon=`echo $lolon | sed 's/-/m/'` |
---|
263 | set res=`$SH_UTIL/floatop $xlolon lt 0.0` |
---|
264 | if ($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 |
---|
271 | awk '{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 |
---|
274 | awk '{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 |
---|
278 | set lcnt=`wc -l $textresult` |
---|
279 | set lcnt=$lcnt[1] |
---|
280 | set ypos=0 |
---|
281 | set cont=-K |
---|
282 | while ($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 |
---|
290 | end |
---|
291 | |
---|
292 | touch $resfile |
---|
293 | |
---|
294 | echo "Latitude : $minlat" >>$resfile |
---|
295 | echo "Longitude : $minlon" >>$resfile |
---|
296 | echo "Depth (km) : $mindep" >>$resfile |
---|
297 | echo "Origin time : $origtime" >>$resfile |
---|
298 | echo "Location method : relative travel times" >>$resfile |
---|
299 | echo "--- End of Phase ---" >>$resfile |
---|
300 | |
---|
301 | if ($?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 | |
---|