View Javadoc

1   /*
2    * The TauP Toolkit: Flexible Seismic Travel-Time and Raypath Utilities.
3    * Copyright (C) 1998-2000 University of South Carolina This program is free
4    * software; you can redistribute it and/or modify it under the terms of the GNU
5    * General Public License as published by the Free Software Foundation; either
6    * version 2 of the License, or (at your option) any later version. This program
7    * is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY;
8    * without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
9    * PARTICULAR PURPOSE. See the GNU General Public License for more details. You
10   * should have received a copy of the GNU General Public License along with this
11   * program; if not, write to the Free Software Foundation, Inc., 59 Temple Place -
12   * Suite 330, Boston, MA 02111-1307, USA. The current version can be found at <A
13   * HREF="www.seis.sc.edu">http://www.seis.sc.edu </A> Bug reports and comments
14   * should be directed to H. Philip Crotwell, crotwell@seis.sc.edu or Tom Owens,
15   * owens@seis.sc.edu
16   */
17  package edu.sc.seis.TauP;
18  
19  import java.io.FileNotFoundException;
20  import java.io.IOException;
21  import java.io.OptionalDataException;
22  import java.io.Serializable;
23  import java.io.StreamCorruptedException;
24  import java.util.Vector;
25  
26  /***
27   * Stores and transforms seismic phase names to and from their corresponding
28   * sequence of branches.
29   * 
30   * @version 1.1.3 Wed Jul 18 15:00:35 GMT 2001
31   * @author H. Philip Crotwell
32   */
33  public class SeismicPhase implements Serializable, Cloneable {
34  
35      /*** Enables debugging output. */
36      public transient boolean DEBUG = false;
37  
38      /*** Enables verbose output. */
39      public transient boolean verbose = false;
40  
41      /*** TauModel to generate phase for. */
42      protected TauModel tMod;
43  
44      /***
45       * Used by addToBranch when the path turns within a segment. We assume that
46       * no ray will turn downward so turning implies turning from downward to
47       * upward, ie U.
48       */
49      protected static final int TURN = 0;
50  
51      /***
52       * Used by addToBranch when the path reflects off the top of the end of a
53       * segment, ie ^.
54       */
55      protected static final int REFLECTTOP = 1;
56  
57      /***
58       * Used by addToBranch when the path reflects off the bottom of the end of a
59       * segment, ie v.
60       */
61      protected static final int REFLECTBOT = 2;
62  
63      /***
64       * Used by addToBranch when the path transmits up through the end of a
65       * segment.
66       */
67      protected static final int TRANSUP = 3;
68  
69      /***
70       * Used by addToBranch when the path transmits down through the end of a
71       * segment.
72       */
73      protected static final int TRANSDOWN = 4;
74  
75      /***
76       * The maximum degrees that a Pn or Sn can refract along the moho. Note this
77       * is not the total distance, only the segment along the moho. The default
78       * is 20 degrees.
79       */
80      protected double maxRefraction = 20;
81  
82      /***
83       * The maximum degrees that a Pdiff or Sdiff can diffract along the CMB.
84       * Note this is not the total distance, only the segment along the CMB. The
85       * default is 60 degrees.
86       */
87      protected double maxDiffraction = 60;
88  
89      /***
90       * The source depth within the TauModel that was used to generate this
91       * phase.
92       */
93      protected double sourceDepth;
94  
95      /***
96       * Array of distances corresponding to the ray parameters stored in
97       * rayParams.
98       */
99      protected double[] dist = new double[0];
100 
101     /***
102      * Array of times corresponding to the ray parameters stored in rayParams.
103      */
104     protected double[] time = new double[0];
105 
106     /*** Array of possible ray parameters for this phase. */
107     protected double[] rayParams = new double[0];
108 
109     /*** Minimum ray parameter that exists for this phase. */
110     protected double minRayParam;
111 
112     /*** Maximum ray parameter that exists for this phase. */
113     protected double maxRayParam;
114 
115     /***
116      * Index within TauModel.rayParams that corresponds to maxRayParam. Note
117      * that maxRayParamIndex < minRayParamIndex as ray parameter decreases with
118      * increasing index.
119      */
120     protected int maxRayParamIndex = -1;
121 
122     /***
123      * Index within TauModel.rayParams that corresponds to minRayParam. Note
124      * that maxRayParamIndex < minRayParamIndex as ray parameter decreases with
125      * increasing index.
126      */
127     protected int minRayParamIndex = -1;
128 
129     /*** The minimum distance that this phase can be theoretically observed. */
130     protected double minDistance = 0.0;
131 
132     /*** The maximum distance that this phase can be theoretically observed. */
133     protected double maxDistance = Double.MAX_VALUE;
134 
135     /***
136      * Vector holding all arrivals for this phase at the given distance. They
137      * are stored as Arrival objects.
138      * 
139      * @see Arrival
140      * @see calcTime(double)
141      * @see calcPierce(TauModel)
142      * @see calcPath(TauModel)
143      */
144     protected Vector arrivals = new Vector();
145 
146     /***
147      * Array of branch numbers for the given phase. Note that this depends upon
148      * both the earth model and the source depth.
149      */
150     protected Vector branchSeq = new Vector();
151 
152     /*** The phase name, ie PKiKP. */
153     protected String name;
154 
155     /***
156      * name with depths corrected to be actuall discontinuities in the model.
157      */
158     protected String puristName;
159 
160     /*** Vector containing Strings for each leg. */
161     protected Vector legs = new Vector();
162 
163     /***
164      * temporary branch number so we know where to start add to the branch
165      * sequence. Used in addToBranch() and parseName().
166      */
167     protected transient int currBranch;
168 
169     /***
170      * temporary end action so we know what we did at the end of the last
171      * section of the branch sequence. Used in addToBranch() and parseName().
172      */
173     protected transient int endAction;
174 
175     /***
176      * records the end action for the current leg. Will be one of
177      * SeismicPhase.TURN, SeismicPhase.TRANSDOWN, SeismicPhase.TRANSUP,
178      * SeismicPhase.REFLECTBOT, or SeismicPhase.REFLECTTOP. This allows a check
179      * to make sure the path is correct. Used in addToBranch() and parseName().
180      */
181     protected Vector legAction = new Vector();
182 
183     /***
184      * true if the current leg of the phase is down going. This allows a check
185      * to make sure the path is correct. Used in addToBranch() and parseName().
186      */
187     protected Vector downGoing = new Vector();
188 
189     /***
190      * Vector of wave types corresponding to each leg of the phase.
191      * 
192      * @see legs
193      */
194     protected Vector waveType = new Vector();
195 
196     public static final boolean PWAVE = true;
197 
198     public static final boolean SWAVE = false;
199 
200     // Methods ----------------------------------------------------------
201     // Constructors
202     /***
203      * @param phaseName
204      *            String containing a name of the phase.
205      * @param tMod
206      *            Tau model to be used to construct the phase.
207      */
208     public SeismicPhase(String name, TauModel tMod) {
209         this.name = name;
210         this.sourceDepth = tMod.sourceDepth;
211         this.tMod = tMod;
212     }
213 
214     // Accessor methods
215     public boolean hasArrivals() {
216         if(arrivals.size() > 0) {
217             return true;
218         } else {
219             return false;
220         }
221     }
222 
223     /*** Returns arrival time array. */
224     public Arrival[] getArrivals() {
225         Arrival[] returnArrivals = new Arrival[arrivals.size()];
226         for(int i = 0; i < arrivals.size(); i++) {
227             returnArrivals[i] = (Arrival)((Arrival)arrivals.elementAt(i));
228         }
229         return returnArrivals;
230     }
231 
232     public void setTauModel(TauModel tMod) throws TauModelException {
233         this.tMod = tMod;
234         init();
235     }
236 
237     public TauModel getTauModel() {
238         return tMod;
239     }
240 
241     public void setDEBUG(boolean DEBUG) {
242         this.DEBUG = DEBUG;
243     }
244 
245     public double getMinDistance() {
246         return minDistance;
247     }
248 
249     public double getMaxDistance() {
250         return maxDistance;
251     }
252 
253     public double getMaxRayParam() {
254         return maxRayParam;
255     }
256 
257     public double getMinRayParam() {
258         return minRayParam;
259     }
260 
261     public int getMaxRayParamIndex() {
262         return maxRayParamIndex;
263     }
264 
265     public int getMinRayParamIndex() {
266         return minRayParamIndex;
267     }
268 
269     public double getMaxDiffraction() {
270         return maxDiffraction;
271     }
272 
273     public void setMaxDiffraction(double maxDiffraction) {
274         this.maxDiffraction = maxDiffraction;
275     }
276 
277     public String getName() {
278         return name;
279     }
280 
281     public String getPuristName() {
282         return name;
283     }
284 
285     public String[] getLegs() {
286         String[] legsArray = new String[legs.size()];
287         for(int i = 0; i < legs.size(); i++) {
288             legsArray[i] = (String)legs.elementAt(i);
289         }
290         return legsArray;
291     }
292 
293     public double[] getRayParams() {
294         return (double[])rayParams.clone();
295     }
296 
297     public double[] getDist() {
298         return (double[])dist.clone();
299     }
300 
301     public double[] getTime() {
302         return (double[])time.clone();
303     }
304 
305     public double[] getTau() {
306         double[] tau = new double[dist.length];
307         for(int i = 0; i < dist.length; i++) {
308             tau[i] = time[i] - rayParams[i] * dist[i];
309         }
310         return tau;
311     }
312 
313     // Normal methods
314     public void init() throws TauModelException {
315         legPuller();
316         createPuristName(tMod);
317         parseName(tMod);
318         sumBranches(tMod);
319     }
320 
321     /*** calculates arrival times for this phase. */
322     public void calcTime(double deg) {
323         double tempDeg = deg;
324         if(tempDeg < 0.0) {
325             tempDeg *= -1.0;
326         } // make sure deg is positive
327         while(tempDeg > 360.0) {
328             tempDeg -= 360.0;
329         } // make sure it is less than 360
330         if(tempDeg > 180.0) {
331             tempDeg = 360.0 - tempDeg;
332         } // make sure less than or equal to 180
333         // now we have 0.0 <= deg <= 180
334         double radDist = tempDeg * Math.PI / 180.0;
335         arrivals.removeAllElements();
336         /*
337          * Search all distances 2n*PI+radDist and 2(n+1)*PI-radDist that are
338          * less than the maximum distance for this phase. This insures that we
339          * get the time for phases that accumulate more than 180 degrees of
340          * distance, for instance PKKKKP might wrap all of the way around. A
341          * special case exists at 180, so we skip the second case if
342          * tempDeg==180.
343          */
344         int n = 0;
345         Arrival newArrival;
346         double searchDist;
347         while(n * 2.0 * Math.PI + radDist <= maxDistance) {
348             /*
349              * Look for arrivals that are radDist + 2nPi, ie rays that have done
350              * more than n laps.
351              */
352             searchDist = n * 2.0 * Math.PI + radDist;
353             for(int rayNum = 0; rayNum < (dist.length - 1); rayNum++) {
354                 if(searchDist == dist[rayNum + 1]
355                         && rayNum + 1 != dist.length - 1) {
356                     /* So we don't get 2 arrivals for the same ray. */
357                     continue;
358                 } else if((dist[rayNum] - searchDist)
359                         * (searchDist - dist[rayNum + 1]) >= 0.0) {
360                     /* look for distances that bracket the search distance */
361                     if((rayParams[rayNum] == rayParams[rayNum + 1])
362                             && rayParams.length > 2) {
363                         /*
364                          * Here we have a shadow zone, so it is not really an
365                          * arrival.
366                          */
367                         continue;
368                     }
369                     if(DEBUG) {
370                         System.err.println("SeismicPhase " + name
371                                 + ", found arrival:\n" + "dist "
372                                 + (float)(180 / Math.PI * dist[rayNum]) + " "
373                                 + (float)(180 / Math.PI * searchDist) + " "
374                                 + (float)(180 / Math.PI * dist[rayNum + 1]));
375                     }
376                     double arrivalTime = (searchDist - dist[rayNum])
377                             / (dist[rayNum + 1] - dist[rayNum])
378                             * (time[rayNum + 1] - time[rayNum]) + time[rayNum];
379                     double arrivalRayParam = (searchDist - dist[rayNum + 1])
380                             * (rayParams[rayNum] - rayParams[rayNum + 1])
381                             / (dist[rayNum] - dist[rayNum + 1])
382                             + rayParams[rayNum + 1];
383                     newArrival = new Arrival(this,
384                                              arrivalTime,
385                                              searchDist,
386                                              arrivalRayParam,
387                                              rayNum,
388                                              name,
389                                              puristName,
390                                              sourceDepth);
391                     arrivals.addElement(newArrival);
392                 }
393             }
394             /*
395              * Look for arrivals that are 2(n+1)Pi-radDist, ie rays that have
396              * done more than one half lap plus some number of whole laps.
397              */
398             searchDist = (n + 1) * 2.0 * Math.PI - radDist;
399             if(tempDeg != 180) {
400                 for(int rayNum = 0; rayNum < (dist.length - 1); rayNum++) {
401                     if(searchDist == dist[rayNum + 1]
402                             && rayNum + 1 != dist.length - 1) {
403                         /* So we don't get 2 arrivals for the same ray. */
404                         continue;
405                     } else if((dist[rayNum] - searchDist)
406                             * (searchDist - dist[rayNum + 1]) >= 0.0) {
407                         if((rayParams[rayNum] == rayParams[rayNum + 1])
408                                 && rayParams.length > 2) {
409                             /*
410                              * Here we have a shadow zone, so it is not really
411                              * an arrival.
412                              */
413                             continue;
414                         }
415                         if(DEBUG) {
416                             System.err.println("SeismicPhase " + name
417                                     + ", found arrival:\n" + "dist "
418                                     + (float)(180 / Math.PI * dist[rayNum])
419                                     + " " + (float)(180 / Math.PI * searchDist)
420                                     + " "
421                                     + (float)(180 / Math.PI * dist[rayNum + 1]));
422                         }
423                         double arrivalTime = (searchDist - dist[rayNum])
424                                 / (dist[rayNum + 1] - dist[rayNum])
425                                 * (time[rayNum + 1] - time[rayNum])
426                                 + time[rayNum];
427                         double arrivalRayParam = (searchDist - dist[rayNum])
428                                 * (rayParams[rayNum + 1] - rayParams[rayNum])
429                                 / (dist[rayNum + 1] - dist[rayNum])
430                                 + rayParams[rayNum];
431                         newArrival = new Arrival(this,
432                                                  arrivalTime,
433                                                  searchDist,
434                                                  arrivalRayParam,
435                                                  rayNum,
436                                                  name,
437                                                  puristName,
438                                                  sourceDepth);
439                         arrivals.addElement(newArrival);
440                     }
441                 }
442             }
443             n++;
444         }
445     }
446 
447     /***
448      * changes maxRayParam and minRayParam whenever there is a phase conversion.
449      * For instance, SKP needs to change the maxRayParam because there are SKS
450      * ray parameters that cannot propagate from the cmb into the mantle as a p
451      * wave.
452      */
453     protected void phaseConversion(TauModel tMod,
454                                    int fromBranch,
455                                    int endAction,
456                                    boolean isPtoS) throws TauModelException {
457         if(endAction == TURN) {
458             // can't phase convert for just a turn point
459             throw new TauModelException("Illegal endAction: endAction="
460                     + endAction
461                     + "\nphase conversion are not allowed at turn points.");
462         } else if(endAction == REFLECTTOP) {
463             maxRayParam = Math.min(maxRayParam, tMod.getTauBranch(fromBranch,
464                                                                   isPtoS)
465                     .getMaxRayParam());
466             maxRayParam = Math.min(maxRayParam, tMod.getTauBranch(fromBranch,
467                                                                   !isPtoS)
468                     .getMaxRayParam());
469         } else if(endAction == REFLECTBOT) {
470             maxRayParam = Math.min(maxRayParam, tMod.getTauBranch(fromBranch,
471                                                                   isPtoS)
472                     .getMinTurnRayParam());
473             maxRayParam = Math.min(maxRayParam, tMod.getTauBranch(fromBranch,
474                                                                   !isPtoS)
475                     .getMinTurnRayParam());
476         } else if(endAction == TRANSUP) {
477             maxRayParam = Math.min(maxRayParam, tMod.getTauBranch(fromBranch,
478                                                                   isPtoS)
479                     .getMaxRayParam());
480             maxRayParam = Math.min(maxRayParam,
481                                    tMod.getTauBranch(fromBranch - 1, !isPtoS)
482                                            .getMinTurnRayParam());
483         } else if(endAction == TRANSDOWN) {
484             maxRayParam = Math.min(maxRayParam, tMod.getTauBranch(fromBranch,
485                                                                   isPtoS)
486                     .getMinRayParam());
487             maxRayParam = Math.min(maxRayParam,
488                                    tMod.getTauBranch(fromBranch + 1, !isPtoS)
489                                            .getMaxRayParam());
490         } else {
491             throw new TauModelException("Illegal endAction: endAction="
492                     + endAction);
493         }
494     }
495 
496     /*
497      * Adds the branch numbers from startBranch to endBranch, inclusive, to
498      * branchSeq, in order. Also, currBranch is set correctly based on the value
499      * of endAction. endAction can be one of TRANSUP, TRANSDOWN, REFLECTTOP,
500      * REFLECTBOT, or TURN.
501      */
502     protected void addToBranch(TauModel tMod,
503                                int startBranch,
504                                int endBranch,
505                                boolean isPWave,
506                                int endAction) throws TauModelException {
507         int endOffset;
508         boolean isDownGoing;
509         this.endAction = endAction;
510         if(DEBUG) {
511             System.out.print("start=" + startBranch + " end=" + endBranch
512                     + " endOffset=");
513             if(endAction == TURN) {
514                 System.out.println("TURN");
515             } else if(endAction == REFLECTTOP) {
516                 System.out.println("REFLECTTOP");
517             } else if(endAction == REFLECTBOT) {
518                 System.out.println("REFLECTBOT");
519             } else if(endAction == TRANSUP) {
520                 System.out.println("TRANSUP");
521             } else if(endAction == TRANSDOWN) {
522                 System.out.println("TRANSDOWN");
523             } else {
524                 System.out.println(endAction);
525             }
526         }
527         if(endAction == TURN) {
528             endOffset = 0;
529             isDownGoing = true;
530             minRayParam = Math.max(minRayParam, tMod.getTauBranch(endBranch,
531                                                                   isPWave)
532                     .getMinTurnRayParam());
533         } else if(endAction == REFLECTTOP) {
534             endOffset = 0;
535             isDownGoing = false;
536             maxRayParam = Math.min(maxRayParam, tMod.getTauBranch(endBranch,
537                                                                   isPWave)
538                     .getMaxRayParam());
539         } else if(endAction == REFLECTBOT) {
540             endOffset = 0;
541             isDownGoing = true;
542             maxRayParam = Math.min(maxRayParam, tMod.getTauBranch(endBranch,
543                                                                   isPWave)
544                     .getMinTurnRayParam());
545         } else if(endAction == TRANSUP) {
546             endOffset = -1;
547             isDownGoing = false;
548             maxRayParam = Math.min(maxRayParam, tMod.getTauBranch(endBranch,
549                                                                   isPWave)
550                     .getMaxRayParam());
551         } else if(endAction == TRANSDOWN) {
552             endOffset = 1;
553             isDownGoing = true;
554             maxRayParam = Math.min(maxRayParam, tMod.getTauBranch(endBranch,
555                                                                   isPWave)
556                     .getMinRayParam());
557         } else {
558             throw new TauModelException("Illegal endAction: endAction="
559                     + endAction);
560         }
561         if(isDownGoing) {
562             /* Must be downgoing, so use i++. */
563             for(int i = startBranch; i <= endBranch; i++) {
564                 branchSeq.addElement(new Integer(i));
565                 downGoing.addElement(new Boolean(isDownGoing));
566                 waveType.addElement(new Boolean(isPWave));
567                 legAction.addElement(new Integer(endAction));
568             }
569             if(DEBUG) {
570                 for(int i = startBranch; i <= endBranch; i++) {
571                     System.out.println("i=" + i + " isDownGoing=" + isDownGoing
572                             + " isPWave=" + isPWave + " startBranch="
573                             + startBranch + " endBranch=" + endBranch + " "
574                             + endAction);
575                 }
576             }
577         } else {
578             /* Must be up going so use i--. */
579             for(int i = startBranch; i >= endBranch; i--) {
580                 branchSeq.addElement(new Integer(i));
581                 downGoing.addElement(new Boolean(isDownGoing));
582                 waveType.addElement(new Boolean(isPWave));
583                 legAction.addElement(new Integer(endAction));
584             }
585             if(DEBUG) {
586                 for(int i = startBranch; i >= endBranch; i--) {
587                     System.out.println("i=" + i + " isDownGoing=" + isDownGoing
588                             + " isPWave=" + isPWave + " startBranch="
589                             + startBranch + " endBranch=" + endBranch + " "
590                             + endAction);
591                 }
592             }
593         }
594         currBranch = endBranch + endOffset;
595     }
596 
597     /***
598      * Finds the closest discontinuity to the given depth that can have
599      * refletions and phase transformations.
600      * 
601      * @returns the branch number with the closest top depth.
602      */
603     public int closestBranchToDepth(TauModel tMod, String depthString) {
604         if(depthString.equals("m")) {
605             return tMod.getMohoBranch();
606         } else if(depthString.equals("c")) {
607             return tMod.getCmbBranch();
608         } else if(depthString.equals("i")) { return tMod.getIocbBranch(); }
609         // nonstandard boundary, given by a number, so we must look for it
610         int disconBranch = -1;
611         double disconMax = Double.MAX_VALUE;
612         double disconDepth = (Double.valueOf(depthString)).doubleValue();
613         TauBranch tBranch;
614         for(int i = 0; i < tMod.getNumBranches(); i++) {
615             tBranch = tMod.getTauBranch(i, PWAVE);
616             if(Math.abs(disconDepth - tBranch.getTopDepth()) < disconMax
617                     && !tMod.isNoDisconDepth(tBranch.getTopDepth())) {
618                 disconBranch = i;
619                 disconMax = Math.abs(disconDepth - tBranch.getTopDepth());
620             }
621         }
622         return disconBranch;
623     }
624 
625     /***
626      * Constructs a branch sequence from the given phase name and tau model.
627      */
628     protected void parseName(TauModel tMod) throws TauModelException {
629         String prevLeg;
630         String currLeg = (String)legs.elementAt(0);
631         String nextLeg = currLeg;
632         branchSeq.removeAllElements();
633         boolean isDownGoing;
634         boolean isPWave = PWAVE;
635         boolean isPWavePrev = isPWave;
636         /*
637          * Deal with surface wave velocities first, since they are a special
638          * case.
639          */
640         if(legs.size() == 2 && currLeg.endsWith("kmps")) { return; }
641         /* Make a check for J legs if the model doesn not allow J */
642         if(name.indexOf('J') != -1 && !tMod.sMod.isAllowInnerCoreS()) { throw new TauModelException("'J' phases were not created for this model: "
643                 + name); }
644         /* set currWave to be the wave type for this leg, 'P' or 'S'. */
645         if(currLeg.equals("p") || currLeg.startsWith("P")
646                 || currLeg.equals("K") || currLeg.equals("I")) {
647             isPWave = PWAVE;
648             isPWavePrev = isPWave;
649         } else if(currLeg.equals("s") || currLeg.startsWith("S")
650                 || currLeg.equals("J")) {
651             isPWave = SWAVE;
652             isPWavePrev = isPWave;
653         }
654         /*
655          * First, decide whether the ray is up going or downgoing from the
656          * source. If it is up going then the first branch number would be
657          * tMod.getSourceBranch()-1 and downgoing would be
658          * tMod.getSourceBranch().
659          */
660         if(currLeg.startsWith("P") || currLeg.startsWith("S")) {
661             // Downgoing from source
662             currBranch = tMod.getSourceBranch();
663             isDownGoing = true;
664             endAction = REFLECTBOT; // treat initial downgoing as if it were a
665             // underside reflection
666         } else if(currLeg.equals("p") || currLeg.equals("s")) {
667             // Up going from source
668             isDownGoing = false;
669             endAction = REFLECTTOP; // treat initial upgoing as if it were a
670             // topside reflection
671             if(tMod.getSourceBranch() != 0) {
672                 currBranch = tMod.getSourceBranch() - 1;
673             } else {
674                 /*
675                  * p and s for zero source depth are only at zero distance and
676                  * then can be called P or S.
677                  */
678                 maxRayParam = -1;
679                 minRayParam = -1;
680                 return;
681             }
682         } else {
683             throw new TauModelException("First phase not recognized: "
684                     + currLeg
685                     + " must be one of P, Pg, Pn, Pdiff, p or the S equivalents");
686         }
687         /*
688          * Set maxRayParam to be a horizontal ray leaving the source and set
689          * minRayParam to be a vertical (p=0) ray.
690          */
691         if(tMod.getSourceBranch() != 0) {
692             maxRayParam = Math.max(tMod.getTauBranch(tMod.getSourceBranch() - 1,
693                                                      isPWave)
694                                            .getMinTurnRayParam(),
695                                    tMod.getTauBranch(tMod.getSourceBranch(),
696                                                      isPWave).getMaxRayParam());
697         } else {
698             maxRayParam = tMod.getTauBranch(tMod.getSourceBranch(), isPWave)
699                     .getMaxRayParam();
700         }
701         minRayParam = 0.0;
702         int disconBranch = 0;
703         double legDepth, nextLegDepth = 0.0;
704         boolean isLegDepth, isNextLegDepth = false;
705         endAction = TRANSDOWN;
706         /*
707          * Now loop over all of the phase legs and construct the proper branch
708          * sequence.
709          */
710         currLeg = "START"; // So the prevLeg isn't wrong on the first pass
711         for(int legNum = 0; legNum < legs.size() - 1; legNum++) {
712             prevLeg = currLeg;
713             currLeg = nextLeg;
714             nextLeg = (String)legs.elementAt(legNum + 1);
715             if(DEBUG) {
716                 System.out.println(legNum + "  " + prevLeg + "  " + currLeg
717                         + "  " + nextLeg);
718             }
719             legDepth = nextLegDepth;
720             isLegDepth = isNextLegDepth;
721             // find out if the next leg represents a phase conversion depth
722             try {
723                 nextLegDepth = (new Double(nextLeg)).doubleValue();
724                 isNextLegDepth = true;
725             } catch(NumberFormatException e) {
726                 nextLegDepth = -1;
727                 isNextLegDepth = false;
728             }
729             /* set currWave to be the wave type for this leg, 'P' or 'S'. */
730             isPWavePrev = isPWave;
731             if(currLeg.equals("p") || currLeg.startsWith("P")
732                     || currLeg.equals("I")) {
733                 isPWave = PWAVE;
734             } else if(currLeg.equals("s") || currLeg.startsWith("S")
735                     || currLeg.equals("J")) {
736                 isPWave = SWAVE;
737             } else if(currLeg.equals("K")) {
738                 /*
739                  * here we want to use whatever isPWave was on the last leg so
740                  * do nothing. This makes sure we us the correct maxRayParam
741                  * from the correct TauBranch within the outer core. In other
742                  * words K has a high slowness zone if it entered the outer core
743                  * as a mantle P wave, but doesn't if it entered as a mantle S
744                  * wave. It shouldn't matter for inner core to outer core type
745                  * legs.
746                  */
747             }
748             // check to see if there has been a phase conversion
749             if(branchSeq.size() > 0 && isPWavePrev != isPWave) {
750                 phaseConversion(tMod,
751                                 ((Integer)branchSeq.elementAt(branchSeq.size() - 1)).intValue(),
752                                 endAction,
753                                 isPWavePrev);
754             }
755             /* Deal with p and s case first. */
756             if(currLeg.equals("p") || currLeg.equals("s")) {
757                 if(nextLeg.startsWith("v")) {
758                     throw new TauModelException("p and s must always be up going "
759                             + " and cannot come immediately before a top-side reflection."
760                             + " currLeg=" + currLeg + " nextLeg=" + nextLeg);
761                 } else if(nextLeg.startsWith("^")) {
762                     disconBranch = closestBranchToDepth(tMod,
763                                                         nextLeg.substring(1));
764                     if(currBranch >= disconBranch) {
765                         addToBranch(tMod,
766                                     currBranch,
767                                     disconBranch,
768                                     isPWave,
769                                     REFLECTTOP);
770                     } else {
771                         throw new TauModelException("Phase not recognized: "
772                                 + currLeg + " followed by " + nextLeg
773                                 + " when currBranch=" + currBranch
774                                 + " > disconBranch=" + disconBranch);
775                     }
776                 } else if(nextLeg.equals("m")
777                         && currBranch >= tMod.getMohoBranch()) {
778                     addToBranch(tMod,
779                                 currBranch,
780                                 tMod.getMohoBranch(),
781                                 isPWave,
782                                 TRANSUP);
783                 } else if(nextLeg.startsWith("P") || nextLeg.startsWith("S")
784                         || nextLeg.equals("END")) {
785                     addToBranch(tMod, currBranch, 0, isPWave, REFLECTTOP);
786                 } else if(isNextLegDepth) {
787                     disconBranch = closestBranchToDepth(tMod, nextLeg);
788                     addToBranch(tMod,
789                                 currBranch,
790                                 disconBranch,
791                                 isPWave,
792                                 TRANSUP);
793                 } else {
794                     throw new TauModelException("Phase not recognized: "
795                             + currLeg + " followed by " + nextLeg);
796                 }
797                 /* Now deal with P and S case. */
798             } else if(currLeg.equals("P") || currLeg.equals("S")) {
799                 if(nextLeg.equals("P") || nextLeg.equals("S")
800                         || nextLeg.equals("Pn") || nextLeg.equals("Sn")
801                         || nextLeg.equals("END")) {
802                     if(endAction == TRANSDOWN || endAction == REFLECTTOP) {
803                         // downgoing, so must first turn in mantle
804                         addToBranch(tMod,
805                                     currBranch,
806                                     tMod.getCmbBranch() - 1,
807                                     isPWave,
808                                     TURN);
809                     }
810                     addToBranch(tMod, currBranch, 0, isPWave, REFLECTTOP);
811                 } else if(nextLeg.startsWith("v")) {
812                     disconBranch = closestBranchToDepth(tMod,
813                                                         nextLeg.substring(1));
814                     if(currBranch <= disconBranch - 1) {
815                         addToBranch(tMod,
816                                     currBranch,
817                                     disconBranch - 1,
818                                     isPWave,
819                                     REFLECTBOT);
820                     } else {
821                         throw new TauModelException("Phase not recognized: "
822                                 + currLeg + " followed by " + nextLeg
823                                 + " when currBranch=" + currBranch
824                                 + " < disconBranch=" + disconBranch);
825                     }
826                 } else if(nextLeg.startsWith("^")) {
827                     disconBranch = closestBranchToDepth(tMod,
828                                                         nextLeg.substring(1));
829                     if(prevLeg.equals("K")) {
830                         addToBranch(tMod,
831                                     currBranch,
832                                     disconBranch,
833                                     isPWave,
834                                     REFLECTTOP);
835                     } else if(prevLeg.startsWith("^") || prevLeg.equals("P")
836                             || prevLeg.equals("S") || prevLeg.equals("p")
837                             || prevLeg.equals("s") || prevLeg.equals("START")) {
838                         addToBranch(tMod,
839                                     currBranch,
840                                     tMod.getCmbBranch() - 1,
841                                     isPWave,
842                                     TURN);
843                         addToBranch(tMod,
844                                     currBranch,
845                                     disconBranch,
846                                     isPWave,
847                                     REFLECTTOP);
848                     } else if((prevLeg.startsWith("v") && disconBranch < closestBranchToDepth(tMod,
849                                                                                               prevLeg.substring(1)))
850                             || (prevLeg.equals("m") && disconBranch < tMod.getMohoBranch())
851                             || (prevLeg.equals("c") && disconBranch < tMod.getCmbBranch())) {
852                         addToBranch(tMod,
853                                     currBranch,
854                                     disconBranch,
855                                     isPWave,
856                                     REFLECTTOP);
857                     } else {
858                         throw new TauModelException("Phase not recognized: "
859                                 + currLeg + " followed by " + nextLeg
860                                 + " when currBranch=" + currBranch
861                                 + " > disconBranch=" + disconBranch);
862                     }
863                 } else if(nextLeg.equals("c")) {
864                     addToBranch(tMod,
865                                 currBranch,
866                                 tMod.getCmbBranch() - 1,
867                                 isPWave,
868                                 REFLECTBOT);
869                 } else if(nextLeg.equals("K")) {
870                     addToBranch(tMod,
871                                 currBranch,
872                                 tMod.getCmbBranch() - 1,
873                                 isPWave,
874                                 TRANSDOWN);
875                 } else if(nextLeg.equals("m")
876                         || (isNextLegDepth && nextLegDepth < tMod.getCmbDepth())) {
877                     // treat the moho in the same wasy as 410 type
878                     // discontinuities
879                     disconBranch = closestBranchToDepth(tMod, nextLeg);
880                     if(DEBUG) {
881                         System.out.println("DisconBranch=" + disconBranch
882                                 + " for " + nextLeg);
883                         System.out.println(tMod.getTauBranch(disconBranch,
884                                                              isPWave)
885                                 .getTopDepth());
886                     }
887                     if(endAction == TURN || endAction == REFLECTBOT
888                             || endAction == TRANSUP) {
889                         // upgoing section
890                         if(disconBranch > currBranch) {
891                             // check for discontinuity below the current branch
892                             // when the ray should be upgoing
893                             throw new TauModelException("Phase not recognized: "
894                                     + currLeg
895                                     + " followed by "
896                                     + nextLeg
897                                     + " when currBranch="
898                                     + currBranch
899                                     + " > disconBranch=" + disconBranch);
900                         }
901                         addToBranch(tMod,
902                                     currBranch,
903                                     disconBranch,
904                                     isPWave,
905                                     TRANSUP);
906                     } else {
907                         // downgoing section, must look at the leg after the
908                         // next
909                         // leg to determine whether to convert on the downgoing
910                         // or
911                         // upgoing part of the path
912                         String nextNextLeg = (String)legs.elementAt(legNum + 2);
913                         if(nextNextLeg.equals("p") || nextNextLeg.equals("s")) {
914                             // convert on upgoing section
915                             addToBranch(tMod,
916                                         currBranch,
917                                         tMod.getCmbBranch() - 1,
918                                         isPWave,
919                                         TURN);
920                             addToBranch(tMod,
921                                         currBranch,
922                                         disconBranch,
923                                         isPWave,
924                                         TRANSUP);
925                         } else if(nextNextLeg.equals("P")
926                                 || nextNextLeg.equals("S")) {
927                             if(disconBranch > currBranch) {
928                                 // discon is below current loc
929                                 addToBranch(tMod,
930                                             currBranch,
931                                             disconBranch - 1,
932                                             isPWave,
933                                             TRANSDOWN);
934                             } else {
935                                 // discon is above current loc, but we have a
936                                 // downgoing ray, so this is an illegal ray for
937                                 // this source depth
938                                 maxRayParam = -1;
939                                 if(DEBUG) {
940                                     System.out.println("Cannot phase convert on the "
941                                             + "downgoing side if the discontinuity is above "
942                                             + "the phase leg starting point, "
943                                             + currLeg
944                                             + " "
945                                             + nextLeg
946                                             + " "
947                                             + nextNextLeg
948                                             + ", so this phase, "
949                                             + getName()
950                                             + " is illegal for this sourceDepth.");
951                                 }
952                                 return;
953                             }
954                         } else {
955                             throw new TauModelException("Phase not recognized: "
956                                     + currLeg
957                                     + " followed by "
958                                     + nextLeg
959                                     + " followed by " + nextNextLeg);
960                         }
961                     }
962                 } else {
963                     throw new TauModelException("Phase not recognized: "
964                             + currLeg + " followed by " + nextLeg);
965                 }
966             } else if(currLeg.startsWith("P") || currLeg.startsWith("S")) {
967                 if(currLeg.equals("Pdiff") || currLeg.equals("Sdiff")) {
968                     /*
969                      * in the diffracted case we trick addToBranch into thinking
970                      * we are turning, but then make the maxRayParam equal to
971                      * minRayParam, which is the deepest turning ray.
972                      */
973                     if(maxRayParam >= tMod.getTauBranch(tMod.getCmbBranch() - 1,
974                                                         isPWave)
975                             .getMinTurnRayParam()
976                             && minRayParam <= tMod.getTauBranch(tMod.getCmbBranch() - 1,
977                                                                 isPWave)
978                                     .getMinTurnRayParam()) {
979                         addToBranch(tMod,
980                                     currBranch,
981                                     tMod.getCmbBranch() - 1,
982                                     isPWave,
983                                     TURN);
984                         maxRayParam = minRayParam;
985                         if(nextLeg.equals("END")) {
986                             addToBranch(tMod,
987                                         currBranch,
988                                         0,
989                                         isPWave,
990                                         REFLECTTOP);
991                         }
992                     } else {
993                         // can't have head wave as ray param is not within range
994                         maxRayParam = -1;
995                         if(DEBUG) {
996                             System.out.println("Cannot have the head wave "
997                                     + currLeg + " within phase " + name
998                                     + " for this sourceDepth and/or path.");
999                         }
1000                         return;
1001                     }
1002                 } else if(currLeg.equals("Pg") || currLeg.equals("Sg")
1003                         || currLeg.equals("Pn") || currLeg.equals("Sn")) {
1004                     if(currBranch >= tMod.getMohoBranch()) {
1005                         /*
1006                          * Pg, Pn, Sg and Sn must be above the moho and so is
1007                          * not valid for rays coming upwards from below,
1008                          * possibly due to the source depth. Setting maxRayParam =
1009                          * -1 effectively disallows this phase.
1010                          */
1011                         maxRayParam = -1;
1012                         if(DEBUG) {
1013                             System.out.println("(currBranch >= tMod.getMohoBranch() "
1014                                     + currBranch
1015                                     + " "
1016                                     + tMod.getMohoBranch()
1017                                     + " so there cannot be a "
1018                                     + currLeg
1019                                     + " phase for this sourceDepth and/or path.");
1020                         }
1021                         return;
1022                     }
1023                     if(currLeg.equals("Pg") || currLeg.equals("Sg")) {
1024                         addToBranch(tMod,
1025                                     currBranch,
1026                                     tMod.getMohoBranch() - 1,
1027                                     isPWave,
1028                                     TURN);
1029                         addToBranch(tMod, currBranch, 0, isPWave, REFLECTTOP);
1030                     } else if(currLeg.equals("Pn") || currLeg.equals("Sn")) {
1031                         /*
1032                          * in the refracted case we trick addToBranch into
1033                          * thinking we are turning below the moho, but then make
1034                          * the minRayParam equal to maxRayParam, which is the
1035                          * head wave ray.
1036                          */
1037                         if(maxRayParam >= tMod.getTauBranch(tMod.getMohoBranch(),
1038                                                             isPWave)
1039                                 .getMaxRayParam()
1040                                 && minRayParam <= tMod.getTauBranch(tMod.getMohoBranch(),
1041                                                                     isPWave)
1042                                         .getMaxRayParam()) {
1043                             addToBranch(tMod,
1044                                         currBranch,
1045                                         tMod.getMohoBranch(),
1046                                         isPWave,
1047                                         TURN);
1048                             addToBranch(tMod,
1049                                         currBranch,
1050                                         tMod.getMohoBranch(),
1051                                         isPWave,
1052                                         TRANSUP);
1053                             minRayParam = maxRayParam;
1054                             if(nextLeg.equals("END") || nextLeg.startsWith("P")
1055                                     || nextLeg.startsWith("S")) {
1056                                 addToBranch(tMod,
1057                                             currBranch,
1058                                             0,
1059                                             isPWave,
1060                                             REFLECTTOP);
1061                             }
1062                         } else {
1063                             // can't have head wave as ray param is not within
1064                             // range
1065                             maxRayParam = -1;
1066                             if(DEBUG) {
1067                                 System.out.println("Cannot have the head wave "
1068                                         + currLeg + " within phase " + name
1069                                         + " for this sourceDepth and/or path.");
1070                             }
1071                             return;
1072                         }
1073                     }
1074                 } else {
1075                     throw new TauModelException("Phase not recognized: "
1076                             + currLeg + " followed by " + nextLeg);
1077                 }
1078             } else if(currLeg.equals("K")) {
1079                 /* Now deal with K. */
1080                 if(nextLeg.equals("P") || nextLeg.equals("S")) {
1081                     if(prevLeg.equals("P") || prevLeg.equals("S")
1082                             || prevLeg.equals("K")) {
1083                         addToBranch(tMod,
1084                                     currBranch,
1085                                     tMod.getIocbBranch() - 1,
1086                                     isPWave,
1087                                     TURN);
1088                     }
1089                     addToBranch(tMod,
1090                                 currBranch,
1091                                 tMod.getCmbBranch(),
1092                                 isPWave,
1093                                 TRANSUP);
1094                 } else if(nextLeg.equals("K")) {
1095                     if(prevLeg.equals("P") || prevLeg.equals("S")
1096                             || prevLeg.equals("K")) {
1097                         addToBranch(tMod,
1098                                     currBranch,
1099                                     tMod.getIocbBranch() - 1,
1100                                     isPWave,
1101                                     TURN);
1102                     }
1103                     addToBranch(tMod,
1104                                 currBranch,
1105                                 tMod.getCmbBranch(),
1106                                 isPWave,
1107                                 REFLECTTOP);
1108                 } else if(nextLeg.equals("I") || nextLeg.equals("J")) {
1109                     addToBranch(tMod,
1110                                 currBranch,
1111                                 tMod.getIocbBranch() - 1,
1112                                 isPWave,
1113                                 TRANSDOWN);
1114                 } else if(nextLeg.equals("i")) {
1115                     addToBranch(tMod,
1116                                 currBranch,
1117                                 tMod.getIocbBranch() - 1,
1118                                 isPWave,
1119                                 REFLECTBOT);
1120                 } else {
1121                     throw new TauModelException("Phase not recognized: "
1122                             + currLeg + " followed by " + nextLeg);
1123                 }
1124             } else if(currLeg.equals("I") || currLeg.equals("J")) {
1125                 /* And now consider inner core, I and J. */
1126                 addToBranch(tMod,
1127                             currBranch,
1128                             tMod.getNumBranches() - 1,
1129                             isPWave,
1130                             TURN);
1131                 if(nextLeg.equals("I") || nextLeg.equals("J")) {
1132                     addToBranch(tMod,
1133                                 currBranch,
1134                                 tMod.getIocbBranch(),
1135                                 isPWave,
1136                                 REFLECTTOP);
1137                 } else if(nextLeg.equals("K")) {
1138                     addToBranch(tMod,
1139                                 currBranch,
1140                                 tMod.getIocbBranch(),
1141                                 isPWave,
1142                                 TRANSUP);
1143                 }
1144             } else if(currLeg.equals("m")) {} else if(currLeg.equals("c")) {} else if(currLeg.equals("i")) {} else if(currLeg.startsWith("^")) {} else if(currLeg.startsWith("v")) {} else if(isLegDepth) {} else {
1145                 throw new TauModelException("Phase not recognized: " + currLeg
1146                         + " followed by " + nextLeg);
1147             }
1148         }
1149     }
1150 
1151     /***
1152      * Tokenizes a phase name into legs, ie PcS becomes 'P'+'c'+'S' while p^410P
1153      * would become 'p'+'^410'+'P'. Once a phase name has been broken into
1154      * tokens we can begin to construct the sequence of branches to which it
1155      * corresponds. Only minor error checking is done at this point, for
1156      * instance PIP generates an exception but ^410 doesn't. It also appends
1157      * "END" as the last leg.
1158      * 
1159      * @throws TauModelException
1160      *             if the phase name cannot be tokenized.
1161      */
1162     protected void legPuller() throws TauModelException {
1163         int offset = 0;
1164         legs.removeAllElements();
1165         /* Special case for surface wave velocity. */
1166         if(name.endsWith("kmps")) {
1167             try {
1168                 double vel = Double.valueOf(name.substring(0, name.length() - 4))
1169                         .doubleValue();
1170                 legs.addElement(name);
1171             } catch(NumberFormatException e) {
1172                 throw new TauModelException("Invalid phase name:\n" + name);
1173             }
1174         } else while(offset < name.length()) {
1175             if(name.charAt(offset) == 'K' || name.charAt(offset) == 'I'
1176                     || name.charAt(offset) == 'J' || name.charAt(offset) == 'p'
1177                     || name.charAt(offset) == 's' || name.charAt(offset) == 'm'
1178                     || name.charAt(offset) == 'c' || name.charAt(offset) == 'i') {
1179                 // Do the easy ones, ie K,I,J,p,s,m,c,i
1180                 legs.addElement(name.substring(offset, offset + 1));
1181                 offset = offset + 1;
1182             } else if(name.charAt(offset) == 'P' || name.charAt(offset) == 'S') {
1183                 /*
1184                  * Now it gets complicated, first see if the next char is part
1185                  * of a different leg or we are at the end.
1186                  */
1187                 if(offset + 1 == name.length()
1188                         || name.charAt(offset + 1) == 'P'
1189                         || name.charAt(offset + 1) == 'S'
1190                         || name.charAt(offset + 1) == 'K'
1191                         || name.charAt(offset + 1) == 'p'
1192                         || name.charAt(offset + 1) == 's'
1193                         || name.charAt(offset + 1) == 'm'
1194                         || name.charAt(offset + 1) == 'c'
1195                         || name.charAt(offset + 1) == '^'
1196                         || name.charAt(offset + 1) == 'v'
1197                         || Character.isDigit(name.charAt(offset + 1))) {
1198                     legs.addElement(name.substring(offset, offset + 1));
1199                     offset++;
1200                 } else if(name.charAt(offset + 1) == 'g'
1201                         || name.charAt(offset + 1) == 'b'
1202                         || name.charAt(offset + 1) == 'n') {
1203                     /* The leg is not described by one letter, check for 2. */
1204                     legs.addElement(name.substring(offset, offset + 2));
1205                     offset = offset + 2;
1206                 } else if(name.length() >= offset + 5
1207                         && (name.substring(offset, offset + 5).equals("Sdiff") || name.substring(offset,
1208                                                                                                  offset + 5)
1209                                 .equals("Pdiff"))) {
1210                     legs.addElement(name.substring(offset, offset + 5));
1211                     offset = offset + 5;
1212                 } else {
1213                     throw new TauModelException("Invalid phase name:\n"
1214                             + name.substring(offset) + " in " + name);
1215                 }
1216             } else if(name.charAt(offset) == '^' || name.charAt(offset) == 'v') {
1217                 /*
1218                  * Top side or bottom side reflections, check for standard
1219                  * boundaries and then check for numerical ones.
1220                  */
1221                 if(name.charAt(offset + 1) == 'm'
1222                         || name.charAt(offset + 1) == 'c'
1223                         || name.charAt(offset + 1) == 'i') {
1224                     legs.addElement(name.substring(offset, offset + 2));
1225                     offset = offset + 2;
1226                 } else if(Character.isDigit(name.charAt(offset + 1))) {
1227                     int numDigits = 1;
1228                     while(offset + numDigits + 1 < name.length()
1229                             && Character.isDigit(name.charAt(offset + numDigits
1230                                     + 1))) {
1231                         numDigits++;
1232                     }
1233                     legs.addElement(name.substring(offset, offset + numDigits
1234                             + 1));
1235                     offset = offset + numDigits + 1;
1236                 } else {
1237                     throw new TauModelException("Invalid phase name:\n"
1238                             + name.substring(offset) + " in " + name);
1239                 }
1240             } else if(Character.isDigit(name.charAt(offset))
1241                     || name.charAt(offset) == '.') {
1242                 String numString = name.substring(offset, offset + 1);
1243                 offset++;
1244                 while(Character.isDigit(name.charAt(offset))
1245                         || name.charAt(offset) == '.') {
1246                     numString += name.substring(offset, offset + 1);
1247                     offset++;
1248                 }
1249                 try {
1250                     Double d = new Double(numString);
1251                     legs.addElement(numString);
1252                 } catch(NumberFormatException e) {
1253                     throw new TauModelException("Invalid phase name: "
1254                             + numString + "\n" + e.getMessage() + " in " + name);
1255                 }
1256             } else {
1257                 throw new TauModelException("Invalid phase name:\n"
1258                         + name.substring(offset) + " in " + name);
1259             }
1260         }
1261         legs.addElement(new String("END"));
1262         if(!phaseValidate()) { throw new TauModelException("Phase failed validation: "
1263                 + name); }
1264     }
1265 
1266     protected void createPuristName(TauModel tMod) {
1267         String currLeg = (String)legs.elementAt(0);
1268         /*
1269          * Deal with surface wave velocities first, since they are a special
1270          * case.
1271          */
1272         if(legs.size() == 2 && currLeg.endsWith("kmps")) {
1273             puristName = name;
1274             return;
1275         }
1276         puristName = "";
1277         double legDepth;
1278         int intLegDepth;
1279         int disconBranch;
1280         // only loop to size()-1 as last leg is always "END"
1281         for(int legNum = 0; legNum < legs.size() - 1; legNum++) {
1282             currLeg = (String)legs.elementAt(legNum);
1283             // find out if the next leg represents a
1284             // phase conversion or reflection depth
1285             if(currLeg.startsWith("v") || currLeg.startsWith("^")) {
1286                 disconBranch = closestBranchToDepth(tMod, currLeg.substring(1));
1287                 legDepth = tMod.getTauBranch(disconBranch, true).getTopDepth();
1288                 puristName += currLeg.substring(0, 1);
1289                 if(legDepth == Math.rint(legDepth)) {
1290                     intLegDepth = (int)legDepth;
1291                     puristName += intLegDepth;
1292                 } else {
1293                     puristName += legDepth;
1294                 }
1295             } else {
1296                 try {
1297                     legDepth = (new Double(currLeg)).doubleValue();
1298                     // only get this far if the currLeg is a number,
1299                     // otherwise exception
1300                     disconBranch = closestBranchToDepth(tMod, currLeg);
1301                     legDepth = tMod.getTauBranch(disconBranch, true)
1302                             .getTopDepth();
1303                     if(legDepth == Math.rint(legDepth)) {
1304                         intLegDepth = (int)legDepth;
1305                         puristName += intLegDepth;
1306                     } else {
1307                         puristName += legDepth;
1308                     }
1309                 } catch(NumberFormatException e) {
1310                     puristName += currLeg;
1311                 }
1312             }
1313         }
1314     }
1315 
1316     /***
1317      * Sums the appropriate branches for this phase.
1318      * 
1319      * @throws TauModelException
1320      *             if the topDepth of the high slowness zone is not contained
1321      *             within the TauModel. This should never happen and would
1322      *             indicate an invalid TauModel.
1323      */
1324     protected void sumBranches(TauModel tMod) throws TauModelException {
1325         if(maxRayParam < 0.0 || minRayParam > maxRayParam) {
1326             /* Phase has no arrivals, possibly due to source depth. */
1327             rayParams = new double[0];
1328             minRayParam = -1;
1329             maxRayParam = -1;
1330             dist = new double[0];
1331             time = new double[0];
1332             maxDistance = -1;
1333             return;
1334         }
1335         /* Special case for surface waves. */
1336         if(name.endsWith("kmps")) {
1337             dist = new double[2];
1338             time = new double[2];
1339             rayParams = new double[2];
1340             dist[0] = 0.0;
1341             time[0] = 0.0;
1342             rayParams[0] = tMod.radiusOfEarth
1343                     / Double.valueOf(name.substring(0, name.length() - 4))
1344                             .doubleValue();
1345             dist[1] = 2 * Math.PI;
1346             time[1] = 2
1347                     * Math.PI
1348                     * tMod.radiusOfEarth
1349                     / Double.valueOf(name.substring(0, name.length() - 4))
1350                             .doubleValue();
1351             rayParams[1] = rayParams[0];
1352             minDistance = 0.0;
1353             maxDistance = 2 * Math.PI;
1354             return;
1355         }
1356         /*
1357          * Find the ray parameter index that corresponds to the minRayParam and
1358          * maxRayParam.
1359          */
1360         for(int i = 0; i < tMod.rayParams.length; i++) {
1361             if(tMod.rayParams[i] >= minRayParam) {
1362                 minRayParamIndex = i;
1363             }
1364             if(tMod.rayParams[i] >= maxRayParam) {
1365                 maxRayParamIndex = i;
1366             }
1367         }
1368         if(maxRayParamIndex == 0
1369                 && minRayParamIndex == tMod.rayParams.length - 1) {
1370             // all ray parameters are valid so just copy
1371             rayParams = new double[tMod.rayParams.length];
1372             System.arraycopy(tMod.rayParams,
1373                              0,
1374                              rayParams,
1375                              0,
1376                              tMod.rayParams.length);
1377         } else if(maxRayParamIndex == minRayParamIndex) {
1378             if(name.indexOf("Sdiff") != -1 || name.indexOf("Pdiff") != -1) {
1379                 rayParams = new double[2];
1380                 rayParams[0] = minRayParam;
1381                 rayParams[1] = minRayParam;
1382             } else if(name.indexOf("Pn") != -1 || name.indexOf("Sn") != -1) {
1383                 rayParams = new double[2];
1384                 rayParams[0] = minRayParam;
1385                 rayParams[1] = minRayParam;
1386             } else if(name.endsWith("kmps")) {
1387                 rayParams = new double[2];
1388                 rayParams[0] = 0;
1389                 rayParams[1] = maxRayParam;
1390             } else {
1391                 rayParams = new double[2];
1392                 rayParams[0] = minRayParam;
1393                 rayParams[1] = minRayParam;
1394             }
1395         } else {
1396             if(DEBUG) {
1397                 System.out.println("maxRayParamIndex=" + maxRayParamIndex
1398                         + " minRayParamIndex=" + minRayParamIndex
1399                         + " tMod.rayParams.length=" + tMod.rayParams.length
1400                         + " tMod.rayParams[0]=" + tMod.rayParams[0]
1401                         + " maxRayParam=" + maxRayParam);
1402             }
1403             // only a subset of ray parameters are valid so only use those
1404             rayParams = new double[minRayParamIndex - maxRayParamIndex + 1];
1405             System.arraycopy(tMod.rayParams,
1406                              maxRayParamIndex,
1407                              rayParams,
1408                              0,
1409                              minRayParamIndex - maxRayParamIndex + 1);
1410         }
1411         dist = new double[rayParams.length];
1412         time = new double[rayParams.length];
1413         /* initialize the counter for each branch to 0. 0 is P and 1 is S. */
1414         int[][] timesBranches = new int[2][tMod.getNumBranches()];
1415         for(int i = 0; i < timesBranches[0].length; i++) {
1416             timesBranches[0][i] = 0;
1417             timesBranches[1][i] = 0;
1418         }
1419         /* Count how many times each branch appears in the path. */
1420         for(int i = 0; i < branchSeq.size(); i++) {
1421             if(((Boolean)waveType.elementAt(i)).booleanValue()) {
1422                 timesBranches[0][((Integer)branchSeq.elementAt(i)).intValue()]++;
1423             } else {
1424                 timesBranches[1][((Integer)branchSeq.elementAt(i)).intValue()]++;
1425             }
1426         }
1427         /* Sum the branches with the appropriate multiplier. */
1428         for(int j = 0; j < tMod.getNumBranches(); j++) {
1429             if(timesBranches[0][j] != 0) {
1430                 for(int i = maxRayParamIndex; i < minRayParamIndex + 1; i++) {
1431                     dist[i - maxRayParamIndex] += timesBranches[0][j]
1432                             * tMod.getTauBranch(j, PWAVE).getDist(i);
1433                     time[i - maxRayParamIndex] += timesBranches[0][j]
1434                             * tMod.getTauBranch(j, PWAVE).time[i];
1435                 }
1436             }
1437             if(timesBranches[1][j] != 0) {
1438                 for(int i = maxRayParamIndex; i < minRayParamIndex + 1; i++) {
1439                     dist[i - maxRayParamIndex] += timesBranches[1][j]
1440                             * tMod.getTauBranch(j, SWAVE).getDist(i);
1441                     time[i - maxRayParamIndex] += timesBranches[1][j]
1442                             * tMod.getTauBranch(j, SWAVE).time[i];
1443                 }
1444             }
1445         }
1446         if(name.indexOf("Sdiff") != -1 || name.indexOf("Pdiff") != -1) {
1447             if(tMod.sMod.depthInHighSlowness(tMod.cmbDepth - 1e-10,
1448                                              minRayParam,
1449                                              (name.charAt(0) == 'P'))) {
1450                 /*
1451                  * No diffraction if there is a high slowness zone at the CMB.
1452                  */
1453                 minRayParam = -1;
1454                 maxRayParam = -1;
1455                 maxDistance = -1;
1456                 dist = new double[0];
1457                 time = new double[0];
1458                 rayParams = new double[0];
1459                 return;
1460             } else {
1461                 dist[1] = dist[0] + maxDiffraction * Math.PI / 180.0;
1462                 time[1] = time[0] + maxDiffraction * Math.PI / 180.0
1463                         * minRayParam;
1464             }
1465         } else if(name.indexOf("Pn") != -1 || name.indexOf("Sn") != -1) {
1466             dist[1] = dist[0] + maxRefraction * Math.PI / 180.0;
1467             time[1] = time[0] + maxRefraction * Math.PI / 180.0 * minRayParam;
1468         } else if(maxRayParamIndex == minRayParamIndex) {
1469             dist[1] = dist[0];
1470             time[1] = time[0];
1471         }
1472         minDistance = Double.MAX_VALUE;
1473         maxDistance = 0.0;
1474         for(int j = 0; j < dist.length; j++) {
1475             if(dist[j] < minDistance) {
1476                 minDistance = dist[j];
1477             }
1478             if(dist[j] > maxDistance) {
1479                 maxDistance = dist[j];
1480             }
1481         }
1482         /*
1483          * Now check to see if our ray parameter range includes any ray
1484          * parameters that are associated with high slowness zones. If so, then
1485          * we will need to insert a "shadow zone" into our time and distance
1486          * arrays. It is represented by a repeated ray parameter.
1487          */
1488         DepthRange[] hsz;
1489         int hSZIndex;
1490         TimeDist tempTD;
1491         int layerNum;
1492         int indexOffset;
1493         boolean foundOverlap = false;
1494         boolean isPWave;
1495         int branchNum;
1496         int dummy;
1497         for(dummy = 0, isPWave = true; dummy < 2; dummy++, isPWave = false) {
1498             hsz = tMod.sMod.getHighSlowness(isPWave);
1499             hSZIndex = 0;
1500             indexOffset = 0;
1501             for(int i = 0; i < hsz.length; i++) {
1502                 if(maxRayParam > hsz[i].rayParam
1503                         && hsz[i].rayParam > minRayParam) {
1504                     /*
1505                      * There is a high slowness zone within our ray parameter
1506                      * range so we might need to add a shadow zone. We need to
1507                      * check to see if this wave type, P or S, is part of the
1508                      * phase at this depth/ray parameter.
1509                      */
1510                     branchNum = tMod.findBranch(hsz[i].topDepth);
1511                     foundOverlap = false;
1512                     for(int legNum = 0; legNum < branchSeq.size(); legNum++) {
1513                         // check for downgoing legs that cross the high slowness
1514                         // zone
1515                         // with the same wave type
1516                         if(((Integer)branchSeq.elementAt(legNum)).intValue() == branchNum
1517                                 && ((Boolean)waveType.elementAt(legNum)).booleanValue() == isPWave
1518                                 && ((Boolean)downGoing.elementAt(legNum)).booleanValue() == true
1519                                 && ((Integer)branchSeq.elementAt(legNum - 1)).intValue() == branchNum - 1
1520                                 && ((Boolean)waveType.elementAt(legNum - 1)).booleanValue() == isPWave
1521                                 && ((Boolean)downGoing.elementAt(legNum - 1)).booleanValue() == true) {
1522                             foundOverlap = true;
1523                             break;
1524                         }
1525                     }
1526                     if(foundOverlap) {
1527                         double[] newdist = new double[dist.length + 1];
1528                         double[] newtime = new double[time.length + 1];
1529                         double[] newrayParams = new double[rayParams.length + 1];
1530                         for(int j = 0; j < rayParams.length; j++) {
1531                             if(rayParams[j] == hsz[i].rayParam) {
1532                                 hSZIndex = j;
1533                                 break;
1534                             }
1535                         }
1536                         System.arraycopy(dist, 0, newdist, 0, hSZIndex);
1537                         System.arraycopy(time, 0, newtime, 0, hSZIndex);
1538                         System.arraycopy(rayParams,
1539                                          0,
1540                                          newrayParams,
1541                                          0,
1542                                          hSZIndex);
1543                         newrayParams[hSZIndex] = hsz[i].rayParam;
1544                         /* Sum the branches with the appropriate multiplier. */
1545                         newdist[hSZIndex] = 0.0;
1546                         newtime[hSZIndex] = 0.0;
1547                         for(int j = 0; j < tMod.getNumBranches(); j++) {
1548                             if(timesBranches[0][j] != 0
1549                                     && tMod.getTauBranch(j, PWAVE)
1550                                             .getTopDepth() < hsz[i].topDepth) {
1551                                 newdist[hSZIndex] += timesBranches[0][j]
1552                                         * tMod.getTauBranch(j, PWAVE).dist[maxRayParamIndex
1553                                                 + hSZIndex - indexOffset];
1554                                 newtime[hSZIndex] += timesBranches[0][j]
1555                                         * tMod.getTauBranch(j, PWAVE).time[maxRayParamIndex
1556                                                 + hSZIndex - indexOffset];
1557                             }
1558                             if(timesBranches[1][j] != 0
1559                                     && tMod.getTauBranch(j, SWAVE)
1560                                             .getTopDepth() < hsz[i].topDepth) {
1561                                 newdist[hSZIndex] += timesBranches[1][j]
1562                                         * tMod.getTauBranch(j, SWAVE).dist[maxRayParamIndex
1563                                                 + hSZIndex - indexOffset];
1564                                 newtime[hSZIndex] += timesBranches[1][j]
1565                                         * tMod.getTauBranch(j, SWAVE).time[maxRayParamIndex
1566                                                 + hSZIndex - indexOffset];
1567                             }
1568                         }
1569                         System.arraycopy(dist,
1570                                          hSZIndex,
1571                                          newdist,
1572                                          hSZIndex + 1,
1573                                          dist.length - hSZIndex);
1574                         System.arraycopy(time,
1575                                          hSZIndex,
1576                                          newtime,
1577                                          hSZIndex + 1,
1578                                          time.length - hSZIndex);
1579                         System.arraycopy(rayParams,
1580                                          hSZIndex,
1581                                          newrayParams,
1582                                          hSZIndex + 1,
1583                                          rayParams.length - hSZIndex);
1584                         indexOffset++;
1585                         dist = newdist;
1586                         time = newtime;
1587                         rayParams = newrayParams;
1588                     }
1589                 }
1590             }
1591         }
1592     }
1593 
1594     /***
1595      * Calculates the "pierce points" for the arrivals stored in arrivals. The
1596      * pierce points are stored within each arrival object.
1597      */
1598     public void calcPierce(TauModel tMod) throws TauModelException {
1599         double rayParamA, rayParamB;
1600         double distRayParam;
1601         double distA, distB;
1602         double timeA, timeB;
1603         double distRatio;
1604         int rayNum = 0;
1605         int branchNum;
1606         boolean isPWave;
1607         int numAdded;
1608         double branchDist, branchTime, branchDepth, prevBranchTime;
1609         double turnDepth;
1610         Arrival currArrival;
1611         for(int arrivalNum = 0; arrivalNum < arrivals.size(); arrivalNum++) {
1612             currArrival = (Arrival)arrivals.elementAt(arrivalNum);
1613             branchDist = 0.0;
1614             branchTime = 0.0;
1615             prevBranchTime = branchTime;
1616             /*
1617              * Find the ray parameter index that corresponds to the arrival ray
1618              * parameter in the TauModel, ie it is between rayNum and rayNum+1,
1619              * We know that it must be <tMod.rayParams.length-1 since the last
1620              * ray parameter sample is 0, at least in a spherical model...
1621              */
1622             for(int i = 0; i < tMod.rayParams.length - 1; i++) {
1623                 if(tMod.rayParams[i] >= currArrival.rayParam) {
1624                     rayNum = i;
1625                 } else {
1626                     break;
1627                 }
1628             }
1629             currArrival.pierce = new TimeDist[branchSeq.size() + 2];
1630             // here we use ray parameter and dist info stored within the
1631             // SeismicPhase so we can use currArrival.rayParamIndex, which
1632             // may not correspond to rayNum (for tMod.rayParams).
1633             rayParamA = rayParams[currArrival.rayParamIndex];
1634             rayParamB = rayParams[currArrival.rayParamIndex + 1];
1635             distA = dist[currArrival.rayParamIndex];
1636             distB = dist[currArrival.rayParamIndex + 1];
1637             distRatio = (currArrival.dist - distA) / (distB - distA);
1638             distRayParam = distRatio * (rayParamB - rayParamA) + rayParamA;
1639             /* First pierce point is always 0 distance at the source depth. */
1640             currArrival.pierce[0] = new TimeDist(distRayParam,
1641                                                  0.0,
1642                                                  0.0,
1643                                                  tMod.sourceDepth);
1644             numAdded = 1;
1645             /*
1646              * Loop from 0 but already done 0, so the pierce point when the ray
1647              * leaves branch i is stored in i+1. Use linear interpolation
1648              * between rays that we know.
1649              */
1650             for(int i = 0; i < branchSeq.size(); i++) {
1651                 branchNum = ((Integer)branchSeq.elementAt(i)).intValue();
1652                 isPWave = ((Boolean)waveType.elementAt(i)).booleanValue();
1653                 if(DEBUG) {
1654                     System.out.println(i + " branchNum =" + branchNum
1655                             + " downGoing=" + (Boolean)downGoing.elementAt(i)
1656                             + "  isPWave=" + isPWave);
1657                 }
1658                 /*
1659                  * Save the turning depths for the ray parameter for both P and
1660                  * S waves. This way we get the depth correct for any rays that
1661                  * turn within a layer. We have to do this on a per branch basis
1662                  * because of converted phases, e.g. SKS.
1663                  */
1664                 try {
1665                     if(distRayParam > tMod.getTauBranch(branchNum, isPWave)
1666                             .getMaxRayParam()) {
1667                         turnDepth = tMod.getTauBranch(branchNum, isPWave)
1668                                 .getTopDepth();
1669                     } else if(distRayParam <= tMod.getTauBranch(branchNum,
1670                                                                 isPWave)
1671                             .getMinRayParam()) {
1672                         turnDepth = tMod.getTauBranch(branchNum, isPWave)
1673                                 .getBotDepth();
1674                     } else {
1675                         if(isPWave
1676                                 || tMod.sMod.depthInFluid((tMod.getTauBranch(branchNum,
1677                                                                              isPWave)
1678                                         .getTopDepth() + tMod.getTauBranch(branchNum,
1679                                                                            isPWave)
1680                                         .getBotDepth()) / 2.0)) {
1681                             turnDepth = tMod.sMod.findDepth(distRayParam,
1682                                                             tMod.getTauBranch(branchNum,
1683                                                                               isPWave)
1684                                                                     .getTopDepth(),
1685                                                             tMod.getTauBranch(branchNum,
1686                                                                               isPWave)
1687                                                                     .getBotDepth(),
1688                                                             PWAVE);
1689                         } else {
1690                             turnDepth = tMod.sMod.findDepth(distRayParam,
1691                                                             tMod.getTauBranch(branchNum,
1692                                                                               isPWave)
1693                                                                     .getTopDepth(),
1694                                                             tMod.getTauBranch(branchNum,
1695                                                                               isPWave)
1696                                                                     .getBotDepth(),
1697                                                             isPWave);
1698                         }
1699                     }
1700                 } catch(SlownessModelException e) {
1701                     // shouldn't happen but...
1702                     System.err.println("SeismicPhase.calcPierce: Caught SlownessModelException: "
1703                             + e.getMessage());
1704                     e.printStackTrace();
1705                     turnDepth = 0.0;
1706                 }
1707                 if(name.indexOf("Pdiff") != -1 || name.indexOf("Pn") != -1
1708                         || name.indexOf("Sdiff") != -1
1709                         || name.indexOf("Sn") != -1) {
1710                     /* head waves and diffracted waves are a special case. */
1711                     distA = tMod.getTauBranch(branchNum, isPWave)
1712                             .getDist(rayNum);
1713                     timeA = tMod.getTauBranch(branchNum, isPWave).time[rayNum];
1714                     distB = tMod.getTauBranch(branchNum, isPWave)
1715                             .getDist(rayNum);
1716                     timeB = tMod.getTauBranch(branchNum, isPWave).time[rayNum];
1717                 } else {
1718                     distA = tMod.getTauBranch(branchNum, isPWave)
1719                             .getDist(rayNum);
1720                     timeA = tMod.getTauBranch(branchNum, isPWave).time[rayNum];
1721                     distB = tMod.getTauBranch(branchNum, isPWave)
1722                             .getDist(rayNum + 1);
1723                     timeB = tMod.getTauBranch(branchNum, isPWave).time[rayNum + 1];
1724                 }
1725                 branchDist += distRatio * (distB - distA) + distA;
1726                 prevBranchTime = branchTime;
1727                 branchTime += distRatio * (timeB - timeA) + timeA;
1728                 if(((Boolean)downGoing.elementAt(i)).booleanValue()) {
1729                     branchDepth = Math.min(tMod.getTauBranch(branchNum, isPWave)
1730                                                    .getBotDepth(),
1731                                            turnDepth);
1732                 } else {
1733                     branchDepth = Math.min(tMod.getTauBranch(branchNum, isPWave)
1734                                                    .getTopDepth(),
1735                                            turnDepth);
1736                 }
1737                 // make sure ray actually propagates in this branch, leave
1738                 // a little room for numerical "chatter"
1739                 if(Math.abs(prevBranchTime - branchTime) > 1e-10) {
1740                     currArrival.pierce[numAdded++] = new TimeDist(distRayParam,
1741                                                                   branchTime,
1742                                                                   branchDist,
1743                                                                   branchDepth);
1744                     if(DEBUG) {
1745                         System.out.println(" branchTime=" + branchTime
1746                                 + " branchDist=" + branchDist + " branchDepth="
1747                                 + branchDepth);
1748                         System.out.println("incrementTime = "
1749                                 + (distRatio * (timeB - timeA)) + " timeB="
1750                                 + timeB + " timeA=" + timeA);
1751                     }
1752                 }
1753             }
1754             /*
1755              * Here we worry about the special case for head and diffracted
1756              * waves. It is assumed that a phase can be a diffracted wave or a
1757              * head wave, but not both. Nor can it be a head wave or diffracted
1758              * wave for both P and S.
1759              */
1760             if(name.indexOf("Pdiff") != -1 || name.indexOf("Pn") != -1
1761                     || name.indexOf("Sdiff") != -1 || name.indexOf("Sn") != -1) {
1762                 TimeDist[] diffTD;
1763                 int i = 0;
1764                 if(name.indexOf("Pdiff") != -1 || name.indexOf("Sdiff") != -1) {
1765                     diffTD = new TimeDist[numAdded + 1];
1766                     while(currArrival.pierce[i].depth != tMod.cmbDepth) {
1767                         diffTD[i] = currArrival.pierce[i];
1768                         i++;
1769                     }
1770                     diffTD[i] = currArrival.pierce[i];
1771                     double diffractDist = currArrival.dist - dist[0];
1772                     double diffractTime = diffractDist * diffTD[i].p;
1773                     diffTD[i + 1] = new TimeDist(diffTD[i].p,
1774                                                  diffTD[i].time + diffractTime,
1775                                                  diffTD[i].dist + diffractDist,
1776                                                  tMod.cmbDepth);
1777                     i++;
1778                     System.arraycopy(currArrival.pierce,
1779                                      i,
1780                                      diffTD,
1781                                      i + 1,
1782                                      numAdded - i);
1783                     for(int j = i + 1; j < diffTD.length; j++) {
1784                         diffTD[j].dist += diffractDist;
1785                         diffTD[j].time += diffractTime;
1786                     }
1787                     numAdded++;
1788                 } else {
1789                     // now deal with Pn and Sn case
1790                     int numFound = 0;
1791                     int indexInString = -1;
1792                     // can't have both Pn and Sn in a phase, so one of these
1793                     // should do nothing
1794                     while((indexInString = name.indexOf("Pn", indexInString + 1)) != -1) {
1795                         numFound++;
1796                     }
1797                     while((indexInString = name.indexOf("Sn", indexInString + 1)) != -1) {
1798                         numFound++;
1799                     }
1800                     double refractDist = currArrival.dist - dist[0];
1801                     diffTD = new TimeDist[numAdded + numFound];
1802                     int j = 0;
1803                     while(j < numFound) {
1804                         while(currArrival.pierce[i].depth != tMod.mohoDepth) {
1805                             diffTD[i + j] = currArrival.pierce[i];
1806                             diffTD[i + j].dist += j * refractDist / numFound;
1807                             i++;
1808                         }
1809                         diffTD[i + j] = currArrival.pierce[i];
1810                         diffTD[i + j].dist += j * refractDist / numFound;
1811                         diffTD[i + j + 1] = new TimeDist(diffTD[i].p,
1812                                                          diffTD[i + j].time
1813                                                                  + refractDist
1814                                                                  / numFound
1815                                                                  * diffTD[i].p,
1816                                                          diffTD[i + j].dist
1817                                                                  + refractDist
1818                                                                  / numFound,
1819                                                          tMod.mohoDepth);
1820                         i++;
1821                         j++;
1822                     }
1823                     System.arraycopy(currArrival.pierce,
1824                                      i,
1825                                      diffTD,
1826                                      i + j,
1827                                      numAdded - i);
1828                     for(int k = i + j; k < diffTD.length; k++) {
1829                         diffTD[k].dist += refractDist;
1830                     }
1831                     numAdded += j;
1832                 }
1833                 currArrival.pierce = diffTD;
1834             } else {
1835                 TimeDist[] temp = new TimeDist[numAdded];
1836                 System.arraycopy(currArrival.pierce, 0, temp, 0, numAdded);
1837                 currArrival.pierce = temp;
1838             }
1839         }
1840     }
1841 
1842     /*** calculates the path this phase takes through the earth model. */
1843     public void calcPath(TauModel tMod) {
1844         TimeDist[] tempTimeDist;
1845         Vector pathList = new Vector(10);
1846         int arraySize;
1847         Arrival currArrival;
1848         int rayNum = 0;
1849         int branchNum;
1850         boolean isPWave;
1851         for(int arrivalNum = 0; arrivalNum < arrivals.size(); arrivalNum++) {
1852             pathList.removeAllElements();
1853             currArrival = (Arrival)arrivals.elementAt(arrivalNum);
1854             /*
1855              * Find the ray parameter index that corresponds to the arrival ray
1856              * parameter in the TauModel, ie it is between rayNum and rayNum+1.
1857              */
1858             for(int i = 0; i < tMod.rayParams.length; i++) {
1859                 if(tMod.rayParams[i] >= currArrival.rayParam) {
1860                     rayNum = i;
1861                 }
1862             }
1863             tempTimeDist = new TimeDist[1];
1864             tempTimeDist[0] = new TimeDist(currArrival.rayParam,
1865                                            0.0,
1866                                            0.0,
1867                                            tMod.sourceDepth);
1868             pathList.addElement(tempTimeDist);
1869             try {
1870                 for(int i = 0; i < branchSeq.size(); i++) {
1871                     branchNum = ((Integer)branchSeq.elementAt(i)).intValue();
1872                     isPWave = ((Boolean)waveType.elementAt(i)).booleanValue();
1873                     if(DEBUG) {
1874                         System.out.println("i="
1875                                 + i
1876                                 + " branchNum="
1877                                 + branchNum
1878                                 + " isPWave="
1879                                 + isPWave
1880                                 + " downgoing="
1881                                 + ((Boolean)downGoing.elementAt(i)).booleanValue());
1882                     }
1883                     tempTimeDist = tMod.getTauBranch(branchNum, isPWave)
1884                             .path(currArrival.rayParam,
1885                                   ((Boolean)downGoing.elementAt(i)).booleanValue(),
1886                                   tMod.sMod);
1887                     if(tempTimeDist != null) {
1888                         pathList.addElement(tempTimeDist);
1889                     }
1890                     /*
1891                      * Here we worry about the special case for head and
1892                      * diffracted waves.
1893                      */
1894                     if(branchNum == tMod.cmbBranch - 1
1895                             && ((Integer)branchSeq.elementAt(i + 1)).intValue() == tMod.cmbBranch - 1
1896                             && (name.indexOf("Pdiff") != -1 || name.indexOf("Sdiff") != -1)) {
1897                         TimeDist[] diffTD = new TimeDist[1];
1898                         diffTD[0] = new TimeDist(currArrival.rayParam,
1899                                                  (currArrival.dist - dist[0])
1900                                                          * currArrival.rayParam,
1901                                                  currArrival.dist - dist[0],
1902                                                  tMod.cmbDepth);
1903                         pathList.addElement(diffTD);
1904                     } else if(branchNum == tMod.mohoBranch - 1
1905                             && ((Integer)branchSeq.elementAt(i + 1)).intValue() == tMod.mohoBranch - 1
1906                             && (name.indexOf("Pn") != -1 || name.indexOf("Sn") != -1)) {
1907                         int numFound = 0;
1908                         int indexInString = -1;
1909                         // can't have both Pn and Sn in a phase, so one of these
1910                         // should do nothing
1911                         while((indexInString = name.indexOf("Pn",
1912                                                             indexInString + 1)) != -1) {
1913                             numFound++;
1914                         }
1915                         while((indexInString = name.indexOf("Sn",
1916                                                             indexInString + 1)) != -1) {
1917                             numFound++;
1918                         }
1919                         TimeDist[] headTD = new TimeDist[1];
1920                         headTD[0] = new TimeDist(currArrival.rayParam,
1921                                                  (currArrival.dist - dist[0])
1922                                                          / numFound
1923                                                          * currArrival.rayParam,
1924                                                  (currArrival.dist - dist[0])
1925                                                          / numFound,
1926                                                  tMod.mohoDepth);
1927                         pathList.addElement(headTD);
1928                     }
1929                 }
1930                 arraySize = 0;
1931                 for(int i = 0; i < pathList.size(); i++) {
1932                     arraySize += ((TimeDist[])pathList.elementAt(i)).length;
1933                 }
1934                 currArrival.path = new TimeDist[arraySize];
1935                 TimeDist cummulative = new TimeDist(currArrival.rayParam,
1936                                                     0.0,
1937                                                     0.0,
1938                                                     currArrival.sourceDepth);
1939                 TimeDist[] branchPath;
1940                 int numAdded = 0;
1941                 for(int i = 0; i < pathList.size(); i++) {
1942                     branchPath = (TimeDist[])pathList.elementAt(i);
1943                     for(int j = 0; j < branchPath.length; j++) {
1944                         cummulative.add(branchPath[j]);
1945                         cummulative.depth = branchPath[j].depth;
1946                         currArrival.path[numAdded] = (TimeDist)cummulative.clone();
1947                         numAdded++;
1948                     }
1949                 }
1950             } catch(SlownessModelException e) {
1951                 // This shouldn't happen
1952                 System.out.println("Caught SlownessModelException in "
1953                         + "SeismicPhase.calcPath(): " + e.getMessage());
1954                 e.printStackTrace();
1955                 System.out.println("Skipping phase " + currArrival.name);
1956             }
1957         }
1958     }
1959 
1960     /***
1961      * Performs consistency checks on the previously tokenized phase name stored
1962      * in legs.
1963      */
1964     public boolean phaseValidate() {
1965         String currToken = (String)legs.elementAt(0);
1966         String prevToken;
1967         boolean prevIsReflect = false;
1968         /* Special cases for diffracted waves. */
1969         if(legs.size() == 2
1970                 && (currToken.equals("Pdiff") || currToken.equals("Sdiff") || currToken.endsWith("kmps"))
1971                 && ((String)legs.elementAt(1)).equals("END")) { return true; }
1972         /* Check first leg. */
1973         if(!(currToken.equals("Pg") || currToken.equals("Pb")
1974                 || currToken.equals("Pn") || currToken.equals("Pdiff")
1975                 || currToken.equals("Sg") || currToken.equals("Sb")
1976                 || currToken.equals("Sn") || currToken.equals("Sdiff")
1977                 || currToken.equals("P") || currToken.equals("S")
1978                 || currToken.equals("p") || currToken.equals("s"))) { return false; }
1979         for(int i = 1; i < legs.size(); i++) {
1980             prevToken = currToken;
1981             currToken = (String)legs.elementAt(i);
1982             /* Check for 2 reflections with no leg between them. */
1983             if(currToken.startsWith("^") || currToken.startsWith("v")
1984                     || currToken.equals("m") || currToken.equals("c")
1985                     || currToken.equals("i")) {
1986                 if(prevIsReflect) {
1987                     return false;
1988                 } else {
1989                     prevIsReflect = true;
1990                 }
1991             } else {
1992                 prevIsReflect = false;
1993             }
1994             /* Check for "END" before the end. */
1995             if(prevToken.equals("END")) { return false; }
1996             /* Check for P or S next to I or J */
1997             if((currToken.startsWith("P") || currToken.startsWith("S")
1998                     || currToken.startsWith("p") || currToken.startsWith("s")
1999                     || currToken.equals("m") || currToken.equals("c"))
2000                     && (prevToken.equals("I") || prevToken.equals("J") || prevToken.equals("i"))) { return false; }
2001             if((prevToken.startsWith("P") || prevToken.startsWith("S")
2002                     || prevToken.startsWith("p") || prevToken.startsWith("s")
2003                     || prevToken.equals("m") || prevToken.equals("c"))
2004                     && (currToken.equals("I") || currToken.equals("J") || currToken.equals("i"))) { return false; }
2005             /* Check for m next to K. */
2006             if(prevToken.equals("m") && currToken.equals("K")) { return false; }
2007             if(currToken.equals("m") && prevToken.equals("K")) { return false; }
2008         }
2009         /* Make sure legs end in "END". */
2010         if(!currToken.equals("END")) { return false; }
2011         return true;
2012     }
2013 
2014     public String toString() {
2015         String desc = name + ": ";
2016         for(int i = 0; i < legs.size(); i++) {
2017             desc += legs.elementAt(i) + " ";
2018         }
2019         desc += "\n";
2020         for(int i = 0; i < branchSeq.size(); i++) {
2021             desc += (Integer)branchSeq.elementAt(i) + " ";
2022         }
2023         desc += "\n";
2024         desc += "minRayParam=" + minRayParam + " maxRayParam=" + maxRayParam;
2025         desc += "\n";
2026         desc += "minDistance=" + (minDistance * 180.0 / Math.PI)
2027                 + " maxDistance=" + (maxDistance * 180.0 / Math.PI);
2028         return desc;
2029     }
2030 
2031     public void dump() {
2032         for(int j = 0; j < dist.length; j++) {
2033             System.out.println(j + "  " + dist[j] + "  " + rayParams[j]);
2034         }
2035     }
2036 
2037     /*** Returns an independent clone of the SeismicPhase. */
2038     public Object clone() {
2039         SeismicPhase newObject;
2040         try {
2041             newObject = (SeismicPhase)super.clone();
2042             newObject.branchSeq = new Vector(branchSeq.size());
2043             for(int i = 0; i < branchSeq.size(); i++) {
2044                 newObject.branchSeq.addElement(new Integer(((Integer)branchSeq.elementAt(i)).intValue()));
2045             }
2046             for(int i = 0; i < legs.size(); i++) {
2047                 newObject.legs.addElement(new Integer(((Integer)legs.elementAt(i)).intValue()));
2048             }
2049             newObject.dist = (double[])dist.clone();
2050             newObject.time = (double[])time.clone();
2051             newObject.rayParams = (double[])rayParams.clone();
2052             for(int i = 0; i < arrivals.size(); i++) {
2053                 newObject.arrivals.addElement((Arrival)((Arrival)arrivals.elementAt(i)).clone());
2054             }
2055             return newObject;
2056         } catch(CloneNotSupportedException e) {
2057             // Can't happen, but...
2058             System.err.println("Caught CloneNotSupportedException: "
2059                     + e.getMessage());
2060             throw new InternalError(e.toString());
2061         }
2062     }
2063 
2064     public static void main(String args[]) {
2065         TauModel tMod;
2066         TauModel tModDepth;
2067         try {
2068             if(args.length < 3) {
2069                 System.out.println("Usage: SeismicPhase modelfile depth phasename [phasename ...]");
2070             }
2071             tMod = TauModel.readModel(args[0]);
2072             tModDepth = tMod.depthCorrect(Double.valueOf(args[1]).doubleValue());
2073             int offset;
2074             for(int i = 2; i < args.length; i++) {
2075                 offset = 0;
2076                 System.out.println("-----");
2077                 SeismicPhase sp = new SeismicPhase(args[i], tModDepth);
2078                 sp.init();
2079                 System.out.println(sp);
2080                 sp.dump();
2081             }
2082             System.out.println("-----");
2083         } catch(FileNotFoundException e) {
2084             System.out.println(e.getMessage());
2085         } catch(OptionalDataException e) {
2086             System.out.println(e.getMessage());
2087         } catch(StreamCorruptedException e) {
2088             System.out.println(e.getMessage());
2089         } catch(IOException e) {
2090             System.out.println(e.getMessage());
2091         } catch(ClassNotFoundException e) {
2092             System.out.println(e.getMessage());
2093         } catch(TauModelException e) {
2094             System.out.println(e.getMessage());
2095             e.printStackTrace();
2096         }
2097     }
2098 }