- #1
atsquishy
- 5
- 0
I want to output an excel file with the results of the trajectories of a two body problem, with initial position and velocity. But my program is not compiling. Any suggestions/problems that you can see?
Code:
#include <iostream>
#include <vector>
#include <fstream>
using namespace std;
void RK4_func(double,double,double,double,double,double,double,double,double,double,double,double);
int main(){
//initialise variables
double time_end,time_step,mass1,x1x,y1y,v1x,v1y,mass2,x2x,y2y,v2x,v2y;
//Get the value
cout << "Enter Simulation Duration in seconds" << endl;
cin >> time_step;
cout << "\n";
cout << "Enter Precision (Time step)" << endl;
cin >> time_end;
cout << "\n";
cout << "Enter mass of body 1 with 10^20kg scale" << endl;
cin >> mass1;
cout << "\n";
cout << "Enter initial body 1 x position" << endl;
cin >> x1x;
cout << "\n";
cout << "Enter initial body 1 y position" << endl;
cin >> y1y;
cout << "\n";
cout << "Enter initial body 1 x velocity" << endl;
cin >> v1x;
cout << "\n";
cout << "Enter initial body 1 x velocity" << endl;
cin >> v1y;
cout << "\n";
cout << "Enter mass of body 2 with 10^20kg scale" << endl;
cin >> mass2;
cout << "\n";
cout << "Enter initial body 2 x position" << endl;
cin >> x2x;
cout << "\n";
cout << "Enter initial body 2 y position" << endl;
cin >> y2y;
cout << "\n";
cout << "Enter initial body 2 x velocity" << endl;
cin >> v2x;
cout << "\n";
cout << "Enter initial body 2 y velocity" << endl;
cin >> v2y;
cout << "\n";
//Pass values int RK4_func
mass1=mass1*10^20
mass2=mass2*10^20
RK4_func(time_end,time_step;mass1;x1x,y1y,v1x,v1y,mass2,x2x,y2y,v2x,v2y);
//Output data
return 0;
}
//Assume function enters data in SI units
void RK4_func (double timeend, double steps, double m1, double x1, double y1, double vx1, double vy1, double m2, double x2, double y2, double vx2,double vy2)
{
//Excel file creation
ofstream outfile1 ("RK4.xls");
outfile1<<"time"<<'\t'<<"x1"<<'\t'<<"y1"<<'\t'<<"vx1"<<'\t'<<"vy1"<<'\t'<<"x2"<<'\t'<<"y2"<<'\t'<<"vx2"<<'\t'<<"vy2"<<endl;
outfile1<<0<<'\t'<<x1<<'\t'<<y1<<'\t'<<vx1<<'\t'<<vy1<<'\t'<<x2<<'\t'<<y2<<'\t'<<vx2<<'\t'<<vy2<<endl;
constant double G = 6.67E-11;
double datalength = 1+(timeend/steps);
double data [8][datalength];
double t=0;
data[0][0]=x1;
data[1][0]=y1;
data[2][0]=x2;
data[3][0]=y2;
data[4][0]=vx1;
data[5][0]=vy1;
data[6][0]=vx2;
data[7][0]=vy2;
for (int i=0; i<datalength; i++)
{
vector<double> interdata(12), k1 (8), k2 (8), k3 (8), k4 (8);
// let 0~3 indicate position x1y1x2y2 4~7 indicate vx1vy1vx2vy2 8~11 indicate ax1ay1ax2ay2
//acceleration x1 y1 x2 y2
interdata[0]=data[0][i];
interdata[1]=data[1][i];
interdata[2]=data[2][i];
interdata[3]=data[3][i];
interdata[4]=data[4][i];
interdata[5]=data[5][i];
interdata[6]=data[6][i];
interdata[7]=data[7][i];
interdata[8] = G*m2*/(((((data[0][i]-data[2][i])^2)+((data[1][i]-data[3][i])^2))^(1/2))^3)*(data[2][i]-data[0][i]);
interdata[9] = G*m2*/(((((data[0][i]-data[2][i])^2)+((data[1][i]-data[3][i])^2))^(1/2))^3)*(data[3][i]-data[1][i]);
interdata[10] = G*m1*/(((((data[0][i]-data[2][i])^2)+((data[1][i]-data[3][i])^2))^(1/2))^3)*(data[0][i]-data[2][i]);
interdata[11] = G*m1*/(((((data[0][i]-data[2][i])^2)+((data[1][i]-data[3][i])^2))^(1/2))^3)*(data[1][i]-data[3][i]);
//k1 is Euler Step
k1[0] = steps*interdata[8];
k1[1] = steps*interdata[9];
k1[2] = steps*interdata[10];
k1[3] = steps*interdata[11];
k1[4] = steps*interdata[4];
k1[5] = steps*interdata[5];
k1[6] = steps*interdata[6];
k1[7] = steps*interdata[7];
interdata[0] = interdata[0]+0.5*k1[4];
interdata[1] = interdata[1]+0.5*k1[5];
interdata[2] = interdata[2]+0.5*k1[6];
interdata[3] = interdata[3]+0.5*k1[7];
interdata[4] = interdata[4]+0.5*k1[0];
interdata[5] = interdata[5]+0.5*k1[1];
interdata[6] = interdata[6]+0.5*k1[2];
interdata[7] = interdata[7]+0.5*k1[3];
interdata[8] = G*m2*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[2]-interdata[0]);
interdata[9] = G*m2*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[3]-interdata[1]);
interdata[10] = G*m1*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[0]-interdata[2]);
interdata[11] = G*m1*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[1]-interdata[3]);
//k2
k2[0] = steps*interdata[8];
k2[1] = steps*interdata[9];
k2[2] = steps*interdata[10];
k2[3] = steps*interdata[11];
k2[4] = steps*interdata[4];
k2[5] = steps*interdata[5];
k2[6] = steps*interdata[6];
k2[7] = steps*interdata[7];
interdata[0] = interdata[0]+0.5*k2[4];
interdata[1] = interdata[1]+0.5*k2[5];
interdata[2] = interdata[2]+0.5*k2[6];
interdata[3] = interdata[3]+0.5*k2[7];
interdata[4] = interdata[4]+0.5*k2[0];
interdata[5] = interdata[5]+0.5*k2[1];
interdata[6] = interdata[6]+0.5*k2[2];
interdata[7] = interdata[7]+0.5*k2[3];
interdata[8] = G*m2*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[2]-interdata[0]);
interdata[9] = G*m2*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[3]-interdata[1]);
interdata[10] = G*m1*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[0]-interdata[2]);
interdata[11] = G*m1*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[1]-interdata[3]);
//k3
k3[0] = steps*interdata[8];
k3[1] = steps*interdata[9];
k3[2] = steps*interdata[10];
k3[3] = steps*interdata[11];
k3[4] = steps*interdata[4];
k3[5] = steps*interdata[5];
k3[6] = steps*interdata[6];
k3[7] = steps*interdata[7];
interdata[0] = interdata[0]+k3[4];
interdata[1] = interdata[1]+k3[5];
interdata[2] = interdata[2]+k3[6];
interdata[3] = interdata[3]+k3[7];
interdata[4] = interdata[4]+k3[0];
interdata[5] = interdata[5]+k3[1];
interdata[6] = interdata[6]+k3[2];
interdata[7] = interdata[7]+k3[3];
interdata[8] = G*m2*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[2]-interdata[0]);
interdata[9] = G*m2*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[3]-interdata[1]);
interdata[10] = G*m1*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[0]-interdata[2]);
interdata[11] = G*m1*/(((((interdata[0]-interdata[2])^2)+((interdata[1]-interdata[3])^2))^(1/2))^3)*(interdata[1]-interdata[3]);
//k4
k4[0] = steps*interdata[8];
k4[1] = steps*interdata[9];
k4[2] = steps*interdata[10];
k4[3] = steps*interdata[11];
k4[4] = steps*interdata[4];
k4[5] = steps*interdata[5];
k4[6] = steps*interdata[6];
k4[7] = steps*interdata[7];
data[0][i+1] = data[0][i]+(1.0/6.0)*(k1[0]+2*k2[0]+2*k3[0]+k4[0]);
data[1][i+1] = data[1][i]+(1.0/6.0)*(k1[1]+2*k2[1]+2*k3[1]+k4[1]);
data[2][i+1] = data[2][i]+(1.0/6.0)*(k1[2]+2*k2[2]+2*k3[2]+k4[2]);
data[3][i+1] = data[3][i]+(1.0/6.0)*(k1[3]+2*k2[3]+2*k3[3]+k4[3]);
data[4][i+1] = data[4][i]+(1.0/6.0)*(k1[4]+2*k2[4]+2*k3[4]+k4[4]);
data[5][i+1] = data[5][i]+(1.0/6.0)*(k1[5]+2*k2[5]+2*k3[5]+k4[5]);
data[6][i+1] = data[6][i]+(1.0/6.0)*(k1[6]+2*k2[6]+2*k3[6]+k4[6]);
data[7][i+1] = data[7][i]+(1.0/6.0)*(k1[7]+2*k2[7]+2*k3[7]+k4[7]);
}
//output new values
for (int i=0; i<datalength; i++){
outfile1<<i*steps<<'\t'<<data[0][i]<<'\t'<<data[1][i]<<'\t'<<data[4][i]<<'\t'<<data[5][i]<<'\t'<<data[2][i]<<'\t'<<data[3][i]<<'\t'<<data[6][i]<<'\t'<<data[7][i]<<endl;
}
outfile1.close();
}