1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28 package edu.sc.seis.TauP;
29
30 import java.io.File;
31 import java.io.FileNotFoundException;
32 import java.io.IOException;
33 import java.io.InputStreamReader;
34 import java.io.OptionalDataException;
35 import java.io.StreamCorruptedException;
36 import java.io.StreamTokenizer;
37 import java.io.Writer;
38
39 /***
40 * Calculates travel time curves
41 * at known slowness samples.
42 *
43 * @version 1.1.3 Wed Jul 18 15:00:35 GMT 2001
44
45
46
47 * @author H. Philip Crotwell
48 *
49 */
50 public class TauP_Curve extends TauP_Time {
51
52 /*** should the output file be a compete script? */
53 protected boolean gmtScript = false;
54
55 /*** should the output times use a reducing velocity? */
56 protected boolean reduceTime = false;
57
58 /*** the reducing velocity to use if reduceTime == true,
59 * in units of radians/second . */
60 protected double reduceVel = .125*Math.PI/180;
61
62 protected TauP_Curve() {
63 super();
64 }
65
66 public TauP_Curve(TauModel tMod) throws TauModelException {
67 super(tMod);
68 }
69
70 public TauP_Curve(String modelName) throws TauModelException {
71 super(modelName);
72 }
73
74 public boolean isGmtScript() {
75 return gmtScript;
76 }
77
78 public void setGmtScript(boolean gmtScript) {
79 this.gmtScript = gmtScript;
80 }
81
82 public boolean isReduceTime() {
83 return reduceTime;
84 }
85
86 public void setReduceTime(boolean reduceTime) {
87 this.reduceTime = reduceTime;
88 }
89
90 /*** @returns reducing velocity in degrees/second. The internal
91 * usage is radians/second. */
92 public double getReduceVelDeg() {
93 return 180.0/Math.PI*reduceVel;
94 }
95
96 /*** set the reducing velocity, in degrees/second. The internal
97 * representation is radians/second. */
98 public void setReduceVelDeg(double reduceVel) {
99 if (reduceVel > 0.0) {
100 this.reduceVel = Math.PI/180.0*reduceVel;
101 }
102 }
103
104 /*** @returns reducing velocity in kilometers/second. The internal
105 * usage is radians/second. */
106 public double getReduceVelKm() {
107 return reduceVel*tMod.getRadiusOfEarth();
108 }
109
110 /*** set the reducing velocity, in kilometers/second. The internal
111 * representation is radians/second. */
112 public void setReduceVelKm(double reduceVel) {
113 if (reduceVel > 0.0) {
114 if (tMod != null) {
115 this.reduceVel = reduceVel/tMod.getRadiusOfEarth();
116 } else {
117 this.reduceVel = reduceVel/ 6371.0;
118 }
119 }
120 }
121
122 public void calculate(double degrees) {
123
124
125
126
127 recalcPhases();
128 }
129
130 public void init() throws IOException {
131
132 super.init();
133
134 if (gmtScript) {
135 String psFile;
136 if (outFile.endsWith(".gmt")) {
137 psFile = outFile.substring(0,outFile.length()-4)+".ps";
138 } else {
139 psFile = outFile+".ps";
140 }
141
142 dos.writeBytes("#!/bin/sh\n");
143 dos.writeBytes(
144 "#\n# This script will plot curves using GMT. If you want to\n"+
145 "#use this as a data file for psxy in another script, delete these"+
146 "\n# first lines, as well as the last line.\n#\n");
147 dos.writeBytes("/bin/rm -f "+psFile+"\n\n");
148 }
149 }
150
151 public void printStdUsage() {
152 String className = this.getClass().getName();
153 className =
154 className.substring(className.lastIndexOf('.')+1,className.length());
155
156 System.out.println("Usage: "+className.toLowerCase()+" [arguments]");
157 System.out.println(" or, for purists, java "+this.getClass().getName()+
158 " [arguments]");
159 System.out.println("\nArguments are:");
160
161 System.out.println(
162 "-ph phase list -- comma separated phase list\n"+
163 "-pf phasefile -- file containing phases\n\n"+
164 "-mod[el] modelname -- use velocity model \"modelname\" for calculations\n"+
165 " Default is iasp91.\n\n"+
166 "-h depth -- source depth in km\n\n");
167 }
168
169 public void printStdUsageTail() {
170 System.out.println(
171 "\n-o outfile -- output is redirected to \"outfile\" instead of taup_curve.gmt\n"+
172 "-debug -- enable debugging output\n"+
173 "-verbose -- enable verbose output\n"+
174 "-version -- print the version\n"+
175 "-help -- print this out, but you already know that!\n\n");
176 }
177
178 public void printUsage() {
179 printStdUsage();
180 System.out.println("-gmt -- outputs curves as a complete GMT script.");
181 System.out.println("-reddeg velocity -- outputs curves with a reducing velocity (deg/sec).");
182 System.out.println("-redkm velocity -- outputs curves with a reducing velocity (km/sec).");
183 printStdUsageTail();
184 }
185
186
187 public void start() throws IOException, TauModelException {
188 double tempDepth;
189
190 if (depth!=-1*Double.MAX_VALUE) {
191
192 depthCorrect( Double.valueOf(
193 toolProps.getProperty("taup.source.depth", "0.0")).doubleValue());
194 calculate(degrees);
195 printResult(dos);
196 } else {
197 StreamTokenizer tokenIn = new StreamTokenizer(
198 new InputStreamReader(System.in));
199 tokenIn.parseNumbers();
200 tokenIn.wordChars(',' , ',');
201 tokenIn.wordChars('_' , '_');
202 System.out.print("Enter Depth: ");
203 tokenIn.nextToken();
204 tempDepth = tokenIn.nval;
205 if (tempDepth < 0.0 || depth > tMod.getRadiusOfEarth()) {
206 System.out.println("Depth must be >= 0.0 and "+
207 "<= tMod.getRadiusOfEarth().\ndepth = "+tempDepth);
208 return;
209 }
210
211 depthCorrect(tempDepth);
212 calculate(degrees);
213 printResult(dos);
214 }
215 }
216
217 public void destroy() throws IOException {
218 if (gmtScript) {
219 dos.writeBytes("END\n");
220 }
221 super.destroy();
222 }
223
224 public void printResult(Writer out) throws IOException {
225
226 SeismicPhase phase;
227 Arrival currArrival;
228 double[] dist, time, rayParams;
229 double arcDistance, timeReduced;
230 Arrival[] phaseArrivals;
231 double maxTime = -1*Double.MAX_VALUE, minTime = Double.MAX_VALUE;
232
233 if (gmtScript) {
234 String scriptStuff = "";
235 String psFile;
236 if (outFile.endsWith(".gmt")) {
237 psFile = outFile.substring(0,outFile.length()-4)+".ps";
238 } else {
239 psFile = outFile+".ps";
240 }
241
242 for (int phaseNum=0;phaseNum<phases.size(); phaseNum++) {
243 phase = (SeismicPhase)phases.elementAt(phaseNum);
244 dist = phase.getDist();
245 time = phase.getTime();
246 if (dist.length > 0) {
247 arcDistance = Math.acos(Math.cos(dist[dist.length-1]));
248 if (reduceTime) {
249 scriptStuff +=
250 (float)(180.0/Math.PI*arcDistance)
251 +" "+ (float)(time[dist.length-1]-arcDistance/reduceVel)+
252 " 10 0 0 9 "+phase.getName()+"\n";
253
254
255 for (int i=0;i<time.length;i++) {
256 timeReduced = time[i]-Math.acos(Math.cos(dist[i]))/reduceVel;
257 if (timeReduced > maxTime) {maxTime = timeReduced;}
258 if (timeReduced < minTime) {minTime = timeReduced;}
259 }
260 } else {
261 scriptStuff +=
262 (float)(180.0/Math.PI*arcDistance)
263 +" "+ (float)time[dist.length-1]+
264 " 10 0 0 1 "+phase.getName()+"\n";
265
266
267 for (int i=0;i<time.length;i++) {
268 if (time[i] > maxTime) {maxTime = time[i];}
269 if (time[i] < minTime) {minTime = time[i];}
270 }
271 }
272 }
273 }
274
275 maxTime = Math.ceil(maxTime/100)*100;
276 minTime = Math.floor(minTime/100)*100;
277 out.write("pstext -JX6 -P -R0/180/"+minTime+"/"+maxTime+
278 " -B20/100/:.'"+modelName+"': -K > "+psFile+" <<END\n");
279 out.write(scriptStuff);
280 out.write("END\n\n");
281
282 out.write(
283 "psxy -JX -R -M -O >> "+psFile+" <<END\n");
284 }
285
286 for (int phaseNum=0;phaseNum<phases.size(); phaseNum++) {
287 phase = (SeismicPhase)phases.elementAt(phaseNum);
288 dist = phase.getDist();
289 time = phase.getTime();
290 rayParams = phase.getRayParams();
291 if (dist.length > 0) {
292 out.write("> "+phase.getName()+" for a source depth of "+depth+
293 " kilometers in the "+modelName+" model\n");
294 }
295 for (int i=0; i<dist.length; i++) {
296
297
298 arcDistance = Math.acos(Math.cos(dist[i]));
299 if (reduceTime) {
300 timeReduced = time[i]-arcDistance/reduceVel;
301 } else {
302 timeReduced = time[i];
303 }
304 out.write((float)(180.0/Math.PI*arcDistance)
305 +" "+ (float)timeReduced+"\n");
306
307 if (i<dist.length-1 && (rayParams[i] == rayParams[i+1]) &&
308 rayParams.length > 2) {
309
310 out.write("> Shadow Zone\n");
311 continue;
312 }
313
314
315
316 if (i<dist.length-1 && Math.sin(dist[i]) > 0 &&
317 Math.sin(dist[i+1]) <0) {
318 phase.calcTime(180.0);
319 phaseArrivals = phase.getArrivals();
320 int j=0;
321 while (j<phaseArrivals.length) {
322 if ((phase.rayParams[i]-phaseArrivals[j].rayParam)*
323 (phaseArrivals[j].rayParam-phase.rayParams[i+1])>0) {
324 if (reduceTime) {
325 out.write("180.0 "+
326 (float)(phaseArrivals[j].time- Math.PI/reduceVel)+"\n");
327 } else {
328 out.write("180.0 "+(float)phaseArrivals[j].time+"\n");
329 }
330 break;
331 }
332 j++;
333 }
334 }
335
336
337
338 if (i<dist.length-1 && Math.sin(dist[i]) < 0 &&
339 Math.sin(dist[i+1]) > 0) {
340 phase.calcTime(0.0);
341 phaseArrivals = phase.getArrivals();
342 int j=0;
343 while (j<phaseArrivals.length) {
344 if ((phase.rayParams[i]-phaseArrivals[j].rayParam)*
345 (phaseArrivals[j].rayParam-phase.rayParams[i+1])>0) {
346 out.write("0.0 "+(float)phaseArrivals[j].time+"\n");
347 break;
348 }
349 j++;
350 }
351 }
352 }
353
354 }
355 }
356
357 public String[] parseCmdLineArgs(String[] args) throws IOException {
358 int i=0;
359 String[] leftOverArgs;
360 int numNoComprendoArgs = 0;
361 File tempFile;
362
363 leftOverArgs = super.parseCmdLineArgs(args);
364 String[] noComprendoArgs = new String[leftOverArgs.length];
365
366 while (i<leftOverArgs.length) {
367 if (leftOverArgs[i].equalsIgnoreCase("-gmt")) {
368 gmtScript = true;
369 } else if (leftOverArgs[i].equals("-reddeg")) {
370 setReduceTime(true);
371 setReduceVelDeg(Double.valueOf(leftOverArgs[i+1]).doubleValue());
372 i++;
373 } else if (leftOverArgs[i].equals("-redkm")) {
374 setReduceTime(true);
375 setReduceVelKm(Double.valueOf(leftOverArgs[i+1]).doubleValue());
376 i++;
377 } else if (leftOverArgs[i].equals("-help")) {
378 noComprendoArgs[numNoComprendoArgs++] = leftOverArgs[i];
379 } else {
380 noComprendoArgs[numNoComprendoArgs++] = leftOverArgs[i];
381 }
382 i++;
383 }
384
385 if (numNoComprendoArgs > 0) {
386 String[] temp = new String[numNoComprendoArgs];
387 System.arraycopy(noComprendoArgs,0,temp,0,numNoComprendoArgs);
388 return temp;
389 } else {
390 return new String[0];
391 }
392 }
393
394 /*** Allows TauP_Curve to run as an application. Creates an instance
395 * of TauP_Curve.
396 * . */
397 public static void main(String[] args)
398 throws FileNotFoundException,
399 IOException,
400 StreamCorruptedException,
401 ClassNotFoundException,
402 OptionalDataException
403 {
404 boolean doInteractive = true;
405 try {
406 TauP_Curve tauPCurve = new TauP_Curve();
407 tauPCurve.outFile = "taup_curve.gmt";
408
409 String[] noComprendoArgs = tauPCurve.parseCmdLineArgs(args);
410 if (noComprendoArgs.length > 0) {
411 for (int i=0;i<noComprendoArgs.length;i++) {
412 if (noComprendoArgs[i].equals("-help") ||
413 noComprendoArgs[i].equals("-version")) {
414 System.exit(0);
415 }
416 }
417 System.out.println("I don't understand the following arguments, continuing:");
418 for (int i=0;i<noComprendoArgs.length;i++) {
419 System.out.print(noComprendoArgs[i]+" ");
420 if (noComprendoArgs[i].equals("-help")) {
421 System.out.println();
422 System.exit(0);
423 }
424 }
425 System.out.println();
426 noComprendoArgs = null;
427 }
428
429 for (int i=0; i<args.length; i++) {
430 if (args[i] == "-h") {
431 doInteractive = false;
432 }
433 }
434
435 if (tauPCurve.DEBUG) {
436 System.out.println("Done reading "+tauPCurve.modelName);
437 }
438
439 tauPCurve.init();
440 if (doInteractive) {
441 tauPCurve.start();
442 } else {
443
444 tauPCurve.depthCorrect( Double.valueOf(
445 tauPCurve.toolProps.getProperty("taup.source.depth", "0.0")
446 ).doubleValue());
447 tauPCurve.calculate(tauPCurve.degrees);
448 tauPCurve.printResult(tauPCurve.dos);
449 }
450
451 tauPCurve.destroy();
452 } catch (TauModelException e) {
453 System.out.println("Caught TauModelException: "+e.getMessage());
454 e.printStackTrace();
455 }
456
457 }
458 }