Main Page   Class Hierarchy   Compound List   File List   Compound Members   File Members  

tlay.c

Go to the documentation of this file.
00001 
00002 /*
00003  *   THIS FILE IS UNDER RCS - DO NOT MODIFY UNLESS YOU HAVE
00004  *   CHECKED IT OUT USING THE COMMAND CHECKOUT.
00005  *
00006  *    $Id: tlay_8c-source.html 2161 2006-05-19 16:55:03Z paulf $
00007  *
00008  *    Revision history:
00009  *     $Log$
00009  *     Revision 1.1  2006/05/19 16:55:02  paulf
00009  *     first inclusion
00009  *
00010  *     Revision 1.1  2000/02/14 18:51:48  lucky
00011  *     Initial revision
00012  *
00013  *
00014  */
00015 
00016 /*
00017  * tlay.c : Travel-time in layered half space (ala Eaton)
00018  *
00019  *$ 95Oct19 LDD Explicitly declared return types for all functions.
00020  *              Added function prototypes to tlay.h
00021  */
00022 /*********************C O P Y R I G H T   N O T I C E ***********************/
00023 /* Copyright 1991 by Carl Johnson.  All rights are reserved. Permission     */
00024 /* is hereby granted for the use of this product for nonprofit, commercial, */
00025 /* or noncommercial publications that contain appropriate acknowledgement   */
00026 /* of the author. Modification of this code is permitted as long as this    */
00027 /* notice is included in each resulting source module.                      */
00028 /****************************************************************************/
00029 #include <math.h>
00030 #include <stdio.h>
00031 #include <stdlib.h>
00032 #include <kom.h>
00033 #include <tlay.h>
00034 
00035 #define Panic(x) (fprintf(stderr,"Panic point %d in tlay.c\n",(x)),exit(-1))
00036 
00037 #define MAXLAY 20
00038 static int nLay = 0;
00039 int iRef;
00040 
00041 /* Function prototypes from other source files*/
00042 void mnbrak(float *, float *, float *, float *, float *, float *, float(*)(float) );
00043 float brent(float, float, float, float (*)(float), float, float *);
00044 
00045 /* Global variables is used by t_fun and t_dis function                 */
00046 static int    iSrc;     /* Source layer index, for direct ray           */
00047 static double zSrc;     /* Source depth, for direct ray calc            */
00048 static double xSrc;     /* Source distance, for direct ray calc         */
00049 
00050 /* Global variables set by t_set and used by t_lay for fast calculations*/
00051 static double zTop[2*MAXLAY];
00052 static double vLay[2*MAXLAY];
00053 
00054 /* Travel time command processing       */
00055 int initMod   = 0;
00056 double PoverS = 1.72;
00057 
00058 /***********************************************************************
00059  *   t_com():  Process all recognized commands.                        *
00060  *             Return 1 on success, 0 if command was not recognized    *
00061  ***********************************************************************/
00062 int t_com( void )
00063 {
00064         double z, v, r, t, dtdr, dtdz;
00065         double tz, tr;
00066         double rmax, rdel;
00067         TPHASE treg[4];
00068         int np;
00069         int i;
00070 
00071         if(k_its("lay")) {
00072                 z = k_val();
00073                 v = k_val();
00074                 t_model(z, v);
00075                 return 1;
00076         }
00077 
00078         if(k_its("psratio")) {
00079                 PoverS = k_val();
00080                 return 1;
00081         }
00082 
00083         if(k_its("ray")) {
00084                 z = k_val();
00085                 r = k_val();
00086                 t = t_lay(r, z, &dtdr, &dtdz);
00087                 fprintf(stdout, "%6.2f : %6.2f %6.2f %6.2f %6.2f\n",
00088                         t, r, z, dtdr, dtdz);
00089                 return 1;
00090         }
00091 
00092         if(k_its("raytst")) {
00093                 t_set();
00094                 z = k_val();
00095                 rmax = k_val();
00096                 rdel = k_val();
00097                 for(r=rdel; r<rmax; r+=rdel) {
00098                         tr = t_lay(r+0.1, z, &dtdr, &dtdz);
00099                         tz = t_lay(r, z+0.1, &dtdr, &dtdz);
00100                         t = t_lay(r, z, &dtdr, &dtdz);
00101                         fprintf(stdout,
00102                                 "%d %6.1f %6.1f %6.2f %6.2f %6.2f %6.2f\n",
00103                                 iRef, z, r, dtdr, dtdz,
00104                                 10.0*(tr-t), 10.0*(tz-t));
00105                 }
00106                 return 1;
00107         }
00108 
00109         if(k_its("region")) {
00110                 t_set();
00111                 z = k_val();
00112                 rmax = k_val();
00113                 rdel = k_val();
00114                 for(r=rdel; r<rmax; r+=rdel) {
00115                         np = t_region(r, z, treg);
00116                         fprintf(stdout, "%6.1f : ", r);
00117                         for(i=0; i<np; i++)
00118                                 fprintf(stdout, "%6.2f ", treg[i].t);
00119                         fprintf(stdout, "\n");
00120                 }
00121                 return 1;
00122         }
00123 
00124         return 0;
00125 }
00126 
00127 /**************************************************************************
00128  * t_model()  Add layer to velocity model.                                *
00129  **************************************************************************/
00130 int t_model(double z, double v)
00131 {
00132         int i;
00133 
00134         if(nLay < MAXLAY) {
00135                 zTop[nLay] = z;
00136                 vLay[nLay] = v;
00137                 nLay++;
00138         }
00139 
00140         for(i=1; i<nLay; i++) {
00141                 zTop[2*nLay-1-i] = 2*zTop[nLay-1] - zTop[i];
00142                 vLay[2*nLay-1-i] = vLay[i-1];
00143         }
00144         zTop[nLay] = zTop[nLay-1] + 0.01;
00145 
00146 /*      printf("nLay = %d\n", nLay);
00147         for(i=0; i<2*nLay-1; i++)
00148                 printf("%d : %6.1f %6.1f\n", i, zTop[i], vLay[i]); */ /*DEBUG*/
00149         return nLay;
00150 }
00151 
00152 /**************************************************************************
00153  *   t_set()  Set up travel time calculations                             *
00154  **************************************************************************/
00155 int t_set( void )
00156 {
00157         if(initMod)
00158                 return nLay;
00159         initMod = 1;
00160         return nLay;
00161 }
00162 
00163 /**************************************************************************
00164  *   t_fun()  Calculate direct travel time from takeoff angle             *
00165  **************************************************************************/
00166 float t_fun(float r)
00167 {
00168         int i;
00169         double q;
00170         double p;
00171         float t;
00172         float td;
00173 
00174         p = atan(r/zSrc);
00175         t = 0.0;
00176         for(i=0; i<=iSrc; i++) {
00177                 if(i == iSrc) {
00178                         td = (float) (fabs(zSrc - zTop[i]) / cos(p) / vLay[i]);
00179                 } else {
00180                         q = asin(vLay[i]*sin(p)/vLay[iSrc]);
00181                         td = (float)((zTop[i+1] - zTop[i]) / cos(q) / vLay[i]);
00182                 }
00183                 t += td;
00184         }
00185         return t;
00186 }
00187 
00188 /**************************************************************************
00189  *   t_dis()  Calculate direct ray distance from takeoff angle            *
00190  **************************************************************************/
00191 float t_dis( float r )
00192 {
00193         int i;
00194         double q;
00195         double p;
00196         float x;
00197         float xd;
00198 
00199         p = atan(r/zSrc);
00200         x = (float)0.0;
00201         for(i=0; i<=iSrc; i++) {
00202                 if(i == iSrc) {
00203                         xd = (float) ((zSrc - zTop[i]) * tan(p));
00204                 } else {
00205                         q = asin(vLay[i]*sin(p)/vLay[iSrc]);
00206                         xd = (float)((zTop[i+1] - zTop[i]) * tan(q));
00207                 }
00208                 x += xd;
00209         }
00210         return (float) ((x-xSrc)*(x-xSrc));
00211 }
00212 
00213 /**************************************************************************
00214  *  t_lay()  Calculate travel times                                       *
00215  **************************************************************************/
00216 double t_lay( double  r,        /* Epicentral distance          */
00217               double  z,        /* Hypocentral depth            */
00218               double *dtdr,     /* dTdR                         */
00219               double *dtdz )    /* dTdZ                         */
00220 {
00221         double tmin;            /* Minimum trav time            */
00222         double t;               /* Travel time                  */
00223         double x;               /* Critical dist                */
00224         double p;               /* Takeoff angle                */
00225         double q;               /* Layer incidence              */
00226         double td;              /* Travel time inc              */
00227         double xd;              /* Crit. dist inc               */
00228         double tdir;
00229         double toa;
00230         double sgn;
00231         float ax, bx, cx;
00232         float fa, fb, fc;
00233         float xmin;
00234         int isrc;
00235         int i, j;
00236 
00237         t_set();
00238         tmin = 10000.0;
00239         sgn = 1.0;
00240         if(z < 0.0) {
00241                 z = -z;
00242                 sgn = -sgn;
00243         }
00244 /* Calculate source layer                       */
00245         isrc = 0;
00246         for(i=1; i<nLay; i++) {
00247                 if(z > zTop[i])
00248                         isrc = i;
00249         }
00250         iSrc = isrc;
00251         zSrc = z;
00252         if ( zSrc < .01 ) zSrc = .01;              /* Added by WMK 2/12/96 */
00253         xSrc = r;
00254 
00255 /* Calculate minimum refraction         */
00256         iRef = 0;
00257         for(i=isrc+1; i<nLay; i++) {
00258                 t = 0.0;
00259                 x = 0.0;
00260                 p = 1.0/vLay[i];
00261                 for(j=0; j<i; j++) {
00262                         q = asin(p * vLay[j]);
00263                         td = (zTop[j+1]-zTop[j]) / cos(q) / vLay[j];
00264                         xd = (zTop[j+1]-zTop[j]) * tan(q);
00265                         t += td;                /* Upgoing ray          */
00266                         x += xd;
00267                         if(j == isrc) { /* Source in layer      */
00268                                 t += td * (zTop[j+1] - z)
00269                                         / (zTop[j+1] - zTop[j]);
00270                                 x += xd * (zTop[j+1] - z)
00271                                         / (zTop[j+1] - zTop[j]);
00272                         }
00273                         if(j > isrc) {  /* Downgoing ray        */
00274                                 t += td;
00275                                 x += xd;
00276                         }
00277                 }
00278                 if(x < r) {
00279                         t += (r-x) / vLay[i];
00280                         if(t < tmin) {
00281                                 iRef = i;
00282                                 tmin = t;
00283                                 toa = 3.1415927 - asin(p * vLay[isrc]);
00284                         }
00285                 }
00286         }
00287 
00288 /* Travel time of direct ray */
00289         ax = (float)(0.8 * xSrc);
00290         bx = (float)(1.5 * xSrc);
00291         mnbrak(&ax, &bx, &cx, &fa, &fb, &fc, t_dis);
00292         brent(ax, bx, cx, t_dis, (float)0.001, &xmin);
00293         tdir = t_fun(xmin);
00294         if(tdir < tmin) {
00295                 iRef = 0;
00296                 tmin = tdir;
00297                 toa = atan(xmin/zSrc);     /* changed from z to zSrc by WMK */
00298                 if(toa < 0.0)
00299                         toa += 3.141592654;
00300         }
00301 
00302         *dtdz = sgn * cos(toa)/vLay[isrc];
00303         *dtdr = sin(toa)/vLay[isrc];
00304         return tmin;
00305 }
00306 
00307 /**************************************************************************
00308  *  t_direct()  Calculate travel time of direct P                         *
00309  **************************************************************************/
00310 double t_direct( double  r,
00311                  double  z,
00312                  double *dtdr,
00313                  double *dtdz )
00314 {
00315         double toa;
00316         double tp;
00317         float ax, bx, cx;
00318         float fa, fb, fc;
00319         float xmin;
00320         int i;
00321 
00322 /* source layer */
00323         iSrc = 0;
00324         for(i=1; i<nLay; i++)
00325                 if(z > zTop[i])
00326                         iSrc = i;
00327         zSrc = z;
00328         xSrc = r;
00329 
00330 /* Calculate travel time, need to take care of Z > 0 by direct calcuation */
00331         ax = (float)(0.8 * xSrc);
00332         bx = (float)(1.5 * xSrc);
00333         mnbrak(&ax, &bx, &cx, &fa, &fb, &fc, t_dis);
00334         brent(ax, bx, cx, t_dis, (float) 0.001, &xmin);
00335         toa = atan(xmin/z);
00336         if(toa < 0.0)
00337                 toa += 3.141592654;
00338         tp = t_fun(xmin);
00339         *dtdz = cos(toa) / vLay[iSrc];
00340         *dtdr = sin(toa) / vLay[iSrc];
00341         return tp;
00342 }
00343 
00344 /**************************************************************************
00345  *  t_pmp()  Calculate travel time of mantle reflection                   *
00346  **************************************************************************/
00347 double t_pmp( double  r,
00348               double  z,
00349               double *dtdr,
00350               double *dtdz )
00351 {
00352         double toa;
00353         double tpmp;
00354         float ax, bx, cx;
00355         float fa, fb, fc;
00356         float xmin;
00357         int nlay;
00358         int i;
00359 
00360 /* source layer */
00361         nlay = nLay;
00362         nLay = 2*nlay - 1;
00363         iSrc = 0;
00364         for(i=1; i<nLay; i++)
00365                 if(z > zTop[i])
00366                         iSrc = i;
00367         zSrc = z;
00368         xSrc = r;
00369 
00370 /* Calcuate travel time of direct ray from image source         */
00371         ax = (float)(0.8 * xSrc);
00372         bx = (float)(1.5 * xSrc);
00373         mnbrak(&ax, &bx, &cx, &fa, &fb, &fc, t_dis);
00374         brent(ax, bx, cx, t_dis, (float)0.001, &xmin);
00375         toa = atan(xmin/z);
00376         if(toa < 0.0)
00377                 toa += 3.141592654;
00378         tpmp = t_fun(xmin);
00379         *dtdz = -cos(toa) / vLay[iSrc];
00380         *dtdr = sin(toa) / vLay[iSrc];
00381         nLay = nlay;
00382         return tpmp;
00383 }
00384 
00385 /**************************************************************************
00386  *  t_phase()  Calculate travel time for a given phase                    *
00387  **************************************************************************/
00388 double t_phase( int     ph,             /* Phase index (0-5)    */
00389                 double  r,              /* Epicentral distance  */
00390                 double  z,              /* Hypocentral depth    */
00391                 double *dtdr,           /* dTdR                 */
00392                 double *dtdz )          /* dTdZ                 */
00393 {
00394         double t;
00395 
00396         t_set();
00397         switch(ph) {
00398         case 0: /* P    */
00399         case 2: /* Pn   */
00400                 t = t_lay(r, z, dtdr, dtdz);
00401                 break;
00402         case 1: /* S    */
00403         case 3: /* Sn   */
00404                 t = PoverS * t_lay(r, z, dtdr, dtdz);
00405                 *dtdr *= PoverS;
00406                 *dtdz *= PoverS;
00407                 break;
00408         case 4: /* Pg   */
00409                 nLay--;
00410                 t = t_lay(r, z, dtdr, dtdz);
00411                 nLay++;
00412                 break;
00413         case 5: /* Sg   */
00414                 nLay--;
00415                 t = PoverS * t_lay(r, z, dtdr, dtdz);
00416                 nLay++;
00417                 *dtdr *= PoverS;
00418                 *dtdz *= PoverS;
00419                 break;
00420         case 6: /* P*, Direct P */
00421                 t = t_direct(r, z, dtdr, dtdz);
00422                 break;
00423         case 7: /* PmP */
00424                 t = t_pmp(r, z, dtdr, dtdz);
00425                 break;
00426         }
00427         return t;
00428 }
00429 
00430 /**************************************************************************
00431  *  t_region()  Calculate regional phase travel times: P, Pg, S, and Sg   *
00432  **************************************************************************/
00433 int t_region( double  r,        /* Epicentral distance                  */
00434               double  z,        /* Hypocentral depth                    */
00435               TPHASE *treg)     /* Travtime and derive for each phase   */
00436 {
00437         double t;
00438         double dtdr;
00439         double dtdz;
00440 
00441         t = t_lay(r, z, &dtdr, &dtdz);
00442         treg[0].phase = 0;
00443         treg[0].t = t;
00444         treg[0].dtdr = dtdr;
00445         treg[0].dtdz = dtdz;
00446         treg[1].phase = 1;
00447         treg[1].t = PoverS * t;
00448         treg[1].dtdr = PoverS * dtdr;
00449         treg[1].dtdz = PoverS * dtdz;
00450         if(nLay < 2)
00451                 return 2;
00452         nLay--;
00453         t = t_lay(r, z, &dtdr, &dtdz);
00454         nLay++;
00455         if(t < treg[0].t + 0.1)
00456                 return 2;
00457         treg[0].phase = 2;
00458         treg[1].phase = 3;
00459         treg[2].phase = 4;
00460         treg[2].t = t;
00461         treg[2].dtdr = dtdr;
00462         treg[2].dtdz = dtdz;
00463         treg[3].phase = 5;
00464         treg[3].t = PoverS * t;
00465         treg[3].dtdr = PoverS * dtdr;
00466         treg[3].dtdz = PoverS * dtdz;
00467         return 4;
00468 }

Generated on Tue May 6 09:16:11 2003 for Earthworm Libs by doxygen1.3-rc3