www.pudn.com > Mie.rar > fernaldxiuzheng1.m, change:2014-03-14,size:3040b


clear; 
[fn,pn]=uigetfile('*.mie','打开Mie散射回波信号'); 
if fn==0 
    return; 
end  
a=transpose(zeros(1,12000)); 
filename=dir([pn,'\*.mie']); 
filelength1=length(filename); 
for i=1:filelength1 
    fid1=fopen([pn,filename(i).name]); 
    data1=textscan(fid1,'%f','HeaderLines',9); 
    data=-data1{1,1}(:,1); 
    a=a+data; 
    fclose(fid1);    
end 
MieB=a./filelength1;%12000个数 
b=MieB; 
zt=30;%zt是探测的最大距离 
dz=0.0025; 
vr=b; 
i=length(vr);%i为回波信号的采样数 
z=dz:dz:zt; 
plot(z,vr); 
pr=(vr/50)/(5.3*100000)/(50*0.001)/0.1;%这步可以不做,直接令pr=vr就行 
plot(z,pr); 
pr=vr; 
noise=sum(pr(((zt-5)/dz):(zt/dz)))/(5/dz+1); %将可探测距离内的最后五公里信号的平均值作为背景噪声 
pk=pr-noise 
% for j=1:i 
%     if pk(j,1)<0 
%         pk(j,1)=0; 
%     end 
% end 
subplot(2,2,1); 
plot(z,pk); 
xlabel('探测距离 km'); 
ylabel('原始信号 v'); 
grid on; 
pk=pk'; 
pz=pk.*z.*z;   
subplot(2,2,2); 
plot(z,pz) 
for j=1:1:10 
    pz(1)=(pz(1)+pz(2)+pz(3)+pz(4)+pz(5))/5; 
     pz(2)=(pz(1)+pz(2)+pz(3)+pz(4)+pz(5))/5; 
for k=3:1:11998 
    pz(k)=(pz(k-2)+pz(k-1)+pz(k)+pz(k+1)+pz(k+2))/5; 
end 
 pz(11999)=(pz(11997)+pz(11998)+pz(11999)+pz(11996)+pz(12000))/5; 
    pz(12000)=(pz(11997)+pz(11998)+pz(11999)+pz(11996)+pz(12000))/5; 
end 
subplot(2,2,3); 
plot(z,pz); 
xlabel('探测距离 km'); 
ylabel('距离较正 v*km^2'); 
grid on; 
sr=log(abs(pz)); 
subplot(2,2,2); 
plot(z,sr); 
xlabel('探测距离 km'); 
ylabel('距离较正'); 
grid on; 
 
 
 
bochang=532; 
b=bochang; 
 
 
z0=8;%标定边界值 
 
 
h2=1.54*10^-3*exp(-z0/7)*(532/b)^4;%h2是标定边界值大气分子后向散射系数 
x1=h2;%x1是气溶胶消光系数 
x2=h2*8*pi/3;%x2是大气分子消光系数 
h1=x1/50;%%%%%%%%2014年1月2号改此步,原为h1=50/x1 
sa=50; 
sm=8/3*pi; 
h2=1.54*10^(-3).*exp(-z/7); 
x2=h2*8/3*pi; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%从这里是那个大公式,需要你帮我看看 
sum=0; 
s=50/(8*pi/3); 
z=dz:dz:zt; 
sum1=ones(1,z0/0.0025); 
for j=1:z0/0.0025-1 
    sum=0; 
    for i=j:z0/0.0025-1 
        sum=sum+(x2(1,i)+x2(1,(i+1)))*0.0025/2; 
    end 
  sum1(1,j)=sum;  
end 
sum1(1,z0/0.0025)=sum1(1,z0/0.0025-1); 
sum2=ones(1,z0/0.0025); 
%%%%%%%%%%从71-79是需要你看看的第一个积分,主要是下限是z(同r),是整个大公式的自变量,会变化的,我用int型试了~不管用 
 
x0=pz.*exp(-14*(s-1)*(x2(1,z0/0.0025)-x2)); 
 
sum=0; 
 
for j=1:z0/0.0025-1 
    sum=0; 
    for i=j:z0/0.0025-1 
        sum=sum+(x0(1,i)+x0(1,(i+1)))*0.0025/2; 
    end 
    sum2(1,j)=sum; 
end 
sum2(1,z0/0.0025)=sum2(1,z0/0.0025-1); 
 
 
sum3=ones(1,z0/0.0025); 
for j=1:z0/0.0025 
    sum3(1,j)=pz(1,j).*exp(2*(s-1).*sum1(1,j)); 
end 
x2a=ones(1,z0/0.0025); 
for j=1:z0/0.0025 
    x2a(1,j)=x2(1,j); 
end 
z=dz:dz:z0; 
 
 
x=-s.*x2a+(sum3./(pz(z0/0.0025)./(x1+s.*x2(z0/0.0025))+2.*sum2));%%%%%%%%%%%%%%%2014年1月4日,X1错写为X2 
 
 
 
z=dz:dz:z0; 
xf=ones(1,z0/0.0025); 
xf=x2a(1:z0/0.0025); 
subplot(2,2,4); 
semilogx(x,z,xf,z,'--'); 
ylabel('探测距离 km'); 
xlabel('消光系数 km^-1'); 
grid on; 
 
 
%%%%%%%%%%%%%光学厚度 
sum=0; 
    for i=1:z0/0.0025-1 
        sum=sum+(x(1,i)+x(1,(i+1)))*0.0025/2; 
    end 
    houdu=sum;