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.DataInputStream;
20  import java.io.DataOutputStream;
21  import java.io.IOException;
22  import java.io.Serializable;
23  
24  /***
25   * provides storage and methods for distance, time and tau increments for a
26   * branch. A branch is a group of layers bounded by discontinuities or reversals
27   * in slowness gradient.
28   * 
29   * @version 1.1.3 Wed Jul 18 15:00:35 GMT 2001
30   * @author H. Philip Crotwell
31   */
32  public class TauBranch implements Serializable, Cloneable {
33  
34      /*** Turns on debugging output. */
35      transient public boolean DEBUG = false;
36  
37      /*** The type of wave for this branch, P or S. */
38      protected boolean isPWave;
39  
40      /*** The minimum depth of this layer. */
41      private double topDepth;
42  
43      /*** The maximum depth of this layer. */
44      private double botDepth;
45  
46      /***
47       * The maximum ray parameter that can penetrate into this branch. Time,
48       * distance, and tau increments are undefined (0.0) for ray parameters
49       * larger than this.
50       */
51      private double maxRayParam;
52  
53      /***
54       * The minimum ray parameter that is turned, but not reflected, in this
55       * branch. If no rays turn, then it is set equal to maxRayParam.
56       */
57      private double minTurnRayParam;
58  
59      /***
60       * The minimum ray parameter that is turned or critically reflected in this
61       * branch. If no rays turn, then it is set equal to maxRayParam.
62       */
63      private double minRayParam;
64  
65      /***
66       * Holds distance (radians) evaluated at the ith ray parameter for this
67       * branch.
68       */
69      protected double[] dist;
70  
71      /***
72       * Holds time (seconds) evaluated at the ith ray parameter for this branch.
73       */
74      protected double[] time;
75  
76      /*** Holds tau evaluated at the ith ray parameter for this branch. */
77      protected double[] tau;
78  
79      // Constructors --------------------------------------------------------
80      public TauBranch(double topDepth, double botDepth, boolean isPWave) {
81          this.topDepth = topDepth;
82          this.botDepth = botDepth;
83          this.isPWave = isPWave;
84      }
85  
86      // Methods -------------------------------------------------------------
87      //Accessor methods
88      /*** @returns the minimum (top) depth of this layer. */
89      public double getTopDepth() {
90          return topDepth;
91      }
92  
93      /*** @returns the maximum (bottom) depth of this layer. */
94      public double getBotDepth() {
95          return botDepth;
96      }
97  
98      /***
99       * @returns the maximum ray parameter that can penetrate into this branch.
100      *          Time, distance, and tau increments are undefined, set to 0.0,
101      *          for ray parameters larger than this.
102      */
103     public double getMaxRayParam() {
104         return maxRayParam;
105     }
106 
107     /***
108      * @returns the minimum ray parameter that is turned, but not reflected, in
109      *          this branch.
110      */
111     public double getMinTurnRayParam() {
112         return minTurnRayParam;
113     }
114 
115     /***
116      * @returns the minimum ray parameter that is turned or reflected in this
117      *          branch.
118      */
119     public double getMinRayParam() {
120         return minRayParam;
121     }
122 
123     /***
124      * @returns an array, cloned, containing distance (radians) evaluated at the
125      *          i_th ray parameter for this branch.
126      */
127     public double[] getDist() {
128         return (double[])dist.clone();
129     }
130 
131     /***
132      * @returns the distance (radians) evaluated at the
133      *          i_th ray parameter for this branch.
134      */
135     public double getDist(int i) {
136         return dist[i];
137     }
138 
139     /***
140      * @returns an array, cloned, containing time (seconds) evaluated at the
141      *          i_th ray parameter for this branch.
142      */
143     public double[] getTime() {
144         return (double[])time.clone();
145     }
146 
147     /***
148      * @returns the time (seconds) evaluated at the
149      *          i_th ray parameter for this branch.
150      */
151     public double getTime(int i) {
152         return time[i];
153     }
154 
155     /***
156      * @returns an array, cloned, containing tau (seconds) evaluated at the i_th
157      *          ray parameter for this branch.
158      */
159     public double[] getTau() {
160         return (double[])tau.clone();
161     }
162 
163     /***
164      * @returns tau (seconds) evaluated at the i_th
165      *          ray parameter for this branch.
166      */
167     public double getTau(int i) {
168         return tau[i];
169     }
170 
171     // normal methods
172     /***
173      * Calculates tau for this branch, between slowness layers topLayerNum and
174      * botLayerNum, inclusive.
175      * 
176      * @exception NoSuchLayerException
177      *                if a needed slowness layer cannot be found.
178      * @exception SlownessModelException
179      *                if there is a problem with the slowness model
180      * @exception TauModelException
181      *                if the slownessmodel and taumodel are not compatible
182      */
183     public void createBranch(SlownessModel sMod,
184                              double minPSoFar,
185                              double rayParams[]) throws NoSuchLayerException,
186             SlownessModelException, TauModelException {
187         int layerNum;
188         TimeDist timeDist;
189         SlownessLayer layer;
190         double p;
191         int topLayerNum = sMod.layerNumberBelow(getTopDepth(), isPWave);
192         int botLayerNum = sMod.layerNumberAbove(getBotDepth(), isPWave);
193         SlownessLayer topSLayer = sMod.getSlownessLayerClone(topLayerNum,
194                                                              isPWave);
195         SlownessLayer botSLayer = sMod.getSlownessLayerClone(botLayerNum,
196                                                              isPWave);
197         if(topSLayer.getTopDepth() != getTopDepth() || botSLayer.getBotDepth() != getBotDepth()) {
198             if(topSLayer.getTopDepth() != getTopDepth()
199                     && Math.abs(topSLayer.getTopDepth() - getTopDepth()) < 0.000001) {
200                 // really close, so just move top
201                 System.err.println("Changing topDepth: " + "\ntopDepth: "
202                         + getTopDepth() + " " + topSLayer.getTopDepth() + "\nbotDepth: "
203                         + getBotDepth() + " " + botSLayer.getBotDepth());
204                 topDepth = topSLayer.getTopDepth();
205             } else if(botSLayer.getBotDepth() != getBotDepth()
206                     && Math.abs(botSLayer.getBotDepth() - getBotDepth()) < 0.000001) {
207                 // really close, so just move bottom
208                 System.err.println("Changing botDepth: " + "\ntopDepth: "
209                         + getTopDepth() + " " + topSLayer.getTopDepth() + "\nbotDepth: "
210                         + getBotDepth() + " " + botSLayer.getBotDepth());
211                 botDepth = botSLayer.getBotDepth();
212             } else {
213                 // bad match, throw exception
214                 throw new TauModelException("createBranch: TauBranch not compatible with slowness sampling:"
215                         + "\ntopDepth: "
216                         + getTopDepth()
217                         + " "
218                         + topSLayer.getTopDepth()
219                         + "\nbotDepth: " + getBotDepth() + " " + botSLayer.getBotDepth());
220             }
221         }
222         /*
223          * Here we set minTurnRayParam to be the ray parameter that turns within
224          * the layer, not including total reflections off of the bottom.
225          * maxRayParam is the largest ray parameter that can penetrate this
226          * branch. minRayParam is the minimum ray parameter that turns or is
227          * totally reflected in this branch.
228          */
229         maxRayParam = minPSoFar;
230         minTurnRayParam = sMod.getMinTurnRayParam(getBotDepth(), isPWave);
231         minRayParam = sMod.getMinRayParam(getBotDepth(), isPWave);
232         tau = new double[rayParams.length];
233         dist = new double[rayParams.length];
234         time = new double[rayParams.length];
235         for(int rayNum = 0; rayNum < rayParams.length; rayNum++) {
236             p = rayParams[rayNum];
237             timeDist = calcTimeDist(sMod, topLayerNum, botLayerNum, p);
238             dist[rayNum] = timeDist.dist;
239             time[rayNum] = timeDist.time;
240             tau[rayNum] = time[rayNum] - p * dist[rayNum];
241             if(DEBUG && (rayNum % ((int)rayParams.length / 10) == 0)) {
242                 System.out.print(rayNum + ", ");
243             }
244         }
245     }
246 
247     /***
248      * calculates the time and distance increments for the given ray parameter.
249      * The topDepth and botDepth must be correct as they determine the bounds on
250      * the integration/summing.
251      * 
252      * @exception NoSuchLayerException
253      *                if topLayerNum or botLayerNum are not in the slowness
254      *                model.
255      * @exception SlownessModelException
256      *                if the ray with ray parameter p turns within a layer
257      *                instead of at the bottom.
258      */
259     public TimeDist calcTimeDist(SlownessModel sMod,
260                                  int topLayerNum,
261                                  int botLayerNum,
262                                  double p) throws NoSuchLayerException,
263             SlownessModelException {
264         int layerNum;
265         TimeDist timeDist = new TimeDist(p);
266         SlownessLayer layer;
267         if(p <= getMaxRayParam()) {
268             layerNum = topLayerNum;
269             layer = sMod.getSlownessLayer(layerNum, isPWave);
270             while(layerNum <= botLayerNum && p <= layer.getTopP() && p <= layer.getBotP()) {
271                 timeDist.add(sMod.layerTimeDist(p, layerNum, isPWave));
272                 layerNum++;
273                 if(layerNum <= botLayerNum) {
274                     layer = sMod.getSlownessLayerClone(layerNum, isPWave);
275                 }
276             }
277             if((layer.getTopP() - p) * (p - layer.getBotP()) > 0) { throw new SlownessModelException("Ray turns in the middle of this"
278                     + " layer. layerNum = "
279                     + layerNum
280                     + " sphericalRayParam "
281                     + p + " layer =" + layer); }
282         }
283         return timeDist;
284     }
285 
286     /***
287      * Inserts the distance, time, and tau increment for the slowness sample
288      * given to the branch. This is used for making the depth correction to a
289      * tau model for a non-surface source.
290      * 
291      * @throws TauModelException
292      *             if the tau branch is not compatable with the slowness
293      *             sampling
294      * @see edu.sc.seis.TauP.TauModel.depthCorrect(double)
295      */
296     protected void insert(double rayParam, SlownessModel sMod, int index)
297             throws NoSuchLayerException, SlownessModelException,
298             TauModelException {
299         int topLayerNum = sMod.layerNumberBelow(getTopDepth(), isPWave);
300         int botLayerNum = sMod.layerNumberAbove(getBotDepth(), isPWave);
301         SlownessLayer topSLayer = sMod.getSlownessLayerClone(topLayerNum,
302                                                              isPWave);
303         SlownessLayer botSLayer = sMod.getSlownessLayerClone(botLayerNum,
304                                                              isPWave);
305         if(topSLayer.getTopDepth() != getTopDepth() || botSLayer.getBotDepth() != getBotDepth()) { throw new TauModelException("insert: TauBranch depths not compatible with slowness sampling:"
306                 + "\ntopDepth: "
307                 + getTopDepth()
308                 + " "
309                 + topSLayer.getTopDepth()
310                 + "\nbotDepth: " + getBotDepth() + " " + botSLayer.getBotDepth()); }
311         TimeDist td = new TimeDist(rayParam, 0.0, 0.0);
312         TimeDist temptd;
313         if(topSLayer.getBotP() >= rayParam && topSLayer.getTopP() >= rayParam) {
314             for(int i = topLayerNum; i <= botLayerNum; i++) {
315                 if(sMod.getSlownessLayer(i, isPWave).getBotP() < rayParam) {
316                     // so we don't sum below the turning depth
317                     break;
318                 } else {
319                     temptd = sMod.layerTimeDist(rayParam, i, isPWave);
320                     td.dist += temptd.dist;
321                     td.time += temptd.time;
322                 }
323             }
324         }
325         shiftBranch(index);
326         dist[index] = td.dist;
327         time[index] = td.time;
328         tau[index] = td.time - rayParam * td.dist;
329     }
330 
331     /***
332      * generates a new tau branch by "subtracting" the given tau branch from
333      * this tau branch. The given tau branch is assumed to by the upper part of
334      * this branch.
335      * 
336      * @argument indexP specifies where a new ray coresponding to a P wave
337      *           sample has been added, it is -1 if no ray parameter has been
338      *           added to topBranch.
339      * @arguement indexS is similar to indexP except for a S wave sample. Note
340      *            that although the ray parameters for indexP and indexS were
341      *            for the P and S waves that turned at the source depth, both
342      *            ray parameters need to be added to both P and S branches.
343      */
344     protected TauBranch difference(TauBranch topBranch,
345                                    int indexP,
346                                    int indexS,
347                                    SlownessModel sMod,
348                                    double minPSoFar,
349                                    double rayParams[])
350             throws NoSuchLayerException, SlownessModelException,
351             TauModelException {
352         if(topBranch.getTopDepth() != getTopDepth() || topBranch.getBotDepth() > getBotDepth()) {
353             if(topBranch.getTopDepth() != getTopDepth()
354                     && Math.abs(topBranch.getTopDepth() - getTopDepth()) < 0.000001) {
355                 // really close, so just move top
356                 topDepth = topBranch.getTopDepth();
357             } else {
358                 // bad match, throw exception
359                 throw new TauModelException("createBranch: TauBranch not compatible with slowness sampling:"
360                         + "\ntopDepth: "
361                         + getTopDepth()
362                         + " "
363                         + topBranch.getTopDepth()
364                         + "\nbotDepth: " + getBotDepth() + " " + topBranch.getBotDepth());
365             }
366         }
367         if(topBranch.isPWave != isPWave) { throw new TauModelException("Can't difference branches: "
368                 + "topBranch.topDepth="
369                 + topBranch.getTopDepth()
370                 + " topDepth="
371                 + getTopDepth()
372                 + " topBranch.botDepth="
373                 + topBranch.getBotDepth()
374                 + " botDepth="
375                 + getBotDepth()
376                 + " waveTypes:"
377                 + topBranch.isPWave
378                 + " " + isPWave); }
379         // find the top and bottom slowness layers of bottom half
380         int topLayerNum = sMod.layerNumberBelow(topBranch.getBotDepth(), isPWave);
381         int botLayerNum = sMod.layerNumberAbove(getBotDepth(), isPWave);
382         SlownessLayer topSLayer = sMod.getSlownessLayer(topLayerNum, isPWave);
383         SlownessLayer botSLayer = sMod.getSlownessLayer(botLayerNum, isPWave);
384         if(topSLayer.getTopDepth() != topBranch.getBotDepth()
385                 || botSLayer.getBotDepth() != getBotDepth()) { throw new TauModelException("difference: TauBranch not compatible with slowness sampling:"
386                 + "\ntopDepth: "
387                 + topBranch.getBotDepth()
388                 + " "
389                 + topSLayer.getTopDepth()
390                 + "\nbotDepth: "
391                 + getBotDepth()
392                 + " "
393                 + botSLayer.getBotDepth() + "\n" + topSLayer + "\n" + botSLayer); }
394         // make sure indexP and indexS really correspond to
395         // new ray parameters at the top of this branch
396         SlownessLayer sLayer = sMod.getSlownessLayer(sMod.layerNumberBelow(topBranch.getBotDepth(),
397                                                                            true),
398                                                      true);
399         if(indexP >= 0 && sLayer.getTopP() != rayParams[indexP]) { throw new TauModelException("difference: P wave index doesn't match top layer.\n "
400                 + sMod.layerNumberBelow(topBranch.getBotDepth(), true)
401                 + "\n"
402                 + rayParams[indexP - 1]
403                 + "\n"
404                 + rayParams[indexP]
405                 + "\n"
406                 + rayParams[indexP + 1]
407                 + "\nP="
408                 + sLayer
409                 + "\nS="
410                 + sMod.getSlownessLayer(sMod.layerNumberBelow(topBranch.getBotDepth(),
411                                                               false),
412                                         false)); }
413         sLayer = sMod.getSlownessLayer(sMod.layerNumberBelow(topBranch.getBotDepth(),
414                                                              false),
415                                        false);
416         if(indexS >= 0 && sLayer.getTopP() != rayParams[indexS]) { throw new TauModelException("difference: S wave index doesn't match top layer. "
417                 + sMod.layerNumberBelow(topBranch.getBotDepth(), false)
418                 + " "
419                 + rayParams[indexS - 1]
420                 + " "
421                 + rayParams[indexS]
422                 + " "
423                 + rayParams[indexS + 1] + "\n" + sLayer); }
424         sLayer = null;
425         // construct the new TauBranch, going from the bottom of
426         // the top half to the bottom of the whole branch
427         TauBranch botBranch = new TauBranch(topBranch.getBotDepth(),
428                                             getBotDepth(),
429                                             isPWave);
430         botBranch.maxRayParam = topBranch.getMinRayParam();
431         botBranch.minTurnRayParam = getMinTurnRayParam();
432         botBranch.minRayParam = getMinRayParam();
433         double PRayParam = -1.0, SRayParam = -1.0;
434         TimeDist timeDistP = new TimeDist();
435         TimeDist timeDistS = new TimeDist();
436         int arrayLength = dist.length;
437         if(indexP != -1) {
438             arrayLength++;
439             PRayParam = rayParams[indexP];
440             timeDistP = botBranch.calcTimeDist(sMod,
441                                                topLayerNum,
442                                                botLayerNum,
443                                                PRayParam);
444         }
445         if(indexS != -1 && indexS != indexP) {
446             arrayLength++;
447             SRayParam = rayParams[indexS];
448             timeDistS = botBranch.calcTimeDist(sMod,
449                                                topLayerNum,
450                                                botLayerNum,
451                                                SRayParam);
452         } else {
453             // in case indexS==indexP then we only need one
454             indexS = -1;
455         }
456         // allocate enough space
457         botBranch.dist = new double[arrayLength];
458         botBranch.time = new double[arrayLength];
459         botBranch.tau = new double[arrayLength];
460         if(indexP == -1) {
461             // then both are -1 so no new ray parameters added
462             for(int i = 0; i < dist.length; i++) {
463                 botBranch.dist[i] = dist[i] - topBranch.dist[i];
464                 botBranch.time[i] = time[i] - topBranch.time[i];
465                 botBranch.tau[i] = tau[i] - topBranch.tau[i];
466             }
467         } else {
468             if(indexS == -1) {
469                 // only indexP != 0
470                 for(int i = 0; i < indexP; i++) {
471                     botBranch.dist[i] = dist[i] - topBranch.dist[i];
472                     botBranch.time[i] = time[i] - topBranch.time[i];
473                     botBranch.tau[i] = tau[i] - topBranch.tau[i];
474                 }
475                 botBranch.dist[indexP] = timeDistP.dist;
476                 botBranch.time[indexP] = timeDistP.time;
477                 botBranch.tau[indexP] = timeDistP.time - PRayParam
478                         * timeDistP.dist;
479                 for(int i = indexP; i < dist.length; i++) {
480                     botBranch.dist[i + 1] = dist[i] - topBranch.dist[i + 1];
481                     botBranch.time[i + 1] = time[i] - topBranch.time[i + 1];
482                     botBranch.tau[i + 1] = tau[i] - topBranch.tau[i + 1];
483                 }
484             } else {
485                 // both indexP and indexS != -1 so we have two new samples
486                 for(int i = 0; i < indexS; i++) {
487                     botBranch.dist[i] = dist[i] - topBranch.dist[i];
488                     botBranch.time[i] = time[i] - topBranch.time[i];
489                     botBranch.tau[i] = tau[i] - topBranch.tau[i];
490                 }
491                 botBranch.dist[indexS] = timeDistS.dist;
492                 botBranch.time[indexS] = timeDistS.time;
493                 botBranch.tau[indexS] = timeDistS.time - SRayParam
494                         * timeDistS.dist;
495                 for(int i = indexS; i < indexP; i++) {
496                     botBranch.dist[i + 1] = dist[i] - topBranch.dist[i + 1];
497                     botBranch.time[i + 1] = time[i] - topBranch.time[i + 1];
498                     botBranch.tau[i + 1] = tau[i] - topBranch.tau[i + 1];
499                 }
500                 // put into indexP+1 as we have already shifted by 1
501                 // due to indexS
502                 botBranch.dist[indexP + 1] = timeDistP.dist;
503                 botBranch.time[indexP + 1] = timeDistP.time;
504                 botBranch.tau[indexP + 1] = timeDistP.time - PRayParam
505                         * timeDistP.dist;
506                 for(int i = indexP; i < dist.length; i++) {
507                     botBranch.dist[i + 2] = dist[i] - topBranch.dist[i + 2];
508                     botBranch.time[i + 2] = time[i] - topBranch.time[i + 2];
509                     botBranch.tau[i + 2] = tau[i] - topBranch.tau[i + 2];
510                 }
511             }
512         }
513         return botBranch;
514     }
515 
516     public void shiftBranch(int index) {
517         double[] newDist = new double[dist.length + 1];
518         System.arraycopy(dist, 0, newDist, 0, index);
519         newDist[index] = 0.0;
520         System.arraycopy(dist, index, newDist, index + 1, dist.length - index);
521         dist = newDist;
522         double[] newTime = new double[time.length + 1];
523         System.arraycopy(time, 0, newTime, 0, index);
524         newTime[index] = 0.0;
525         System.arraycopy(time, index, newTime, index + 1, time.length - index);
526         time = newTime;
527         double[] newTau = new double[tau.length + 1];
528         System.arraycopy(tau, 0, newTau, 0, index);
529         newTau[index] = 0.0;
530         System.arraycopy(tau, index, newTau, index + 1, tau.length - index);
531         tau = newTau;
532     }
533 
534     public TimeDist[] path(double rayParam,
535                            boolean downgoing,
536                            SlownessModel sMod) throws SlownessModelException {
537         if(rayParam > getMaxRayParam()) { return null; }
538         Assert.isTrue(rayParam >= 0.0, "ray parameter must not be negative.");
539         int topLayerNum;
540         int botLayerNum;
541         try {
542             topLayerNum = sMod.layerNumberBelow(getTopDepth(), isPWave);
543             botLayerNum = sMod.layerNumberAbove(getBotDepth(), isPWave);
544         } catch(NoSuchLayerException e) {
545             throw new SlownessModelException("Caught NoSuchLayerException. This likely means the"
546                     + "SlownessModel and TauModel are out of sync. "
547                     + e.getMessage());
548         }
549         TimeDist[] thePath = new TimeDist[botLayerNum - topLayerNum + 1];
550         int sLayerNum;
551         int pathIndex = 0;
552         double turnDepth;
553         SlownessLayer sLayer, turnSLayer;
554         /*** check to make sure layers and branches are compatable. */
555         sLayer = sMod.getSlownessLayer(topLayerNum, isPWave);
556         if(sLayer.getTopDepth() != getTopDepth()) { throw new SlownessModelException("Branch and Slowness model are not compatible! "
557                 + sLayer.getTopDepth() + " != " + getTopDepth() + "=topDepth"); }
558         sLayer = sMod.getSlownessLayer(botLayerNum, isPWave);
559         if(sLayer.getBotDepth() != getBotDepth()) { throw new SlownessModelException("Branch and Slowness model are not compatible! "
560                 + sLayer.getBotDepth() + " != " + getBotDepth() + "=botDepth"); }
561         if(downgoing) {
562             sLayerNum = topLayerNum;
563             sLayer = sMod.getSlownessLayer(sLayerNum, isPWave);
564             while(sLayer.getBotP() >= rayParam && sLayerNum <= botLayerNum) {
565                 if(!sLayer.isZeroThickness()) {
566                     thePath[pathIndex] = sMod.layerTimeDist(rayParam,
567                                                             sLayerNum,
568                                                             isPWave);
569                     thePath[pathIndex].depth = sLayer.getBotDepth();
570                     pathIndex++;
571                 }
572                 sLayerNum++;
573                 if(sLayerNum <= botLayerNum) {
574                     sLayer = sMod.getSlownessLayer(sLayerNum, isPWave);
575                 }
576             }
577             if(sLayerNum <= botLayerNum && !sLayer.isZeroThickness()) {
578                 turnDepth = sLayer.bullenDepthFor(rayParam,
579                                                   sMod.getRadiusOfEarth());
580                 turnSLayer = new SlownessLayer(sLayer.getTopP(),
581                                                sLayer.getTopDepth(),
582                                                rayParam,
583                                                turnDepth);
584                 thePath[pathIndex] = turnSLayer.bullenRadialSlowness(rayParam,
585                                                                      sMod.getRadiusOfEarth());
586                 thePath[pathIndex].depth = turnSLayer.getBotDepth();
587                 pathIndex++;
588             }
589         } else {
590             // upgoing
591             sLayerNum = botLayerNum;
592             sLayer = sMod.getSlownessLayer(sLayerNum, isPWave);
593             while(sLayer.getTopP() <= rayParam || sLayer.isZeroThickness()) {
594                 sLayerNum--;
595                 sLayer = sMod.getSlownessLayer(sLayerNum, isPWave);
596             }
597             if(sLayer.getBotP() < rayParam) {
598                 turnDepth = sLayer.bullenDepthFor(rayParam,
599                                                   sMod.getRadiusOfEarth());
600                 turnSLayer = new SlownessLayer(sLayer.getTopP(),
601                                                sLayer.getTopDepth(),
602                                                rayParam,
603                                                turnDepth);
604                 thePath[pathIndex] = turnSLayer.bullenRadialSlowness(rayParam,
605                                                                      sMod.getRadiusOfEarth());
606                 thePath[pathIndex].depth = turnSLayer.getTopDepth();
607                 pathIndex++;
608                 sLayerNum--;
609                 if(sLayerNum >= topLayerNum) {
610                     sLayer = sMod.getSlownessLayer(sLayerNum, isPWave);
611                 }
612             }
613             while(sLayerNum >= topLayerNum) {
614                 if(!sLayer.isZeroThickness()) {
615                     thePath[pathIndex] = sMod.layerTimeDist(rayParam,
616                                                             sLayerNum,
617                                                             isPWave);
618                     thePath[pathIndex].depth = sLayer.getTopDepth();
619                     pathIndex++;
620                 }
621                 sLayerNum--;
622                 if(sLayerNum >= topLayerNum) {
623                     sLayer = sMod.getSlownessLayer(sLayerNum, isPWave);
624                 }
625             }
626         }
627         TimeDist[] tempPath = new TimeDist[pathIndex];
628         System.arraycopy(thePath, 0, tempPath, 0, pathIndex);
629         return tempPath;
630     }
631 
632     public void writeToStream(DataOutputStream dos) throws IOException {
633         dos.writeInt(getClass().getName().length());
634         dos.writeBytes(getClass().getName());
635         dos.writeDouble(getTopDepth());
636         dos.writeDouble(getBotDepth());
637         dos.writeDouble(getMaxRayParam());
638         dos.writeDouble(getMinRayParam());
639         dos.writeDouble(getMinTurnRayParam());
640         dos.writeInt(dist.length);
641         for(int i = 0; i < dist.length; i++) {
642             dos.writeDouble(dist[i]);
643         }
644         dos.writeInt(time.length);
645         for(int i = 0; i < time.length; i++) {
646             dos.writeDouble(time[i]);
647         }
648         dos.writeInt(tau.length);
649         for(int i = 0; i < tau.length; i++) {
650             dos.writeDouble(tau[i]);
651         }
652     }
653 
654     public static TauBranch readFromStream(DataInputStream dis)
655             throws IOException, ClassNotFoundException, InstantiationException,
656             IllegalAccessException {
657         int length;
658         byte[] classString = new byte[dis.readInt()];
659         dis.read(classString);
660         Class tBranchClass = Class.forName(new String(classString));
661         TauBranch tBranch = (TauBranch)tBranchClass.newInstance();
662         tBranch.topDepth = dis.readDouble();
663         tBranch.botDepth = dis.readDouble();
664         tBranch.maxRayParam = dis.readDouble();
665         tBranch.minRayParam = dis.readDouble();
666         tBranch.minTurnRayParam = dis.readDouble();
667         length = dis.readInt();
668         tBranch.dist = new double[length];
669         for(int i = 0; i < tBranch.dist.length; i++) {
670             tBranch.dist[i] = dis.readDouble();
671         }
672         length = dis.readInt();
673         tBranch.time = new double[length];
674         for(int i = 0; i < tBranch.time.length; i++) {
675             tBranch.time[i] = dis.readDouble();
676         }
677         length = dis.readInt();
678         tBranch.tau = new double[length];
679         for(int i = 0; i < tBranch.tau.length; i++) {
680             tBranch.tau[i] = dis.readDouble();
681         }
682         return tBranch;
683     }
684 
685     /***
686      * Returns a clone of this TauBranch object. Note that super.clone() handles
687      * all normal variables while the arrays need to be cloned separately to
688      * generate a new array as opposed to a new reference to the old array.
689      * 
690      * @see Cloneable
691      */
692     public Object clone() {
693         TauBranch newObject;
694         try {
695             newObject = (TauBranch)super.clone();
696             newObject.dist = (double[])dist.clone();
697             newObject.time = (double[])time.clone();
698             newObject.tau = (double[])tau.clone();
699             return newObject;
700         } catch(CloneNotSupportedException e) {
701             // Can't happen, but...
702             System.err.println("Caught CloneNotSupportedException: "
703                     + e.getMessage());
704             throw new InternalError(e.toString());
705         }
706     }
707 
708     public String toString() {
709         String desc = "Tau Branch\n";
710         desc += " topDepth = " + getTopDepth() + "\n";
711         desc += " botDepth = " + getBotDepth() + "\n";
712         desc += " maxRayParam=" + getMaxRayParam() + " minTurnRayParam="
713                 + getMinTurnRayParam();
714         desc += " minRayParam=" + getMinRayParam() + "\n";
715         /*
716          * for (int i=0;i <tau.length;i++) { desc += "i = "+i+" time="+time[i]+"
717          * dist="+dist[i]+" tau="+tau[i]+"\n"; }
718          */
719         return desc;
720     }
721 }