Hiatus (?)

2:45:00 AM Leo Cahya D 0 Comments

Hello Digital World,


Lama nih gak ngeblog, karena memang lagi sibuk simulasi numerik kampus gitu deh, mana ditambah ngerjain paper + poster hagi. Lalu jobseeking gitu deh, nyobain pertamina, schlumberger, haliburton, dll. Susah juga ya cari kerjaan? hahaha

So, untuk self-challenge selanjutnya saya mencoba membuat perhitungan waktu tiba (traveltime) 3D dengan menggunakan Group Marching Method untuk update value-nya. Berarti eikonal solver gitu mas? Nope, saya mencoba menggunakan Biliniear Interpolation untuk interpolasi waktu tiba tiap grid-nya. Susah sih, makanya lama nih belum selesei2..
Niatnya kalau bisa jadi, mau saya aplikasikan ke Tomografi.

Well, yang jelas untuk postingan mendatang saya akan membahas sekitar itu..
Lama gak ngepost = emang lama bikin programnya (hehe)


REVISI :
ternyata saya keterima kerja. 6 bulan pelatihan, pas sudah diangkat malah banyak game-game bagus rilis gitu deh. so.. :p

see ya,

L

0 comments:

Komputasi Geofisika,

Crosshole Traveltime Tomography without Ray metode LSQR dan extra Tikonov orde 1

6:52:00 AM Leo Cahya D 0 Comments


Hello digital world :)



Lama juga ya saya tidak nge-blog. Seperti biasa saya disibukkan riset di kampus dan sedikit sekali waktu buat eksperimen komputasi, sampai-sampai saya mencuri waktu disela tugas sambil membuat program ini. Awalnya saya tertarik gara-gara paper berjudul "Traveltime tomography of crosshole radar data without ray tracing". Point yang menarik dari paper ini adalah :

- Forward untuk tomography tidak melulu menggunakan ray-tracing. Bisa pakai eikonal, bisa selesaikan rumus propagasi gelombang lalu dipick.  Untuk kasus seismik saya coba menggunakan forward engine dari skripsi saya yang sudah menggunakan kartu grafis untuk komputasi paralel (CUDA C). Program forward 2D saya berbasis  FDTD elastik Virieux dan CPML komastich.

- Matriks Jacobi tidak perlu dihitung tiap iterasi! Ya! :) Saya juga kaget ternyata ada caranya untuk estimasi jacobi dari iterasi awal dengan menggunakan metode Broyden. Dengan ini waktu inversi jadi lebih singkat.
 

Oke, sekarang kita coba hajar, inversi yang saya gunakan adalah LSQR standar tanpa roughening matrix (cuma matriks identitas) dan pakai tikhonov orde 1. Loh kok gak pakai roughening matrix juga? iyasih untuk parameter skala besar sering tidak aman, tapi ya dicari hasil yang aman saja hehe (ubah-ubah parameter awal inversi) :p

DESAIN SURVEY
Oke jadi semisal saya punya daerah survey dengan estimasi perambatan gelombang seperti ini :

Ray Coverage

Testing forward dulu~


Berikut hasilnya untuk tiap model-model unik yang saya ujikan :

MODEL 1

True Model 1
No smoothing
  
Tikhonov 2nd

Well, hasil di atas menunjukkan lebih cantik hasil inversi yang menggunakan tikhonov regulation. 


MODEL 2


True Model 2

No Smoothing

Tikhonov 2nd (gak memuaskan)

Yup, model yang agak kompleks hasilnya kurang memuaskan jika diberi smoothing/roughening matrix. Mungkin cukup susah dalam mencari nilai parameter yang smooth pada model kompleks seperti itu.

Well that's it. Susah juga ya optimisasi skala besar. Ini baru dua dimensi, belum 3 dimensi hehe.

See ya,

L



0 comments:

Programming,

Inversi 1D Occam Magnetotellurik

8:06:00 PM Leo Cahya D 12 Comments

Hello Digital World,

Beberapa waktu lalu saya sempat posting tentang nonlinearitas dari fungsi forward dalam magnetotelluric yang hasil inversinya menunjukkan multi solusi hanya dari sebuah data pengukuran. Well, kali ini saya akan mencoba membuat jenis inversi lain untuk fungsi forward magnetotelluric ini yang namanya Occam[1].



