View Javadoc

1   /*
2     The TauP Toolkit: Flexible Seismic Travel-Time and Raypath Utilities.
3     Copyright (C) 1998-2000 University of South Carolina
4   
5     This program is free software; you can redistribute it and/or
6     modify it under the terms of the GNU General Public License
7     as published by the Free Software Foundation; either version 2
8     of the License, or (at your option) any later version.
9   
10    This program is distributed in the hope that it will be useful,
11    but WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13    GNU General Public License for more details.
14  
15    You should have received a copy of the GNU General Public License
16    along with this program; if not, write to the Free Software
17    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
18  
19    The current version can be found at 
20    <A HREF="www.seis.sc.edu">http://www.seis.sc.edu</A>
21  
22    Bug reports and comments should be directed to 
23    H. Philip Crotwell, crotwell@seis.sc.edu or
24    Tom Owens, owens@seis.sc.edu
25  
26  */
27  
28  package edu.sc.seis.TauP;
29  
30  import java.io.File;
31  import java.io.IOException;
32  import java.io.Serializable;
33  
34  /***
35    * This class provides storage and methods for generating slowness-depth
36    * pairs in a spherical earth model.
37    *
38    * @version 1.1.3 Wed Jul 18 15:00:35 GMT 2001
39  
40  
41  
42    * @author H. Philip Crotwell
43    *
44    */
45  public class SphericalSModel extends SlownessModel
46        implements Serializable, Cloneable {
47  
48  // METHODS ----------------------------------------------------------------
49  
50        /*** Just for debugging purposes. */
51     public static void main(String[] args) {
52  	   System.out.println("Starting main");
53  		
54        VelocityModel vMod = new VelocityModel();
55        SphericalSModel sMod = new SphericalSModel();
56        String modelFilename;
57        if (args.length == 1) {
58           modelFilename = args[0];
59        } else {
60        	vMod.setFileType("tvel");
61           modelFilename = File.separator+"MacintoshHD"+File.separator+"Philip"+File.separator+"TauP"+File.separator+"VModels"+File.separator+"iasp91.tvel";
62        }
63        boolean DEBUG = true;
64  
65        try {
66           DEBUG = true;
67           vMod.readVelocityFile(modelFilename);
68           System.out.println("Done reading.");
69           if (DEBUG) System.out.println(vMod);
70           DEBUG = true;
71   
72           sMod.DEBUG = true;
73  
74           sMod.createSample(vMod);
75  
76           if (sMod.DEBUG) System.out.println(sMod);
77           sMod.validate();
78        } catch (IOException e) {
79           System.out.println("Tried to read!\n Caught IOException "
80                              + e.getMessage()+"\n"+e.getClass().getName());
81  					e.printStackTrace();
82        } catch (VelocityModelException e) {
83           System.out.println("Tried to read!\n Caught VelocityModelException "
84                              + e.getMessage());
85  					e.printStackTrace();
86        } catch (SlownessModelException e) {
87           System.out.println(vMod);
88           System.out.println(sMod);
89           System.out.println("Tried to create slowness!\n Caught SlownessModelException "
90                              + e.getMessage());
91  					e.printStackTrace();
92        } finally {
93           System.out.println("Done!\n");
94        }
95  
96     }
97     
98        /*** Returns the slowness for a velocity at a depth. 
99  		 *  @exception SlownessModelException if velocity is zero. */
100    public double toSlowness(double velocity, double depth) throws SlownessModelException {
101 		if (velocity == 0.0) {
102 			throw new SlownessModelException("Divide by zero in toSlowness()"+
103 				"\ndepth = "+depth+"\nThis likely has to do with using S velocities in the outer core");
104 		}
105       return (radiusOfEarth-depth)/velocity;
106    }
107    
108       /*** Returns the velocity for a slowness at a depth. 
109 		 *  @exception SlownessModelException if slowness is zero. */
110    public double toVelocity(double slowness, double depth) throws SlownessModelException {
111 		if (slowness == 0.0) {
112 			throw new SlownessModelException("Divide by zero in toVelocity()"+
113 				"\ndepth = "+depth+"\nPossibly this is due to depth at center of the earth?");
114 		}
115       return (radiusOfEarth-depth)/slowness;
116    }
117 
118       /*** Converts a velocity layer into a slowness layer. 
119        *  @exception SlownessModelException if velocity layer is malformed. */
120    public SlownessLayer toSlownessLayer(VelocityLayer vLayer, boolean isPWave) 
121 	throws SlownessModelException {
122        return new SlownessLayer(vLayer, true, radiusOfEarth, isPWave);
123    }
124 
125       /*** Returns the depth for a slowness given a velocity gradient. 
126 		 *  @exception SlownessModelException if the velocity gradient
127 		 *     exactly balances the spherical decrease in slowness. */
128    public double interpolate(double p, double topVelocity, double topDepth,
129          double slope) throws SlownessModelException {
130       double depth;
131       double denominator = p * slope + 1.0;
132 
133       if (denominator != 0.0) {
134          depth = (radiusOfEarth +
135             p* (topDepth*slope-topVelocity))/ denominator;
136       } else {
137             /* Uh oh, this is a neg velocity gradient that just
138              * balances the slowness gradient of the spherical
139              * slowness. In this case we should equally space the
140              * depths. ????
141              * This probably won't happen, but...
142              */ 
143          depth = Double.MAX_VALUE;
144          throw new SlownessModelException("Neg velocity gradient "+
145             "just balances the earth flattening!"+
146             " What should I do?!?!?!? topDepth= "+topDepth);
147       }
148       return depth;
149    }
150 
151       /*** Calculates the time and distance increments accumulated by a
152        *  ray of spherical ray parameter p when passing through layer layerNum.
153        *  for the easy cases of zero ray parameter, the center of the earth,
154        *  and constant velocity layers.
155        *  Note that this gives 1/2 of the true range and time increments since
156        *  there will be both an up going and a downgoing path.
157        *
158        *  @exception SlownessModelException occurs if the ray with the given
159        *     spherical ray parameter cannot propagate within this layer, or
160        *     if the ray turns within this layer but not at the bottom.
161        */
162    public TimeDist layerTimeDist(double sphericalRayParam, int layerNum, 
163    boolean isPWave) throws SlownessModelException {
164 
165       double swapDouble;
166       double b;        // temporary variable makes the calculations less ugly.
167  
168          // To hold the return values.
169       TimeDist timedist = new TimeDist(sphericalRayParam);
170 
171       SlownessLayer sphericalLayer = getSlownessLayer(layerNum, isPWave);
172       double topRadius = radiusOfEarth-sphericalLayer.getTopDepth(); // radius to top
173       double botRadius = radiusOfEarth-sphericalLayer.getBotDepth(); // radius to bot
174 
175          /* First we make sure that a ray with this ray parameter can propagate
176           * within this layer and doesn't turn in the middle of the layer.
177           * If not, then throw an exception. */
178       if (sphericalRayParam > 
179           Math.max(sphericalLayer.getTopP(),sphericalLayer.getBotP())) {
180          throw new SlownessModelException("Ray cannot propagate within this"+
181             " layer. layerNum = "+ layerNum+
182             " sphericalRayParam="+sphericalRayParam+"\n"+sphericalLayer);
183       }
184       if (sphericalRayParam < 0.0) {
185          throw new SlownessModelException("Ray Parameter is negative!!! "+
186             sphericalRayParam);
187       } 
188       if (sphericalRayParam > 
189             Math.min(sphericalLayer.getTopP(),sphericalLayer.getBotP())) {
190          if (DEBUG) {
191             System.out.println("Ray Turns in layer, velocities: "+
192                topRadius/sphericalRayParam+" "+
193                topRadius/sphericalLayer.getTopP()+" "+
194                botRadius/sphericalLayer.getBotP());
195             System.out.println("depths        top "+sphericalLayer.getTopDepth()+
196                "  bot "+sphericalLayer.getBotDepth());
197          }
198          throw new SlownessModelException("Ray turns in the middle of this"+
199             " layer. \nlayerNum = "+ layerNum+
200             " sphericalRayParam "+sphericalRayParam+
201             " sphericalLayer =  "+sphericalLayer+"\n");
202       }
203 
204          /* Check to see if this layer has zero thickness, if so then it is
205           * from a critically reflected slowness sample. So we should just
206           * return 0.0 for time and distance increments. */
207       if (sphericalLayer.getTopDepth() == sphericalLayer.getBotDepth()) {
208          timedist.time = 0.0;
209          timedist.dist = 0.0;
210          return timedist;
211       }
212  
213          /* Check to see if this layer contains the center of the earth. 
214           * If so then
215           * the spherical ray parameter should be 0.0 and we calculate 
216           * the range and
217           * time increments using a constant velocity layer (sphere). 
218           * See eq 43 and 44
219           * of Buland and Chapman, although we implement them slightly 
220           * differently. 
221           * Note that the distance and time increments are 
222           * for just downgoing or just up going, ie top of the layer 
223           * to the center of
224           * the earth or vice versa but not both. This is in keeping with 
225           * the convention
226           * that these are one way distance and time increments. We will
227           * multiply the result by 2 at the end, or if we are doing a 
228           * 1.5D model, the
229           * other direction may be different.
230           * The time increment for a ray of zero ray parameter passing
231           * half way through a sphere of constant velocity is just the spherical
232           * slowness at the top of the sphere. An amazingly simple result!
233           */
234       if (sphericalRayParam == 0.0 && sphericalLayer.getBotDepth() == radiusOfEarth){
235          if (layerNum != getNumLayers(isPWave)-1) throw new SlownessModelException(
236              "There are layers deeper than the center of the earth!");
237              
238          timedist.dist = Math.PI/2.0;
239          timedist.time = sphericalLayer.getTopP();
240         
241          if (DEBUG) {
242             System.out.println("Center of Earth: dist "+timedist.dist+
243                " time "+timedist.time);
244          }
245          if (timedist.dist<0.0 || timedist.time<0.0 ||
246              Double.isNaN(timedist.time) || Double.isNaN(timedist.dist)) {
247             throw new SlownessModelException("CoE timedist <0.0 or NaN: "+
248                "sphericalRayParam= "+sphericalRayParam+
249                " botDepth = "+sphericalLayer.getBotDepth()+
250                " dist="+timedist.dist+
251                " time="+timedist.time);
252          }
253          return timedist;
254       }
255        
256          /* Now we check to see if this is a constant velocity layer and
257           * if so than we can do a simple triangle calculation to get 
258           * the range and time increments. 
259           * To get the time increment we first calculate the path length
260           * through the layer using law of cosines, noting that the angle
261           * at the top of the layer can be obtained from the spherical Snell's
262           * Law. The time increment is just the path length divided by 
263           * the velocity. To get the distance we first find the angular 
264           * distance traveled, using law of sines.
265           */
266       if (Math.abs(topRadius/sphericalLayer.getTopP() -
267                    botRadius/sphericalLayer.getBotP()) < slownessTolerance ) {
268 
269             // temp variables
270          double vel = botRadius / sphericalLayer.getBotP();              // velocity
271             /* In cases of a ray turning at the bottom of the layer numerical
272              * roundoff can cause botTerm to be very small (1e-9) but negative
273              * which causes the sqrt to return NaN. We check for values that
274              * are within the numerical chatter of zero and just set them to
275              * zero. */
276          double topTerm, botTerm;
277          topTerm = topRadius*topRadius-
278             sphericalRayParam*sphericalRayParam*vel*vel;
279          if (Math.abs(topTerm) < slownessTolerance) {topTerm = 0.0;}
280 
281          if (sphericalRayParam == sphericalLayer.getBotP()) {
282                /* In this case the ray turns at the bottom of this layer
283                 * so sphericalRayParam*vel == botRadius and botTerm should
284                 * be zero. We check for this case specially because numerical
285                 * chatter can cause small round offs that lead to botTerm being
286                 * negative, causing a sqrt(-1) error.  */
287             botTerm = 0.0;
288          } else {
289             botTerm = botRadius*botRadius-
290                sphericalRayParam*sphericalRayParam*vel*vel;
291          }
292          
293             // Use b for temp storage of the length of the ray path.
294          b = Math.sqrt(topTerm) - Math.sqrt(botTerm);
295 
296          timedist.time = b/vel;
297 
298          timedist.dist=Math.asin(b*sphericalRayParam*vel/(topRadius*botRadius));
299 
300          if (timedist.dist<0.0 || timedist.time<0.0 ||
301              Double.isNaN(timedist.time) || Double.isNaN(timedist.dist)) {
302             throw new SlownessModelException("CVL timedist <0.0 or NaN: "+
303                "\nsphericalRayParam= "+sphericalRayParam+
304                "\n botDepth = "+sphericalLayer.getBotDepth()+
305                "\n topDepth = "+sphericalLayer.getTopDepth()+
306                "\n topRadius="+topRadius+" botRadius="+botRadius+
307                "\n dist="+timedist.dist+
308                "\n time="+timedist.time+
309                "\n b="+b+
310                "\n topTerm="+topTerm+
311                "\n botTerm="+botTerm+
312                "\n vel    ="+vel+"\n"+
313                "\n bR^2   ="+(botRadius*botRadius)+
314                "\n p^2v^2 ="+sphericalRayParam*sphericalRayParam*vel*vel+
315                "\n tR^2   ="+(topRadius*topRadius)+
316                "\n p^2v^2 ="+sphericalRayParam*sphericalRayParam*vel*vel);
317 
318          }
319          return timedist;
320       }
321 
322          /* OK, the layer is not a constant velocity layer or 
323           * the center of the earth and p is not zero
324           * so we have to do it the hard way...
325           *
326           */
327       return sphericalLayer.bullenRadialSlowness(sphericalRayParam, 
328          radiusOfEarth);
329    }
330     
331       /*** Performs consistency check on the velocity model.
332        *  @return true if successful, throws SlownessModelException otherwise.
333        *  @exception SlownessModelException if any check fails
334        */
335    public boolean validate() 
336       throws SlownessModelException 
337    {
338       boolean isOK = super.validate();
339    
340       double prevDepth= 0.0;
341       DepthRange highSZoneDepth, fluidZone;
342       SlownessLayer sLayer;
343 
344       boolean isPWave = true;
345       for (int j=0;j<2; j++, isPWave = false) {
346          for (int i=0;i<getNumLayers(isPWave);i++) {
347             sLayer = getSlownessLayer(i,isPWave);
348             prevDepth = sLayer.getBotDepth();
349          
350          /* No slowness layer should have a depth greater than radiusOfEarth.*/
351             if (prevDepth > radiusOfEarth) {
352                isOK = false;
353                throw new SlownessModelException(
354                   "Slowness layer has a depth larger than the radius of "+
355                   "the earth in a spherical model. max depth = "+prevDepth+
356                   " radiusOfEarth = "+radiusOfEarth);
357             } else {
358                isOK |= true;
359             }
360          }
361       }
362       
363 
364          /* Everything checks out OK so return true. */
365       return isOK;
366    }
367 
368       /* Returns a clone of this slowness model. All fields are correctly
369        * copied so modifications to the clone do not affect the original. */
370    public Object clone() {
371       SphericalSModel newObject = (SphericalSModel)super.clone();
372 
373       return newObject;
374    }
375 
376    public String toString() {
377       int topCriticalLayerNum;
378       int botCriticalLayerNum;
379 
380       String desc = "spherical model:\n"+super.toString();
381 
382       return desc;
383    }
384 
385 }
386