1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41
| function [uchap,d_square] = Sph(y,Q,M,C)
ro=y*inv(M);
n=4;
T(n)=C;
S=ro;
i=n;
L(i)=sqrt(T(i)/Q(i,i))+S(i);
u(i)=-(sqrt(T(i)/Q(i,i))+S(i))-1;
while(1)
u(i)=u(i)+1;
if u(i)>L(i)
if (i==n)
break;
else
i=i+1;
end
else
if (i>1)
epsilon(i)=ro(i)-u(i);
T(i-1)=T(i)-Q(i,i)*(S(i)-u(i))^2;
S(i-1)=ro(i-1)+sum(Q(i-1,i:n)*epsilon(i:n));
i=i-1;
L(i)=sqrt(T(i)/Q(i,i))+S(i);
u(i)=-(sqrt(T(i)/Q(i,i))+S(i))-1;
else
d_square_chap=T(n)-T(1)+Q(1,1)*(S(1)-u(1))^2;
if (d_square_chap < d_square)
uchap=u;
d_square=d_square_chap;
T(n)=d_square_chap;
i=n;
L(i)=sqrt(T(i)/Q(i,i))+S(i);
u(i)=-(sqrt(T(i)/Q(i,i))+S(i))-1;
end
end
end
end
end |
Partager