c++ code

//C++ code for calculating the motion of a double pendulum.
#include <iostream>
#include <fstream>
#include <cmath>
using namespace std;

// function prototypes here
double thetaup(double&, double);
double thetalo(double&, double);
double omegaup(double&, double);
double omegalo(double, double&, double);

//variables
double period = (2), periodmax = 10, r = 1.0, g = 1.0, anought = 0.0, t = 0.0, dt = 0.1, tmax = 100;
double p = ((g/r)*dt); //variables.. could tweak for user input
double nmax = ((tmax/dt));


//starting main
int main()
{
	double thetaupvar = 0.1, thetalovar = 0.0, omegaupvar = 0.0, omegalovar = 0.0;
	ofstream results("results.txt");	//set up file output and below setting up column headings
	ofstream chaosgraph("graph.txt");
	results << "Time t \t" << "theta (upper) \t" << "omega (upper) \t" << "theta (lower) \t" << "omega (lower)" << endl;
	chaosgraph << "Period \t anought" << endl;

	for(int z=0; z<1300000; z++)		// For 200 period 'lines' times 200 anought 'lines'
	{

	for(int n =0; n <=nmax; n++)		// loop here, for t to 100
	{
		t += dt;						// add to t, then do functions that that value of t
	 	//results << t << "\t" << thetaup(thetaupvar, omegaupvar) << "\t" << omegaup(omegaupvar, thetaupvar) << "\t" << thetalo(thetalovar, omegalovar) << "\t" << omegalo(t, omegalovar, thetalovar) << endl;
		thetaup(thetaupvar, omegaupvar);
		omegaup(omegaupvar, thetaupvar);	// commented out above to speed up code;
		thetalo(thetalovar, omegalovar);	// calling functions here
		omegalo(t, omegalovar, thetalovar);
	}

	if(thetalovar<0)			// CHECK: is theta negative?
	{
		chaosgraph << period << "\t" << anought << endl;	//yes - so plot, then reset parameters for next run;
		thetaupvar = 0.1, thetalovar = 0.0, omegaupvar = 0.0, omegalovar = 0.0, t = 0.0;
	}
	else													//no - so just reset the parameters.
	{
		thetaupvar = 0.1, thetalovar = 0.0, omegaupvar = 0.0, omegalovar = 0.0, t = 0.0;
	}

	if(period>=periodmax)		// CHECK: is the period at the end of a line?
	{
		anought += 0.0002;		// yes - reset the period and increment anought
		period = 2;
	}
	else						// no - keep going with the line
	{
		period += 0.005;
	}

	if(anought>=0.25 && period>=periodmax)	// quit when top right corner
	{
		exit(99);
	}

	}
	return 0;


}

// begin called functions
double thetaup(double& thetaupvar, double omegaupvar)			// displacement of upper bob
{
	thetaupvar += omegaupvar * dt;
	return thetaupvar;
}

double omegaup(double& omegaupvar, double thetaupvar)			// angular velocity of upper bob
{
	omegaupvar -= p * sin(thetaupvar);
	return omegaupvar;
}

double thetalo(double& thetalovar, double omegalovar)			// displacement of lower bob
{
	thetalovar += omegalovar * dt;
	return thetalovar;
}

double omegalo(double t, double& omegalovar, double thetalovar)			// angular velocity of lower bob (works out angular acceleration inside function)
{
	double x;
	x = (((2*3.1415926)*t)/(period));									// x is what sin processes below = 2 pi time t over period
	double a;
	a = anought * sin(x);												// a is angular acceleration is worked here
	omegalovar -= ((1/r)*((g*sin(thetalovar)) + (a*cos(thetalovar))))*dt;	// new value for omega (lower)
	return omegalovar;
}


//end

//have a cookie, code


Now energy

//C++ code for calculating the motion of a double pendulum.
#include <iostream>
#include <fstream>
#include <cmath>
using namespace std;

// function prototypes here
double thetaup(double&, double);
double thetalo(double&, double);
double omegaup(double&, double);
double omegalo(double, double&, double);
double kenergy(double&, double, double, double, double);
double penergy(double&, double, double);
double tenergy(double&, double, double);

//variables
double period = (7), periodmax = 9, r = 1.0, g = 1.0, anought = 0.21, t = 0.0, dt = 0.1, tmax = 100;
//taking relative masses mass upper = 20, mass lower = 1.
double p = ((g/r)*dt); //variables.. could tweak for user input
double nmax = ((tmax/dt));


//starting main
int main()
{
	double thetaupvar = 0.1, thetalovar = 0.0, omegaupvar = 0.0, omegalovar = 0.0;
	double kinetic = 0.0, potential = 0.0, total = 0.0;
	ofstream results("results.txt");	//set up file output and below setting up column headings
	ofstream extend("energies.txt");
	ofstream test("this is only a test.txt");
	results << "Time t \t" << "theta (upper) \t" << "omega (upper) \t" << "theta (lower) \t" << "omega (lower)" << endl;
	extend << "Kinetic \t Potential \t Total" << endl;
	test << "time t \t potential" << endl;


		for(int n =0; n <=nmax; n++)		// loop here, for t to 100
		{
			t += dt;						// add to t, then do functions that that value of t
	 		thetaup(thetaupvar, omegaupvar);
			omegaup(omegaupvar, thetaupvar);
			thetalo(thetalovar, omegalovar);	// calling functions here
			omegalo(t, omegalovar, thetalovar);
			extend << kenergy(kinetic, omegaupvar, omegalovar, thetaupvar, thetalovar) << "\t" << penergy(potential, thetaupvar, thetalovar) << "\t" << tenergy(total, kinetic, potential) << endl;
			test << t << "\t" << penergy(potential, thetaupvar, thetalovar) << endl;
		}

	return 0;


}

// begin called functions
double thetaup(double& thetaupvar, double omegaupvar)			// displacement of upper bob
{
	thetaupvar += omegaupvar * dt;
	return thetaupvar;
}

double omegaup(double& omegaupvar, double thetaupvar)			// angular velocity of upper bob
{
	omegaupvar -= p * sin(thetaupvar);
	return omegaupvar;
}

double thetalo(double& thetalovar, double omegalovar)			// displacement of lower bob
{
	thetalovar += omegalovar * dt;
	return thetalovar;
}

double omegalo(double t, double& omegalovar, double thetalovar)			// angular velocity of lower bob (works out angular acceleration inside function)
{
	double x;
	x = (((2*3.1415926)*t)/(period));									// x is what sin processes below = 2 pi time t over period
	double a;
	a = anought * sin(x);												// a is angular acceleration is worked here
	omegalovar -= ((1/r)*((g*sin(thetalovar)) + (a*cos(thetalovar))))*dt;	// new value for omega (lower)
	return omegalovar;
}

double kenergy(double& kinetic, double omegaupvar, double omegalovar, double thetaupvar, double thetalovar)
{
	kinetic = (((10.5)*(omegaupvar)*(omegaupvar))+((0.5)*(omegalovar)*(omegalovar))+((omegaupvar)*(omegalovar)*cos((thetaupvar)-(thetalovar))));
	return kinetic;
}

double penergy(double& potential, double thetaupvar, double thetalovar)
{
	potential = (((21)*(1-cos(thetaupvar)))-(1-cos(thetalovar)));
	return potential;
}

double tenergy(double& total, double kinetic, double potential)
{
	total = ((kinetic) + (potential));
	return total;
}


//end

//have a cookie, code

//