修正後のプログラム

#include<iostream>
#include<math.h>
#include<stdio.h>
#define     PI      3.141592653589793238463
#define		G		 9.80665
using namespace std;

const double rigid=1000/1e-4; //0.1mmで100kg重
const double S=PI*1.5e-2*1.5e-2; //直径3cmのピストン
const double L0=10e-2;//ストローク長10cm
const double P0=4100e2/60*66; //60psi=4100HPa
const double Patm=1024e2;//1atm
const double g=500;//1m/sの速度の時50kg重のダンピング、結構適当
const double dT=1e-3;//1msごとに計算
const double m1=2;
const double m2=40;
const double velocity=20e3/3600; //20km/s
double Fs(double l, double l_){
	if (l<=0) {
		return(rigid*l);
	}else if (l>L0) {
		return(rigid*(l-L0));
	}else{
		double P=P0*L0/(L0-l);
		double tmp=(P-Patm)*S+g*(l-l_)/dT;
		return(tmp);
	}
}

double y(double t){
	double R1=0.02;
	double R2=26.0*0.0254/2.0;
	double w_=sqrt(R1*R1+2.0*R1*R2);
	double x0=1.0;
	double x=velocity*t-x0;
	double height=0.0;
	if( -w_ < x && x < w_) {
		height=sqrt((R1+R2)*(R1+R2)-x*x)-R2;
	}
	return(height);
}
		
void main(void){
	long int it;
	int i;
	double y1[3];
	double y2[3];
	for(i=0;i<=2;i++) y1[i]=0.0;
	for(i=0;i<=2;i++) y2[i]=L0-0.03;

	for(it=0;it*dT<=2;it++){
		double l=L0-(y2[1]-y1[1]);
		double l_=L0-(y2[0]-y1[0]);
		double Fe=m1*G+Fs(l,l_)+m1*(y((it+1)*dT)+y((it-1)*dT)-2.0*y((it)*dT))/dT/dT;
		if(y1[1]>y(it*dT)) Fe=0.0;
		if(Fe<0) Fe=0.0;
	
		y1[2]=2.0*y1[1]-y1[0]+(-G-Fs(l,l_)/m1+Fe/m1)*dT*dT;
		y2[2]=2.0*y2[1]-y2[0]+(-G+Fs(l,l_)/m2)*dT*dT;
		if(y1[2]<y((it+1)*dT)) y1[2]=y((it+1)*dT); //地面の下に突っ込まないように条件式を入れた。


		if (it%1==0)
			cout <<it*dT<<" "<<velocity*it*dT<<" "<<y(it*dT)<<" "<<y1[1]<<" "<<y2[1]<<" "<<Fe/G<<" "<<l<<"\n";
		y1[0]=y1[1];
		y1[1]=y1[2];
		y2[0]=y2[1];
		y2[1]=y2[2];

	}
}