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  /***
29   * ReflTransCoefficient.java
30   *
31   *  Reflection and transmission coefficients for body waves. Methods
32   *  for calculating coefficients for each of the possible interactions are
33   *  provided. Calculations are done using the
34   *  com.visualnumerics.javagrande.Complex class from VisualNumerics.
35   *  It is further assume that the incoming ray is coming from the "top".
36   *  If the ray is actually coming from the bottom, the flip the velocities.
37   *
38   * @see "Aki and Richards page 144-151"
39   * @see "Lay and Wallace page 98 "
40   * @see <A HREF="http://math.nist.gov/javanumerics/">Java Numerics</A>
41   *
42   * Created: Wed Feb 17 12:25:27 1999
43   *
44   * @author Philip Crotwell
45   * @version 1.1.3 Wed Jul 18 15:00:35 GMT 2001
46  
47  
48  
49   */
50  
51  package edu.sc.seis.TauP;
52  
53  import java.io.Serializable;
54  //import com.visualnumerics.javagrande.*;
55  
56  
57  
58  public class ReflTransCoefficient implements Serializable, Cloneable {
59  
60      // IMPORTANT!!!!
61      // Where ever "CX" appears in this class, it is used as a shorthand for
62      // the Complex class, so CX.times() is the same as Complex.times, but
63      // the code is, IMHO, less cluttered.
64  
65      /*** just to avoid having Complex all over the place. */
66      private static final Complex CX = new Complex();
67  
68      protected double topVp;
69      protected double topVs;
70      protected double topDensity;
71      protected double botVp;
72      protected double botVs;
73      protected double botDensity;
74  
75      // "flat earth" ray parameter
76      protected double rp;
77  
78      // temp variables to make calculations less ugly
79      // first 3 lines follow both Aki and Richards and Lay and Wallace
80      protected double a, b, c, d;
81      protected Complex det, E, F, G, H;
82      protected Complex fsA; // used only in free surface calculations
83      // store the vertical slownesses for both the top and bottom halfspaces
84      // for both P and S waves
85      protected Complex topVertSlownessP, topVertSlownessS;
86      protected Complex botVertSlownessP, botVertSlownessS;
87      // we need the squared terms so often that it is worthwhile to store them
88      protected double sqBotVs;    // botVs squared
89      protected double sqTopVs;       // topVs squared
90      protected double sqBotVp;    // botVp squared
91      protected double sqTopVp;       // topVp squared
92      protected double sqRP; // rp squared
93  
94      // remember last calculated ray param and wave type to avoid repeating
95      protected double lastRayParam = -1.0;
96      protected boolean lastInIsPWave = true;
97  
98      public ReflTransCoefficient(double topVp,
99                  double topVs,
100                 double topDensity,
101                 double botVp,
102                 double botVs,
103                 double botDensity) {
104     this.topVp = topVp;
105     this.topVs = topVs;
106     this.topDensity = topDensity;
107     this.botVp = botVp;
108     this.botVs = botVs;
109     this.botDensity = botDensity;
110     }
111 
112     protected void calcTempVars(double rayParam, boolean inIsPWave) {
113 
114     this.rp = rayParam;   // ray parameter
115 
116     if (rayParam != lastRayParam && inIsPWave == lastInIsPWave) {
117         lastRayParam = -1.0;  // in case of failure in method
118     
119         sqBotVs = botVs * botVs;    // botVs squared
120         sqTopVs = topVs * topVs;       // topVs squared
121         sqBotVp = botVp * botVp;    // botVp squared
122         sqTopVp = topVp * topVp;       // topVp squared
123         sqRP = rp * rp; // rp squared
124 
125         topVertSlownessP = Complex.sqrt(new Complex(1.0 / sqTopVp - sqRP));
126         topVertSlownessS = Complex.sqrt(new Complex(1.0 / sqTopVs - sqRP));
127         botVertSlownessP = Complex.sqrt(new Complex(1.0 / sqBotVp - sqRP));
128         botVertSlownessS = Complex.sqrt(new Complex(1.0 / sqBotVs - sqRP));
129 
130         a = botDensity * (1.0 - 2 * sqBotVs * sqRP) -
131         topDensity  * (1.0 - 2 * sqTopVs  * sqRP);
132         b = botDensity * (1.0 - 2 * sqBotVs * sqRP) +
133         2 * topDensity * sqTopVs * sqRP;
134         c = topDensity * (1.0 - 2 * sqTopVs * sqRP) +
135         2 * botDensity * sqBotVs * sqRP;
136         d = 2 * (botDensity * sqBotVs - topDensity * sqTopVs);
137 
138         // math with complex objects is hard to read, so we give
139         // the formulas as comments
140 
141         // E = b * topVertSlownessP + c * botVertSlownessP
142         E = CX.plus(topVertSlownessP.times( b ),
143             botVertSlownessP.times( c ));
144 
145         // F = b * topVertSlownessS + c * botVertSlownessS
146         F = CX.plus(topVertSlownessS.times( b ),
147             botVertSlownessS.times( c ));
148 
149         // G = a - d * topVertSlownessP * botVertSlownessS
150         G = CX.minus(new Complex( a ),
151              CX.times(d,
152                   CX.times(topVertSlownessP,
153                        botVertSlownessS)));
154 
155         // H = a - d * botVertSlownessP * topVertSlownessS
156         H = CX.minus(new Complex( a ),
157              CX.times(d,
158                   CX.times(botVertSlownessP,
159                        topVertSlownessS)));
160 
161         // det = E * F + G * H * sqRP
162         det = CX.plus(CX.times( E, F),
163               CX.times(G, H).times(sqRP));
164 
165         // fsA = ((1/sqBotVs) - 2 * sqRP)^2 +
166         //       4 * sqRP * botVertSlownessP * botVertSlownessS
167         fsA = CX.plus(new Complex(((1/sqBotVs) - 2 * sqRP) *
168                       ((1/sqBotVs) - 2 * sqRP)),
169               CX.times(botVertSlownessP,
170                    botVertSlownessS).times(4*sqRP));
171         lastRayParam = rayParam;
172         lastInIsPWave = inIsPWave;
173     }
174     }
175 
176     // FREE SURFACE
177 
178     /*** Calculates incident P wave to reflected
179     P wave complex coefficient at free surface.
180     Only topVp, topVs, and topDensity are used, the bottom
181     values are ignored. This is a little strange as free surface rays are
182     always upgoing, but it mantains consistency with the solid-solid
183     calculations.
184     <P>
185     = (-1*((1/sqTopVs) - 2 * sqRP)^2 + <BR>
186     4 * sqRP * topVertSlownessP * topVertSlownessS) / A
187         
188     */
189     public  Complex getComplexFreePtoPRefl(double rayParam) {
190     calcTempVars(rayParam, true);
191     Complex numerator = CX.plus(new Complex( -1.0 *
192                          ((1/sqTopVs) - 2 * sqRP) *
193                          ((1/sqTopVs) - 2 * sqRP)),
194                     CX.times(topVertSlownessP,
195                          topVertSlownessS).times(4*sqRP));
196     return CX.over( numerator, fsA);
197     }
198     
199     /*** Calculates incident P wave to reflected
200     P wave coefficient at free surface.
201     @see #getComplexFreePtoPRefl(double)
202     */
203     public  double getFreePtoPRefl(double rayParam) {
204     return CX.abs(getComplexFreePtoPRefl(rayParam));
205     }
206        
207  
208     /***
209        Calculates incident P wave to reflected SV wave
210        complex coefficient at free surface.
211 
212        = (4 * (topVp/topVs) * rp * topVertSlownessP *
213        ((1/sqTopVs) - 2 * sqRP)) / fsA
214     */
215     public  Complex getComplexFreePtoSVRefl(double rayParam) {
216         calcTempVars(rayParam, true);
217     double realNumerator = 4 * (topVp/topVs) * rp *
218         ((1/sqTopVs) - 2 * sqRP);
219     return CX.over(CX.times(realNumerator, topVertSlownessP),
220                fsA);
221     }
222 
223     /***
224        Calculates incident P wave to reflected SV wave
225        coefficient at free surface.
226        @see #getComplexFreePtoSVRefl(double)
227     */
228     public  double getFreePtoSVRefl(double rayParam) {
229     return CX.abs(getComplexFreePtoSVRefl(rayParam));
230     }
231 
232 
233     /*** Calculates incident SV wave to reflected P wave
234     complex coefficient at free surface.
235     <P>
236     = (4 * (topVs/topVp) * rp * topVertSlownessS *<BR>
237     ((1/sqTopVs) - 2 * sqRP)) / fsA
238     */
239     public Complex  getComplexFreeSVtoPRefl(double rayParam) {
240     calcTempVars(rayParam, false);
241     double realNumerator = 4 * (topVs/topVp) * rp *
242         ((1/sqTopVs) - 2 * sqRP);
243     return CX.over(CX.times(realNumerator, topVertSlownessS),
244                fsA);
245     }
246 
247     /*** Calculates incident SV wave to reflected P wave
248     coefficient at free surface.
249     @see #getComplexFreeSVtoPRefl(double)
250     */
251     public  double getFreeSVtoPRefl(double rayParam) {
252         return CX.abs(getComplexFreeSVtoPRefl(rayParam));
253     }
254         
255     /*** Calculates incident SV wave to reflected SV wave
256         complex coefficient at free surface.
257     <P>
258 
259     = (-1 * ((1/sqTopVs) - 2 * sqRP)^2 + <BR>
260     4 * sqRP * topVertSlownessP * topVertSlownessS) / fsA
261 */
262     public Complex  getComplexFreeSVtoSVRefl(double rayParam) {
263     calcTempVars(rayParam, false);
264     double realNumerator = -1 * (((1/sqTopVs) - 2 * sqRP) *
265                      ((1/sqTopVs) - 2 * sqRP));
266     Complex numerator =  CX.plus(realNumerator,
267                      CX.times(4 * sqRP,
268                           CX.times(topVertSlownessP,
269                                topVertSlownessS)));
270     return CX.over(numerator, fsA);
271     }
272 
273     /*** Calculates incident SV wave to reflected SV wave
274        coefficient at free surface. */
275     public  double getFreeSVtoSVRefl(double rayParam) {
276         return CX.abs(getComplexFreeSVtoSVRefl(rayParam));
277     }
278         
279     /*** Calculates incident SH wave to reflected SH wave
280     complex coefficient at free surface. Free surface SH is always 1.*/
281     public Complex  getComplexFreeSHtoSHRefl(double rayParam) {
282     return  new Complex(1);
283     }
284         
285     /*** Calculates incident SH wave to reflected SH wave
286     coefficient at free surface. */
287     public  double getFreeSHtoSHRefl(double rayParam) {
288     return 1;
289     }
290 
291     // Solid-Solid interface
292 
293     /*** Calculates incident P wave to reflected
294     P wave Complex coefficient.
295     <P>
296     = ((b*topVertSlownessP - c*botVertSlownessP)*F - <BR>
297     (a + d*topVertSlownessP * botVertSlownessS)*H*sqRP) / det
298     */
299     public Complex  getComplexPtoPRefl(double rayParam) {
300     calcTempVars(rayParam, true);
301     Complex FTerm = CX.times(CX.minus(CX.times(b, topVertSlownessP),
302                       CX.times(c, botVertSlownessP)),
303                  F);
304     Complex HTerm = CX.times(CX.plus(a,
305                      CX.times(d,
306                           CX.times(topVertSlownessP,
307                                botVertSlownessS))),
308                  CX.times(H,
309                       sqRP));
310     return CX.over( CX.minus(FTerm, HTerm),
311             det);
312     }
313     
314     /*** Calculates incident P wave to reflected
315     P wave coefficient. */
316     public  double getPtoPRefl(double rayParam) {
317     return CX.abs(getComplexPtoPRefl(rayParam));
318     }
319     
320     /***
321        Calculates incident P wave to reflected SV wave
322        Complex coefficient.
323        <P>
324        = -2 * topVertSlownessP * <BR>
325        (a * b + c * d *botVertSlownessP *botVertSlownessS) * <BR>
326        rp * (topVp/topVs)) / det
327 
328     */
329     public Complex  getComplexPtoSVRefl(double rayParam) {
330         calcTempVars(rayParam, true);
331     double realNumerator = -2 * rp * (topVp/topVs);
332     Complex middleTerm = CX.plus(a * b,
333                      CX.times(c * d,
334                           CX.times( botVertSlownessP,
335                             botVertSlownessS)));
336     Complex numerator = CX.times(CX.times(realNumerator,
337                      topVertSlownessP),
338                      middleTerm);
339     return CX.over(numerator, det);
340     }
341         
342     /***
343        Calculates incident P wave to reflected SV wave
344        coefficient. */
345     public  double getPtoSVRefl(double rayParam) {
346         return CX.abs(getComplexPtoSVRefl(rayParam));
347     }
348     
349     /*** Calculates incident P wave to transmitted P wave
350     Complex coefficient.
351     <P>
352     = ( 2 * topDensity * topVertSlownessP * F * <BR>
353     (topVp / botVp)) / det
354 
355     */
356     public Complex  getComplexPtoPTrans(double rayParam) {
357     calcTempVars(rayParam, true);
358 
359     double realNumerator = 2 * topDensity * (topVp / botVp);
360     return CX.over( CX.times(realNumerator,
361                  CX.times(topVertSlownessP,
362                       F)),
363             det);
364     }
365         
366     /*** Calculates incident P wave to transmitted P wave
367     coefficient. */
368     public  double getPtoPTrans(double rayParam) {
369         return CX.abs(getComplexPtoPTrans(rayParam));
370     }
371     
372     /***
373        Calculates incident P wave to transmitted SV wave
374        Complex coefficient.
375        <P>
376        = (2 * topDensity * topVertSlownessP * H * rp * (topVp / botVs)) / <BR>
377        det
378     */
379     public Complex getComplexPtoSVTrans(double rayParam) {
380     calcTempVars(rayParam, true);
381 
382     double realNumerator = 2 * topDensity * rp * (topVp / botVs);
383     Complex numerator = CX.times(realNumerator,
384                      CX.times(topVertSlownessP,
385                           H));
386     return CX.over(numerator,
387                det);
388     }
389     
390     /***
391        Calculates incident P wave to transmitted SV wave
392        coefficient. */
393     public  double getPtoSVTrans(double rayParam) {
394         return CX.abs(getComplexPtoSVTrans(rayParam));
395     }
396 
397     /*** Calculates incident SV wave to reflected P wave
398     Complex coefficient.
399     <P>
400     = (-2 * topVertSlownessS * <BR>
401     (a * b + c * d *  botVertSlownessP * botVertSlownessS) * <BR>
402     rp * (topVs / topVp)) / <BR>
403     det
404     */
405     public Complex  getComplexSVtoPRefl(double rayParam) {
406         calcTempVars(rayParam, false);
407     double realNumerator = -2 * rp * (topVs / topVp);
408     //double realNumerator = -2 * rp * (topVs / topVp);
409     Complex middleTerm = CX.plus( a * b,
410                       CX.times(c * d,
411                            CX.times(botVertSlownessP,
412                             botVertSlownessS)));
413     Complex numerator = CX.times(realNumerator,
414                      CX.times(topVertSlownessS,
415                           middleTerm));
416     return CX.over(numerator, det);
417     }
418     
419     /*** Calculates incident SV wave to reflected P wave
420     coefficient. */
421     public  double getSVtoPRefl(double rayParam) {
422         return CX.abs(getComplexSVtoPRefl(rayParam));
423     }
424     
425     /*** Calculates incident SV wave to reflected SV wave
426        Complex coefficient.
427        <P>
428        = -1 * ((b * topVertSlownessS - c * botVertSlownessS) * E - <BR>
429        (a + b * botVertSlownessP * topVertSlownessS) * G * sqRP) / <BR>
430        det
431 
432     */
433     public Complex  getComplexSVtoSVRefl(double rayParam) {
434         calcTempVars(rayParam, false);
435     Complex abNumerator =
436         CX.times( CX.plus(a,
437                   CX.times(d,
438                        CX.times(botVertSlownessP,
439                         topVertSlownessS))),
440               CX.times( G, sqRP));
441     Complex bcNumerator = CX.times( CX.plus( CX.times( b,
442                                topVertSlownessS),
443                          CX.times( c,
444                                botVertSlownessS)),
445                     E);
446     return CX.over( CX.minus(abNumerator, bcNumerator),
447             det);
448     }
449         
450     /*** Calculates incident SV wave to reflected SV wave
451        coefficient. */
452     public  double getSVtoSVRefl(double rayParam) {
453         return CX.abs(getComplexSVtoSVRefl(rayParam));
454     }
455     
456     /*** Calculates incident SV wave to transmitted P wave
457     Complex coefficient.
458     <P>
459     = -2 * topDensity * topVertSlownessS * G * rp * (topVs / botVp) / <BR>
460     det
461     */
462     public Complex getComplexSVtoPTrans(double rayParam) {
463     calcTempVars(rayParam, false);
464     double realNumerator = -2 * topDensity * rp * (topVs / botVp);
465     Complex numerator = CX.times(realNumerator,
466                      CX.times(topVertSlownessS,
467                           G ));
468     return CX.over(numerator, det);
469     }
470         
471     /*** Calculates incident SV wave to transmitted P wave
472     coefficient.
473     */
474     public  double getSVtoPTrans(double rayParam) {
475         return CX.abs(getComplexSVtoPTrans(rayParam));
476     }
477     
478     /***
479        Calculates incident SV wave to transmitted SV wave
480        Complex coefficient.
481        <P>
482        = 2 * topDensity * topVertSlownessS * E * (topVs / botVs) / <BR>
483        det
484     */
485     public Complex getComplexSVtoSVTrans(double rayParam) {
486     calcTempVars(rayParam, false);
487     double realNumerator = 2 * topDensity * rp * (topVs / botVs);
488     Complex numerator = CX.times(realNumerator,
489                      CX.times(topVertSlownessS,
490                           E ));
491     return CX.over(numerator, det);
492     }
493          
494     /***
495        Calculates incident SV wave to transmitted SV wave
496        coefficient.
497     */
498     public  double getSVtoSVTrans(double rayParam) {
499         return CX.abs(getComplexSVtoSVTrans(rayParam));
500     }
501        
502     // SH waves
503 
504 
505     /*** Calculates incident SH wave to reflected SH wave
506     Complex coefficient.
507     <P>
508     mu = Vs * Vs * density
509     <P>
510     = (topMu * topVertSlownessS - botMu * botVertSlownessS) / <BR>
511       (topMu * topVertSlownessS + botMu * botVertSlownessS)
512     */
513     public Complex  getComplexSHtoSHRefl(double rayParam) {
514         calcTempVars(rayParam, false);
515     double topMu = topVs * topVs * topDensity;
516     double botMu = botVs * botVs * botDensity;
517     Complex topTerm = CX.times(topMu, topVertSlownessS);
518     Complex botTerm = CX.times(botMu, botVertSlownessS);
519 
520     return CX.over( CX.minus(topTerm, botTerm),
521             CX.plus(topTerm, botTerm));
522     }
523 
524     /*** Calculates incident SH wave to reflected SH wave
525     coefficient. */
526     public  double getSHtoSHRefl(double rayParam) {
527         return CX.abs(getComplexSHtoSHRefl(rayParam));
528     }
529 
530     /*** Calculates incident SH wave to transmitted SH wave
531     Complex coefficient.
532     <P>
533     mu = Vs * Vs * density
534     <P>
535     = 2 * topMu * topVertSlownessS / <BR>
536     (topMu * topVertSlownessS + botMu * botVertSlownessS)
537     */
538     public Complex  getComplexSHtoSHTrans(double rayParam) {
539         calcTempVars(rayParam, false);
540     double topMu = topVs * topVs * topDensity;
541     double botMu = botVs * botVs * botDensity;
542     Complex topTerm = CX.times(topMu, topVertSlownessS);
543     Complex botTerm = CX.times(botMu, botVertSlownessS);
544 
545     return CX.over( CX.times(topTerm, 2 ),
546             CX.plus(topTerm, botTerm));
547     }
548 
549     /*** Calculates incident SH wave to transmitted SH wave
550     coefficient. */
551     public  double getSHtoSHTrans(double rayParam) {
552         return CX.abs(getComplexSHtoSHTrans(rayParam));
553     }
554 
555     public static void main(String[] args) {
556     double topVp = 4.98;
557     double topVs = 2.9;
558     double topDensity = 2.667;
559     double botVp = 8.0;
560     double botVs = 4.6;
561     double botDensity = 3.38;
562     double depth;
563     double radiusOfEarth;
564 
565     double DtoR = Math.PI / 180.0;
566     double RtoD = 180.0 / Math.PI;
567 
568     double[] RPP = new double[91];
569     double[] RPS = new double[91];
570     double[] RSP = new double[91];
571     double[] RSS = new double[91];
572 
573     double[] TPP = new double[91];
574     double[] TPS = new double[91];
575     double[] TSP = new double[91];
576     double[] TSS = new double[91];
577 
578     ReflTransCoefficient coeff =
579         new ReflTransCoefficient(topVp,
580                      topVs,
581                      topDensity,
582                      botVp,
583                      botVs,
584                      botDensity);
585     double rayParam;
586     for (int i=0; i<= 90; i++) {
587         rayParam = Math.sin(DtoR * i) / topVp;
588         RPP[i] = coeff.getPtoPRefl(rayParam);
589         RPS[i] = coeff.getPtoSVRefl(rayParam);
590         TPP[i] = coeff.getPtoPTrans(rayParam);
591         TPS[i] = coeff.getPtoSVTrans(rayParam);
592 
593         rayParam = Math.sin(DtoR * i) / topVs;
594         RSP[i] = coeff.getSVtoPRefl(rayParam);
595         RSS[i] = coeff.getSVtoSVRefl(rayParam);
596         TSP[i] = coeff.getSVtoPTrans(rayParam);
597         TSS[i] = coeff.getSVtoSVTrans(rayParam);
598     }
599 
600     try {
601         java.io.Writer out =
602         new java.io.BufferedWriter(
603             new java.io.FileWriter("refltrans.gmt"));
604         out.write("#!/bin/sh\n\n");
605         out.write("/bin/rm -f refltrans.ps\n\n");
606         out.write("psbasemap -K -P -JX6 -R0/90/0/2  -B10/1 > refltrans.ps\n");
607 
608 //      out.write("psxy -K -O -JX -R -M -W1/255/0/0 >> refltrans.ps <<END\n");
609 //      for (int i=0; i<=90; i++) {
610 //      out.write(i+" "+RPP[i]+"\n");
611 //      }
612 //      out.write("END\n");
613 
614 //      out.write("psxy -K -O -JX -R -M -W1/0/255/0 >> refltrans.ps <<END\n");
615 //      for (int i=0; i<=90; i++) {
616 //      out.write(i+" "+RPS[i]+"\n");
617 //      }
618 //      out.write("END\n");
619 //      out.write("psxy -K -O -JX -R -M -W1/0/0/255 >> refltrans.ps <<END\n");
620 //      for (int i=0; i<=90; i++) {
621 //      out.write(i+" "+TPP[i]+"\n");
622 //      }
623 //      out.write("END\n");
624 //      out.write("psxy -K -O -JX -R -M -W1/255/255/0 >> refltrans.ps <<END\n");
625 //      for (int i=0; i<=90; i++) {
626 //      out.write(i+" "+TPS[i]+"\n");
627 //      }
628 //      out.write("END\n");
629         out.write("psxy -K -O -JX -R -M -W1/255/0/0 >> refltrans.ps <<END\n");
630         for (int i=0; i<=90; i++) {
631         out.write(i+" "+RSP[i]+"\n");
632         }
633         out.write("END\n");
634         out.write("psxy -K -O -JX -R -M -W1/0/255/0 >> refltrans.ps <<END\n");
635         for (int i=0; i<=90; i++) {
636         out.write(i+" "+RSS[i]+"\n");
637         }
638         out.write("END\n");
639         out.write("psxy -K -O -JX -R -M -W1/0/0/255 >> refltrans.ps <<END\n");
640         for (int i=0; i<=90; i++) {
641         out.write(i+" "+TSP[i]+"\n");
642         }
643         out.write("END\n");
644         out.write("psxy -O -JX -R -M -W1/255/0/255 >> refltrans.ps <<END\n");
645         for (int i=0; i<=90; i++) {
646         out.write(i+" "+TSS[i]+"\n");
647         }
648         out.write("END\n");
649         out.close();
650     } catch (java.io.IOException e) {
651         System.err.println(e);
652     }
653     }
654 
655 } // ReflTransCoefficient