%%My attempt to do FDTD with a staggered grid

%Boundary Limits
IPLM=220;           JPLM=220;                               %Total Number of Points
IR=210;             JU=210;                                 %Right Boundary
IL=10;              JB=10;                                  %Left Boundary

[i,j]=meshgrid(1:JPLM,1:IPLM);

%Left-Bottom        Middle-Middle       right-top              total-total
il=i(:,1:IL);       im=i(:,IL+1:IR);    ir=i(:,IR+1:IPLM);     it=[il,im,ir];
jb=j(1:JB,:);       jm=j(JB+1:JU,:);    ju=j(JU+1:IPLM,:);     jt=[jb',jm',ju']';

ihi=2:IPLM;        ilo=1:IPLM-1;
jhi=2:JPLM;        jlo=1:JPLM-1;

%Velocity Field                      %Pressure Filed
ux=zeros(JPLM,IPLM);                 px=zeros(JPLM,IPLM);        
uy=zeros(JPLM,IPLM);                 py=zeros(JPLM,IPLM);
p=zeros(JPLM,IPLM);

%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(.01);   

%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;
alphat=alpha_max*((ju-JU)./(JPLM-JU)).^2;                   gammat=alphat/k;
alphab=alpha_max*((JB+1-jb)./(JB)).^2;                      gammab=alphab/k;


%Coefficients for artificial attenuation
e1=ones(JPLM,IPLM);                         e1(:,[IR+1:IPLM])=exp(-alphar.*h./k);
                                            e1(:,[1:IL])=exp(-alphal.*h./k);
e2=(h/(delta*k)).*ones(JPLM,IPLM);          e2(:,[IR+1:IPLM])=(1-exp(-alphar.*h./k))./(delta.*alphar);
                                            e2(:,[1:IL])=(1-exp(-alphal.*h./k))./(delta.*alphal);        
e3=ones(JPLM,IPLM);                         e3([JU+1:JPLM],:)=exp(-gammat.*h);
                                            e3([1:JB],:)=exp(-gammab.*h);
e4=(h/(delta*pho)).*ones(JPLM,IPLM);        e4([JU+1:JPLM],:)=(1-exp(-gammat.*h))./(delta.*pho.*gammat);
                                            e4([1:JB],:)=(1-exp(-gammab.*h))./(delta.*pho.*gammab);
                                            

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

for t=1:nsteps 
   
    px(:,ihi)=e1(:,ihi).*px(:,ihi)-e2(:,ihi).*(ux(:,ihi)-ux(:,ilo));
    py(jhi,:)=e1(jhi,:).*py(jhi,:)-e2(jhi,:).*(uy(jhi,:)-uy(jlo,:));
    ux(:,ilo)=e3(:,ilo).*ux(:,ilo)-e4(:,ilo).*(px(:,ihi)-px(:,ilo)+py(:,ihi)-py(:,ilo));
    uy(jlo,:)=e3(jlo,:).*uy(jlo,:)-e4(jlo,:).*(px(jhi,:)-px(jlo,:)+py(jhi,:)-py(jlo,:));
    px(JPLM/2,IPLM/2)=p1source(t);
    py(JPLM/2,IPLM/2)=p1source(t);
    
    
   % px(:,1)=0;      px(1,:)=0;   px(JPLM,:)=0;   px(:,IPLM)=0;
   % py(:,1)=0;      py(1,:)=0;   py(JPLM,:)=0;   py(:,IPLM)=0;  
    p=px+py;
    max(max(p))
    max(max(ux))
    imagesc(p);
    colorbar;
    pause(.05);
    
end