%%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