#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;
const double S=PI*1.5e-2*1.5e-2;
const double L0=10e-2;
const double P0=4100e2/60*66;
const double Patm=1024e2;
const double g=500;
const double dT=1e-3;
const double m1=2;
const double m2=40;
const double velocity=20e3/3600;
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];
}
}