ElasticWaveFDTDSimTransientEng.m

%*******************************************************************
% FDTD based Seismic Wave simulation written for Earthquake Shield
% Application. More detailed information about Earthquake Shield
% Application, see article "SeismicCrystals And Earthquake Shield
% Application" written by B. Baykant ALAGZ.
% 
% It uses two images file in BMP format. One describes density map and the
% other describes bulk modules map. Color codes are used in these pictures
% to design simulated area of land such that,
% Black color was reserved for host material
% White color was reserved for host material
% Red color was reserved for wave blocks (Refraction walls)
% Blue color is reserved for wave source
% Bulk modules map and density map are usually the same. Same picture
% can be set following variables,
% YogunlukDosya='DepremEliptic.bmp';
% BulkDosya='DepremEliptic.bmp'; 
% 
% If you use article or these codes in your study, please cite 
% B.B. Alagoz, Seismic Crystals And Earthquake Shield Application, 
% OncuBilim Algorithm And Systems Labs. Vol.09, Art.No:02,(2009)
%
% Programmer: B. Baykant ALAGZ
% Ver.1.0 date: 08.2008
%******************************************************************

clear
%--------Aada grlen simlasyon parametrelerini ayarlayabilirsiniz-----
%*************************************************************************

%Normalize Ortam parametreleri 
% normalize K: Normalize bulk modulu K=Kscat/Khost
% Burada Kscat: Datc bulk modulu
% Khost=Yayc Bulk modulu
% normalize q: Normalize yogunluk q=qhost/qscat
% Burada qscat: Datc yogunluk
% qhost=Yayc yogunluk
% Ortam Haritas Paint programnda izilip bitmap (.bmp) olarak kayt
% edilmelidir.

YogunlukDosya='DepremEliptic.bmp';
BulkDosya='DepremEliptic.bmp';

% Paint dosyasnda kullanlan Renk Kodlamas:
% Siyah: Yaylma Ortam, Beyaz : Datc Bloklar
% Krmz : Yansma Duvar  Mavi : Dalga Kayna Yeil: Belirlenmedi.

NormYogunlukSiyah=1;%Dalga yaylma ortam normalize yogunluk
NormYogunlukBeyaz=10;%Datc ortam normalize yogunluk
NormYogunlukKirmizi=0;%Engelleyici duvar normalize yogunluk
NormYogunlukMavi=1;%Engelleyici duvar normalize yogunluk
NormYogunlukYesil=1;

NormBulkSiyah=1;%Dalga yaylma ortam normalize bulk
NormBulkBeyaz=0.1;%Datc ortam normalize bulk
NormBulkKirmizi=0;%Engelleyici duvar normalize bulk
NormBulkMavi=1;%Engelleyici duvar normalize bulk
NormBulkYesil=1;

%Dalga simlasyon ortam genislik ve yukseklik pixel(nokta) olarak.
%Resimden ortam (K,q) parametreleri ykleniyorsa otomatik olarak
%Resim boyutlarna eitlenir.
AlanGenislik=300;
AlanYukseklik=300;

%Simlasyon Adm says SimulasyonSure ile belirlenir.
SimulasyonSure=1600; %simulasyon sresi adim olarak.
%Durma Periyodu sfrdan farkl seilirse simulasyonu DurmaPeriyodu admda
%bir durdurup o anki basn dalga yaylm resmini ekrana izer ve
%ilerlemek iin bir tua baslmasn bekler.
DurmaPeriyodu=50; %10,20,50 periyotlar seilebilir.

% Dalga yaylma ortamn Kesilen blmden sonras iin grntler.
% (Sadece Sonik kristal yapnn arkasn grntlemek iin tasarlanmtr)

KesilenBolumXBas=250;
KesilenBolumXBit=650;
KesilenBolumYBas=150;
KesilenBolumYBit=900;


%Host materyal dalga hizi m/sn.
Co=3000;

