Tracing out the motion of a double pendulum in gnuplot (and gif conversion)?

856 views Asked by At

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). 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.

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)

/* Header Files */
#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') {

        printf("Please enter a length l1:\n");
        scanf("%lf", &l1);
        printf("Please enter a length l2:\n");
        scanf("%lf", &l2);
        printf("Please enter a mass m1:\n");
        scanf("%lf", &m1);
        printf("Please enter a mass m2:\n");
        scanf("%lf", &m2);
        printf("Please enter an angle theta1:\n");
        scanf("%lf", &th1);
        printf("Please enter an angle theta2:\n");
        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;
}
0

There are 0 answers