Occam ini salah satu metode inversi Linearized aka Least Square aka kelompok metode inversi ribet yang harus ngitung matriks jacobi mana sering singular nyebelin bikin inversi gagal.

Occam Inversion
Dasarnya sih mirip dengan metode Levenberg-Marquardt aka Damping Gauss-Newton. Hanya saja kita tambahkan parameter delta untuk smoothing (regulasi tikonov orde 1)[3] :


dan update parameternya :


dimana m merupakan parameter yang akan diestimasikan, J matriks jacobi, alpha adalah damping parameter dan d adalah selisih data estimasi dan data pengukuran.

Bagaimana cara membuatnya?
- Pertama kita perlu membuat forward operator-nya yang dapat dilihat di appendix[1]. Ternyata simpel ya?

- Lalu, membuat matriks Jacobi aka Sensitivitas-nya. Apa itu matriks Jacobi? ya matriks dari fungsi forward yang seperti ini :


Ada opsi 3 metode dalam membuat matriks Jacobi [2] :
1. Adjoint equation,
2. Sensitivity equation (Appendix setelah forward MT [1]),
3. Pertubation aka Numerical method.




Saya sih lebih suka menggunakan yang metode numeris (karena malas membuat program analitiknya yang lebih ribet dari forwardnya haha~).
Metode numeris yang saya gunakan adalah forward difference :
 
Rumus dan Skemanya
- Selesaikan rumus update parameter Occam secara iterative. Paling bagus setiap iterasi dibuat multiple damping parameter dan dipilih RMS error yang lebih kecil.


Kira-kira contoh hasilnya seperti ini untuk data sintetik :



Well, kenapa ya harus dibuat smooth ? Setelah program ini saya buat jadi sedikit dapat pencerahan maksud dari Om Constable membuat Occam. Pernah mendengar Occam's Razor ?



Dari pernyataan Occam tersebut, maksud Om Constable membuat metode Inversi Occam ini (dan mungkin asal penamaan metode inversinya juga) adalah mencari solusi paling simpel lewat bentuk model yang smooth (meskipun tidak super fit dengan data pengukuran) untuk menghindari kasus overparameterized seperti pada postingan saya sebelumnya. Maksud saya disini bukan berarti metode lain salah, tetapi metode ini merupakan salah satu opsi untuk menghindari kasus "ke-lebay-an" hasil inversi yang sedang kita tidak harapkan pada suatu pengolahan

So that's it. Lain kali saya update dengan programnya, masih belum rapi nih.

See ya,
-L

Referensi:
[1] Constable, S.C., R.L. Parker, and C.G. Constable, 1987: Occam's Inversion: a practical algorithm for generating smooth models from EM sounding data, Geophysics, 52, pp. 289-300.
[2]  P.  R. McGILLIVRAY and D. W. OLDENBURG. 1990. Methods for calculating Fréchet derivatives and sensitivities for the non-linear inverse problem: a comparative study. Geophysical Prospecting.
[3]  Pei, Donghong, Ph.D. 2007. Modeling and inversion of dispersion curves of surface waves in shallow site investigations. UNIVERSITY OF NEVADA, RENO.

12 comments:

Event

Asisten IESO 2014 (lagi..)

6:51:00 AM Leo Cahya D 1 Comments

Hello digital world,




Kemarin tanggal 13-18 Maret 2014, Saya, Theo, dan Putri Portugal menjadi asisten praktek lagi untuk pelatihan olimpiade SMA yang rencananya akan maju di IESO (International Earth Science Olympiad) 2014 di US kalau tidak salah. Seperti tahun sebelumnya, lokasinya di Hotel University, bedanya perwakilan dosen dari geofisika cuma Pak Afif dan Pak Eddy karena Pak Herlan sibuk ke Jerman rumornya.

Jadwal kegiatan untuk sesi Geofisika seperti tahun sebelumnya, 3 hari kelas, 1 hari lapangan, dan 1 hari ujian seleksi. Saya kali ini dapat bagian mengajarkan metode refraksi dan geolistrik. Hal-hal yang diajarkan pun masih simpel, baru asumsi lapisan datar untuk refraksi, dipole-dipole untuk geolistrik. 

Sesi lapangan sesuai ide Pak Eddy, saya ajak anak-anak tersebut profiling 2D menggunakan dipole-dipole di depan Lab Geofisika yang teduh dan kebetulan Pohon Taloknya sedang berbuah hehe. Lintasan yang mereka buat hanya 20 meter dengan a=1 m dan mereka udah bilang capek, hehe yasudahlah namanya cuma hanya untuk pemahaman.

