function [shotmig,tmig,xmig,ymig,xgrid,ygrid]=kirk_shot3d(shotrec,t,x,y,xshot,yshot,velmod,tv,xv,yv,params)
% KIRK_SHOT3D: Kirchhoff prestack shot-record time migration in 3D   
% 
% [shotmig,tmig,xmig,ymig,xgrid,ygrid]=kirk_shot3d(shotrec,t,x,y,xshot,yshot,velmod,tv,xv,yv,params)
%
% KIRK_SHOT3D migrates a single 3D shot record using prestack time
% migration. The 3D shot record is stored in a 2D matrix and it is not
% necessary that the traces fall on a regular spatial grid.
% The algorithm is simple traveltime path summation with a few
% options as described in the parameters vector. Traveltime from source to
% scatterpoint (i.e the image point) is approximated by a Dix equation
% using the rms velocity from the model at the lateral position halfway
% between source and receiver and at the vertical traveltime of the 
% scatterpoint. Similarly, from the scatterpoint to a receiver, a Dix 
% equation using the rms velocity at halfway between scatterpoint and 
% receiver is used. There is no topographic compensation at present. The
% source and all receivers are assumed to be on the same horizontal plane.
% The Kirchhoff integral weighting is done as described in Bleistein, Cohen,
% and Stockwell p247 eqn 5.2.22. Although this derivation is for constant
% velocity, the adjustments for time migration are reasonably obvious. The
% algorithm works entirely in the time domain including the sample
% interpolation.
%
% shotrec ... 2D matrix containing the shot record. One trace per column.
% t ... time coordinate vector for shotrec. 
%           Requirement: length(t)=size(shotrec,1);
% x ... row vector of the first space coordinate for shotrec. Let
%   [nt,ntr]=size(shotrec), then we must have ntr=length(x);
%           Requirement: length(x)=size(shotrec,2);
% y ... row vector of the first space coordinate for shotrec. Let
%   [nt,ntr]=size(shotrec), then we must have ntr=length(y);
%           Requirement: length(y)=size(shotrec,2)
% velmod ... 3D velocity model. This is a 3D matrix of RMS velocities as a
%           function of lateral position and time. For a constant velocity
%           migration this may be input as a single scalar. For a v(t)
%           migration, this may be input as a column vector (time series). 
%           In these latter two cases, it is internally expanded into a 3D volume.
% tv ... time coordinate vector for velmod
%        ignored if velmod input as scalar.
%           Requirement: length(tv)=size(velmod,1);
% xv ... space coordinate vector for velmod
%        ignored if velmod input as scalar or vector
%           Requirement: length(xv)=size(velmod,2);
% yv ... space coordinate vector for velmod
%        ignored if velmod input as scalar or vector
%           Requirement: length(yv)=size(velmod,3);
% REQUIREMENT: the span of tv, xv and yv must equal or exceed that of t, x and y.
% NOTE: While the data may be spatially irregular, the velocity model must
% be regular.
% NOTE: when a velocity function is required at a location not directly
% given in velmod, nearest-neighbor interpolation is used.
% xshot ... x coordinate of the shot. It must be in the same coordinate
%       system as vector x. 
% yshot ... y coordinate of the shot. It must be in the same coordinate
%       system as vector y. 
%
% params ... vector of migration parameters. An entry of nan gets the
% default.
%
%    params(1--3) : migration aperture and its taper
%    params(1) ... physical aperture in meters. This is the largest lateral
%               distance a trace may have from the output location and still be
%               allowed to contribute.
%	          default is the largest diagonal of the velocity model
%    params(2) ... width of the aperture taper
%                 default is 0.05*params(1)
%    params(3) ... = 0, linear taper
%                     = 1, cosine taper   
%                 default is 1 (cosine taper)
%    params(4-5) : scattering angle limit (degrees)
%    params(4) ... maximum scattering angle limit (degrees)
%                 default = 60
%    params(5) ... width of angle limit taper (linear in cosine)
%                 default = min([0.15*params(4), 90-params(4)])
% NOTE: The angle limit described by params(4-5) is applied to both the
% source angle (angle w.r.t vertical of the vector from imagepoint to source)
% and the receiver angle (similarly defined). If either angle exceeds this
% limit, then the response is either damped or attenuated altogether.
%    params(6) ... not used
%    params(7) : determine how seismic amplitudes are interpolated
%    params(7) ... = 1, linear interpolation
%                     = 2, cubic interpolation
%                     = 3, spline interpolation
%                     = 4, sinc interpolation
%                 default = 1
%    params(8--11) : defines migration target window
%	 params(8) ... tmin of migration target window
%	          default = min(tv) 
%	 params(9) ... tmax of migration target window
%	          default = max(tv)
%    params(10) ... xmin of target window
%	          default = min(xv)
%    params(11) ... xmax of migration target window
%	          default = max(xv) 
%    params(12) ... ymin of target window
%	          default = min(yv)
%    params(13) ... ymax of migration target window
%	          default = max(yv) 
%    params(14) ... dx = spatial grid size of migrated volume. Same in x
%               and y.
%             default = min([mean(diff(x)) mean(diff(y))])
%
%    params(15) : box-car anti-aliasing filter
%    params(15) ... = 0, no box-car filter used;
%                      = 1, box-car filter will be used.
%                 default is 0.
%
% OUTPUT arguments
%
%    shotmig ... 2D matrix of the output migrated shot record
%    tmig ... t coordinates of migrated data (depends on params(6)
%    xmig ... x coordinates of migrated data (length(xmig)=size(shotmig,2))
%    ymig ... y coordinates of migrated data (length(ymig)=size(shotmig,2))
%    xgrid ... the unique x coordinates in the output grid
%    ygrid ... the unique y coordinates in the output grid
%
% By G.F. Margrave
% CREWES Project, U of Calgary, 2007
%
% 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



%tstart=clock;
[nsamp,ntr]=size(shotrec);
[nvsamp,nvx,nvy]=size(velmod);

% check the validity input arguments

if(length(t)~=nsamp)
    error('Incorrect time specification')
end
t=t(:);
dt=t(2)-t(1);

if(length(x)~=ntr)
    error('Incorrect x specification for input data')
end

if(length(y)~=ntr)
    error('Incorrect y specification for input data')
end



%  ---- test velocity info ----
%test for alternate velocity input
if(nvsamp*nvx*nvy==1)
    %scalar velocity input
    velmod=velmod*ones(nt,2,2);
    xv=[min(x) max(v)];
    yv=[min(y) max(y)];
    nvx=2;nvy=2;
    tv=t;
end
if(nvsamp~=1 && nvx*nvy==1)
    %v(t) input
    tmp=velmod;
    velmod=zeros(length(tmp),2,2);
    velmod(:,1,1)=tmp;
    velmod(:,1,2)=tmp;
    velmod(:,2,1)=tmp;
    velmod(:,2,2)=tmp;
    xv=[min(x) max(x)];
    yv=[min(y) max(y)];
    nvx=2;nvy=2;
end
        
if(length(tv)~=nvsamp)
    error('Time vector for velocity model is incorrect')
end
if(length(xv)~=nvx)
    error('Space coordinate vector xv for velocity model is incorrect')
end
if(length(yv)~=nvy)
    error('Space coordinate vector yv for velocity model is incorrect')
end
if(min(tv)>min(t) || max(tv)<max(t))
    error('Time vector for velocity model must span that for data')
end
if(min(xv)>min(x) || max(xv)<max(x))
    error('Space vector xv for velocity model must span vector x for data')
end
if(min(yv)>min(y) || max(yv)<max(y))
    error('Space vector yv for velocity model must span vector y for data')
end

%  ---- examine parameters ----
nparams=15; 				% number of defined parameters
					
if(nargin<11) 				% no parameters inputted
	params= nan*ones(1,nparams); 
end 	 	

if(length(params)<nparams) 
	params = [params nan*ones(1,nparams-length(params))];
end

%assign parameter defaults

if( isnan(params(1)) ) 
	aper = sqrt((max(xv)-min(xv))^2+(max(yv)-min(yv))^2);
else
    aper = params(1);
end

if( isnan(params(2)) )
    width1 = aper/20;
else
    width1 = params(2);
end

if( isnan(params(3)) )
    itaper1 = 1;
else
    itaper1 = params(3);
end

if( isnan(params(4)) )
	ang_limit = pi/3;
else
	ang_limit = params(4)*pi/180;
end

if( isnan(params(5)) )
	width2 = 0.15*ang_limit;
    if(width2>pi/2-ang_limit)
        width2=pi/2-ang_limit;
    end
else
	width2 = params(5)*pi/180;
end
anglemax = ang_limit + width2;


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(tv);
else
		tmig1 = params(8);
end

if( isnan(params(9)) ) 
		tmig2 = max(tv);
else
		tmig2 = params(9);
end
if tmig2 < tmig1
	error(['the target time window start time should be smaller than the end time !'...
            'i.e. paraams(8) < params(9)']);