1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17 package edu.sc.seis.TauP;
18
19 import java.io.File;
20 import java.io.FileNotFoundException;
21 import java.io.IOException;
22 import java.io.OptionalDataException;
23 import java.io.StreamCorruptedException;
24 import java.io.Writer;
25
26 /***
27 * Calculate travel paths for different phases using a linear interpolated ray
28 * parameter between known slowness samples.
29 *
30 * @version 1.1.3 Wed Jul 18 15:00:35 GMT 2001
31 * @author H. Philip Crotwell
32 */
33 public class TauP_Path extends TauP_Pierce {
34
35 protected float mapWidth = (float)6.0;
36
37 protected boolean gmtScript = false;
38
39 protected double maxPathInc = 1.0;
40
41 protected static Format float8_4 = new Format("%8.4f");
42
43 protected TauP_Path() {
44 super();
45 outFile = null;
46 }
47
48 public TauP_Path(TauModel tMod) throws TauModelException {
49 super(tMod);
50 outFile = null;
51 }
52
53 public TauP_Path(String modelName) throws TauModelException {
54 super(modelName);
55 outFile = null;
56 }
57
58 public TauP_Path(TauModel tMod, String outFileBase)
59 throws TauModelException {
60 super(tMod);
61 setOutFileBase(outFileBase);
62 }
63
64 public TauP_Path(String modelName, String outFileBase)
65 throws TauModelException {
66 super(modelName);
67 setOutFileBase(outFileBase);
68 }
69
70 /*** Sets the output file base, appending ".gmt" for the filename. */
71 public void setOutFileBase(String outFileBase) {
72 if(outFileBase != null && outFileBase.length() != 0) {
73 outFile = outFileBase + ".gmt";
74 } else {
75 outFile = "taup_path.gmt";
76 }
77 }
78
79 /***
80 * Sets the gmt map width to be used with the output script and for creating
81 * the circles for each discontinuity. Default is 6 inches.
82 */
83 public void setMapWidth() {
84 }
85
86 /***
87 * Gets the gmt map width to be used with the output script and for creating
88 * the circles for each discontinuity.
89 */
90 public float getMapWidth() {
91 return mapWidth;
92 }
93
94 public boolean getGmtScript() {
95 return gmtScript;
96 }
97
98 public void setGmtScript(boolean gmtScript) {
99 this.gmtScript = gmtScript;
100 }
101
102 public double getMaxPathInc() {
103 return maxPathInc;
104 }
105
106 public void setMaxPathInc(double maxPathInc) {
107 this.maxPathInc = maxPathInc;
108 }
109
110 public void calculate(double degrees) throws TauModelException {
111 depthCorrect(getSourceDepth());
112 recalcPhases();
113 clearArrivals();
114 calcPath(degrees);
115 }
116
117 public void calcPath(double degrees) {
118 this.degrees = degrees;
119 SeismicPhase phase;
120 Arrival[] phaseArrivals;
121 double calcTime, calcDist;
122 for(int phaseNum = 0; phaseNum < phases.size(); phaseNum++) {
123 phase = (SeismicPhase)phases.elementAt(phaseNum);
124 phase.setDEBUG(DEBUG);
125 phase.calcTime(degrees);
126 if(phase.hasArrivals()) {
127 phase.calcPath(tModDepth);
128 phaseArrivals = phase.getArrivals();
129 for(int i = 0; i < phaseArrivals.length; i++) {
130 arrivals.addElement(phaseArrivals[i]);
131 }
132 }
133 }
134 }
135
136 public void printResult(Writer out) throws IOException {
137 double calcTime, calcDist, prevDepth, calcDepth;
138 double lat, lon;
139 Arrival currArrival;
140 int interpNum, maxInterpNum;
141 double interpInc;
142 double radiusOfEarth = tModDepth.getRadiusOfEarth();
143 boolean longWayRound;
144 Format float12_5 = new Format("%12.5f");
145 for(int i = 0; i < arrivals.size(); i++) {
146 currArrival = (Arrival)arrivals.elementAt(i);
147 out.write("> " + currArrival.name + " at "
148 + outForms.formatDistance(currArrival.getDistDeg())
149 + " degrees for a "
150 + outForms.formatDepth(currArrival.sourceDepth)
151 + " km deep source in the " + modelName + " model arriving at "
152 + outForms.formatTime(currArrival.getTime())
153 + " s with rayParam "+outForms.formatRayParam(Math.PI/180*currArrival.getRayParam())+" s/deg.\n");
154 longWayRound = false;
155 if((currArrival.dist * 180 / Math.PI) % 360 > 180) {
156 longWayRound = true;
157 }
158 calcTime = 0.0;
159 calcDist = 0.0;
160 calcDepth = currArrival.sourceDepth;
161 for(int j = 0; j < currArrival.path.length; j++) {
162 calcTime = currArrival.path[j].time;
163 calcDepth = currArrival.path[j].depth;
164 prevDepth = calcDepth;
165 calcDist = currArrival.path[j].dist * 180.0 / Math.PI;
166 if(longWayRound && calcDist != 0.0) {
167 calcDist = -1.0 * calcDist;
168 }
169 out.write(outForms.formatDistance(calcDist) + " "
170 + outForms.formatDepth(radiusOfEarth - calcDepth));
171 if(!gmtScript) {
172 printLatLon(out, calcDist);
173 }
174 out.write("\n");
175 if(j < currArrival.path.length - 1
176 && (currArrival.rayParam != 0.0 && 180.0
177 / Math.PI
178 * (currArrival.path[j + 1].dist - currArrival.path[j].dist) > maxPathInc)) {
179
180
181 maxInterpNum = (int)Math.ceil((currArrival.path[j + 1].dist - currArrival.path[j].dist)
182 * 180.0 / Math.PI / maxPathInc);
183 for(interpNum = 1; interpNum < maxInterpNum; interpNum++) {
184 calcTime += (currArrival.path[j + 1].time - currArrival.path[j].time)
185 / maxInterpNum;
186 if(longWayRound) {
187 calcDist -= (currArrival.path[j + 1].dist - currArrival.path[j].dist)
188 / maxInterpNum * 180.0 / Math.PI;
189 } else {
190 calcDist += (currArrival.path[j + 1].dist - currArrival.path[j].dist)
191 / maxInterpNum * 180.0 / Math.PI;
192 }
193 calcDepth = prevDepth + interpNum
194 * (currArrival.path[j + 1].depth - prevDepth)
195 / maxInterpNum;
196 out.write(outForms.formatDistance(calcDist)
197 + " "
198 + outForms.formatDepth(radiusOfEarth
199 - calcDepth));
200 if(!gmtScript) {
201 printLatLon(out, calcDist);
202 }
203 out.write("\n");
204 }
205 }
206 prevDepth = currArrival.path[j].depth;
207 }
208 }
209 }
210
211 protected void printLatLon(Writer out, double calcDist) throws IOException {
212 double lat, lon;
213 if(eventLat != Double.MAX_VALUE && eventLon != Double.MAX_VALUE
214 && azimuth != Double.MAX_VALUE) {
215 lat = SphericalCoords.latFor(eventLat, eventLon, calcDist, azimuth);
216 lon = SphericalCoords.lonFor(eventLat, eventLon, calcDist, azimuth);
217 out.write(" " + outForms.formatLatLon(lat) + " "
218 + outForms.formatLatLon(lon));
219 } else if(stationLat != Double.MAX_VALUE
220 && stationLon != Double.MAX_VALUE
221 && backAzimuth != Double.MAX_VALUE) {
222 lat = SphericalCoords.latFor(stationLat, stationLon, degrees
223 - calcDist, backAzimuth);
224 lon = SphericalCoords.lonFor(stationLat, stationLon, degrees
225 - calcDist, backAzimuth);
226 out.write(" " + outForms.formatLatLon(lat) + " "
227 + outForms.formatLatLon(lon));
228 } else if(stationLat != Double.MAX_VALUE
229 && stationLon != Double.MAX_VALUE
230 && eventLat != Double.MAX_VALUE && eventLon != Double.MAX_VALUE) {
231 azimuth = SphericalCoords.azimuth(eventLat,
232 eventLon,
233 stationLat,
234 stationLon);
235 backAzimuth = SphericalCoords.azimuth(stationLat,
236 stationLon,
237 eventLat,
238 eventLon);
239 lat = SphericalCoords.latFor(eventLat, eventLon, calcDist, azimuth);
240 lon = SphericalCoords.lonFor(eventLat, eventLon, calcDist, azimuth);
241 out.write(" " + outForms.formatLatLon(lat) + " "
242 + outForms.formatLatLon(lon));
243 }
244 }
245
246 public void init() throws IOException {
247 super.init();
248 if(gmtScript) {
249 String psFile;
250 if(outFile == null) {
251 outFile = "taup_path.gmt";
252 psFile = "taup_path.ps";
253 } else if(outFile.endsWith(".gmt")) {
254 psFile = outFile.substring(0, outFile.length() - 4) + ".ps";
255 } else {
256 psFile = outFile + ".ps";
257 }
258 dos.writeBytes("#!/bin/sh\n");
259 dos.writeBytes("#\n# This script will plot ray paths using GMT. If you want to\n"
260 + "#use this as a data file for psxy in another script, delete these"
261 + "\n# first lines, to the last psxy, as well as the last line.\n#\n");
262 dos.writeBytes("/bin/rm -f " + psFile + "\n\n");
263 dos.writeBytes("# draw surface and label distances.\n"
264 + "psbasemap -K -P -R0/360/0/6371 -JP" + mapWidth
265 + " -B30p/500N > " + psFile + "\n\n");
266 dos.writeBytes("# draw circles for branches, note these are scaled for a \n"
267 + "# map using -JP"
268 + mapWidth
269 + "\n"
270 + "psxy -K -O -P -R -JP -Sc -A >> "
271 + psFile
272 + " <<ENDLAYERS\n");
273 for(int j = 0; j < 2; j++) {
274 dos.writeBytes("0.0 0.0 "
275 + (float)(tMod.getRadiusOfEarth() * mapWidth / tMod.getRadiusOfEarth())
276 + "\n");
277 for(int i = 0; i < tMod.tauBranches[j].length; i++) {
278 dos.writeBytes("0.0 0.0 "
279 + (float)((tMod.getRadiusOfEarth() - tMod.tauBranches[j][i].getBotDepth())
280 * mapWidth / tMod.getRadiusOfEarth())
281 + "\n");
282 }
283 }
284 dos.writeBytes("ENDLAYERS\n\n");
285 dos.writeBytes("# draw paths\n");
286 dos.writeBytes("psxy -P -R -O -JP -M -A >> " + psFile + " <<END\n");
287 }
288 }
289
290 public void printUsage() {
291 printStdUsage();
292 System.out.println("-gmt -- outputs path as a complete GMT script.");
293 printStdUsageTail();
294 }
295
296 public String[] parseCmdLineArgs(String[] args) throws IOException {
297 int i = 0;
298 String[] leftOverArgs;
299 int numNoComprendoArgs = 0;
300 File tempFile;
301 leftOverArgs = super.parseCmdLineArgs(args);
302 String[] noComprendoArgs = new String[leftOverArgs.length];
303 while(i < leftOverArgs.length) {
304 if(leftOverArgs[i].equalsIgnoreCase("-gmt")) {
305 gmtScript = true;
306 } else if(leftOverArgs[i].equals("-help")) {
307 noComprendoArgs[numNoComprendoArgs++] = leftOverArgs[i];
308 } else {
309 noComprendoArgs[numNoComprendoArgs++] = leftOverArgs[i];
310 }
311 i++;
312 }
313 if(numNoComprendoArgs > 0) {
314 String[] temp = new String[numNoComprendoArgs];
315 System.arraycopy(noComprendoArgs, 0, temp, 0, numNoComprendoArgs);
316 return temp;
317 } else {
318 return new String[0];
319 }
320 }
321
322 public void start() throws IOException, TauModelException, TauPException {
323 super.start();
324 }
325
326 public void destroy() throws IOException {
327 if(gmtScript) {
328 dos.writeBytes("END\n");
329 }
330 super.destroy();
331 }
332
333 /***
334 * Allows TauP_Path to run as an application. Creates an instance of
335 * TauP_Path and calls TauP_Path.init() and TauP_Path.start().
336 */
337 public static void main(String[] args) throws FileNotFoundException,
338 IOException, StreamCorruptedException, ClassNotFoundException,
339 OptionalDataException {
340 try {
341 TauP_Path tauPPath = new TauP_Path();
342 tauPPath.setOutFileBase("taup_path");
343 String[] noComprendoArgs = tauPPath.parseCmdLineArgs(args);
344 if(noComprendoArgs.length > 0) {
345 for(int i = 0; i < noComprendoArgs.length; i++) {
346 if(noComprendoArgs[i].equals("-help")
347 || noComprendoArgs[i].equals("-version")) {
348 System.exit(0);
349 }
350 }
351 System.out.println("I don't understand the following arguments, continuing:");
352 for(int i = 0; i < noComprendoArgs.length; i++) {
353 System.out.print(noComprendoArgs[i] + " ");
354 if(noComprendoArgs[i].equals("-help")
355 || noComprendoArgs[i].equals("-version")) {
356 System.out.println();
357 System.exit(0);
358 }
359 }
360 System.out.println();
361 noComprendoArgs = null;
362 }
363 tauPPath.init();
364 if(tauPPath.DEBUG) {
365 System.out.println("Done reading " + tauPPath.modelName);
366 }
367 tauPPath.start();
368 tauPPath.destroy();
369 } catch(TauModelException e) {
370 System.out.println("Caught TauModelException: " + e.getMessage());
371 e.printStackTrace();
372 } catch(TauPException e) {
373 System.out.println("Caught TauPException: " + e.getMessage());
374 e.printStackTrace();
375 }
376 }
377 }