Selain kegiatan belajar, menurut saya makan malam merupakan acara yang menarik juga, pada saat itu saya sempat berkenalan dengan dosen dan mahasiswa Astronomi ITB, ramah-ramah orangnya.




Sepertinya itu saja yang bisa diceritakan. Well, semoga sukses semua deh buat mereka, Aamiin :)

See ya,

L

1 comments:

Programming

Inversi 1D Magnetotelurrik menggunakan metode VFSA

6:24:00 AM Leo Cahya D 7 Comments

Hello Digital World,


First of all, nonlinearity is a b*tch.
Yap, fungsi-fungsi asumsi dalam metode-metode geofisika itu memang nonlinear, saya baru tau maksud dari bercandaan "kalau orang geofisika ditanya 1+1 pasti bakal dijawab, 'Anda maunya berapa?' " setelah saya membuat beberapa program inversi seperti gravitasi dan MT1D ini.  

Jadi ceritanya saya merampungkan project mt1d dengan teman saya theo yang sebenarnya untuk suatu lomba dua tahun silam (tapi karena ada kegiatan lapangan batal deh). Acuan yang saya pakai menggunakan metode forward yang sama dengan buatan Pak Hendra Grandis*, sedangkan untuk inversinya saya menggunakan metode Very Fast Simulated Annealing (VFSA) besutan Ingber (1989).

Kenapa saya menggunakan VFSA?
Ya menurut saya dalam menyelesaikan masalah inversi perlu dilihat dulu berapa jumlah parameter model yang diinversikan dan lama proses perhitungan forwardnya. Jika parameter < 200 jenis (dipost selanjutnya mungkin akan saya jelaskan kenapa saya berpendapat seperti ini melalui experimen dengan problem simpel) dan waktu untuk forward sebuah model kurang dari 30 detik metode random walk semacam monte carlo dan simulated annealing opsi yang baik, dimana anda tak perlu repot-repot menghitung matriks jacobi dan stress apabila terjadi singularitas.

Hasilnya?
Error kecil tapi model beda dari hasil forward awal. Yup, super ultimate mega nonlinear deh. Bisa dilihat hasil 3 kali running untuk suatu data sintetik dengan model yang sama tetapi hasilnya berbeda.
(Jumlah parameter = 21 (11 nilai rho dan 10 ketebalan lapisan) )

Test 1

Test 2

Test 3

Baru ini, belum seismik haha.. Man gotta dream like never seen obstacles~

Berikut program matlabnya jika ingin dipakai :
forward: mt1dfor2.p
inversi : MT_1D_Inversill.m

*Referensi : H. Grandis, M. Menvielle, M. Roussignol. Bayesian inversion with Markov chains—I. The magnetotelluricone‐dimensional case

see ya
- L

7 comments:

Programming

Forward Modelling Ketinggian Geoid dari model Kontras Densitas

8:12:00 PM Leo Cahya D 0 Comments

Hello Digital World,

Sebulan yang lalu saya dapat sebuah tantangan dari salah seorang klien (ceileh!) untuk membuat forward modelling ketinggian geoid dari model dua dimensi kontras densitas. Jrengjreng.. berbekal dari Michael E. Chapman yang berjudul Techniques interpretation geoid anomalies, saya pun terima quest tersebut karena fulus-nya lumayan.

Oke, yang ada di kepala saya sekarang membuat semacam forward data gravitasi tetapi outputnya N (geoid height) seperti post saya sebelumnya. Berarti saya harus mencari rumus asumsi persegi (yang kemudian kalau memang bisa saya asumsikan untuk tiap grid nanti) dan eng ing eng.. rumus yang saya gunakan pada paper tersebut adalah :


dimana y adalah jarak horizontal dan z adalah kedalaman. Nilai x untuk lebar tidak ada karena asumsi lebar tak terhingga.


Dan memang bisa ternyata :))

Asumsi -200 dan +200 kg/m^3  kontras densitas


Still, meskipun asumsi pengukuran masih pada tempat yang datar alias z=0, yang penting klien sudah puas hehe.

See ya,

L

0 comments:

Programming,

Prosesing Data Seismik menggunakan Matlab (Part 1) - Load data, otak-atik dan cek geometri

10:55:00 PM Leo Cahya D 1 Comments

Hello, digital world..

