kirk_mig2.m

Univ. of Calgary CREWS的免费地震研究软件

if itaper2 ~= 1 & itaper2 ~= 0
	error('the angle limit taper type: params(6) should be 0 and 1 !');
end
if( isnan(params(7)) )
	interp_type = 1;
else
	interp_type = params(7);
end
if interp_type < 1 | interp_type > 4
	error('the interpolation indexx paarams(7) should be 1, 2, 3 and 4 !');
end
if( isnan(params(8)) ) 
		tmig1 = min(t);
else
		tmig1 = params(8);
end
if( isnan(params(9)) ) 
		tmig2 = max(t);
else
		tmig2 = params(9);
end
if tmig2 < tmig1
	error('the target time window start time should be smaller than the end time !');
	error('i.e. paraams(8) < params(9)');
end
if( isnan(params(10)) ) 
		xmig1 = min(x);
else
		xmig1 = params(10);
end
if( isnan(params(11)) ) 
		xmig2 = max(x);
else
		xmig2 = params(11);
end
if xmig2 < xmig1
 	error('the start location of target trace range should be prerior to the end location');
	error('i.e. params(10) < params(11)');
end
if( isnan(params(12)) )
		ibcfilter = 0;
else
		ibcfilter = params(12);
end
if ibcfilter 
	% get a cumulative array from aryin
	arycum=cumsum(aryin);
end
	
%aperture in traces and the taper coefficient
traper0 = .5*aper/dx;
traper1 = width1/dx;
traper = round(traper0+traper1);
if itaper1 == 0
	coef1 = lin_taper(traper0,traper0+traper1);
else
	coef1 = cos_taper(traper0,traper0+traper1);
end
%one way time
dt1=.5*dt;
t1=t/2;
t2= t1.^2;
%compute maximum time needed
vmin=min(min(aryvel));
vmax=max(max(aryvel));
tmax=sqrt( .25*tmig2^2 + ((.5*aper+width1)/vmin)^2);
%pad input to tmaxin
npad=ceil(tmax/dt1)-nsamp+5;
if( npad > 0)
	aryin= [aryin; zeros(npad,ntr)];
	t1 = [t1',(nsamp+1:nsamp+npad)*dt1]';
	if ibcfilter
		for j=1:npad
			arycum=[arycum; arycum(nsamp,:)];
		end
	end
end
%output samples targeted
samptarget=near(t,tmig1,tmig2);
tmig=t(samptarget);
%output traces desired
trtarget= near(x,xmig1,xmig2);
xmig=x(trtarget);
%initialize output array
arymig=zeros(length(samptarget),length(trtarget));
%loop over migrated traces
%
kmig=0;
disp([' ']);
disp([' --- Total number of traces to be migrated : ' int2str(length(trtarget)) ' ---']);
disp([' ']);
%flops1=flops;
clock1=clock;
for ktr=trtarget 			% ktr--the location of output trace
	kmig=kmig+1; 			% numerator
	%determine traces in aperture
	n1=max([1 ktr-traper]);
	n2=min([ntr ktr+traper]);
	truse=n1:n2;
	%offsets and velocity 
	offset2=((truse-ktr)*dx).^2;
	v2 = aryvel(:,ktr).^2;
	
	% loop over traces in aperture
	disp([' --- Trace No.' ,int2str(kmig) ,' The traces in aperture : ' ,int2str(length(truse))]);
	for kaper=1:length(truse)
	
		% offset times
		t_aper = sqrt( offset2(kaper)./v2(samptarget) + t2(samptarget) );
		%cosine theta amplitude correction
		if truse(kaper) == ktr
			costheta = ones(size(samptarget'));
			tanalpha = zeros(size(samptarget'));
		else
			costheta = 0.5*tmig./t_aper;
			tanalpha = sqrt(1-costheta.^2);
					
			%angle limit and the taper
		
			ind = find( costheta < cos(angle1) );
			i1 = ind(length(ind));
			ind = find( costheta < cos(ang_limit) );
			i2 = ind(length(ind));
		
			if i1 < i2
			        if itaper2  ==  0
					coef2 = lin_taper(i2,i1);
				else
					coef2 = cos_taper(i2,i1);
				end
				costheta(1:i1) = zeros(i1,1);
				costheta(i1+1:i2) = coef2(i2-i1:-1:1)'.*costheta(i1+1:i2);
			end
		end
	
		% boxcar anti-aliasing filter
		if ibcfilter
			lt0=(dx*tanalpha./aryvel(samptarget,ktr)/dt1);
			indt = (t_aper/dt1)+1;
			lentr = nsamp+npad;
			lt = ones(lentr,1)*max(lt0);
			lt(indt)=lt0;
			lt(max(indt)+1:lentr) = ones(lentr-max(indt),1)*min(lt0);
			it = (1:lentr)';
			l1=it-lt-1;
			l2=it+lt;
			ind = find(l1 < 1);
			l1(ind) = ones(length(ind),1);
			ind = find(l2> lentr);
			l2(ind)=ones(length(ind),1)*lentr;
			tmp0=t1;
			tmp0(1) = arycum(1,truse(kaper));
			ind = 2:lentr;
			tmp0(ind) = (arycum(l2(ind),truse(kaper))-arycum(l1(ind),truse(kaper)))...
				    ./(l2(ind)-l1(ind));
		else
			tmp0 = aryin(:,truse(kaper));
		end	
		%interpolation
		% Linear
		if interp_type == 1
			tnumber = t_aper/dt1;
			it0 = floor( tnumber ) + 1;
			it1 = it0+1; 
			xt0 = tnumber - it0 + 1;
			xt1 = it0-tnumber;
			tmp = xt1.*tmp0(it0)+xt0.*tmp0(it1);
		end
		% Spline
		if interp_type == 2
			tmp = interp1(t1,tmp0,t_aper,'spline');
		end
		% Cubic
		if interp_type == 3
			tmp = interp1(t1,tmp0,t_aper,'cubic');
		end
		% Sinc
		if interp_type == 4
			tmp = sinci(tmp0,t1,t_aper);
		end
		% aperture taper
		ccoef = 1. ;
		if abs(truse(kaper)-ktr)*dx > 0.5*aper
			ccoef = coef1( round(abs(truse(kaper)-ktr)-traper0) );
		end
		if abs(1-ccoef) > 0.05
			tmp = tmp .* ccoef;
		end
		
		ind = find( costheta < 0.999);
		costheta(ind) = sqrt(costheta(ind).^3);
		tmp(ind) = tmp(ind) .* costheta(ind);
		arymig(:,kmig)= arymig(:,kmig)+tmp;
		
	end
	
	% scaling and 45 degree phase shift
	scalemig = aryvel(samptarget,kmig).*sqrt(pi.*(tmig+0.0001)) ;
	arymig(:,kmig) = arymig(:,kmig)./scalemig ;
end
% 45 degree phase shift
arymig = conv45(arymig);
%flops2=flops;
runtime=etime(clock,clock1);
%disp(['Total floating operation2 --' int2str(flops2-flops1)]);
disp(['Total run time --' int2str(runtime)]);
>>>>>>> 1.4