symsb_sd_hhc1d_vvklTHETA=geotiffread('theta1111.tif');HH=geotiffread('hh1111.tif');VV=geotiffread('vv1111.tif');p=nan(size(THETA));i1=inversion1();i2=d_hh-i1;i3=[c1,d_hh,...
syms b_s d_hh c1 d_vv k l
THETA=geotiffread('theta1111.tif');
HH=geotiffread('hh1111.tif');
VV=geotiffread('vv1111.tif');
p=nan(size(THETA));
i1=inversion1();
i2=d_hh-i1;
i3=[c1,d_hh,d_vv];
for n=1:size(THETA,1)
for o=1:size(THETA,2)
theta=THETA(n,o);
hh=HH(n,o);
vv=VV(n,o);
i4=[theta,hh,vv];
i=subs(i2,i3,i4);
j1=solve(i(1,1),b_s);
j2=solve(i(2,1),b_s);
j=[real(j1);real(j2)];
if ~isnan(theta) && ~isnan(vv) && ~isnan(hh)
if ~isempty(j1) && ~isempty(j2)
l1=solve(j(1,1)==1+1.5995*k+1.861*k^3,k);
l2=solve(j(2,1)==1+1.5995*k+1.861*k^3,k);
l=[l1;l2];
elseif ~isempty(j1) && isempty(j2)
l=solve(j==1+1.5995*k+1.861*k^3,k);
elseif isempty(j1) && ~isempty(j2)
l=solve(j==1+1.5995*k+1.861*k^3,k);
end
for m=1:numel(l)
if isreal(l(m,1)) && l(m,1)>=0.1 && l(m,1)<0.45
p(n,o)=l(m,1);
disp(num2str([n,o]));
end
end
end
end
end
HH和VV必须都算出来,才能进行下一步。
展开