%*************************************************************
% ZO Modeling and Postack Migration Code
% Modeling is exploding reflector (smearing a point reflectivity along
% hyperbola in data space) , migration is smearing data
% along circles in model space. I have killed geometrical spreading.
%
% dx       - input- spatial sampling interval (assume dx=dz)
% c        - input- velocity m/s
% option   - input- =0 syncline model
%                   =1 point scatterer model
% (nx,nz)  - input- grid dimensions of model, trace at each surface pt
%  app     - input- migration aperture
%                   app=nx means full aperture of traces
%                   app=nx/3 means 1/3 of traces are used
% ioption  - input- 0  syncline model
%                   1  6-pt scatterer model
%                   2  1-pt scatterer model
%(refl,MIG)-output- (ZO traces, Poststack Migrated section)
%**************************************************************
clear all
dx=25;dz=dx;c=1500; nx=75;nz=75;ntrace=nx;dt=.5*0.008;time_table=1;
ntime=round(1.2*sqrt(nx^2+nz^2)*dx/c/dt);app=round(nx);option=2;
%
%**************************************************************
% Create Source Wavelet, Reflectivity Model MIG, and Synthetic Data refl
%**************************************************************
%
    rick=ricker(100,15,dt);tshift=1;     
%
%**************************************************************
% Create Reflectivity Syncline Model MIG, and Synthetic Data refl
% Make sure app=nx
%**************************************************************
%
   if option==0;
     MIG=zeros(nx,nz); 
     for i=1:nx/3;MIG(i,round(nz/3))=1;end;
     for i=round(nx*2/3):nx;MIG(i,round(nz/3))=1;end;
     for i=round(nx/3)+1: round(nx*2/3-1);iz=15*sin((i-nx/3+1)*1.9*pi*3/nx/2);
     MIG(i,round(iz+nz/3))=1;end;MOD=MIG;
   end
%
%**************************************************************
% Pt 6-pt Scatterer Model;    
%**************************************************************
   if option==1;
     MIG=zeros(nx,nz);
     for i=2:4;z0=5+(i-1)*20;MIG(round(nx/2)-2,z0)=1; 
     MIG(round(nx/2)+2,z0)=1; end
     MOD=MIG;
   end

%
%**************************************************************
% Pt 1-pt Scatterer Model;    
%**************************************************************
   if option==2;
     MIG=zeros(nx,nz);
     for i=4:4;z0=5+(i-1)*20;MIG(round(nx/2)-2,z0)=1; end
     MOD=MIG;
   end


%
    [cdp1]= model1(MIG,time_table,nx,nx,nz,ntrace,ntime,dt,app,rick,tshift,dx,c);    
    [nx,ntime]=size(cdp1); refl=cdp1;n=nx;
    q=0; migi=zeros(nx,nz); 
%
%**************************************************************
% Migrate Data
% Set app= # of traces you want migrated about image point
%**************************************************************
%
app=round(nx/2);
refl1=[diff(refl')' zeros(nx,1)];
refl2=refl1;subsample=7;
for ii=1:subsample:nx;refl2(ii,:)=0;end;
%refl1=refl2-refl1;
[MIG]=migrate1(refl1,time_table,nx,nx,nz,ntrace,ntime,dt,app,rick,tshift,dx,c,MOD);
%
%**************************************************************
% Plot Migration Section
%**************************************************************
%
   plotit(MIG,MOD,nx,nz,dx,dz,dt,refl1,ntime)
print -dpsc fig15a.ps

figure(2)
surf([1:nx]*dx,[1:nx]*dx,MIG)
xlabel('Z (m)','Fontsize',22)
zlabel('Migration Amplitude','Fontsize',22)
ylabel('X (m)','Fontsize',22)
title('Side View of Migration Section','Fontsize',22)
print -dpsc fig15b.ps


