yuan8503
2009-03-30, 22:02
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')
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')