/*
* str.c
* (c) Neil Gershenfeld 3/2/03
* integrate the equations of motion for a damped string
*/

#include <math.h>
#include <stdio.h>
#include <X11/Xlib.h>

#define MAX_NODES 10000
#define WIN_SIZE 500
#define WORLD_MIN -3.0
#define WORLD_MAX 3.0

int i,n,nnodes,node,npts;
float d2_1,d2_0,Spring,Dissipation,Force,
	Y0[MAX_NODES],Y1[MAX_NODES],Y2[MAX_NODES];
XPoint Y2Buf[MAX_NODES],YoldBuf[MAX_NODES];
int S;
Display *D;
Window W;
GC Gc,GcRev;
XEvent Event;

main() {
	printf("Number of nodes?  "); scanf("%d",&nnodes);
	printf("Spring constant?  "); scanf("%f",&Spring);
	printf("Bending dissipation?  "); scanf("%f",&Dissipation);
	D = XOpenDisplay("");
	S = DefaultScreen(D);
	W = XCreateSimpleWindow(D,DefaultRootWindow(D),0,0,
		 WIN_SIZE,WIN_SIZE,1,WhitePixel(D,S),WhitePixel(D,S));
	XSelectInput(D,W,ButtonPressMask);
	XStoreName(D,W,"str");
	XMapRaised(D,W);
	Gc = XCreateGC (D,W,0L,(XGCValues *) 0);
	XSetForeground(D,Gc,BlackPixel(D,S));
	XSetBackground(D,Gc,WhitePixel(D,S));
	XSetLineAttributes(D,Gc,0,LineSolid,CapButt,JoinMiter);
	GcRev = XCreateGC(D,W,0L,(XGCValues *) 0);
	XSetForeground(D,GcRev,WhitePixel(D,S));
	XSetBackground(D,GcRev,BlackPixel(D,S));
	XSetLineAttributes(D,GcRev,0,LineSolid,CapButt,JoinMiter);
	for (i = 0; i <= nnodes; ++i) {
		Y2[i] = Y1[i] = Y0[i] = 0;
		YoldBuf[i].y = Y2Buf[i].y = (int) ((WIN_SIZE-1.0) * 0.5);
		YoldBuf[i].x = Y2Buf[i].x = (i * (WIN_SIZE-1.0)) / nnodes;
	}
	while (1) {
		for (i = 1; i < nnodes; ++i) {
			Y0[i] = Y1[i];
			Y1[i] = Y2[i];
		}
		for (i = 1; i < nnodes; ++i) {
			d2_1 = Y1[i+1] + Y1[i-1] - 2 * Y1[i];
			d2_0 = Y0[i+1] + Y0[i-1] - 2 * Y0[i];
			Y2[i] = Spring * d2_1 + Dissipation * (d2_1 - d2_0) +
					  2 * Y1[i] - Y0[i];
		}
		if (XCheckWindowEvent(D,W,ButtonPressMask,&Event) == True) {
			 node =  (int) (nnodes * Event.xbutton.x / (WIN_SIZE-1.0));
			if ((node > 0) && (node < nnodes))
				Y2[node] = 3.0 - (6.0 * Event.xbutton.y / (WIN_SIZE-1.0));
		}
		for (i = 1; i < nnodes; ++i) {
			YoldBuf[i].y = Y2Buf[i].y;
			Y2Buf[i].y = (int) ((WIN_SIZE-1.0) * (3.0-Y2[i]) / 6.0);
		}
		XDrawLines(D,W,GcRev,YoldBuf,(nnodes+1),CoordModeOrigin);
		XDrawLines(D,W,Gc,Y2Buf,(nnodes+1),CoordModeOrigin);
	  XFlush(D);
	  usleep(10);		
	 }
}