c - 在gnuplot(和gif转换)中追踪双摆的运动?

c - 在gnuplot(和gif转换)中追踪双摆的运动?,第1张

我已经编写了一个C程序来追踪双摆的运动,但我很难获得gnuplot(从我的class="superseo">c程序控制)来追踪质量的路径(example)。 到目前为止,我已经创建了程序,使得它在每个间隔产生许多png图像(使用runge kutta方法),但是我想将其作为gif输出,因此一条线实时地描绘出质量的路径。

在下面的代码中,我假设问题出现在for循环中滚出gnuplot(为了避免浪费你的时间来筛选它)

/* 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 个答案:

没有答案
本文经用户投稿或网站收集转载,如有侵权请联系本站。

发表评论

0条回复