[yuan8503]

function liu

load('g_liu');
G=6.67*10^(-3);
q=0.45;
X=1:2000:100001;
[m,n]=size(X);
x=X;
h=(g_liu(:,2))';
i=1;
j=1;
ggg=zeros(n,n-1);

while j<=n %
while i<=n-1 %
if x(i)==x(j)
T(i)=pi/2;
else
T(i)=atan(h(i)/(x(i)-x(j)));
end
f1=(x(i)-x(j))*h(i+1)-h(i)*(x(i+1)-x(j));
f2=(h(i+1)-h(i)).^2+(x(i+1)-x(i)).^2;
f3=(x(i+1)-x(i))*(T(i)-T(i+1)); %这一行出错了，不知道原因，请教高手指点，非常感谢。

f4=h(i+1)-h(i);
f5=((x(i+1)-X(j)).^2+h(i+1).^2)/((x(i)-X(j)).^2+h(i).^2);
f6=log(f5);
ggg(j,i)=2*G*q*f1*(f3+f4*f6)/f2;
i=i+1;
end
j=j+1;
i=1;
end

gg=-sum(ggg,2);
g_lilun=gg';
g_shice=(g_shice(:,1))';

disp('?í???ì???????????ì???????ó????:')
g_err=sum((g_shice-g_lilun).^2)

figure(1)
plot(x,g_lilun,'r')
title('?í???????ì?????ú??')
xlabel('?à??(m)')
ylabel('?????ì????(g.u.)')
grid on
figure(2)
plot(x,g_shice)
title('?????????ì?????ú??')
xlabel('?à??(m)')
ylabel('?????ì????(g.u.)')
grid on
figure(3)
plotyy(x,g_lilun,x,g_shice,'plot','plot')