kalo pake ini tetep gak bisa dibuka…
behhhh susahhhhhhhh
ini dia
% PROGRAM REKURSIF
clear
clc
disp(‘=========== PROGRAM REKURSIF MAGNETOTELLURIK ============’);
disp(‘========== Oleh :AGUNG WAHYU SAPUTRO/ 12405004 =============’);
disp(‘===========================================================’);
disp(”)
disp(‘Petunjuk :’);
disp(‘1.Jumlah lapisan dan frekuensi ditentukan pada sesi input’);
disp(‘2.Input pada tiap lapisan diinputkan nilai rho dan ketebalan’);
%masukan
disp(”)
disp(‘INPUT RHO LAPISAN N+1,FREKUENSI,dan JUMLAH LAPISAN’);
pt=input(‘masukkan rho pada lapisan n+1 = ‘ );
fa=input(‘masukkan batas bawah frekuensi ‘);
fb=input(‘masukkan batas atas frekuensi ‘);
nf=input(‘tentukan jumlah sampling perdekade ‘);
% frekuensi
a=exp(log(10)/nf); %faktor pengali
nd=log10(fb)-log10(fa); % jumlah dekade
nft=nf*nd; %jumlah frek total
for i=1:1:nft
f(i)=fa*a;
fa=f(i);
end
L=input(‘masukkan jumlah lapisan = ‘);
disp(‘===========================================================’);
disp(”)
disp(‘INPUT RHO DAN KETEBALAN TIAP LAPISAN’);
for i=1:1:L
pa(i)=input([‘masukkan harga rho lapisan ke ‘,num2str(i),’= ‘])
h(i)=input([‘masukkan tebal lapisan ke ‘,num2str(i),’= ‘])
end
%——————————————————————
n=length(f);
% iterasi besar
for m=1:1:n;
Zn=complex(sqrt(pi*f(m)*4e-07*pi*pt),sqrt(pi*f(m)*4e-07*pi*pt));
Z=complex(sqrt(pa.*pi*f(m)*4e-07*pi),sqrt(pa.*pi*f(m)*4e-07*pi));
Zt=Zn;
c=sqrt(pi*f(m)*4e-07*pi./pa);
k=complex(c,c);
for a=L:-1:1 % pengulangan kecil
% perhitungan Rj dan Zj
Rj=(Z(a)-Zt)/(Z(a)+Zt);
Zj=Z(a)*(1-Rj*exp(-2*k(a)*h(a)))/(1+Rj*exp(-2*k(a)*h(a)));
Zt=Zj;
end
% akhir pengulangan kecil
Zo(m)= Zj;
end % akhir iterasi besar
% mencari perioda
T=1./f;
figure(1)
plot(log10(f),real(Zo),’*’),title(‘plot Z0 dengan beberapa frekuensi’), hold on
xlabel(‘frekuensi(10^x)’)
legend(‘komponen real’),hold on
plot(log10(f),imag(Zo),’*r’),title(‘plot Z0 dengan beberapa frekuensi’),grid on
xlabel(‘frekuensi(10^x)’)
legend(‘komponen real’,’komponen imaginer’)
ylabel(‘Zo(pada permukaan)’)
figure(2)
plot(log10(T),real(Zo),’*’),title(‘plot Z0 dengan beberapa frekuensi’), hold on
xlabel(‘perioda’)
legend(‘komponen real’),hold on
plot(log10(T),imag(Zo),’*r’),title(‘plot Z0 dengan beberapa frekuensi’),grid on
xlabel(‘perioda / (1/f) ‘)
legend(‘komponen real’,’komponen imaginer’)
ylabel(‘Zo(pada permukaan)’)
% rho apparent and fasa
for m=1:1:n
rho(m)=1/(2*pi*f(m)*4e-07*pi)*(abs(Zo(m)))^2;
fasa(m)=atan(imag(Zo(m))/real(Zo(m)))*180/pi;
end
figure(3)
loglog(f,rho,’*’),title(‘hasil plot apparent resistivity dengan frekuensi berbeda’)
axis equal
xlabel(‘frekuensi’),ylabel(‘apparent resistivity’),grid on
figure(4)
loglog(T,rho,’*’),title(‘hasil plot apparent resistivity dengan frekuensi berbeda’)
axis equal
xlabel(‘perioda / (1/f)’),ylabel(‘apparent resistivity’),grid on
figure(5)
plot(log10(T),fasa,’o’),title(‘plot fasa VS frekuensi’)
xlabel(‘frekuensi(10^x)’),ylabel(‘fasa’),grid on
kukirim juga attacment nya dibawah