1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28 /***
29 * TauP_WKBJ.java
30 *
31 *
32 *
33 *
34 * @author Philip Crotwell
35 * @version 1.1.3 Wed Jul 18 15:00:35 GMT 2001
36
37
38
39 */
40
41 package edu.sc.seis.TauP;
42
43 import java.io.FileNotFoundException;
44 import java.io.IOException;
45 import java.io.OptionalDataException;
46 import java.io.StreamCorruptedException;
47 import edu.sc.seis.seisFile.sac.SacTimeSeries;
48
49 public class TauP_WKBJ extends TauP_Time {
50
51 /*** deltaT of the seismogram, default is .05 which gives
52 20 sps. */
53 protected double deltaT = .05;
54
55 /*** number of samples in the seismogram. Default is 100. */
56 protected int numSamples = 1000;
57
58 /*** start time of the seismogram relative to the origin time.
59 default is 0. */
60 protected double startTime = 0;
61
62 public TauP_WKBJ() {
63 super();
64 }
65
66 public TauP_WKBJ(TauModel tMod) throws TauModelException {
67 super(tMod);
68 }
69
70 public TauP_WKBJ(String modelName) throws TauModelException {
71 super(modelName);
72 }
73
74 /***
75 * Get the value of deltaT.
76 * @return Value of deltaT.
77 */
78 public double getDeltaT() {return deltaT;}
79
80 /***
81 * Set the value of deltaT.
82 * @param v Value to assign to deltaT.
83 */
84 public void setDeltaT(double v) {this.deltaT = v;}
85
86 /***
87 * Get the value of numSamples.
88 * @return Value of numSamples.
89 */
90 public int getNumSamples() {return numSamples;}
91
92 /***
93 * Set the value of numSamples.
94 * @param v Value to assign to numSamples.
95 */
96 public void setNumSamples(int v) {this.numSamples = v;}
97
98 /***
99 * Get the value of startTime.
100 * @return Value of startTime.
101 */
102 public double getStartTime() {return startTime;}
103
104 /***
105 * Set the value of startTime.
106 * @param v Value to assign to startTime.
107 */
108 public void setStartTime(double v) {this.startTime = v;}
109
110
111 public void calculate(double degrees) throws TauModelException {
112 recalcPhases();
113 clearArrivals();
114 calcWKBJ(degrees);
115 }
116
117 public void calcWKBJ(double degrees) throws TauModelException {
118 this.degrees = degrees;
119 System.out.println("In calcWKBJ for "+degrees+" degrees.");
120
121 SeismicPhase phase;
122 Arrival[] phaseArrivals;
123 double calcTime, calcDist;
124 Theta thetaAtX;
125 double rayParam, nextRayParam, minRayParam;
126 double theta, nextTheta;
127
128 ReflTransCoefficient rtCoeff;
129
130 for (int phaseNum=0;phaseNum<phases.size(); phaseNum++) {
131 System.out.println("Phase "+((SeismicPhase)phases.elementAt(phaseNum)).getName()+".");
132 phase = (SeismicPhase)phases.elementAt(phaseNum);
133 phase.setDEBUG(DEBUG);
134 phase.calcTime(degrees);
135 minRayParam = phase.getMinRayParam();
136
137 if (phase.hasArrivals()) {
138
139
140 phaseArrivals = phase.getArrivals();
141
142 for (int i=0;i<phaseArrivals.length;i++) {
143 System.out.println("Arrival number "+i);
144 thetaAtX = new Theta(phase, phaseArrivals[i].getDist());
145 System.out.println("Got Theta");
146 float[] seismogramPoints = new float[numSamples];
147 rayParam = thetaAtX.getMaxRayParam();
148 System.out.println("Got ray param");
149 theta = thetaAtX.getTheta(rayParam);
150 System.out.println("Got theta for ray param");
151 setStartTime(320);
152 nextRayParam = thetaAtX.getStepRayParam(rayParam,
153 getDeltaT());
154 nextTheta = thetaAtX.getTheta(nextRayParam);
155 int n = 0;
156 try {
157 while (nextRayParam >= minRayParam) {
158
159 n = (int)Math.round((theta-getStartTime()) /
160 getDeltaT());
161 if (n >= 0 && n < seismogramPoints.length) {
162
163
164
165
166 System.out.println(n+" "+seismogramPoints[n]);
167
168 }
169 rayParam = nextRayParam;
170 theta = nextTheta;
171 nextRayParam = thetaAtX.getStepRayParam(rayParam,
172 getDeltaT());
173 nextTheta = thetaAtX.getTheta(nextRayParam);
174 }
175 } catch (ArrayIndexOutOfBoundsException e) {
176
177 System.out.println("ArrayIndexOutOfBoundsException: "+e);
178 }
179
180
181 SacTimeSeries sac = new SacTimeSeries();
182 sac.y = seismogramPoints;
183 sac.npts = seismogramPoints.length;
184 sac.leven = SacTimeSeries.TRUE;
185 sac.iftype = SacTimeSeries.ITIME;
186 sac.delta = (float)getDeltaT();
187 sac.t0 = (float)phaseArrivals[i].getTime();
188 sac.kt0 = phaseArrivals[i].getName();
189 sac.o = 0;
190 sac.b = 320;
191 sac.e = sac.b + (sac.npts-1)* sac.delta;
192
193 try {
194 sac.write("tempsacfile");
195 } catch (IOException e) {
196 }
197
198 }
199 }
200 }
201 }
202
203 /*** Allows TauP_Time to run as an application. Creates an instance
204 * of TauP_Time.
205 * . */
206 public static void main(String[] args)
207 throws FileNotFoundException,
208 IOException,
209 StreamCorruptedException,
210 ClassNotFoundException,
211 OptionalDataException
212 {
213 try {
214 long prevTime = 0;
215 long currTime;
216
217 prevTime = System.currentTimeMillis();
218 TauP_Time tauPTime = new TauP_WKBJ();
219 String[] noComprendoArgs = tauPTime.parseCmdLineArgs(args);
220 if (noComprendoArgs.length > 0) {
221 for (int i=0;i<noComprendoArgs.length;i++) {
222 if (noComprendoArgs[i].equals("-help") ||
223 noComprendoArgs[i].equals("-version")) {
224 System.exit(0);
225 }
226 }
227 System.out.println("I don't understand the following arguments, continuing:");
228 for (int i=0;i<noComprendoArgs.length;i++) {
229 System.out.print(noComprendoArgs[i]+" ");
230 }
231 System.out.println();
232 noComprendoArgs = null;
233 }
234 currTime = System.currentTimeMillis();
235
236 prevTime = System.currentTimeMillis();
237 tauPTime.init();
238 currTime = System.currentTimeMillis();
239 if (tauPTime.DEBUG) {
240 System.out.println("taup model read time="+(currTime-prevTime));
241 }
242
243 tauPTime.start();
244 tauPTime.destroy();
245
246 } catch (TauModelException e) {
247 System.out.println("Caught: "+e.getMessage());
248 e.printStackTrace();
249 } catch (TauPException e) {
250 System.out.println("Caught: "+e.getMessage());
251 e.printStackTrace();
252 }
253
254 }
255
256 }