%DeltaT:her simlasyon admnn gerek zaman karl yani Simulasyon zaman znrlkl
%Her simlasyon adm DeltaT saniye.Toplam simulasyon suresi gerek zamanda SimulasyonSure*DeltaT
%DeltaX=0.005;% Herbir pixelin gerekte ka metreye karlk gelecei

%DeltaT=0.000001;

% DeltaX ve DeltaY yaylma ortam znrl. Bir pixel(nokta) gerek uzayda ka metre
%DeltaX=0.01;
%DeltaY=0.01;
%DeltaU=Co*DeltaT; bir simlasyon admnda dalgann kateddii gerek yol

%zmn nmerik kararll iin Courant Koulu salanmaldr.(Kaynak[40])
%Courant Koulu: Rx=Ry < 1/(sqrt(2)*Mpk) seilmelidir. Burada Mpk=Max(K,q)
%Mpk: K ve q ortam parametrelerinin maksimum deeridir.

%-------Otomatik kararl Rx ve Ry belirleme iin--------
%Mpk=Max(Kcelhav,qcelhav);
%KaralilikSinir=1/(sqrt(2)*Mpk);
%Rx=Yakinlik*KaralilikSinir;
%Ry=Rx;%Karesel nokta(pixel) dalm
%--------------------------------------

Rx=0.350;%Courant Says
Ry=0.350;%Courant Says
DeltaX=40;
DeltaT=DeltaX*Rx/Co;
%DeltaT=0.000001;
%Sinizoidal dalga reteci: Genlik*sin(2*pi*f*it);
% f:frekans ve Genlik:sins genlii girilir.
% iaretin periyodu 1/f adet simulasyon zaman. Buda DeltaT/f gerek zaman yapar.
fg=10;
f=fg*DeltaT;
%f=0.0180;
Genlik=1;

%--------Buradan Sonras Simulasyon Kodlar Deitirilmesi Tavsiye edilmez--
%**************************************************************************
%
disp('FDTD Simlasyonu Balamtr... Ltfen Simlasyonun sonlanmasn bekleyiniz...');
disp('-------------------------------------------');

DeltaU=Co*DeltaT;
DeltaX=DeltaU/Rx;
DeltaY=DeltaU/Ry;



%Ortam Normalize Yogunluk Haritasnn karlmas
YogunlukSiyah=NormYogunlukSiyah;%Yaylma ortam datc yok iken normalize yogunluk (qhost/qhost) 1 dir.
YogunlukBeyaz=NormYogunlukBeyaz;%Datc olduu blgede normalize yogunluk (qhost/qscat) dk olmal. 0.1 dir 
YogunlukKirmizi=NormYogunlukKirmizi;
YogunlukMavi=NormYogunlukMavi;
YogunlukYesil=NormYogunlukYesil;


ResimSiyahEsik=130;


DuzlemEni=300
DuzlemBoyu=300
Yogunluk = imread(YogunlukDosya,'bmp');
[DuzlemEniy,DuzlemBoyuy,rgb]=size(Yogunluk);
YogunlukHamHarita=double(Yogunluk);
figure(1)
colormap('gray')
if rgb==3
for i=1:DuzlemEniy
    for j=1:DuzlemBoyuy
          YogunlukHaritasi(i,j)=(YogunlukHamHarita(i,j,1)+YogunlukHamHarita(i,j,2)+YogunlukHamHarita(i,j,3))/3;
    end
end
else
   YogunlukHaritasi=HamHarita; 
end

