%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2-4 FD Soln to ACOUSTIC WAVE EQUATION (cnst density)
% (nx,nz,nt)      - input- (Horizontal,Vertical) gridpt dimens. of vel 
%                            model & # Time Steps
% FRE             - input- Peak frequency of Ricker wavelet
% BVEL            - input- NXxNZ matrix of background velocity model
% (dx,dt)         - input- (space, time) sample intervals
% (xs,zs)         - input- (x,z) coordinates of line source
% RICKER(nt)      - input- nt vector of source time histories
% (p2,p1,p0)      -calcul- (future,present,past) NXxNZ matrices of 
%                           modeled pressure field
% (p0,p1)         -output- Old and present pressure panels at time nt.
% REALDATA(nx,nt) -output- CSG seismograms at z=2
% Han, July. 10, 2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
c   = 4.0; %km/s
FRE = 20; 
nx  = 300; 
nz  = nx/2; % Define model size, src. freq., velocity
dx  = c/FRE/20;  % Set spatial sampling dx
dz  = dx;
nt  = 800; % Define pt src coord.

dt  = .6*(sqrt(dx*dx+dz*dz)/2.0)/c;   % Set dt for stability condition
nxm = nx-1;
nzm = nz-1; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%show the Ricker wavelet 
t       =   [0:1:nt-1]*dt - 0.95/FRE;
RICKER  =   zeros(length(t));
RICKER  =   (1 - t.*t * FRE^2 *pi^2  ) .*exp(- t.^2 * pi^2 * FRE^2 );

figure;
plot  ( [0:nt-1]*dt, RICKER );
title ( 'Ricker Wavelet' );
xlabel( 'Time (s)' );
ylabel( 'Amplitude' );
pause ( 1.5 );
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%constructing the velocity model
BVEL  =  ones(nz,nx)*c;
BVEL(nz-round(nz/2):nz, :) = BVEL(nz-round(nz/2):nz, :)*1.4; % 2-Layer velocity model
% velocity in the 2nd layer is 1.4 time of the first layer above it
REALDATA  = zeros(nx,nt); % Define shot gather matrix
p0 = zeros(nz,nx);
p1 = p0; 
p2 = p0;
cns = (dt/dx*BVEL).^2; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%LOOP OVER TIME STEPS, forward modeling
zs   =  2;                 % source depth, zs=1 is free surface, so zs=2
xs   =  round(nx/2.0);     % source horizontal 
zr   =  2;                 % Receiver Depth


tic;                   % A pair of tic and toc is for timing.
for i = 12:12          % this for loop can be used for a single shot
%parfor i = 1:24       % this is for parallel, one core calc one shot 
    xs = (i-1)*(nx/25) + (nx/25);  % reset source horizontal position
    p0 = zeros([nz nx]);
    p1 = p0;
    p2 = p0; 
    p0(zs, xs) = RICKER(1); 
    p1(zs, xs) = RICKER(2);    
    
    RICKER   = (1 - t.*t * FRE^2 *pi^2  ) .*exp(- t.^2 * pi^2 * FRE^2 ); 
    REALDATA = zeros(nx,nt);
    REALDATA(:,1) = p0(zr,:)';
    REALDATA(:,2) = p1(zr,:)';
    laplace_p1 = zeros([nz nx]);
    alpha = dt/dx;
    BVEL  = ones(nz,nx)*c;
    BVEL(nz-round(nz/2):nz, :) = BVEL(nz-round(nz/2):nz, :)*1.4;        
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for it = 3:1:nt                   
        laplace_p1 = Laplacian(p1,nz,nx,p0,alpha,BVEL); %4th order in space
        p2 = 2*p1 - p0 + cns.*laplace_p1; %p2  = 2*p1 - p0 + 4*cns.*del2(p1); % 2nd order FD of Acoustic WE
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%% A POINT SOURCE               
        p2(zs,xs) = p2(zs,xs) + RICKER(it); % Add bodypoint src term, adding pulses continuously          
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %  ABC BOUNDARY CONDITIONS
        % if p2(1, :) == 0, then it is a free boundary surface, row = 1,         
        %p2(1, :) = p1(1,:) - alpha*BVEL(1,:).*(p1(1,:) - p1(2,:));          % Top ABC, just perpendicular to boundary       
        p2(1, :) = 0;                                                        % free surface
        p2(nz, :) =  p1(nz, :) - alpha*BVEL(nz, :).*(p1(nz,:)  - p1(nzm,:)); % Bottom ABC
        %p2(nz, :) =  p1(nz, :);                                             % Bottom reflecting
        p2(:, nx) =  p1(:, nx) - alpha*BVEL(:, nx).*(p1(:,nx) - p1(:,nxm));  % RHS ABC                
        p2(:, 1) =  p1(:, 1) + alpha*BVEL(:, 1).*(p1(:, 2) - p1(:, 1));      % LHS ABC
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%             
        REALDATA(:, it) = (p2(zr, :))'; % nx receivers at the surface
        if ( round(it/10)*10 == it )    % update the picture every 10 iterations             
            plt;pause(0.01);
        end % PLOT OUT RESULTS        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        p0 = p1; 
        p1 = p2;  % UPDATE PRESSURE FIELDS
        wavefield_path = strcat('./fd_data/shot', num2str(xs), 'iter',num2str(it), '.dat');
        save_survey(wavefield_path,p2); %cannot save file explicitly in parfor           
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    end
    obs_seismo_path = strcat('./fd_data/seismo_shot', num2str(xs), '.dat');
    save_survey(obs_seismo_path, REALDATA);
end
toc;

