clc
clear
ARpath='E:\work\lxy\发表文章\左旗电阻率\';
% Rainpath='G:\GPM\';
datac=[ARpath,'rainz_hn2011to15.txt'];%rain
datat0=[ARpath,'t_180.txt'];
datat1=[ARpath,'t1_180.txt'];
datat2=[ARpath,'t2_180.txt'];
% datah=[ARpath,'exar_1113.txt'];
dataa=[ARpath,'A_rain1113.txt'];
t=(load(datat0));
t1=(load(datat1));
t2=(load(datat2));
datat=load(datac);
data=datat(:,2);
[mi,ni,]=find(data==999999.0)
[m,n]=find(data~=999999.0)
yi=interp1(m,data(m),mi,'linear');%线性插值
data(mi)=yi;
L=length(data);
%data=-data;
n=0;m=0;k=0;r=0;
N=180;%滞后时间
for i=1:L
c(i,2)=i;
m=m+1; r=r+1;
if i==1
c(i,3)=data(1);
c(i,4)=0;
c(i,5)=0;
c(i,6)=0;
c(i,7)=0;
c(i,8)=0;
elseif i==2
c(i,3)=data(1)+data(2);
c(i,4)=data(1);
c(i,5)=data(1);
c(i,6)=data(1);
c(i,7)=data(1);
c(i,8)=data(1);
elseif i>2 && i<=N
c(i,3)=data(i-1)+data(i);
c(i,4)=data(i-1);
c(i,5)=sum(data(1:i-1));
for j=1:i-1
n=n+1;
W(n)=data(j);
T(n)=j;
T1(n)=1/j;
T2(n)=1/j^2;
% a=sort(T);
a=sort(T,'descend');
a1=sort(T1);
a2=sort(T2);
A1=W.*a;
A4=W.*a1;
A5=W.*a2;
end
B(m)=sum(sum(A1));%c3
B4(m)=sum(sum(A4));%c4
B5(m)=sum(sum(A5));%c5
D=B';
D4=B4';
D5=B5';
n=0;
elseif i>N
c(i,3)=data(i-1)+data(i);
c(i,4)=data(i-1);
c(i,5)=sum(data(1:i-1));
for j1=i-N:i-1
k=k+1;
W1(k)=data(j1);
% A1=W1.*t;
% A2=W1.*t1;
% A3=W1.*t2;
end
B1(r)=sum(sum(W1.*t));
B2(r)=sum(sum(W1.*t1));
B3(r)=sum(sum(W1.*t2));
D1=B1';
D2=B2';
D3=B3';
k=0;
end
end
c(1:N,6)=D(1:N);
c(1:N,7)=D4(1:N);
c(1:N,8)=D5(1:N);
c(N+1:L,6)=D1(N+1:L);
c(N+1:L,7)=D2(N+1:L);
c(N+1:L,8)=D3(N+1:L);
c(:,1)=1;
% C=load(datac);
% H1=load(datah);
% H=H1(:,2);
% %
% C1=c.';%转置
% C2=C1*c;
% C3=inv(C2);%逆矩阵
% C4=C3*C1;
% A=C4*H;
A=load(dataa);
ps(:,2)=c*A;
ps(:,1)=datat(:,1);
