raysimp.m

这是matlab在地球物理数据处理方面的源码

function [traceoff,timebase,csg] = raysimp(h,v)

% Program to calculate using simple Snell's law ray tracing through a stack of 
% n flat-lying layers defined by their velocity and thickness 
% The program also plots the expected traces and supeimposes on this the 
% traveltimes expected for the model using the short offset approximation.
%
% To run    [traceoff,timebase,csg] = raysimp(h,v);
%
% h = thickness of n layers in metres
% v = corresponding velocity of the layers in m/s.
% csg = resulting output set of seismograms observed, the 'pulse' is a hanning window
% traceoff is a vector containing the corresponding trace offsets.
% timebase is a vector which is the corresponding traces sample times.
% Current implementation set to 2 ms sampling, 1 second record, offsets @ 10 m spacing
% from 0 to 1 km.

numlayers = length(h);

maxcritical = (asin(v(1)*v.^(-1)));  % Find the maximum critical takeoff angle for each layer.
maximum_critical_angle = maxcritical*180/pi

maxtheta = 80*pi/180; deltheta = 1*pi/180;
maxoffset = 1000; delx = 10; delt = 0.002; maxtime = 1; 
traceoff = 0:delx:maxoffset;
timebase = 0:delt:maxtime;
csg = zeros(length(timebase),length(traceoff));  % Commons Shot Gather.
[numsamps,numtraces] = size(csg);

[takeoff,H] = meshgrid([0:deltheta:maxtheta],h);  % Assign takeoff angles
[takeoff,V] = meshgrid([0:deltheta:maxtheta],v); 
for i = 1:numlayers
   maxcritical(i) = min([maxcritical(i) maxtheta])
end


[m,n] = size(takeoff);

theta = asin(V.*sin(takeoff)/v(1));  % Calculate angles within each layer., will go complex
% imagesc(theta), colorbar
x = 2*H.*tan(theta);  % Offset in each layer
offset = cumsum(x);

L = 2*H./cos(theta);
traveltime = cumsum(L./V);

for i = 1:numlayers
   ind(i) = floor(maxcritical(i)/deltheta)-1
   arrtimeind = round(spline(x(i,1:ind(i)),traveltime(i,1:ind(i)),traceoff)/delt);
   
   for k = 1:numtraces
      if arrtimeind(k) < numsamps
         csg(arrtimeind(k),k) = 1;
      end
   end
end

waveletsamps = 21;
for i = 1:numtraces
   temp(:,i) = conv(hanning(waveletsamps),csg(:,i));
end

csg = temp(round(waveletsamps/2):numsamps+round(waveletsamps/2)-1,:);
figure
imagesc(traceoff,timebase,csg), xlabel('Offset metres'), ylabel('Time seconds')
title('Comparison of Simple Snell Law Raytracing vs Short Offset Approx (white lines)')
hold on

twoway = 2*h./v
vrms = sqrt(cumsum(v.^2.*twoway)./cumsum(twoway))
[Vrms,X] = meshgrid(vrms,traceoff);
[Twoway,X] = meshgrid(cumsum(twoway),traceoff);

times = sqrt(Twoway.^2 + X.^2./Vrms.^2);
plot(traceoff,times','w')
hold off