%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%一維大地電磁測深模擬退火法反演
%by?xiejingtao
%2016.05.21
%模擬退火法的優點：質量高；初值魯棒性強；簡單、通用、易實現
%缺點：由于要求較高的初始溫度，較慢的降溫速率、較低的終止溫度，以及各溫度下足夠多次的抽樣率，因此優化過程較長。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc
clear?all
close?all

%初始化參數
u=4*pi*1e-7;
f=[1500.000
1500.000
1000.000
800.000????
560.000
480.000
400.000
320.000
240.000
160.000
120.000
80.0000
60.0000
40.0000
30.0000
20.0000
15.0000
10.0000
7.50000
5.00000
3.75000
2.50000
1.87500
1.25000
.937500
.625000
.468800
.312500
.234400
.156300
.117200
.781300E-01
.585900E-01
.390600E-01
.293000E-01
.195300E-01
.146500E-01
.977000E-02
.733000E-02
.488000E-02
.366000E-02
.244000E-02
.183000E-02
.122000E-02
.920000E-03
.610000E-03
.460000E-03];
f=f‘;
nfreq=length（f）;


%給定初始模型
m0=[300;200;100;100;100];

%給定模型取值范圍
area=[01000;01000;01000;01000;01000];

for?j=1:5
????while?m0（j1）area（j2）
????????if?m0（j1）????????????m0（j1）=2*area（j1）-m0（j1）;
????????elseif?m0（j1）>area（j2）
????????????????m0（j1）=2*area（j2）-m0（j1）;
????????end
????end
end

%求原始數據視電阻率
????
%數據參數
m1=[100;20;300];
rho1=m1（1:2:）;
h1=m1（3:）;
n1=length（rho1）;
????
%求視電阻率
for?ifreq=1:nfreq
????w=2*pi*f（ifreq）;
????k=sqrt（i*w*u/rho1（n1））;
????Z=w*u/k;
????for?j=n1:-1:2
????????k=sqrt（i*w*u/rho1（j-1））;
????????Z0=w*u/k;
????????e1=1-exp（2*i*k*h1（j-1））;
????????e2=1+exp（2*i*k*h1（j-1））;
????????Z=Z0*（Z0*e1+Z*e2）/（Z0*e2+Z*e1）;
????end
????rhoas（ifreq）=（abs（Z））^2/（w*u）;
end

for?iter=1:20000
????
????%m0模型視電阻率
????rho0=m0（1:3:）;
????h0=m0（4:5:）;
????n0=length（rho0）;
????
????%求視電阻率
????for?ifreq=1:nfreq
????????w=2*pi*f（ifreq）;
????????k=sqrt（i*w*u/rho0（n0））;
????????Z=w*u/k;
????????for?j=n0:-1:2
????????????k=sqrt（i*w*u/rho0（j-1））;
????????????Z0=w*u/k;
????????????e1=1-exp（2*i*k*h0（j-1））;
????????????e2=1+exp（2*i*k*h0（j-1））;
????????????Z=Z0*（Z0*e1+Z*e2）/（Z0*e