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   * Theta.java
29   *
30   *
31   * Created: Mon Feb 15 13:48:06 1999
32   *
33   * @author Philip Crotwell
34   * @version 1.1.3 Wed Jul 18 15:00:35 GMT 2001
35  
36  
37  
38   */
39  package edu.sc.seis.TauP;
40  
41  
42  public class Theta  {
43      
44      protected double radians;
45      protected double[] thetaAtX;
46      protected double[] rayParams;
47  
48      public Theta(SeismicPhase phase, double radians) {
49  	this.radians = radians;
50  
51  	int minRayParamIndex = phase.getMinRayParamIndex();
52  	int maxRayParamIndex = phase.getMaxRayParamIndex();
53  
54  	rayParams = phase.getRayParams();
55  	thetaAtX = phase.getTau();
56  	TauModel tMod = phase.getTauModel();
57  	// change tau to theta, theta = tau + p*x0
58  	for (int i=0; i<thetaAtX.length; i++) {
59  	    thetaAtX[i] += rayParams[i]*radians;
60  	}
61      }
62      
63      /***
64         * Get the value of radians.
65         * @return Value of radians.
66         */
67      public double getRadians() {return radians;}
68      
69      public double getMaxRayParam() {
70  	return rayParams[0];
71      }
72  
73      public double getStepRayParam(double rayParam, double timeStep) {
74  	double thetaStart = getTheta(rayParam);
75  
76  	// loop until we find a theta s.t. abs(thetaStart-theta) == timeStep
77  	// or we fall off the end of the array, ie ArrayIndexOutOfBounds
78  	boolean found = false;
79  	int i= getThetaIndex(rayParam);
80  	while (Math.abs(thetaAtX[i+1]-thetaStart) <= timeStep) {
81  	    i++;
82  	}
83  
84  	double newTheta;
85  	if (thetaStart <thetaAtX[i+1]) {
86  	    newTheta = thetaStart + timeStep;
87  	} else {
88  	    newTheta = thetaStart - timeStep;
89  	}
90  	// return the interpolated ray parameter.
91  	return linInterp(thetaAtX[i], thetaAtX[i+1],
92  			 rayParams[i], rayParams[i+1],
93  			 newTheta);
94  
95      }
96  
97      public double getTheta(double rayParam) {
98  	if (rayParam > rayParams[0] ||
99  	    rayParam < rayParams[rayParams.length-1]) {
100 	    throw new 
101 		ArrayIndexOutOfBoundsException(rayParam+
102 					       " not in range "+
103 					       rayParams[0]+" to "+
104 					       rayParams[rayParams.length-1]);
105 	}
106  
107 	int currentNum = getThetaIndex(rayParam);
108 	//find theta at given rayParam
109 	double thetaStart = 
110 	    linInterp(rayParams[currentNum], rayParams[currentNum+1],
111 		      thetaAtX[currentNum], thetaAtX[currentNum+1],
112 		      rayParam);
113 	return thetaStart;
114     }
115 
116 
117     protected int getThetaIndex(double rayParam) {
118 	// find index containing rayParam	
119 	int tooSmallNum = 0;
120 	int tooLargeNum = rayParams.length-1;
121 	int currentNum = 0;
122 	boolean found = false;
123 	while ( ! found) {
124 	    currentNum = (int)Math.floor((tooSmallNum + tooLargeNum) / 2.0f);
125 
126 	    if (rayParams[currentNum] >= rayParam &&
127 		rayParams[currentNum+1] < rayParam) {
128 		found = true;
129 	    } else if (rayParams[currentNum] > rayParam) {
130 		tooSmallNum = currentNum;
131 	    } else if (rayParams[currentNum] <= rayParam) {
132 		tooLargeNum = currentNum;
133 	    }
134 	}
135 	return currentNum;
136     }
137 
138     public static double linInterp(double xa, double xb, 
139 				   double ya, double yb, 
140 				   double x) {
141 	return (yb-ya)*(x-xa)/(xb-xa) + ya;
142     }
143 
144     
145 } // Theta