#include <fstream>
#include <iostream>

using namespace std;


#define DT 0.001
//k=0.1,c=0.1,r=1.0
double k=50;    //ばね定数
double c=0.002;
double r=1.0;   //空気抵抗

double va,vb;   //vの初期値
double xa,xb;   //xの初期値
double T_MAX=20;

double fxa(double t,double xa,double va);

double fva(double t,double xa,double va);
double fxb(double t,double xb,double vb);

double fvb(double t,double xb,double vb);


int main(){
    double XA,XB,VA,VB;

    XA=30.0;XB=15.0;    //xの初期値 

    VA=0.0;VB=0.0;  //vの初期値

    xa=XA;xb=XB;
    va=VA;vb=VB;

    double t=0.0;
    double xa1,xa2,xa3,xa4,va1,va2,va3,va4;
    double xb1,xb2,xb3,xb4,vb1,vb2,vb3,vb4;

    //----------------------------------//

    ofstream out1("out/sindo9a.txt");
    if(!out1){
        cout<<"Open Error"<<endl;
        exit(0);
    }
    out1 <<"# r="<<r<<endl
         <<"# k="<<k<<endl
         <<"# c="<<c<<endl
         <<"# XA="<<XA<<endl
         <<"# XB="<<XB<<endl
         <<"# VA="<<VA<<endl
         <<"# VB="<<VB<<endl;
         
    out1 <<endl;

    //----------ルンゲクッタ-----------//

    for(t=0.0;t<T_MAX;t+=DT) {
          
        xa1=fxa(t,xa,va)*DT;
        va1=fva(t,xa,va)*DT;
    
        xa2=fxa(t+0.5*DT,xa+0.5*xa1,va+0.5*va1)*DT;
        va2=fva(t+0.5*DT,xa+0.5*xa1,va+0.5*va1)*DT;
    
        xa3=fxa(t+0.5*DT,xa+0.5*xa2,va+0.5*va2)*DT;
        va3=fva(t+0.5*DT,xa+0.5*xa2,va+0.5*va2)*DT;
    
        xa4=fxa(t+DT,xa+xa3,va+va3)*DT;
        va4=fva(t+DT,xa+xa3,va+va3)*DT;
    
        xa+=(xa1+2.0*xa2+2.0*xa3+xa4)/6;
        va+=(va1+2.0*va2+2.0*va3+va4)/6;

        xb1=fxb(t,xb,vb)*DT;
        vb1=fvb(t,xb,vb)*DT;
    
        xb2=fxb(t+0.5*DT,xb+0.5*xb1,vb+0.5*vb1)*DT;
        vb2=fvb(t+0.5*DT,xb+0.5*xb1,vb+0.5*vb1)*DT;
    
        xb3=fxb(t+0.5*DT,xb+0.5*xb2,vb+0.5*vb2)*DT;
        vb3=fvb(t+0.5*DT,xb+0.5*xb2,vb+0.5*vb2)*DT;
    
        xb4=fxb(t+DT,xb+xb3,vb+vb3)*DT;
        vb4=fvb(t+DT,xb+xb3,vb+vb3)*DT;
    
        xb+=(xb1+2.0*xb2+2.0*xb3+xb4)/6;
        vb+=(vb1+2.0*vb2+2.0*vb3+vb4)/6;
        
        out1 << t << ' ' << xa <<endl;
  }

    out1.close();

    //-----------------------------------//
    xa=XA;xb=XB;
    va=VA;vb=VB;

    ofstream out2("out/sindo9b.txt");
    if(!out2){
        cout<<"Open Error"<<endl;
        exit(0);
    }
    out2 <<"# r="<<r<<endl
         <<"# k="<<k<<endl
         <<"# c="<<c<<endl
         <<"# XA="<<XA<<endl
         <<"# XB="<<XB<<endl
         <<"# VA="<<VA<<endl
         <<"# VB="<<VB<<endl;
    out2 <<endl;

    //----------ルンゲクッタ-----------//

    for(t=0.0;t<T_MAX;t+=DT) {
          
        xa1=fxa(t,xa,va)*DT;
        va1=fva(t,xa,va)*DT;
    
        xa2=fxa(t+0.5*DT,xa+0.5*xa1,va+0.5*va1)*DT;
        va2=fva(t+0.5*DT,xa+0.5*xa1,va+0.5*va1)*DT;
    
        xa3=fxa(t+0.5*DT,xa+0.5*xa2,va+0.5*va2)*DT;
        va3=fva(t+0.5*DT,xa+0.5*xa2,va+0.5*va2)*DT;
    
        xa4=fxa(t+DT,xa+xa3,va+va3)*DT;
        va4=fva(t+DT,xa+xa3,va+va3)*DT;
    
        xa+=(xa1+2.0*xa2+2.0*xa3+xa4)/6;
        va+=(va1+2.0*va2+2.0*va3+va4)/6;

        xb1=fxb(t,xb,vb)*DT;
        vb1=fvb(t,xb,vb)*DT;
    
        xb2=fxb(t+0.5*DT,xb+0.5*xb1,vb+0.5*vb1)*DT;
        vb2=fvb(t+0.5*DT,xb+0.5*xb1,vb+0.5*vb1)*DT;
    
        xb3=fxb(t+0.5*DT,xb+0.5*xb2,vb+0.5*vb2)*DT;
        vb3=fvb(t+0.5*DT,xb+0.5*xb2,vb+0.5*vb2)*DT;
    
        xb4=fxb(t+DT,xb+xb3,vb+vb3)*DT;
        vb4=fvb(t+DT,xb+xb3,vb+vb3)*DT;
    
        xb+=(xb1+2.0*xb2+2.0*xb3+xb4)/6;
        vb+=(vb1+2.0*vb2+2.0*vb3+vb4)/6;
        
        out2 << t << ' ' << xb <<endl;
  }

    out2.close();

    //----------ルンゲクッタ-----------//
    xa=XA;xb=XB;
    va=VA;vb=VB;

    ofstream out3("out/sindo9c.txt");
    if(!out3){
        cout<<"Open Error"<<endl;
        exit(0);
    }
    out3 <<"# r="<<r<<endl
         <<"# k="<<k<<endl
         <<"# c="<<c<<endl
         <<"# XA="<<XA<<endl
         <<"# XB="<<XB<<endl
         <<"# VA="<<VA<<endl
         <<"# VB="<<VB<<endl;
    out3 <<endl;

    //----------ルンゲクッタ-----------//

    for(t=0.0;t<T_MAX;t+=DT) {
          
        xa1=fxa(t,xa,va)*DT;
        va1=fva(t,xa,va)*DT;
    
        xa2=fxa(t+0.5*DT,xa+0.5*xa1,va+0.5*va1)*DT;
        va2=fva(t+0.5*DT,xa+0.5*xa1,va+0.5*va1)*DT;
    
        xa3=fxa(t+0.5*DT,xa+0.5*xa2,va+0.5*va2)*DT;
        va3=fva(t+0.5*DT,xa+0.5*xa2,va+0.5*va2)*DT;
    
        xa4=fxa(t+DT,xa+xa3,va+va3)*DT;
        va4=fva(t+DT,xa+xa3,va+va3)*DT;
    
        xa+=(xa1+2.0*xa2+2.0*xa3+xa4)/6;
        va+=(va1+2.0*va2+2.0*va3+va4)/6;

        xb1=fxb(t,xb,vb)*DT;
        vb1=fvb(t,xb,vb)*DT;
    
        xb2=fxb(t+0.5*DT,xb+0.5*xb1,vb+0.5*vb1)*DT;
        vb2=fvb(t+0.5*DT,xb+0.5*xb1,vb+0.5*vb1)*DT;
    
        xb3=fxb(t+0.5*DT,xb+0.5*xb2,vb+0.5*vb2)*DT;
        vb3=fvb(t+0.5*DT,xb+0.5*xb2,vb+0.5*vb2)*DT;
    
        xb4=fxb(t+DT,xb+xb3,vb+vb3)*DT;
        vb4=fvb(t+DT,xb+xb3,vb+vb3)*DT;
    
        xb+=(xb1+2.0*xb2+2.0*xb3+xb4)/6;
        vb+=(vb1+2.0*vb2+2.0*vb3+vb4)/6;
        
        out3 << t << ' ' << (xa+xb)/2.0 <<endl;
  }

    out3.close();

}

double fxa(double t,double xa,double va) { return va; }


double fva(double t,double xa,double va) { return -k*xa-r*va+c*(xb-xa); }

double fxb(double t,double xb,double vb) { return vb; }


double fvb(double t,double xb,double vb) { return -k*xb-r*vb+c*(xa-xb); }