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

Click For Summary
SUMMARY

The discussion centers on a C program designed to simulate the motion of a double pendulum and output the results as a GIF using gnuplot. The user initially faced challenges with piping data to gnuplot within a loop, which hindered real-time path tracing. The solution involved opening the gnuplot pipe outside the loop and executing the plotting commands within the loop to generate frames, ultimately allowing for successful GIF output. The key takeaway is the importance of managing the gnuplot session correctly to avoid errors related to the multiplot command.

PREREQUISITES
  • Proficiency in C programming, particularly with file I/O and function definitions.
  • Understanding of numerical methods, specifically the Runge-Kutta method for solving differential equations.
  • Familiarity with gnuplot for data visualization and animation generation.
  • Knowledge of physics principles related to pendulum motion, including gravitational effects.
NEXT STEPS
  • Research "C programming file I/O" to enhance understanding of data handling in C.
  • Study "Runge-Kutta method implementation in C" for better numerical simulation techniques.
  • Explore "gnuplot scripting for animations" to optimize the visualization process.
  • Investigate "double pendulum dynamics" to deepen knowledge of the underlying physics.
USEFUL FOR

Developers and researchers in physics simulations, computer graphics, and numerical methods who are looking to implement real-time data visualization and animation in C using gnuplot.

billyp245
Messages
3
Reaction score
0
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:
/* 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;
}
 
Technology news on Phys.org
Is popen() returning a valid pointer? If it's returning a null pointer or -1, the stream isn't valid.
 
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.
 

Similar threads

  • · Replies 16 ·
Replies
16
Views
3K
  • · Replies 1 ·
Replies
1
Views
3K
  • · Replies 15 ·
Replies
15
Views
3K
  • · Replies 5 ·
Replies
5
Views
2K
Replies
3
Views
6K
  • · Replies 2 ·
Replies
2
Views
2K
  • · Replies 13 ·
Replies
13
Views
5K
  • · Replies 8 ·
Replies
8
Views
4K
  • · Replies 2 ·
Replies
2
Views
2K
  • · Replies 11 ·
Replies
11
Views
4K