Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

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

  1. Jun 14, 2015 #1
    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):
    /* 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;
    }
     
  2. jcsd
  3. Jun 14, 2015 #2

    Mark44

    Staff: Mentor

    Is popen() returning a valid pointer? If it's returning a null pointer or -1, the stream isn't valid.
     
  4. Jun 20, 2015 #3
    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.
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: Double Pendulum in C (Trace Motion Real-Time + gif output)?
Loading...