Hari ini mau saya mau iseng-iseng mengolah data seismik menggunakan Matlab setelah melihat buku karangan Wail A Mousa & co., Processing of Seismic Reflection Data Using MATLAB. Kode-kode dan data untuk latihan sudah dibuat beliau dan dapat diunduh di link berikut.

Oke setelah mendapatkan data bernama Book_Seismic_Data.mat tersebut, kita buka dengan perintah :

>>load Book_Seismic_Data.mat

Pada workspace matlab dapat dilihat ada dua buah variabel dari file tersebut, yaitu D untuk matriks rekaman data dengan panjang trace 594 dan panjang sampel rekaman 150, lalu variabel H yang berisi header dari data tiap trace tersebut.

Menampilkan Rekaman
Untuk menampilkan data, pertama-tama harus kita ketahui waktu cuplik (sampling time) rekaman terlebih dahulu dari variabel H dengan menggunakan perintah :

>>dt=H(1).dt/1000/1000;

kenapa perlu dibagi 1e-6 ? karena umumnya header dt itu dalam satuan µsec.
Selanjutnya untuk menampilkan data anda dapat secara manual menggunakan perintah :
>> x=1:size(D,2); z=0:dt:(size(D,1)-1)*dt; 
>> pcolor(x,z,-1*D); shading interp; axis ij; colormap(gray); colorbar;
>> xlabel('trace number'); ylabel('time (s)');
Plot data dengan colormap gray, dapat diketahui jumlah shotnya ada 18

Jika ingin menggunakan warna umum RGB rekaman seismik :
>> x=1:size(D,2); z=0:dt:(size(D,1)-1)*dt; 
>> pcolor(x,z,-1*D); shading interp; axis ij; colormap(seis_colors); colorbar;
>> xlabel('trace number'); ylabel('time (s)');


Dan juga dari program milik Mousa dapat menampilkan wigglenya dengan perintah
>> x=1:size(D,2); z=0:dt:(size(D,1)-1)*dt; 
>>mwigb(D,1,x,z); axis([0 600 0 3]);
>>xlabel('trace number'); ylabel('time (s)') 


Dari gambar di atas dapat diketahui bahwa jumlah shot yang nya ada 18 buah (bisa juga dilihat nomor fldr pada variabel H dengan perintah "H(end).fldr").

Menampilkan Desain Survei Posisi Shot, Receiver, dan Stacking Chart
Pertama-tama kita ektrak terlebih dahulu nilai geometri dari header, anda bisa menggunakan fungsi extract_geometry.m milik Mousa atau run perintah berikut

>>sx=[H.sx]; sy=[H.sy]; gx=[H.gx]; gy=[H.gy]; 
sz=[H.selev]; gz=[H.gelev]; sg=[H.fldr];
count=1; l=1; k=2;
for i=1:length(sg)
    if k<=length(sg) && sg(k-1)==sg(k)
        count=count+1;
    else
        num_trace_per_sg(l)=count;
        num_sg(l)=l;
        count=1;
        l=l+1;
    end
    k=k+1;
end

dimana variabel tersebut merupakan :
sx: sources x-axis locations
sy: sources y-axis locations
gx: receivers x-axis locations
gy: receivers y-axis locations
num_sg: number of shot gathers
num_trace_per_sg: number of traces/shot gather
sz: elevations of the sources
gz: elevations of the receivers

Untuk melihat jumlah trace tiap shot, dapat dilihat menggunakan perintah :
>>figure,stem(shot_gathers,num_trace_per_sg)
xlabel('Shot gather numbers','FontSize',14)
ylabel('Number of traces/shot gather','FontSize',14)
axis([0,max(shot_gathers)+1,0,max(num_trace_per_sg)+2])
set(gca,'YMinorGrid','on')

dapat dilihat bahwa tiap shot memiliki 33 trace. Selanjutnya posisi ketinggian receiver dan shot dapat dilihat dengan perintah :

>>figure,plot(sz,'*-')
xlabel('Number of traces','FontSize',14)
ylabel('Sources elevation (ft)','FontSize',14)
grid

>>figure,plot(gz,'*-')
xlabel('Number of traces','FontSize',14)
ylabel('Receivers elevation (ft)','FontSize',14)
grid

dan stacking chart dapat kita lihat dengan perintah :
>>num_shots=length(shot_gathers);
stacking_chart(sx,gx,num_shots,num_trace_per_sg)


