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 }