%%My attempt to do FDTD with a staggered grid

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

il=[1:IL];            im=[I+1:IR];        it=[IR+1,IPLM];    it=[il,im,it];
ih=[2:IPLM];
il=[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/(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(-
e3=ones(1,IPLM);                     e3(iplm)=exp(-gamma.*h);
e4=(h/(delta*pho)).*ones(1,IPLM);    e4(iplm)=(1-exp(-gamma.*h))./(delta.*pho.*gamma);


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

for t=1:nsteps 
   
    p(ih)=e1(ih).*p(ih)-e2(ih).*(u(ih)-u(il));
    p1(ih)=e1(ih).*p1(ih)-e2(ih).*(u1(ih)-u1(il));
    u(il)=e3(il).*u(il)-e4(il).*(p(ih)-p(il));
    u1(il)=e3(il).*u1(il)-e4(il).*(p1(ih)-p1(il));
    %p1(1)=10*exp(-(t-30.)*(t-30)/100);
    tau=t*fmax;
    %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(1)=p1source(t);
    p1(IPLM)=0;
    %plot(it,p,it,u,it,p1);
    plot(it,p1);
    set(gca,'ylim',[-20,20]);
    pause(.005);
    
end