for i=1:DuzlemEniy
    for j=1:DuzlemBoyuy
          if YogunlukHamHarita(i,j,1)<ResimSiyahEsik && YogunlukHamHarita(i,j,2)<ResimSiyahEsik && YogunlukHamHarita(i,j,3)<ResimSiyahEsik
              q(i,j)=YogunlukSiyah;
              Kaynak(i,j)=0;
              Kayit(i,j)=0;
          elseif  YogunlukHamHarita(i,j,1)>ResimSiyahEsik && YogunlukHamHarita(i,j,2)<ResimSiyahEsik && YogunlukHamHarita(i,j,3)<ResimSiyahEsik
              q(i,j)=YogunlukKirmizi;
              Kaynak(i,j)=0;
              Kayit(i,j)=0;
          elseif  YogunlukHamHarita(i,j,1)<ResimSiyahEsik && YogunlukHamHarita(i,j,2)<ResimSiyahEsik && YogunlukHamHarita(i,j,3)>ResimSiyahEsik
              q(i,j)=YogunlukMavi;
              Kaynak(i,j)=1;
              Kayit(i,j)=0;
          elseif  YogunlukHamHarita(i,j,1)<ResimSiyahEsik && YogunlukHamHarita(i,j,2)>ResimSiyahEsik && YogunlukHamHarita(i,j,3)<ResimSiyahEsik
               q(i,j)=YogunlukYesil;
               Kaynak(i,j)=0;
               Kayit(i,j)=1;
          else 
              q(i,j)=YogunlukBeyaz;
              Kaynak(i,j)=0;
              Kayit(i,j)=0;
          end
    end
end

%Ortam Normalize Bulk Modulu Haritasnn karlmas
BulkSiyah=NormBulkSiyah;%Yaylma ortam datc yok iken normalize Bulk Modulu (Khost/Khost) 1 dir.
BulkBeyaz=NormBulkBeyaz;%Yaylma ortam datc yok iken normalize Bulk Modulu (Kscat/Khost) 1.5 (Buyuk olmal)dir.
BulkKirmizi=NormBulkKirmizi;
BulkMavi=NormBulkMavi;
BulkYesil=NormBulkYesil;

DuzlemEni=300
DuzlemBoyu=300
BulkModule = imread(BulkDosya,'bmp');
[DuzlemEnib,DuzlemBoyub,rgb]=size(BulkModule);
BulkModuleHamHarita=double(BulkModule);

if rgb==3
for i=1:DuzlemEnib
    for j=1:DuzlemBoyub
          BulkModuleHarita(i,j)=(BulkModuleHamHarita(i,j,1)+BulkModuleHamHarita(i,j,2)+BulkModuleHamHarita(i,j,3))/3;
    end
end
else
   BulkModuleHarita=BulkModuleHamHarita; 
end

for i=1:DuzlemEniy
    for j=1:DuzlemBoyuy
          if YogunlukHamHarita(i,j,1)<ResimSiyahEsik && YogunlukHamHarita(i,j,2)<ResimSiyahEsik && YogunlukHamHarita(i,j,3)<ResimSiyahEsik
              K(i,j)=BulkSiyah;
              Kaynak(i,j)=0;
          elseif  YogunlukHamHarita(i,j,1)>ResimSiyahEsik && YogunlukHamHarita(i,j,2)<ResimSiyahEsik && YogunlukHamHarita(i,j,3)<ResimSiyahEsik
              K(i,j)=BulkKirmizi;
              Kaynak(i,j)=0;
          elseif  YogunlukHamHarita(i,j,1)<ResimSiyahEsik && YogunlukHamHarita(i,j,2)<ResimSiyahEsik && YogunlukHamHarita(i,j,3)>ResimSiyahEsik
              K(i,j)=BulkMavi;
              Kaynak(i,j)=1;
          elseif  YogunlukHamHarita(i,j,1)<ResimSiyahEsik && YogunlukHamHarita(i,j,2)>ResimSiyahEsik && YogunlukHamHarita(i,j,3)<ResimSiyahEsik
               K(i,j)=BulkYesil;
               Kaynak(i,j)=0;
          else 
              K(i,j)=BulkBeyaz;
              Kaynak(i,j)=0;
          end
    end
end


%Similasyon alan yogunluk ve bulk resmi alanlarna gre
%ayarlanyor. Her iki alann iine girebilecek en byk alan seiliyor.
AlanGenislik=min(DuzlemEnib,DuzlemEniy);
AlanYukseklik=min(DuzlemBoyub,DuzlemBoyuy);

