%%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']';

iall=1:IPLM;
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
k=4.44e-6;
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(.1);   

%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
e1x=ones(JPLM,IPLM);                        e1x(:,[IR+1:IPLM])=exp(-alphar.*h./k);
                                            e1x(:,[1:IL])=exp(-alphal.*h./k);
e1y=ones(JPLM,IPLM);                        e1y([JU+1:JPLM],:)=exp(-alphat.*h./k);
                                            e1y([1:JB],:)=exp(-alphab.*h./k);
e2x=(h/(delta*k)).*ones(JPLM,IPLM);         e2x(:,[IR+1:IPLM])=(1-exp(-alphar.*h./k))./(delta.*alphar);
                                            e2x(:,[1:IL])=(1-exp(-alphal.*h./k))./(delta.*alphal);        
e2y=(h/(delta*k)).*ones(JPLM,IPLM);         e2y([JU+1:JPLM],:)=(1-exp(-alphat.*h./k))./(delta.*alphat);
                                            e2y([1:JB],:)=(1-exp(-alphab.*h./k))./(delta.*alphab);        
e3x=ones(JPLM,IPLM);                        e3x(:,[IR+1:IPLM])=exp(-gammar.*h);
                                            e3x(:,[1:IL])=exp(-gammal.*h);
e3y=ones(JPLM,IPLM);                        e3y([JU+1:JPLM],:)=exp(-gammat.*h);
                                            e3y([1:JB],:)=exp(-gammab.*h);
e4y=(h/(delta*pho)).*ones(JPLM,IPLM);       e4y([JU+1:JPLM],:)=(1-exp(-gammat.*h))./(delta.*pho.*gammat);
                                            e4y([1:JB],:)=(1-exp(-gammab.*h))./(delta.*pho.*gammab);
e4x=(h/(delta*pho)).*ones(JPLM,IPLM);       e4x(:,[IR+1:IPLM])=(1-exp(-gammar.*h))./(delta.*pho.*gammar);
                                            e4x(:,[1: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)=10*blackmanwin;

for t=1:nsteps 
   
    px(:,ihi)=e1x(:,ihi).*px(:,ihi)-e2x(:,ihi).*(ux(:,ihi)-ux(:,ilo));
    py(jhi,:)=e1y(jhi,:).*py(jhi,:)-e2y(jhi,:).*(uy(jhi,:)-uy(jlo,:));
    ux(:,ilo)=e3x(:,ilo).*ux(:,ilo)-e4x(:,ilo).*(px(:,ihi)-px(:,ilo)+py(:,ihi)-py(:,ilo));
    uy(jlo,:)=e3y(jlo,:).*uy(jlo,:)-e4y(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))
    subplot(2,1,1);
    plot(iall,p(:,IPLM/2),iall,p(JPLM/2,:));
    set(gca,'ylim',[-20,20]);
 
    %imagesc(p);
    colorbar;
    set(gca,'ylim',[-20,20])
    pause(.005);
    
end