请教!程序错误,能否帮忙修改和调试`~谢谢~~~ 急!!!!!!!!

[tommyfool]

降落伞拉直阶段力计算:


#define L 12 /*伞伞绳总长度, 可调*/
#define H 0.01 /*步长, 可调*/
#define D0 11.284 /*伞衣名义半径，为伞绳长加伞衣的名义半径*/
#include <stdio.h>
#include <math.h>
double t[10000]; /*定义全局变量, 可调后面的个数*/
double thta[10000];
double xd[10000];
double yd[10000];
double l[10000];
double vys[10000];
double vw[10000];
double s[10000];
double u[10000];
double n=36; /*伞绳的根数*/
double E=10; /*伞绳的理论弹性强度*/
double e[5]; /*龙格－库塔法中的系数*/
double g=9.8; /*当地重力加速度，计算过程中作常量考虑*/
double rou=0.819; /*当地的密度*/
double Ays=1; /*引导伞的面积*/
double Cdys=0.5; /*引导伞的阻力系数*/
double Ad=0.01; /*伞袋的横截面积*/
double Cys=0.5;
int i,j;
main()
{
FILE *ft,*fthta,*fxd,*fyd,*fl,*fvys,*fvw,*fs,*fu,*fFlt;
double Dys;/*引导伞的阻力*/
double Dd; /*伞袋的阻力*/
double Fsh=0; /*伞绳的拉出阻力*/
double Dw=0; /*回收物的气动阻力*/
double De=0; /*主伞系统已拉出部分的气动阻力*/
double m; /* m′拉动的伞系统单位长度的质量*/ /*以下的在整个程序中为常量，赋以初值*/
double mw=110; /*回收物的质量*/
double me; /*已拉出主伞的质量*/
double mys=5; /*引导伞、伞袋的质量*/
double mv; /*未拉出主伞的质量*/
double Fl=0; /*伞绳内张力*/
double msh1=0.18; /*伞绳的单位长度的质量*/
double msh=2.16; /*伞绳的总质量*/
double msy=3.8; /*伞衣的总质量*/
double msy1=1.284; /*伞衣底边研伞绳方向的单位长度的质量*/
double b=5.921; /*伞衣的底边宽度*/
double Flt[10000]; /*定义一拉直力与时间关系的数组*/
double f1(int i,double Dys,double Dd, double Fsh,double Dw,double me,double mys,double mv,double mw,double m,double De,double Fl);
double f2(int i,double Dys,double Dd, double Fsh,double Dw,double me,double mys,double mv,double mw,double m,double De,double Fl);
double f3(int i,double Dys,double Dd, double Fsh,double Dw,double me,double mys,double mv,double mw,double m,double De,double Fl);
double f4(int i,double Dys,double Dd, double Fsh,double Dw,double me,double mys,double mv,double mw,double m,double De,double Fl);
double f5(int i,double Dys,double Dd, double Fsh,double Dw,double me,double mys,double mv,double mw,double m,double De,double Fl);
double f6(int i,double Dys,double Dd, double Fsh,double Dw,double me,double mys,double mv,double mw,double m,double De,double Fl);
double f7(int i,double Dys,double Dd, double Fsh,double Dw,double me,double mys,double mv,double mw,double m,double De,double Fl);
double f8(int i,double Dys,double Dd, double Fsh,double Dw,double me,double mys,double mv,double mw,double m,double De,double Fl);
double f9(int i,double Dys,double Dd, double Fsh,double Dw,double me,double mys,double mv,double mw,double m,double De,double Fl); /*为解决f8中分母为零而外加的一个函数*/

[tommyfool]

e[2]=e[3]=1;
e[0]=e[1]=e[4]=0.5;
t[0]=0; /*赋初值, 根据伞系统的开始量，可调*/
thta[0]=-0.189;
xd[0]=99.907;
yd[0]=291.768;
l[0]=0;
vys[0]=59.949;
vw[0]=59.949;
s[0]=0;
u[0]=0;
for(i=1;;i++)
{
t[i]=t[i-1]+H;
thta[i]=thta[i-1];
xd[i]=xd[i-1];
yd[i]=yd[i-1];
l[i]=l[i-1];
vys[i]=vys[i-1];
vw[i]=vw[i-1];
s[i]=s[i-1];
u[i]=u[i-1];

Dys=0.5*Cdys*rou*Ays*vys[i]*vys[i]; /*以下根据公式代入,改*/
Dd=0.5*1*rou*Ad*vys[i]*vys[i]; /*伞袋的阻力系数取1*/
if(0<=l[i]<=L) me=msh1*l[i]; /*me是一个阶段性变量*/
else me=msh1*L+msy1*b+(msy1-msy1*b)*(1-4*(l[i]-L)/D0);
mv=msh+msy-me;
if(0<=l[i]<=L) /*m′是一个阶段性量*/
m=msh1;
else m=msy1;
Fl=m*u[i]*u[i]+Fsh;

for(j=1;j<=4;j++) /*此循环用于计算点t[i]对应的各个参数的值，并存入数组中i的位置*/
{

thta[i+j]=H*f1(i+j-1,Dys,Dd,Fsh,Dw,me,mys,mv,mw,m,De,Fl); /*计算Kj,校核ij*/
xd[i+j]=H*f2(i+j-1,Dys,Dd,Fsh,Dw,me,mys,mv,mw,m,De,Fl);
yd[i+j]=H*f3(i+j-1,Dys,Dd,Fsh,Dw,me,mys,mv,mw,m,De,Fl);
l[i+j]=H*f4(i+j-1,Dys,Dd,Fsh,Dw,me,mys,mv,mw,m,De,Fl);
vys[i+j]=H*f5(i+j-1,Dys,Dd,Fsh,Dw,me,mys,mv,mw,m,De,Fl);
vw[i+j]=H*f6(i+j-1,Dys,Dd,Fsh,Dw,me,mys,mv,mw,m,De,Fl);
s[i+j]=H*f7(i+j-1,Dys,Dd,Fsh,Dw,me,mys,mv,mw,m,De,Fl);
if((i+j-1)!=1)
u[i+j]=H*f8(i+j-1,Dys,Dd,Fsh,Dw,me,mys,mv,mw,m,De,Fl);
else
u[i+j]=H*f9(i+j-1,Dys,Dd,Fsh,Dw,me,mys,mv,mw,m,De,Fl);

thta[i]=thta[i]+0.33333*thta[i+j]*e[j]; /*计算总和，每一步自动迭加,校核ij*/
xd[i]=xd[i+j]+0.33333*xd[i+j]*e[j];
yd[i]=yd[i+j]+0.33333*yd[i+j]*e[j];
l[i]=l[i+j]+0.33333*l[i+j]*e[j];
vys[i]=vys[i+j]+0.33333*vys[i+j]*e[j];
vw[i]=vw[i+j]+0.33333*vw[i+j]*e[j];
s[i]=s[i+j]+0.33333*s[i+j]*e[j];
u[i]=u[i+j]+0.33333*u[i+j]*e[j];

t[i+j]=t[i+j-1]+H*e[j-1] ; /*运用古典的龙格—库塔形式，为下一步做好准备,校核e[j-1]*/
thta[i+j]=thta[i+j-1]+H*e[j-1];
xd[i+j]=xd[i+j-1]+H*e[j-1];
yd[i+j]=yd[i+j-1]+H*e[j-1];
l[i+j]=l[i+j-1]+H*e[j-1];
vys[i+j]=vys[i+j-1]+H*e[j-1];
vw[i+j]=vw[i+j-1]+H*e[j-1];
s[i+j]=s[i+j-1]+H*e[j-1];
u[i+j]=u[i+j-1]+H*e[j-1];

}

if(l[i]>=(L+D0/2)) break; /*用长度法，控制循环的结束*/
j=i;
}
for(i=0;i<=j;i++)
{
Flt[i]=m*u[i]*u[i]+Fsh;
}

/*数据存储*/
for(i=0;i<=j;i++)
{
ft=fopen("t.txt","a");
fprintf(ft,"%9.5f\t",t[i]);
fclose(ft);
}
for(i=0;i<=j;i++)
{
fthta=fopen("thta.txt","a");
fprintf(fthta,"%9.5f\t",thta[i]);
fclose(fthta);
}
for(i=0;i<=j;i++)
{
fxd=fopen("xd.txt","a");
fprintf(fxd,"%9.5f\t",xd[i]);
fclose(fxd);
}
for(i=0;i<=j;i++)
{
fyd=fopen("yd.txt","a");
fprintf(fyd,"%9.5f\t",yd[i]);
fclose(fyd);

}
for(i=0;i<=j;i++)
{
fl=fopen("l.txt","a");
fprintf(fl,"%9.5f\t",l[i]);
fclose(fl);

}
for(i=0;i<=j;i++)
{
fvys=fopen("vys.txt","a");
fprintf(fvys,"%9.5f\t",vys[i]);
fclose(fvys);
}


for(i=0;i<=j;i++)
{
fvw=fopen("vw.txt","a");
fprintf(fvw,"%9.5f\t",vw[i]);
fclose(fvw);
}
for(i=0;i<=j;i++)
{
fs=fopen("s.txt","a");
fprintf(fs,"%9.5f\t",s[i]);
fclose(fs);

}
for(i=0;i<=j;i++)
{
fu=fopen("u.txt","a");
fprintf(fu,"%9.5f\t",u[i]);
fclose(fu);

}
for(i=0;i<=j;i++)
{
fFlt=fopen("Flt.txt","a");
fprintf(fFlt,"%9.5f\t",Flt[i]);
fclose(fFlt);

}
}
/*各个子函数*/
double f1(int i,double Dys,double Dd, double Fsh,double Dw,double me,double mys,double mv,double mw,double m,double De,double Fl) /*θ函数,i传的是上一个数值*/
{
double t;
t=-g*cos(thta[i])/vw[i];
return t;
}
double f2(int i,double Dys,double Dd, double Fsh,double Dw,double me,double mys,double mv,double mw,double m,double De,double Fl) /*Xd函数*/
{
double t;
t=vw[i]*cos(thta[i]);
return t;
}
double f3(int i,double Dys,double Dd, double Fsh,double Dw,double me,double mys,double mv,double mw,double m,double De,double Fl) /*Yd函数*/
{
double t;
t=vw[i]*sin(thta[i]);
return t;
}
double f4(int i,double Dys,double Dd, double Fsh,double Dw,double me,double mys,double mv,double mw,double m,double De,double Fl) /*L函数*/
{
double t;
t=sqrt((Fl-Fsh)/m);
return t;
}
double f5(int i,double Dys,double Dd, double Fsh,double Dw,double me,double mys,double mv,double mw,double m,double De,double Fl) /*Vys函数*/
{
double t;
t=-g*sin(thta[i])-(Dys+Dd-Fsh)/(mys+mv);
return t;
}
double f6(int i,double Dys,double Dd, double Fsh,double Dw,double me,double mys,double mv,double mw,double m,double De,double Fl) /*Vw函数*/
{
double t;
t=-g*sin(thta[i])-(Dw+De+m*(vw[i]-vys[i])*(vw[i]-vys[i])+Fsh)/(mw+me);
return t;
}
double f7(int i,double Dys,double Dd, double Fsh,double Dw,double me,double mys,double mv,double mw,double m,double De,double Fl) /*S函数*/
{
double t;
t=vw-vys;
return t;
}
double f8(int i,double Dys,double Dd, double Fsh,double Dw,double me,double mys,double mv,double mw,double m,double De,double Fl) /*u函数*/
{
double t;
t=0.5*n*E*(s[i]*l[i]-u[i]*s[i])/pow(l[i],2.0);
return t;
}
double f9(int i,double Dys,double Dd, double Fsh,double Dw,double me,double mys,double mv,double mw,double m,double De,double Fl) /*u函数,分母为零*/
{
double t;
t=0.5*n*E*u[0]*(vw[i+1]-vys[i+1])/(4*m*pow(u[0],3.0)+n*E*(vw[0]-vys[0]));
return t;
}

[code8238]

编绎器说哪里有错误，有什么错误啊？太长了吧！

[jjgoooooole]

没有错误哦

只有3个warning而已！
不过在main()前加一个void就解决了哈

[jjgoooooole]

运行起好象没什么意义哦！

[tommyfool]

就是输出时有些问题.没有输出~~

[tommyfool]

拉直力计算方程：
dθ/dt=-gcosθ/VW
dxd/dt=VWcosθ
dyd/dt=VWsinθ
dL/dt=[(FL-Fsh)/m’]1/2
du/dt=(nE/2m’u)1/L2(S’L-uS)
dVW/dt=-gsinθ-(DW+FL)/(mW+me)
dVys/dt=-gsinθ-(Dys+Fsh)/(mys+mV)
dS/dt=VW-Vys
利用上述方程用数值积分求出各瞬间的各个参数，然后将计算所得u(即dL/dt)代入式
FL=(dm/dt)(dL/dt)+Fsh=m’u2+Fsh
即可求得伞绳内的张力，而最大张力即为拉直力。

[zhouqingyuan]

好象没有什么问题啊？

[tommyfool]

问题是没有输出

[tommyfool]

自己顶

[subtop]

if(0<=l[i]<=L) me=msh1*l[i]; /*me是一个阶段性变量*/
else me=msh1*L+msy1*b+(msy1-msy1*b)*(1-4*(l[i]-L)/D0);
mv=msh+msy-me;
if(0<=l[i]<=L) /*m′是一个阶段性量*/
m=msh1;



0<=l[i]<=L 是什么意思？比较大小可以这样比较么？！
应该是 if(0<=l[i] && l[i] <=L)