conn@stratus.UUCP (Avery Shealey) (10/16/87)
Thanks for the responces that have been given to me in the past few
weeks. Makes me glad I bought an amiga. This is the graphs and runge kutta
solver that people asked for. Sorry to have to post it, but I am having a
lot of trouble with the mailer. Only about one person I have mailed to
has gotten it. Soory about this.
Delete everything above the cut line and run through /bin/sh to
extract makefile, swing.c, rk.c, graphs.c. It will definitely
compile under Manx 3.4a with all patches. If you have any
questions, comments, or suggestions about the program or code, please
let me know.
Avery Shealey
School of Information & Computer Science, Georgia Tech, Atlanta GA 30332
Internet: conn@stratus.gatech.edu CSNet: conn%stratus@gatech
UUCP: ...!{akgua,allegra,amd,hplabs,seismo,ihnp4}!gatech!stratus!conn
#cut line======================================================
cat <<SHAR_EOF >makefile
CFLAGS = +l
swing: swing.o rk.o graphs.o
ln swing.o rk.o graphs.o -lm32 -lc32
SHAR_EOF
cat <<SHAR_EOF >swing.c
/* Copyright 1987, by Avery Shealey -- this code may be freely
distributed as long as this notice remains intact --
the code may not be sold without express written of
the author */
/* this is the main driver for the differential equation solver -
it is included as an example of how to use the solver and the
graphics - the user must declare the variable numeq and supply
the function derivative to get the runge kutta solver to work -
to use the graphics, simply call the initialization first
and cleanup when you leave -
if you have any suggestions, problems or comments, feel free
to contact me */
/* the program computes the solution to the simple pendulum
problem and plots the phase portrait of the solution -
the pendulum is of length 2 with the force of gravity defined to
be 32 - play with the world values - if you want to see how accurate
the precision is, start the conditions out at about
y[0] = 3.14 -- about pi -- and y[1] = 0 -- no velocity --
- WARNING - WARNING - WARNING - WARNING - WARNING -
unless you have a pendulum fetish, this display will not be very
exciting - it was written only to demonstrate the runge kutta
solver and the graphs interface */
#include <intuition/intuition.h>
#include <exec/types.h>
#include <math.h>
#define NUMEQ 2
int numeq = 2;
double g = 32, len = 2, t, h, y[NUMEQ], dy[NUMEQ];
double xmin, xmax, ymin, ymax;
main()
{
void init(), initgraphics(), cleanupgraphics(), rk();
int GetWinStat();
init();
initgraphics(320, 200, xmin, xmax, ymin, ymax);
moveto(y[0], y[1]);
while (GetWinStat() == TRUE) {
rk(t, h, y, dy);
drawto(y[0], y[1]);
t += h;
}
cleanupgraphics();
}
void derivative(t, y, dy)
double t, y[], dy[];
{
double sin();
dy[0] = y[1];
dy[1] = (-g / len) * sin(y[0]);
}
void init()
{
char *str, *malloc(), *gets();
double atof();
str = malloc(20);
printf("Enter initial t ");
gets(str);
t = atof(str);
printf("Enter initial h ");
gets(str);
h = atof(str);
printf("Enter y[0] ");
gets(str);
y[0] = atof(str);
printf("Enter y[1] ");
gets(str);
y[1] = atof(str);
printf("Enter xmin ");
gets(str);
xmin = atof(str);
printf("Enter xmax ");
gets(str);
xmax = atof(str);
printf("Enter ymin ");
gets(str);
ymin = atof(str);
printf("Enter ymax ");
gets(str);
ymax = atof(str);
}
SHAR_EOF
cat <<SHAR_EOF >rk.c
/* Copyright 1987, by Avery Shealey -- this code may be freely
distributed as long as this notice remains intact --
the code may not be sold without express written of
the author */
/* this module is the runge kutta differential equation solver --
it will solve a system of first order differential equations --
the external variable numeq should contain the number of equations
-- the parameters are :
t - the independent variable
h - the stepsize (generally 0.01 is a good value)
y - the array containing the solution at the last step --
should initially contain the starting values for eqch equation
dy - used to hold temporary values in the computation - this is
not declared locally to avoid the cost of allocating stack
space at each invocation - this could be placed at the start
of this module, but it was declared in the main module - do
what you like with it
*/
extern int numeq;
rk(t, h, y, dy)
double t, h, y[], dy[];
{
void derivative();
double k1[10], k2[10], k3[10], k4[10], temp[10];
register int i;
derivative(t, y, dy);
for (i = 0; i < numeq; i++)
k1[i] = h * dy[i];
for (i = 0; i < numeq; i++)
temp[i] = y[i] + 0.5*k1[i];
derivative(t + 0.5*h, temp, dy);
for (i = 0; i < numeq; i++)
k2[i] = h * dy[i];
for (i = 0; i < numeq; i++)
temp[i] = y[i] + 0.5*k2[i];
derivative(t + 0.5*h, temp, dy);
for (i = 0; i < numeq; i++)
k3[i] = h * dy[i];
for (i = 0; i < numeq; i++)
temp[i] = y[i] + k3[i];
derivative(t + h, temp, dy);
for (i = 0; i < numeq; i++)
k4[i] = h * dy[i];
for (i = 0; i < numeq; i++)
y[i] = y[i] + (1.0 / 6.0) * (k1[i] + 2*k2[i] + 2*k3[i] + k4[i]);
}
SHAR_EOF
cat <<SHAR_EOF >graphs.c
/* Copyright 1987, by Avery Shealey -- this code may be freely
distributed as long as this notice remains intact --
the code may not be sold without express written consent of
the author */
/* this is the graphics module -- it may not be the most efficient,
but is was designed to be easy to use and easily extendable --
the user does not have to worry with any of the normal details --
simply call initgraphics with the necessary startup conditions --
then call getwinstat to see if the window has been closed -- the
module is easily extendible without a lot if trouble -- if you
have any trouble or suggestions, please mail */
#include <intuition/intuition.h>
#include <intuition/intuitionbase.h>
#include <exec/types.h>
/* window for the graphics */
struct NewWindow newwindow = {
0,0,320,200,
-1,-1,
CLOSEWINDOW,
WINDOWDRAG | WINDOWCLOSE | WINDOWDEPTH | SMART_REFRESH,
NULL, NULL,
(UBYTE *)"Graph",
NULL, NULL,
10, 10, 640, 400,
WBENCHSCREEN };
struct IntuitionBase *IntuitionBase;
struct GfxBase *GfxBase;
struct Window *window;
struct RastPort *rp;
int XMAX, YMAX;
double xmin, ymin, xmax, ymax;
double xscale, yscale;
/* cleanup and release resources */
void cleanupgraphics()
{
if (IntuitionBase)
CloseLibrary(IntuitionBase);
if (GfxBase)
CloseLibrary(GfxBase);
if (window)
CloseWindow(window);
}
/* open the libraries and graphics window */
void opengraphics()
{
void *OpenLibrary();
struct Window *OpenWindow();
if ((IntuitionBase = (struct IntuitionBase *)
OpenLibrary("intuition.library", 0)) == NULL) {
printf("unable to open intuition.\n");
cleanupgraphics();
exit(1);
}
if ((GfxBase = (struct GfxBase *)
OpenLibrary("graphics.library",0)) == NULL) {
printf("unable to open graphics.\n");
cleanupgraphics();
exit(2);
}
newwindow.Width = XMAX;
newwindow.Height = YMAX;
if ((window = OpenWindow(&newwindow)) == NULL) {
printf("unable to open window.\n");
cleanupgraphics();
exit(3);
}
rp = window->RPort;
}
/* initialize the graphics -- set up the scale factors
and draw the grid */
/* the parameters are :
width - width of the window
height - height of the window
wxmin, wxmax - min and max x world coordinates
wymin, wymax - min and max y world coordinates
*/
void initgraphics(width, height, wxmin, wxmax, wymin, wymax)
int width, height;
double wxmin, wxmax, wymin, wymax;
{
double i, j;
void moveto(), drawto();
XMAX = width - 1;
YMAX = height - 1;
xmin = wxmin;
xmax = wxmax;
ymin = wymin;
ymax = wymax;
xscale = XMAX / (xmax - xmin);
yscale = YMAX / (ymax - ymin);
opengraphics();
SetAPen(rp, 2);
for (i = xmin; i <= xmax; i += (xmax - xmin)/4) {
moveto(i, ymin);
drawto(i, ymax);
}
for (j = ymin; j <= ymax; j += (ymax - ymin)/4) {
moveto(xmin, j);
drawto(xmax, j);
}
SetAPen(rp, 1);
}
/* plot a point */
void plot(x, y)
double x, y;
{
int xcoord, ycoord;
double floor();
xcoord = (int) floor((xscale * (x - xmin)) + 0.5);
ycoord = YMAX - (int) floor((yscale * (y - ymin)) + 0.5);
Move(rp, xcoord, ycoord);
Draw(rp, xcoord, ycoord);
}
/* draw from the last pen position to the new coordinate */
void drawto(x, y)
double x, y;
{
int xcoord, ycoord;
double floor();
xcoord = (int) floor((xscale * (x - xmin)) + 0.5);
ycoord = YMAX - (int) floor((yscale * (y - ymin)) + 0.5);
Draw(rp, xcoord, ycoord);
}
/* move to a new pen position */
void moveto(x, y)
double x, y;
{
int xcoord, ycoord;
double floor();
xcoord = (int) floor((xscale * (x - xmin)) + 0.5);
ycoord = YMAX - (int) floor((yscale * (y - ymin)) + 0.5);
Move(rp, xcoord, ycoord);
}
/* see if the window close gadget has been selected -- returns
FALSE if window gadget has been selected, TRUE otherwise */
int GetWinStat()
{
struct Message *GetMsg();
struct IntuiMessage *msg;
ULONG class;
int value = TRUE;
while (msg = (struct IntuiMessage *)GetMsg(window->UserPort)) {
class = msg->Class;
ReplyMsg(msg);
if (class == CLOSEWINDOW)
value = FALSE;
}
return(valuembehing.
*.