C code that solves the Chaotic Pendulum Equation
Written by Jon J.V. Maestri
#include
#include
#include
#define PI 3.1415926
/**************************************************************************
* Gathering the input *
**************************************************************************/
int i, r, s, j, tv;
double y[5], h, N, n, l, t1[5], k4[5], t, h2, Tmax, m;
double t2[5], t3[5], t4[5], k1[5], k2[5], k3[5];
double k, p, tk, yt, w, g, d, a, b;
double freq, xx,yy, yo, xo, yp ,x, tt,xi,vi;
main(int argc, char *argv[])
{
double f(int, double, double *, double, double, double);
double myrk4(double *, double, double);
FILE *out1, *out2;
d = 9.8;
g = 9.8;
b=strtod(argv[1], NULL);
a=strtod(argv[2], NULL);
w=strtod(argv[3], NULL);
Tmax = strtod(argv[4], NULL);
tt = strtod(argv[5], NULL);
h = tt;
h2 = h/2.;
j=0;
/* y[1] = 8.; initial position*/
/* y[2] = .1234; initial velocity*/
/* Creating Data Files for two-dimensional Plots */
out1 = fopen("/users/webman/htdocs/animgifs/pend.dat", "wb");
out2 = fopen("/users/webman/htdocs/animgifs/pend2.dat", "wb");
for(t=0;t<=Tmax;t=t+tt)
{
/* Declare the image */
gdImagePtr im;
/* Declare an output file */
FILE *out;
char s[]="/users/webman/htdocs/animgifs/run0001.gif";
/* Declare color indexes */
int black;
int white;
int red;
/* Points of polygon */
gdPoint points[100];
/* Allocate the image: 64 pixels across by 64 pixels tall */
im = gdImageCreate(240, 240);
/* Allocate the color white (red, green and blue all minimum).
Since this is the first color in a new image, it will
be the background color. */
white = gdImageColorAllocate(im, 255, 255, 255);
/
/* Allocate the color black (red, green and blue all maximum). */
black = gdImageColorAllocate(im, 0, 0, 0);
red = gdImageColorAllocate(im, 255, 0, 0);
myrk4(y, t, h);
xx=d*10.*sin(y[1]);
yy=d*10.*cos(y[1]);
printf("%e %e %e %e\n", y[1], y[2], xx, yy);
fprintf(out1, "%e %e\n", t, y[1]);
fprintf(out2, "%e %e\n", y[1], y[2]);
/* Draw a line from the upper left to the lower right,
using white color index. */
gdImageLine(im, 120, 120, 120.+xx, 120.+yy, black);
/* define a circle */
r=10.;
yo=120+yy;
xo=120+xx;
N=100;
yp=yo+r;
/* Draw a Polygone. */
for(i=0;i<=N-1;i=i+2){
x = sqrt(r*r-(yp-yo)*(yp-yo)) + xo;
points[i].x = (int)x;
points[i].y = (int)yp;
points[i+1].x = -sqrt(r*r-(yp-yo)*(yp-yo)) + xo;
points[i+1].y = yp;
yp=yp-(2.*r/50.);
}
/* Paint it in red */
gdImageFilledPolygon(im, points, 100, red);
/* Outline it in red; must be done second */
gdImagePolygon(im, points, 100, red);
/* Open a file for writing. "wb" means "write binary", important
under MSDOS, harmless under Unix. */
tv=j;
if(tv<10){
s[36]= tv+48;
}
if(tv<100 && tv>9){
s[35]=(tv/10)+48;
s[36]= (tv%10)+48;
}
if(tv<1000 && tv>99){
s[34]=(tv/100)+48;
s[35]= (((tv%100))/10)+48;
s[36]= (tv%10)+48;
}
if(tv<10000 && tv>999){
s[33]=(tv/1000)+48;
s[34]=(((tv%1000))/100)+48;
s[35]= (tv%100)/10+48;
s[36]= (tv%10)+48;
}
out=fopen(s,"wb");
/* Output the image to the disk file. */
gdImageGif(im, out);
/* Close the file. */
fclose(out);
/* Destroy the image in memory. */
gdImageDestroy(im);
j=j+1;
}
fclose(out1);
fclose(out2);
}
/**************************************************************************
* Subroutine derivatives *
**************************************************************************/
double f(int i, double t, double y[5], double d, double a, double b)
{
if(i==1) return y[2];
if(i == 2) return (-1*(g/d)*sin(y[1]) - a*y[2] + b*cos(w*t));
}
/**************************************************************************
* Subroutine runge-kutta *
**************************************************************************/
double myrk4(double *y, double t, double h)
{
int i;
for(i=1;i<=2;i++) t1[i]=y[i]+(.5)*(k1[i]=h*f(i,t, y, d, a, b));
for(i=1;i<=2;i++) t2[i]=y[i]+.5*(k2[i]=h*f(i,t+h2, t1, d, a, b));
for(i=1;i<=2;i++) t3[i]=y[i]+ (k3[i]=h*f(i,t+h2, t2, d, a, b));
for(i=1;i<=2;i++) k4[i] = h*f(i,t+h, t3, d, a, b);
for(i=1;i<=2;i++) y[i]+=(k1[i]+2*(k2[i]+k3[i])+k4[i])/6.;
}