untuk selanjutnya, akan kita bahas tentang editing dan QC data, kill trace seismik yang rekamannya kacau dan sebagainya..

see ya,

L

1 comments:

Simulasi Numerik,

NUMERICAL MODELLING : GELOMBANG AKUSTIK (Tanpa batas penyerap)

7:38:00 PM Leo Cahya D 3 Comments



NUMERICAL MODELLING : GELOMBANG AKUSTIK


Di sini kita akan menggunakan metode numerik centered difference staggered grid.

Persamaan di atas disebut dengan centered finite difference time domain 2D orde satu dengan akurasi spasial orde dua. Lalu, persamaan penjalaran gelombang akustik untuk medium 2D yang akan kita selesaikan dengan metode numeris di atas:
(2)
P di atas merupakan pressure, vx dan vz merupakan kecepatan partikel ke arah x dan z. ρ merupakan densitas medium dan ca merupakan kecepatan gelombang P. Nilai dx, dz, dan dt merupakan sampling spasial arah x, z, dan sampling waktu.  Kemudian  persamaan (2) diselesaikan menggunakan pers (1) dengan memisalkan medium kita berukuran Z*X, menjadi :
P(i,j) = P(i,j)-dt*ρ(i,j)*ca(i,j)^2*((vx(i+1,j)-vx(i,j))/dx+(vz(i,j+1)-vz(i,j))/dy)

vx(i,j)= vz(i,j)-dt*(P(i,j)-P(i-1,j))/dx/ρ(i,j)

vz(i,j)= vx(i,j)-dt*(P(i,j)-P(i,j-1))/dy/ρ(i,j)

(i=1,2….,Z;  j=1,2,…,X.)                                                                                                                                                   (3)
Dari persamaan di atas kita ketahui bahwa untuk menyelesaikan P, vx dan vz, kita membutuhkan parameter rho dan kecepatan P (ca) terlebih dahulu.

Membuat Parameter Model untuk simulasi
Untuk membuat matriks parameter model lapisan datar tentu sangat gampang, cukup dengan perulangan menggunakan perintah looping seperti for..end.
misal :
Nx = 201;
Ny = 101;
for i=1:Nx
for j=1:floor(Ny/2);
rho(i,j) = 1800;
c(i,j)   = 2000;
end
for j=floor(Ny/2):Ny;
rho(i,j) = 2300;
c(i,j)   = 2600;
K(i,j)=rho(i,j)*c^2;
end
end

Bagaimana dengan model seperti ini?

Tentu susah bukan jika menggunakan perulangan biasa? untuk itu kita perlu akali dengan cara lain seperti berikut :
1.       Pertama, gambar sendiri model menggunakan image editor software seperti Paint, photoshop, atau Paint.NET (recommended).


2.       Set ukurannya sesuai model kita, misal untuk tugas kita disuruh membuat model 2000 x 1000 m dengan jarak spasial dx/dz=10 m. Maka ukuran matriksnya (ukuran model dibagi jarak spasial) sebesar 200x100.

dan matikan anti-aliasing.


3.       Kemudian gambar model anda dengan warna grayscale (warna dimana nilai R, G, dan B-nya sama dari 0 sampai 255).
misal lapisan 1 R: 233 B:233 G:233, lapisan 2 : R: 213 B:213 G:213.

Kemudian save dengan format BMP dengan setting bit-depth = 8 bit.

4.       Kemudian buka di matlab dengan perintah :
data=imread('namafile.bmp');

Sekarang variabel ‘data’ sudah terisi tanda yang menunjukkan masing-masing lapisan, misal untuk model di atas, lapisan satu pada matriks ditandai dengan nilai 0, sedangkan lapisan 2 bernilai 1.
Untuk pembuatan model, tinggal dikondisikan saja sesuai ketentuan :
misal lapisan 1 vp=2300 & rho=2000, lapisan 2 vp=2500 & rho=2200
Nx = 200;
Nz = 100;
for i=1:Nz
for j=1:Nx
if data(i,j)==0
rho(i,j) = 2000;
c(i,j)   = 2300;
else if data(i,j)==1
rho(i,j) = 2200;
c(i,j)   = 2500;
end
end
end
end

lalu simpan data tersebut. Salah satu contohnya dengan menggunakan perintah
>> dlmwrite('c.txt',c,'delimiter','\t')
>> dlmwrite('rho.txt',rho,'delimiter','\t')


