%%My attempt to do FDTD with a staggered grid

IR=210;
IPLM=220;
IL=10;

il=[1:IL];            im=[IL+1:IR];        ir=[IR+1:IPLM];    it=[il,im,ir];
ihi=[2:IPLM];
ilo=[1:IPLM-1];

%Velocity Field                      Pressure Field
u=zeros(1,IPLM);                     p=zeros(1,IPLM);       %Field for Gaussian Wave
u1=zeros(1,IPLM);                    p1=zeros(1,IPLM);      %Field Blackman window
%Constants and Values
c=1500;                             %Speed                   c=1500m/s
pho=1;                              %Density                 pho=1 gm/cm^3
k=1/(c^2*pho);                      %compressibility         k=4.44*10^-6 N^-1 cm^2
fmax=1e6;                           %Max frequency           fmax=1 Mhz
lambda_min=c/fmax;                  %Minimum wavelength      lambda=
delta=lambda_min/10;                %Grid spacing            m
h=delta/(2*c);                       %Time step               t
alpha_max=-(c*k/(delta))*log(.99);   

%Perfectly absorbing boundary condition
alphar=alpha_max*((ir-IR)./(IPLM-IR)).^2;           gammar=alphar/k;
alphal=alpha_max*((IL+1-il)./(IL)).^2;              gammal=alphal/k;


%Coefficients for artificial attenuation
e1=ones(1,IPLM);                     e1(ir)=exp(-alphar.*h./k);
                                     e1(il)=exp(-alphal.*h./k);
e2=(h/(delta*k)).*ones(1,IPLM);      e2(ir)=(1-exp(-alphar.*h./k))./(delta.*alphar);
                                     e2(il)=(1-exp(-alphal.*h./k))./(delta.*alphal);
e3=ones(1,IPLM);                     e3(ir)=exp(-gammar.*h);
                                     e3(il)=exp(-gammal.*h);
e4=(h/(delta*pho)).*ones(1,IPLM);    e4(ir)=(1-exp(-gammar.*h))./(delta.*pho.*gammar);
                                     e4(il)=(1-exp(-gammal.*h))./(delta.*pho.*gammal);


nsteps=1000;
set(gcf,'doublebuffer','on');
blackmanlength=50;
blackmanwin=blackman(blackmanlength);
p1source=zeros(nsteps,1);
p1source(1:blackmanlength)=20*blackmanwin;

for t=1:nsteps 
   
    %p1(ihi)=e1(ihi).*p1(ihi)-e2(ihi).*(u1(ihi)-u1(ilo));
    %u1(ilo)=e3(ilo).*u1(ilo)-e4(ilo).*(p1(ihi)-p1(ilo));
    p1(ilo)=e1(ilo).*p1(ilo)-e2(ilo).*(u1(ihi)-u1(ilo));
    u1(ihi)=e3(ihi).*u1(ihi)-e4(ihi).*(p1(ihi)-p1(ilo));
    %p1(1)=10*exp(-(t-30.)*(t-30)/100);
    %p1(1)=(fmax/(pi*tau))*cos(tau/pi)*sin(tau/pi)*(.84-4*tau^2/(4*tau^2-1)+...
    %                                            .16*tau^2/(tau^2-1));
    p1(IPLM/2)=p1source(t);
    p1(IPLM)=0;
    p1(1)=0;        
    plot(it,p1);
    set(gca,'ylim',[-20,20]);
    pause(.005);
    
end