function [seismigd,seismigx,zmig,migtime,ttables,seisout,shotout,tout]=migfocishot(seis,x,t,v,xv,dz,xshot,zshot,wlet,tw,fmin,fmax,nfor,ninv,nop,pow,params)
% migfocishot ... depth migration of a single 2D shot record by the FOCI method
%
% [seismigd,seismigx,zmig,migtime,ttables,seisout,shotout,tout]=migfocishot(seis,x,t,v,xv,dz,xshot,zshot,wavelet,tw,fmin,fmax,nfor,ninv,nop,pow,params)
%
% Shot record migration by the FOCI method with spatial resampling to allow shorter
% operator lengths. The algorithm includes a used specified temporal pad
% but applies no spatial (i.e. x) pad. The user must do this extermally.
% Typical pad might be zero trace spans on to the left and right of the
% shot whose size is about 1/2 the shot, thus doubling the shot size.
%
% seis ... input shot record, regularly sampled in x and t
% x ... x coordinates of the columns of seis
% t ... t coordintes of the rows of seis
% v ... velocity matrix. One row for each depth step to be taken.
% xv ... x coordinates for v
% NOTE: seis and v must be assigned x oordinates with the same sample
% interval and origin. The velocity matrix must span at least as large an
% x range as seis. The vertical coordinate of v is assumed to be (k-1)*dz where
% k is row number.
% dz ... depth step size
% xshot ... x coordinate of the shot (may be an array)
% zshot ... z coordinate of shot (same size array as xshot) Should be < dz for
% best results. NOT FUNCTIONAL.
% wlet ... source waveform
% tw ... time coordinate vector for wlet.
% fmin ... minimum frequency (Hz) to migrate
% fmax ... maximum frequency (Hz) to migrate
% nfor ... the forward operator length (in samples)
% ninv ... the inverse operator length (in samples)
% nop ... the final operator length, 0 gets nfor+ninv-1 but no windowing.
% if nop<nfor+ninv-1 then the operator is designed as nfor+ninv-1
% in length and then windowed (nicely) to length nop.
% pow ... vector of length 1 or 2 giving the "pow" values (see
% extrap_design) for the first and possibly second operator
% table. The first operator table is invoked most of the time while
% the second gets invoked as controlled by params(2). This can be
% used to do evanescent filtering only every so often. Setting
% pow=[.01 1] will do this. Setting pow equal to a single number
% causes only a single operator to get built.
% NOTE: The params vector controls a lot of options. To get the default
% behavior, either do not provide params on input, or provide it as a vector
% of nan's. For example, setting params=nan(ones(1,10)) will get the
% default behavior. Setting params=nan*ones(1,10); params(1)=.5; gets
% default behavior in all but the first parameter.
% params(1) ... size in seconds of tpad
% ********* default max(t) **************
% params(2) ... operator table 2 is invoked every this many steps
% ********** default 10 ****************
% params(3) ... rezero the temporal pad every this many steps
% NOT PRESENTLY FUNCTIONAL
% ********* default 10 **************
% params(4) ... eta, a new frequency chunk is defined whenever the
% evanescent wavenumber falls below eta*nyquist
% ********* default .5 ********
% params(5) ... beta, a new frequency chunk is resampled such that its highest
% evanescent wavenumber is beta*nyquist
% ********* default .9 *********
% params(6) ... stab, the stability constant for the deconvolution imaging
% condition
% ********* default .01 ********
% params(7) ... makeTable, if 0, the operator tables from the previous call
% to this function (if any) will be used. If 1, then new tables
% are made regardless. If -1, then the existing tables are used
% (if any) but cleared at the function termination. NOTE: If
% existing tables are used, no effort is made to ensure that
% these tables are consistent with the input parameters. Tables
% should be remade whenever dx,dz,v,fmin,fmax,nfor,ninv,nop,or pow are
% changed.
% ********* default 0 *********
% params(8) ... vcrit. This is the critical velocity that governs spatial
% resampling. If 0 or nan, then the critical velocity is computed as the
% minimum value in the velocity model. Programming this ensures
% that all shot records get the same resampling.
% ********** default 0 *********
% params(9) ... source modelling options. If 1, then the source is computed
% as the theoretically correct 2D Greens function with spatial
% antialias filtering. if 2, then similar to 1 except that there
% is not antialias filtering. If 3, then the source is a simple
% impulse (delta function. Only 1 is theoretically correct, the
% others are errors. This is included to demonstrate these common
% errors.
% ********** default 1 *********
% params(10) ... specifies depth at which the forward modelled shot and the
% downward continued data are desired to be output for analysis.
% A value of nan means no additional output is created.
% If programmed as the negative of the desired depth (e.g. -1000)
% then the migration is terminated at this depth and not
% completed. If programmed, then abs(params(10)) must be >= 2*dz
% ********** default nan **********
%
% seismigd ... migrated shot imaged with the deconvolution imaging
% condition
% seismigx ... migrated shot imaged with the crosscorrelation imaging
% condition
% zmig ... depth coordinate for the above two products
% NOTE: The migrated records will be almost the same size and coordinates as v
% except that the depth coordinate will have one extra sample.
% migtime ... computation time required (excluding time to build tables)
% ttables ... time required to build tables
% seisout ... output downward continued data at the depth specified by
% params(10)
% shotout ... output forward modelled shot at the depth specified by
% params(10)
% tout ... time coordinate vector for seisout and shotout (includes pad)
%
% G.F. Margrave, CREWES/POTSI 2004
%
% NOTE: This SOFTWARE may be used by any individual or corporation for any purpose
% with the exception of re-selling or re-distributing the SOFTWARE.
% 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: [email protected]
%
% 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) This SOFTWARE may be used by any individual or corporation for any purpose
% with the exception of re-selling or re-distributing the SOFTWARE.
%
% 2) The AUTHOR and CREWES must be acknowledged in any resulting publications or
% presentations
%
% 3) This SOFTWARE is provided "as is" with no warranty of any kind
% either expressed or implied. CREWES makes no warranties or representation
% as to its accuracy, completeness, or fitness for any purpose. CREWES
% is under no obligation to provide support of any kind for this SOFTWARE.
%
% 4) CREWES periodically adds, changes, improves or updates this SOFTWARE without
% notice. New versions will be made available at www.crewes.org .
%
% 5) Use this SOFTWARE at your own risk.
%
% END TERMS OF USE LICENSE
if((min(xv)>min(x))||(max(xv)<max(x)))
error('velocity model must be at leas