風の抵抗を考慮した放物運動
風の抵抗を考慮した放物運動のプログラムを作ったのですが、無限ループが終了しません。どうすればよいでしょうか？以下、プログラムです。
#include<stdio.h>
#include<math.h>
#define MAX 2
double g=9.8; /*重力加速度*/
double k=0.0; /*抵抗係数*/
void main()
{
int rhs(double,double*,double*);
double x,y[MAX];
double h=0.01;
double x_old,y_old[MAX],f[MAX];
double k1[MAX],k2[MAX],k3[MAX],k4[MAX],k_work[MAX];
double z;
int i,j;
x=0;y[0]=28.3;y[1]=0.0;y[2]=28.3;y[3]=0.0;
printf("%f %f %f\n",x,y[1],y[3]);
for(i=1; ;i++){ /*ここは無限ループに*/
if(y[3]>0.0){
break;
} /*地面に落ちたら計算終了*/
x_old=x;
for(j=0;j<MAX;j++)y_old[j]=y[j];
x=i*h;
rhs(x_old,y_old,f);
for(j=0;j<MAX;j++){k1[j]=h*f[j];k_work[j]=y_old[j]+k1[j]/2.;}
rhs(x_old+h/2,k_work,f);
for(j=0;j<MAX;j++){k2[j]=h*f[j];k_work[j]=y_old[j]+k2[j]/2.;}
rhs(x_old+h/2,k_work,f);
for(j=0;j<MAX;j++){k3[j]=h*f[j];k_work[j]=y_old[j]+k3[j];}
rhs(x_old+h,k_work,f);
for(j=0;j<MAX;j++) k4[j]=h*f[j];
for(j=0;j<MAX;j++)
y[j]=y_old[j]+(k1[j]+2*k2[j]+2*k3[j]+k4[j])/6;
if(i%1==0)printf("%f %f %f\n",x,y[1],y[3]);
}
}
int rhs(double x,double y[],double f[])
{
f[0]=y[1];
f[1]=y[3];
f[2]=-k*pow(f[0]*f[0]+f[1]*f[1],0.5)*f[0];
f[3]=-k*pow(f[0]*f[0]+f[1]*f[1],0.5)*f[1]-g;
return 0;
}
お礼
なるほど！わかりやすい説明ありがとうございました！