%%writefile double_pendulum.cpp
#include <cmath>
#include <ctime>
#include <iostream>
#include <iomanip>
#include <string>
#include <unistd.h>
#include <valarray>
#define TIME 50000
#define N 2
using std::string;
using std::to_string;
using std::valarray;
const double PI = 3.141592653589793;
// angle function
// d2(ang[0])/dt2 = f1
// d2(ang[1])/dt2 = f2
// ↓
// dang0 := d(ang[0])/dt = fd1
// dang1 := d(ang[1])/dt = fd2
// d(dang0)/dt = f1
// d(dang1)/dt = f2
void credit(void){
std::cout << "Double Pendulum Calculation" << std::endl;
std::cout << "Runge-Kutte with 4-th Order Accuarcy" << std::endl;
std::cout << "Created : 2022.06.09" << std::endl;
std::cout << "Credit : Nuclear Theory Group, Department of Physics, Graduate School of Science, Tohoku University" << std::endl;
std::cout << std::endl;
}
// time evolution functions
double f1(valarray<double> ang, valarray<double> l, valarray<double> m, double g){
double A = m[1]*l[1]*sin(ang[0]-ang[1])*pow(ang[3],2)+(m[0]+m[1])*g*sin(ang[0]);
double B = -l[0]*sin(ang[0]-ang[1])*pow(ang[2],2)+g*sin(ang[1]);
return (
A-B*m[1]*cos(ang[0]-ang[1])
)/(
(m[1]*pow(cos(ang[0]-ang[1]),2)-(m[0]+m[1]))*l[0]
);
}
double f2(valarray<double> ang, valarray<double> l, valarray<double> m, double g){
double A = m[1]*l[1]*sin(ang[0]-ang[1])*pow(ang[3],2)+(m[0]+m[1])*g*sin(ang[0]);
double B = -l[0]*sin(ang[0]-ang[1])*pow(ang[2],2)+g*sin(ang[1]);
return (
A*cos(ang[0]-ang[1])-B*(m[0]+m[1])
)/(
((m[0]+m[1])-m[1]*pow(cos(ang[0]-ang[1]),2))*l[1]
);
}
double fd1(valarray<double> ang, valarray<double> l, valarray<double> m, double g){
return ang[2];
}
double fd2(valarray<double> ang, valarray<double> l, valarray<double> m, double g){
return ang[3];
}
// to enable to access the functions like array
double (* const angleFuncArr[])(valarray<double> ang, valarray<double> l, valarray<double> m, double g){
fd1,fd2,f1,f2,
};
// convert the angle to euclid postion
// the return is string-formated-data
string euclid_position(valarray<double> ang, valarray<double> l){
valarray<double> x,y;
x.resize(N+1), y.resize(N+1);
x[0] = 0.0, y[0] = 0.0;
string text = "";
for(int i = 0; i < N; i++){
text += to_string(x[i]) + "," + to_string(y[i]) + ",";
x[i+1] = x[i] + l[i]*sin(ang[i]);
y[i+1] = y[i] - l[i]*cos(ang[i]);
}
text += to_string(x[N]) + "," + to_string(y[N]);
return text;
}
// generate random string
string gen_random(const int len) {
std::srand( time(NULL) );
static const char alphanum[] =
"0123456789"
"ABCDEFGHIJKLMNOPQRSTUVWXYZ"
"abcdefghijklmnopqrstuvwxyz";
string tmp_s;
tmp_s.reserve(len);
for (int i = 0; i < len; ++i) {
tmp_s += alphanum[rand() % (sizeof(alphanum) - 1)];
}
return tmp_s;
}
int main(void){
credit();
// generate file id
string ID = gen_random(16);
// open file
FILE *fp;
if( N == 2 ){
string filename = "double_pendulum_" + ID + ".txt";
fp = fopen(filename.c_str(), "w");
std::cout << "> The data file is " << filename << std::endl;
}else{
std::cout << "> Ops! 2 != N case is currently not supported!" << std::endl;
return 0;
}
// global parameters
const double g = 9.8;
const double dt = 0.001;
int i, j, t;
// set angles
// angle[TIME][2*number of angles]
valarray<valarray<double>> angle;
angle.resize(TIME);
for(t = 0; t < TIME; t++){ angle[t].resize(2*N); }
// initial angle
double theta0, phi0;
std::cout << "> Enter the init angle theta_1[deg] : ";
std::cin >> theta0;
angle[0][0] = PI*theta0/180.0;
std::cout << "> Enter the init angle theta_1[deg] : ";
std::cin >> phi0;
angle[0][1] = PI*phi0/180.0;
angle[0][2] = 0.0, angle[0][3] = 0.0; // velocity
// length of arms and mass
// taking all identical
valarray <double> l(N);
l = 1.0;
valarray <double> m(N);
m = 1.0;
fprintf(
fp, "{'l1':%.5f,'l2':%.5f,'m1':%.5f,'m2':%.5f,'dt':%.5f,'TIME':%d,'ang0':%.15f,'dang0':%.15f,'ang1':%.15f,'dang1':%.15f}\n"
,l[0],l[1],m[0],m[1],dt,TIME,theta0,phi0,0.0,0.0
);
// RK4 scheme
// k1~k4 + k0
// k0: dummy k ( to calculate ks with single function )
valarray<valarray<double>> k;
k.resize(5);
for(i = 0; i < 5; i++){ k[i].resize(2*N); }
k[0] = 0.0;
// weight
valarray <double> w = {1.0,0.5,0.5,1.0};
for(t = 0; t < TIME-1; t++){ //for each time
for(i = 1; i <= 4; i++){ //calculate k1~k4
for(j = 0; j < 2*N; j++){ //calculate for each angle freedom
k[i][j] = dt*angleFuncArr[j](k[i-1]*w[i-1]+angle[t],l,m,g);
}
}
// time evoluation
angle[t+1] = angle[t] + (k[1]+2.0*k[2]+2.0*k[3]+k[4])/6.0;
// convert angle -> euclidian position
string dat = euclid_position(angle[t],l);
// write to the file
fprintf(fp,"%s\n",dat.c_str());
}
return 0;
}