%K=1*ones(AlanGenislik,AlanYukseklik);
%q=0.5*ones(AlanGenislik,AlanYukseklik);
Vx=zeros(AlanGenislik,AlanYukseklik);
Vy=zeros(AlanGenislik,AlanYukseklik);
P=zeros(AlanGenislik,AlanYukseklik);
Pmax=zeros(AlanGenislik,AlanYukseklik);
Vmax=zeros(AlanGenislik,AlanYukseklik);
Vxs=zeros(AlanGenislik,AlanYukseklik);
Vys=zeros(AlanGenislik,AlanYukseklik);
Ps=zeros(AlanGenislik,AlanYukseklik);
Guc=Ps;
kayitIndex=0;
VbilGecmis=zeros(AlanGenislik-1,AlanYukseklik-1);
DiffV=zeros(AlanGenislik-1,AlanYukseklik-1);
AccMax=zeros(AlanGenislik-1,AlanYukseklik-1);
FMax=zeros(AlanGenislik-1,AlanYukseklik-1);
Acceleration=zeros(AlanGenislik-1,AlanYukseklik-1);
Force=zeros(AlanGenislik-1,AlanYukseklik-1);

for it=1:SimulasyonSure
    fprintf('Tamamlanan Kisim:%f\n',it/SimulasyonSure)
 if (mod(it,DurmaPeriyodu)==0 && DurmaPeriyodu~=0)    
     imagesc(P)
     pause(1);
 end
    %Sinizoidal dalga reteci
     for i=1:AlanGenislik-1
      for j=1:AlanYukseklik-1
       % Eklenir kaynak P=P+aret
       %P(kaynakX,kaynakY)=P(kaynakX,kaynakY)+Genlik*sin(2*pi*f*it);
       % kat kaynak P=aret
        if Kaynak(i,j)==1 %&& it<300
           %P(i,j)=Genlik*sin(2*pi*f*it);
           P(i,j)=Genlik*sin(2*pi*0.5*DeltaT*it)+Genlik*sin(2*pi*1*DeltaT*it)+Genlik*sin(2*pi*2*DeltaT*it)+Genlik*sin(2*pi*3*DeltaT*it)+Genlik*sin(2*pi*4*DeltaT*it)+Genlik*sin(2*pi*5*DeltaT*it)+Genlik*sin(2*pi*10*DeltaT*it);
           %tmax=100;Varyans=200;P(i,j)=Genlik*exp(-((it-tmax)^2)/Varyans);
           DepremMerkezDalga(it)=P(i,j);
        end
        if Kayit(i,j)==1
            kayitIndex=kayitIndex+1;
            GenlikKayitP(kayitIndex)=P(i,j);
            GenlikKayitVx(kayitIndex)=Vx(i,j);
            GenlikKayitVy(kayitIndex)=Vy(i,j);
        end
    % Eklenebilir gauss dalgas exp(-(t-tmax)/Varyans)
    % tmax:Tepe olduu zaman Varyans:Sapmas
    %
    %tmax=100;Varyans=200;P(kaynakX,kaynakY)=Genlik*exp(-((it-tmax)^2)/Varyans);
    %tmax=100;Varyans=200;P(kaynakX,kaynakY)=P(kaynakX,kaynakY)+Genlik*exp(-(it-tmax)^2/Varyans)
      end
    end

%Yatay ve dey hz hesaplamalar
for i=1:AlanGenislik-1
    for j=1:AlanYukseklik-1
        Vx(i,j)=Vx(i,j)-q(i,j)*Rx*(P(i+1,j)-P(i,j));
        Vy(i,j)=Vy(i,j)-q(i,j)*Ry*(P(i,j+1)-P(i,j));
    end
end

  
%Basn hesaplamalar
for i=2:AlanGenislik
    for j=2:AlanYukseklik
        P(i,j)=P(i,j)-K(i,j)*Rx*(Vx(i,j)-Vx(i-1,j))-K(i,j)*Ry*(Vy(i,j)-Vy(i,j-1));
    end
end

