# Double Pendulum in C (Trace Motion Real-Time + gif output)?

Tags:
1. Jun 14, 2015

### billyp245

I have written a C program to trace out the motion of a double pendulum, but am having difficulties in getting gnuplot (controlled from my c program) to trace out the paths of the masses (example video below). Thus far I have created the program such that it produces a number of png images at each interval (using runge kutta method), however I want to output it as a gif instead so a line traces out the path of the masses in real time.

Been stuck for quite a while so any help is greatly appreciated!! In the code below I am assuming the problem occurs with piping out to gnuplot from within the for loop (to save you from wasting your time sifting through it)

Code (Text):
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <assert.h>

// Definitions
#define GRAVITY 9.8
#define INCREMENT 0.0175

// Declerations of functions
double th1_xder(double t1d);
double th1_yder(double t2d);
double th2_xder(double t1d, double th1, double t2d, double th2, double l1, double l2, double m1, double m2);
double th2_yder(double t1d, double th1, double t2d, double th2, double l1, double l2, double m1, double m2);

int main (int argc, char * argv[]) {

double x1=0, y1=0, x2=0, y2=0; //coordinates
double l1=0, l2=0;  //lengths
double m1=0, m2=0;  //masses
double th1=0, th2=0; //angles
double t1d=0, t2d=0; //zero velocity initially

//second order Runge-Kutta equations
double k1=0, k2=0, k3=0, k4=0; //for x1
double q1=0, q2=0, q3=0, q4=0; //for y1
double r1=0, r2=0, r3=0, r4=0; //for x2
double s1=0, s2=0, s3=0, s4=0; //for y2

int i = 0;
double x0=0, y0=0;

printf("Enter '1' to input your own data or '0' to use the preset data\n");
char dummy = 'a';
scanf("%c", &dummy);
assert((dummy == '1') || (dummy == '0'));

if (dummy == '1') {

scanf("%lf", &l1);
scanf("%lf", &l2);
scanf("%lf", &m1);
scanf("%lf", &m2);
scanf("%lf", &th1);
scanf("%lf", &th2);
} else {

l1 = 1;
l2 = 1;
m1 = 1;
m2 = 1;
th1= 90;
th2= 0;
}

th1 = th1*(M_PI/180);
th2 = th2*(M_PI/180);

FILE *fsp;

if((fsp=fopen("origin.dat", "w"))==NULL) {
fprintf(stdout, "cannot open origin.dat\n");
exit (EXIT_FAILURE);
}

fprintf(fsp, "0\t0");
fclose(fsp);

FILE *fout;

if((fout=fopen("testout.dat", "w"))==NULL) {
fprintf(stdout, "cannot open testout.dat\n");
exit (EXIT_FAILURE);
}

printf("%f\t%f\t%f\t%f\n", x1, y1, x2, y2);
fprintf(fout, "%f\t%f\t%f\t%f\n", x1, y1, x2, y2);

for(i = 0; i < 250; i++) {

if ((fout=fopen("testout.dat", "w"))==NULL) {
fprintf(stdout, "cannot open testout.dat\n");
exit (EXIT_FAILURE);
}

k1 = th1_xder(t1d);
q1 = th1_yder(t2d);
r1 = th2_xder(t1d, th1, t2d, th2, l1, l2, m1, m2);
s1 = th2_yder(t1d, th1, t2d, th2, l1, l2, m1, m2);

k2 = th1_xder(t1d + (r1/2));
q2 = th1_yder(t2d + (s1/2));
r2 = th2_xder((t1d + (r1/2)), (th1 + (k1/2)), (t2d +(s1/2)), (th2 + (q1/2)), l1, l2, m1, m2);
s2 = th2_yder((t1d + (r1/2)), (th1 + (k1/2)), (t2d +(s1/2)), (th2 + (q1/2)), l1, l2, m1, m2);

k3 = th1_xder(t1d + (r2/2));
q3 = th1_yder(t2d + (s2/2));
r3 = th2_xder((t1d + (r2/2)), (th1 + (k2/2)), (t2d +(s2/2)), (th2 + (q2/2)), l1, l2, m1, m2);
s3 = th2_yder((t1d + (r2/2)), (th1 + (k2/2)), (t2d +(s2/2)), (th2 + (q2/2)), l1, l2, m1, m2);

k4 = th1_xder(t1d + r3);
q4 = th1_yder(t2d + s3);
r4 = th2_xder((t1d + r3), (th1 + k3), (t2d + s3), (th2 + q3), l1, l2, m1, m2);
s4 = th2_yder((t1d + r3), (th1 + k3), (t2d + s3), (th2 + q3), l1, l2, m1, m2);

t1d = t1d + (r1 + 2*r2 + 2*r3 + r4)/6;
t2d = t2d + (s1 + 2*s2 + 2*s3 + s4)/6;
th1 = th1 + (k1 + 2*k2 + 2*k3 + k4)/6;
th2 = th2 + (q1 + 2*q2 + 2*q3 + q4)/6;

x1 = l1*sin(th1);
y1 = -l1*cos(th1);
x2 = x1 + l2*sin(th2);
y2 = y1 - l2*cos(th2);

printf("%f\t%f\t%f\t%f\n", x1, y1, x2, y2);
fprintf(fout, "%f\t%f\t%f\t%f\n", x1, y1, x2, y2);

fclose(fout);

FILE *gnuplotPipe = popen("gnuplot -persist","w");
if (gnuplotPipe) {
fprintf(gnuplotPipe, "set style data lines\n");
fprintf(gnuplotPipe, "set terminal png nocrop enhanced size 1280,720; set output 'yyy%d.png'\n", i);
fprintf(gnuplotPipe, "set title 'frame%d'\n", i);
fprintf(gnuplotPipe, "set multiplot\n");
fprintf(gnuplotPipe, "set xrange [-2.5:2.5]; set yrange [-2.5:2]\n");
fprintf(gnuplotPipe, "unset key; unset ytics; unset xtics\n");
fprintf(gnuplotPipe, "plot 'testout.dat' using 3:4\n");
fprintf(gnuplotPipe, "plot '-' with lines lw 2 lc rgb 'black', 'testout.dat' u 1:2 w points pt 7 ps 2, 'testout.dat' u 3:4 w points pt 7 ps 2, 'origin.dat' u 1:2 w points pt 7 ps 2 lc 0\n");
fprintf(gnuplotPipe, "%f %f\n", x0, y0);
fprintf(gnuplotPipe, "%f %f\n", x1, y1);
fprintf(gnuplotPipe, "%f %f\n", x2, y2);
fprintf(gnuplotPipe, "e\n");
fprintf(gnuplotPipe, "\n");
fprintf(gnuplotPipe, "set nomultiplot\n");
fflush(gnuplotPipe);
fprintf(gnuplotPipe,"exit \n");
pclose(gnuplotPipe);
}
}

return EXIT_SUCCESS;
}

