基于matlab二维射线追踪程序地震声波正演源程序

function [t,p,L]=traceray_ps(vp,zp,vs,zs,zsrc,zr,zd,x,xcap,pfan,itermax,optflag,pflag,dflag) % TRACERAY_PS: traces a P-S (or S-P) reflection for v(z) % % [t,p,L]=traceray_ps(vp,zp,vs,zs,zsrc,zr,zd,x,xcap,pfan,itermax,optflag,pflag,dflag) % % TRACERAY_PS uses the modified bisection algorithm % to trace a p-s reflection % in a stratified medium. % % vp,zp = velocity and depth vectors giving the p wave model. Depths are % considered to be the tops of homogeneous layers. % vs,zs = velocity and depth vectors giving the s wave model. % for both models, depths are considered to be the tops of homogeneous layers. Thus, % the first velocity applies from the first depth to the second. A constant velocity % would be specified like vp=1500;zp=0; etc. % zsrc = depth of the source (scalar) % zr = depth of receiver (scalar) % zd = depth of the reflection (scalar) % x = vector of desired source-receiver offsets (horizontal) % MUST be non-negative numbers % xcap = (scalar) capture radius % pfan = vector of ray parameters to use for the initial fan of rays % By default, the routine generates the fan using straight rays. % Setting pfan to -1 or omitting it activates default behavior. % Setting pfan to -2 causes the routine to use the pfan used in the % last call to this program. (pfan of -1 and -2 are identical on the % first call.) % itermax = maximum number of iterations allowed for ray capture % =========================== Default = 4 ================================ % optflag = if nonzero then refine captured rays by linear interpolation % =========================== Default = 1 ================================ % pflag = if nonzero, then print information about all failed rays % =========================== Default = 0 ================================ % dflag = 0 no action % = 1 then a new figure window will be opened and the raypaths plotted % = 2 raypaths will be plotted in the current window % =========================== Default = 0 ================================ % NOTE: All z variables must be specified relative to a common datum % % t = vector of traveltimes for each x % p = vector of ray parameters for each x % L = vector of geometrical spreading factors for each x % (if L is not asked for, it will not be calculated) % % G.F. Margrave, CREWES Project, June 1995 % % NOTE: It is illegal for you to use this software for a purpose other % than non-profit education or research UNLESS you are employed by a CREWES % Project sponsor. By using this software, you are agreeing to the terms % detailed in this software's Matlab source file. % BEGIN TERMS OF USE LICENSE % % This SOFTWARE is maintained by the CREWES Project at the Department % of Geology and Geophysics of the University of Calgary, Calgary, % Alberta, Canada. The copyright and ownership is jointly held by % its author (identified above) and the CREWES Project. The CREWES % project may be contacted via email at: crewesinfo@crewes.org % % The term 'SOFTWARE' refers to the Matlab source code, translations to % any other computer language, or object code % % Terms of use of this SOFTWARE % % 1) Use of this SOFTWARE by any for-profit commercial organization is % expressly forbidden unless said organization is a CREWES Project % Sponsor. % % 2) A CREWES Project sponsor may use this SOFTWARE under the terms of the % CREWES Project Sponsorship agreement. % % 3) A student or employee of a non-profit educational institution may % use this SOFTWARE subject to the following terms and conditions: % - this SOFTWARE is for teaching or research purposes only. % - this SOFTWARE may be distributed to other students or researchers % provided that these license terms are included. % - reselling the SOFTWARE, or including it or any portion of it, in any % software that will be resold is expressly forbidden. % - transfering the SOFTWARE in any form to a commercial firm or any % other for-profit organization is expressly forbidden. % % END TERMS OF USE LICENSE if(nargin<14) dflag=0; end if(nargin<13) pflag=0; end if(nargin<12) optflag=1; end if(nargin<11) itermax=4; end if( nargin<10) pfan=-1; end if(~prod(double(size(vp)==size(zp)))) error('vp and zp must be the same size'); end if(~prod(double(size(vs)==size(zs)))) error('vp and zp must be the same size'); end if(length(zsrc)~=1 | length(zr)~=1 | length(zd)~=1 | length(xcap)~=1 ) error(' zsrc,zr,zd and xcap must be scalars ') end ind=find(x<0); if(~isempty(ind)); error('offsets must be nonnegative'); end %make sure zsrc < Zd and zr<zd if(zsrc > zd ) error(' zsrc must be less than zd'); end if(zr > zd ) error(' zr must be less than zd'); end %adjust zsrc,zr,and zd and z2 so that they won't be exactly on layer boundaries zsrc=zsrc+100000*eps; zr=zr+100000*eps; zd=zd-100000*eps; if( zsrc < zp(1) | zsrc < zs(1)) error(' source depth outside model range'); elseif(zr < zp(1) | zr < zs(1) ) error('receiver depth outside model range'); end %determine the layers we propagate through %down leg ind=find(zp>zsrc); if(isempty(ind)) ibegd=length(zp); else ibegd=ind(1)-1; end ind=find(zp>zd); if(isempty(ind)) iendd=length(zp); else iendd=ind(1)-1; end %test for sensible layer indices if(iendd<1 | iendd > length(zp) | ibegd <1 | iendd > length(zp)) %somethings wrong. return nans t=inf*ones(size(x)); p=nan*ones(size(x)); return; end %create v and z arrays for down leg vp=vp(:);zp=zp(:); vpd=[vp(ibegd:iendd);vp(iendd)];%last v is irrelevent. zpd=[zsrc;zp(ibegd+1:iendd);zd]; %up leg ind=find(zs>zr); if(isempty(ind)) ibegu=length(zs); else ibegu=ind(1)-1; end ind=find(zs>zd); if(isempty(ind)) iendu=length(zs); else iendu=ind(1)-1; end %test for sensible layer indices if(iendu<1 | iendu > length(zs) | ibegu <1 | iendu > length(zs)) %somethings wrong. return nans t=inf*ones(size(x)); p=nan*ones(size(x)); return; end %create v and z arrays for up leg vs=vs(:);zs=zs(:); vsu=[vs(ibegu:iendu);vs(iendu)];%last v is irrelevent. zsu=[zr;zs(ibegu+1:iendu);zd]; % combine into one model. We shoot oneway rays through an % equivalent model consiting of the down-leg model followed by % the upleg model (upside down). vmod=[vpd(1:end-1);flipud(vsu(1:end-1))]; tmp=flipud(zsu); zmod=[zpd;cumsum(abs(diff(tmp)))+zpd(end)]; z1=zsrc;%start depth z2=max(zmod)+100000*eps;%end depth %determine pmax vmax=max([vmod]); nearlyone=.999; pmax=nearlyone/vmax; %the meaning of pmax is that it is the maximum ray parameter which will not %be critically refracted anywhere in vp or vs. if(pfan==-2) global PFAN if(isempty(PFAN)) pfan=-1; else pfan=PFAN; PFAN=-2; end else PFAN=0; end if(pfan==-1) %shoot a fan of rays Iterate once if the fan does not span the xrange %guess angle range nrays=max([3*length(x),10]); pfan=linspace(0,1/max(vmod),nrays); % xmaxf=max(x); % xminf=min(x); % if(xminf~=xmaxf) % nrays=3*length(x); % xfan=linspace(xminf,xmaxf,nrays); % %th=atan(xfan/(z2-z1)); % th = pi*linspace(0,90,nrays)/180; % %determine the fan ray parameters % % pfan=sin(th)/vmod(1); % % %put in a max and min p % if(min(pfan)>0) % pfan=[0 pfan]; % xfan=[0 xfan]; % nrays=nrays+1; % end % % if(max(pfan)>pmax) % ind=find(pfan<pmax); % pfan=pfan(ind); % xfan=xfan(ind); % nrays=length(ind); % end % % if(max(pfan)<pmax) % pfan=[pfan pmax]; % xfan=[xfan (z2-z1)*nearlyone/sqrt(1-nearlyone*nearlyone)]; % nrays=nrays+1; % end % % xminf=2.0*min(xfan); % xmaxf=2.0*max(xfan); % else % %3 rays 1