%maksimum genlik kayd
for i=1:AlanGenislik-1
    for j=1:AlanYukseklik-1
        if Pmax(i,j)<abs(P(i,j))
            Pmax(i,j)=abs(P(i,j));
            Guc(i,j)=P(i,j)^2;
        end
        Vbil(i,j)=sqrt(Vx(i,j)^2+Vy(i,j)^2);
        if Vmax(i,j)<Vbil(i,j)
            Vmax(i,j)=Vbil(i,j);
        end        
    end
end

% force applied on unit volume of host material
% F=q*DeltaV 

for i=1:AlanGenislik-1
    for j=1:AlanYukseklik-1
        Acceleration(i,j)=(Vbil(i,j)-VbilGecmis(i,j));
        Force(i,j)=abs((1/q(i,j))*(Vbil(i,j)-VbilGecmis(i,j)));
        if AccMax(i,j) < Acceleration(i,j)
            AccMax(i,j)=Acceleration(i,j);
        end
        if FMax(i,j) < Force(i,j)
            FMax(i,j)=Force(i,j);
        end        
    end
end
VbilGecmis=Vbil;

end


fprintf('---Simlasyon Sonlanmtr---\n\n');
fprintf('---Simlasyon Parametreleri---\n');
fprintf('-----------------------------------------------------------\n');
fprintf('DeltaT(Zaman znrl (sn)):%f \n',DeltaT);
fprintf('DeltaU(Dalgann ilerleme znrl (m)):%f \n',DeltaU);
fprintf('DeltaX(Yaylma ortam dey znrl (m)):%f \n',DeltaX);
fprintf('DeltaY(Yaylma ortam yatay znrl (m)):%f \n',DeltaY);
fprintf('Relatif Bulk Modulu:%f \n',BulkBeyaz);
fprintf('Relatif Younluk:%f \n',YogunlukBeyaz);
fprintf('Similasyon Sresi:%f (adm) \n',SimulasyonSure);
fprintf('Similasyon Sresi:%f (sn) \n',SimulasyonSure*DeltaT);

[x y]=size(P)
XP=1:x;
YP=1:y;
XP=XP*DeltaX/1000;
YP=YP*DeltaY/1000;

colormap('default')

kenar = edge(q,'sobel');
Pr=P;
maksimum=0;
for i=1:DuzlemEniy
    for j=1:DuzlemBoyuy
        if (Pr(i,j)>maksimum)
            maksimum=Pr(i,j);
        end
    end
end
            
for i=1:DuzlemEniy
    for j=1:DuzlemBoyuy
        %if (kenar(i,j)==1)
        if q(i,j)==YogunlukBeyaz
            Pr(i,j)=maksimum;
        end
    end
end

figure(1)
imagesc(XP,YP,Pr),hold on,colorbar,hold off;
xlabel('[km]')
ylabel('[km]')
title('Pressure Wave (p)')


Vmax1=Vmax(KesilenBolumXBas:KesilenBolumXBit,KesilenBolumYBas:KesilenBolumYBit);       
XPk=KesilenBolumXBas:KesilenBolumXBit;
YPk=KesilenBolumYBas:KesilenBolumYBit;
XPk=XPk*DeltaX/1000;
YPk=YPk*DeltaY/1000;

MRichter=log10(Vmax1*1000000);
figure(2)
imagesc(YPk,XPk,MRichter),hold on,colorbar,hold off;
%pixval on; pixel value figrde grlr
%pixval off; pixel value figrde grlmez
xlabel('[km]')
ylabel('[km]')
title('Vibration in Rihter Scale [dB]')

XP1=1:x-1;
YP1=1:y-1;
XP1=XP1*DeltaX/1000;
YP1=YP1*DeltaY/1000;

figure(3)
imagesc(YP1,XP1,AccMax);
xlabel('[km]')
ylabel('[km]')
title('Acceleration On Unit Area of Ground')

figure(4)
imagesc(YP1,XP1,FMax);
xlabel('[km]')
ylabel('[km]')
title('Force On Unit Area of Ground')