double th1_xder(double t1d) {

double k = t1d*INCREMENT;
return k;
}

double th1_yder(double t2d) {

double m = t2d*INCREMENT;
return m;
}

double th2_xder(double t1d, double th1, double t2d, double th2, double l1, double l2, double m1, double m2) {

double l = INCREMENT*((GRAVITY/l1)*((m2/(m1 + m2))*sin(th2)*cos(th1-th2)-sin(th1))-(m2/(m1 + m2))*sin(th1-th2)*((l2/l1)*t2d*t2d + t1d*t1d*cos(th1-th2)))/(1-((m2/(m1 + m2))*cos(th1-th2)*cos(th1-th2)));
return l;
}

double th2_yder(double t1d, double th1, double t2d, double th2, double l1, double l2, double m1, double m2) {

double p = INCREMENT*((GRAVITY/l2)*(sin(th1)*cos(th1-th2)-sin(th1)) + sin(th1-th2)*((l1/l2)*t1d*t1d + (m2/(m1 + m2))*t2d*t2d*cos(th1-th2)))/(1-((m2/(m1 + m2))*cos(th1-th2)*cos(th1-th2)));
return p;
}

2. Jun 14, 2015

### Staff: Mentor

Is popen() returning a valid pointer? If it's returning a null pointer or -1, the stream isn't valid.

3. Jun 20, 2015

### billyp245

Thanks for the suggestion however the stream was valid. I didn't manage to get the tracing of the paths working (as in the above video) however I did get the output of the gif animation. I simply had to open the pipe into gnuplot outside of my loop, do the plotting inside the loop to create all the frames, and then close the pipe once the loop was finished. I had tried that earlier but hadn't placed the 'set multiplot' command outside the loop which was giving me the error.