Program utama simulasi gelombang seismik
Kita awali dengan membuat m-file baru, setelah itu seperti biasa kita deklarasikan penghapusan memori dan penutupan figure-figure.
close all;
clear all;
clc

Setelah itu, kita deklarasikan parameter-parameter awal :
%input
Nx = 200; %panjang grid model
Ny = 100; %kedalaman grid model
nstep=300; %jumlah iterasi
f0=20; %frekuensi source dominan ricker
dx = 10; %jarak spasial horizontal
dy = 10; %jarak spasial vertikal
c=dlmread('c.txt'); %matriks vp model
rho=dlmread('rho.txt'); %matriks rho model

Lalu, kita buat posisi koordinat untuk penampilan simulasi dan penghitungan parameter. Parameter yang diload tadi kita transpose. Kenapa dilakukan transpose? format pengindeksan matlab  untuk array misal a(i,j) dimana i merupakan menunjukkan posisi indeks vertikal kadang sering tertukar oleh penulisan rumus yang terbiasa menulis i pada koordinat bagian kiri sebagai posisi indeks horizontal.
%penghitungan parameter model dan posisi untuk koordinat surf
c=c';
rho=rho';
K=rho.*c.^2;

is=1;
x1 = (0:Nx-2)*dx+dx/2;
y1 = (0:Ny-2)*dy+dy/2;

dt=1./sqrt(1./dx^2+1./dy^2)/max(max(c));
dt







Kemudian,  kita pesan memori awal dan menghitung koordinat untuk tampilan surf simulasi kita.
%pemesanan memori awal
p=zeros(Nx-1,Ny-1); %pressure
u=zeros(Nx,Ny-1); %vx
v=zeros(Nx-1,Ny); %vz

x=zeros(Nx-1,Ny-1);
y=zeros(Nx-1,Ny-1);

rec=zeros(nstep,Nx);

for j=1:Ny-1
    x(:,j)=x1(:);
end

for i=1:Nx-1
    y(i,:)=y1(:);
end

surf(x,y,p); view(2); axis equal; shading interp;
time=0.0;

Setelah itu, perhitungan utama simulasi kita muat dalam program. Tiap perhitungan matriks P, u, dan v, iterasi ditambahkan gangguan berupa source ricker wavelet dan hasil perhitungan pada elemen matriks komponen v pada posisi permukaan model disimpan ke dalam variabel rec sebagai receiver. Untuk tampilan animasi, tiap 10 iterasi kita kondisikan agar data hasil terakhir perhitungan kita surf.

for time_steps = 1: nstep

    for j=1:Ny-1
        for i=1:Nx-1;
            p(i,j)=p(i,j)-dt*K(i,j)*((u(i+1,j)-u(i,j))/dx+(v(i,j+1)-v(i,j))/dy);
        end
    end

   
    for j=1:Ny-1
        for i=2:Nx-1;
            u(i,j)=u(i,j)-dt*(p(i,j)-p(i-1,j))/dx/rho(i,j);
        end
    end

  for j=2:Ny-1
        for i=1:Nx-1;
            v(i,j)=v(i,j)-dt*(p(i,j)-p(i,j-1))/dy/rho(i,j);
        end
  end
 
  %sumber seismik
    a=pi^2*f0^2;
    f=1e8*2.0*a*(time-1.2/f0)*exp(-a*(time-1.2/f0)^2); %ricker wavelet
    v(is,Ny)=v(is,Ny)+dt*f/rho(is,1); %source diletakkan di pojok kiri model

    rec(time_steps,:)=u(:,Ny-2); %receiver diletakkan di permukaan model
   
    time=time+dt;
   
  if (mod(time_steps,10) == 0)
      figure(gcf)
      surf(x,y,p); view(2); axis equal; shading interp;
      fprintf(2,' time = %8.5f\n', time);
      pause(0.1)
  end 
end


Setelah program tersebut dijalankan, rekaman seismik dapat ditampilkan dengan perintah :
clim=[-500,500];
imagesc(x1,0:dt:nstep*dt,rec,clim);
colormap(flipud(gray));
ylabel('waktu (s)');  xlabel('jarak (m)');
Sekian :)


tips : karena tidak diberi kondisi batas penyerap, sebaiknya iterasi (pada program atas disebut nstep) dikurangi jangan sampai gelombang langsung pada permukaan memantul kembali dari bagian kiri model. 

3 comments: