/*
xrel.c -- a simple relativistic motion simulator
Copyright 1996, Steve VanDevender <stevev@hexadecimal.uoregon.edu>

Permission is granted to redistribute source code for this program if
and only if all of the following conditions are met:

1.  This copyright notice must not be removed or modified.
2.  Any changes you make to this source code must be documented.
3.  Binary versions must be distributed with complete source code.
4.  No money will be charged for distribution of binaries or source code.

Compile with something like:

cc -O xrel.c -lXt -lX11 -lm -o xrel

*/

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <time.h>
#include <sys/time.h>
#include <X11/Intrinsic.h>
#include <X11/StringDefs.h>
#include <X11/Shell.h>

#define TEMP_MIN      2000	/* minimum star temperature */
#define TEMP_MAX      10000	/* maximum star temperature */
#define TEMP_RANGE    (TEMP_MAX - TEMP_MIN)
#define VIS_TEMP_MIN  2000	/* range of temperatures easily distinguishable */
#define VIS_TEMP_MAX  7000	/* by color */
#define VIS_TEMP_RANGE (VIS_TEMP_MAX - VIS_TEMP_MIN)
#define ABS_MAG_MIN   1000	/* minimum (dimmest) absolute magnitude of stars */
#define ABS_MAG_MAX   -200	/* maximum (brightest) absolute magnitude of stars */
#define ABS_MAG_RANGE (-(ABS_MAG_MAX - ABS_MAG_MIN))
#define VIS_MAG_MIN   600	/* minimum (dimmest) visible magnitude */
#define VIS_MAG_MAX   -200	/* maximum (brightest) visible magnitude */
#define VIS_MAG_RANGE (-(VIS_MAG_MAX - VIS_MAG_MIN))

#define DAY   86400		/* day (s) */
#define YR    31556925.9747	/* year (s) */
#define C     2.99792458e8	/* speed of light (m/s) */
#define H     6.6256e-34	/* Planck's constant (J*s) */
#define K     1.38054e-23	/* Boltzmann's constant (J/K) */
#define PC    3.08568e16	/* parsec (m) */
#define RED   650.0e-9		/* wavelength of red light (m) */
#define GREEN 530.0e-9		/* wavelength of green light (m) */
#define BLUE  425.0e-9		/* wavelength of blue light (m) */

#define B(l, T) ((2 * H * C * C / ((l) * (l) * (l) * (l) * (l))) / expm1(H * C / ((l) * (T) * K)))
				/* returns intensity at wavelength l of blackbody with temperature T */
#define KNOWN_SPACE (100.0 * PC)	/* radius of "known space" */

/* application resource name strings */
#define Nowncmap "owncmap"
#define COwncmap "Owncmap"
#define Ncolortest "colortest"
#define CColortest "Colortest"
#define Nmax "max"
#define CMax "Max"
#define Ncolors "colors"
#define CColors "Colors"
#define Ngamma "gamma"
#define CGamma "Gamma"
#define Nsize "size"
#define CSize "Size"
#define Nviewangle "viewangle"
#define CViewangle "Viewangle"
#define Nfps "fps"
#define CFps "Fps"
#define Nstars "stars"
#define CStars "Stars"
#define Nchunk "chunk"
#define CChunk "Chunk"
#define Nlattice "lattice"
#define CLattice "Lattice"
#define Nbounce "bounce"
#define CBounce "Bounce"
#define Nseed "seed"
#define CSeed "Seed"
#define Nvelocity "velocity"
#define CVelocity "Velocity"
#define Naccel "accel"
#define CAccel "Accel"
#define Nlimit "limit"
#define CLimit "Limit"
#define Ntimescale "timescale"
#define CTimescale "Timescale"

/* application resource data structure */
typedef struct {
	Font font;
	Bool owncmap;		/* use private colormap? */
	Bool colortest;		/* display color test pattern */
	Bool max;		/* resize window to screen maximum */
	int colors;		/* number of colors to use on pseudocolor display */
	String gamma;		/* gamma correction factor of display */
	Dimension size;		/* size (radius) of stars drawn */
	String viewangle;	/* viewing angle (degrees) */
	int fps;		/* frames per second to update */
	int stars;		/* number of stars to display */
	int seed;		/* random number seed */
	int chunk;		/* number of stars to erase/redraw at once */
	String lattice;		/* generate stars in lattice pattern */
	Bool bounce;		/* reverse acceleration at speed limit? */
	String velocity;	/* ship velocity */
	String accel;		/* ship acceleration (m/s) */
	String limit;		/* maximum speed (fraction of c) */
	String timescale;	/* ratio of simulation time to real time */
} ResData;

/* parameters for a star */
typedef struct {
	double x, y, z;		/* location (m) in frame */
	unsigned short T;	/* surface temperature (K) */
	short M;		/* absolute magnitude * 100 */
	short sx, sy;		/* last position plotted on screen */
} star;

/* application resource table */
XtResource resources[] = {
	{ XtNfont, XtCFont,
	  XtRFont, sizeof(Font), XtOffsetOf(ResData, font),
	  XtRString, "variable" },
	{ Nowncmap, COwncmap,
	  XtRBool, sizeof(Bool), XtOffsetOf(ResData, owncmap),
	  XtRString, "False" },
	{ Ncolortest, CColortest,
	  XtRBool, sizeof(Bool), XtOffsetOf(ResData, colortest),
	  XtRString, "False" },
	{ Nmax, CMax,
	  XtRBool, sizeof(Bool), XtOffsetOf(ResData, max),
	  XtRString, "False" },
	{ Ncolors, CColors,
	  XtRInt, sizeof(int), XtOffsetOf(ResData, colors),
	  XtRString, "128" },
	{ Ngamma, CGamma,
	  XtRString, sizeof(String), XtOffsetOf(ResData, gamma),
	  XtRString, "0.5" },
	{ Nsize, CSize,
	  XtRDimension, sizeof(Dimension), XtOffsetOf(ResData, size),
	  XtRString, "1" },
	{ Nviewangle, CViewangle,
	  XtRString, sizeof(String), XtOffsetOf(ResData, viewangle),
	  XtRString, "90.0" },
	{ Nfps, CFps,
	  XtRInt, sizeof(int), XtOffsetOf(ResData, fps),
	  XtRString, "0" },
	{ Nstars, CStars,
	  XtRInt, sizeof(int), XtOffsetOf(ResData, stars),
	  XtRString, "4096" },
	{ Nchunk, CChunk,
	  XtRInt, sizeof(int), XtOffsetOf(ResData, chunk),
	  XtRString, "100" },
	{ Nseed, CSeed,
	  XtRInt, sizeof(int), XtOffsetOf(ResData, seed),
	  XtRString, "0" },
	{ Nlattice, CLattice,
	  XtRString, sizeof(String), XtOffsetOf(ResData, lattice),
	  XtRString, "0.0" },
	{ Nbounce, CBounce,
	  XtRBool, sizeof(Bool), XtOffsetOf(ResData, bounce),
	  XtRString, "False" },
	{ Nvelocity, CVelocity,
	  XtRString, sizeof(String), XtOffsetOf(ResData, velocity),
	  XtRString, "0.0" },
	{ Naccel, CAccel,
	  XtRString, sizeof(String), XtOffsetOf(ResData, accel),
	  XtRString, "9.8" },
	{ Nlimit, CLimit,
	  XtRString, sizeof(String), XtOffsetOf(ResData, limit),
	  XtRString, "0.95" },
	{ Ntimescale, CTimescale,
	  XtRString, sizeof(String), XtOffsetOf(ResData, timescale),
	  XtRString, "86400" },
};

/* default values for critical resources */
String fallback_resources[] = {
	"*drawing.width: 640",
	"*drawing.height: 480",
	"*font: variable",
	(String) 0
};

/* command line options */
XrmOptionDescRec options[] = {
	{ "-" Nowncmap, "*" Nowncmap, XrmoptionNoArg, "True" },
	{ "-" Ncolortest, "*" Ncolortest, XrmoptionNoArg, "True" },
	{ "-" Nmax, "*" Nmax, XrmoptionNoArg, "True" },
	{ "-" Ncolors, "*" Ncolors, XrmoptionSepArg, 0 },
	{ "-" Ngamma, "*" Ngamma, XrmoptionSepArg, 0 },
	{ "-" Nsize, "*" Nsize, XrmoptionSepArg, 0 },
	{ "-" Nviewangle, "*" Nviewangle, XrmoptionSepArg, 0 },
	{ "-" Nfps, "*" Nfps, XrmoptionSepArg, 0 },
	{ "-" Nstars, "*" Nstars, XrmoptionSepArg, 0 },
	{ "-" Nseed, "*" Nseed, XrmoptionSepArg, 0 },
	{ "-" Nchunk, "*" Nchunk, XrmoptionSepArg, 0 },
	{ "-" Nlattice, "*" Nlattice, XrmoptionSepArg, 0 },
	{ "-" Nbounce, "*" Nbounce, XrmoptionNoArg, "True" },
	{ "-" Nvelocity, "*" Nvelocity, XrmoptionSepArg, 0 },
	{ "-" Naccel, "*" Naccel, XrmoptionSepArg, 0 },
	{ "-" Nlimit, "*" Nlimit, XrmoptionSepArg, 0 },
	{ "-" Ntimescale, "*" Ntimescale, XrmoptionSepArg, 0 },
};

XtAppContext app_context;	/* Xt application context */
Widget top_level, drawing;	/* global widgets */
ResData resdata;		/* application data */

XtIntervalId update_timeout_id;
XtWorkProcId update_workproc_id;

Display *display;		/* display of drawing widget */
int screen;			/* screen of drawing widget */
Window window;			/* window of drawing widget */
Colormap colormap;		/* colormap of drawing widget (set by init_colors()) */
XFontStruct *fontinfo;		/* info about label font */
GC gc;				/* drawing GC */
int width, height;		/* width and height of drawing window */

unsigned long black, white;	/* black and white colors */
unsigned long *colors;		/* array of pixel values for preallocated colors */
int ntemps;			/* number of surface temperatures in the colors array */
int nbrights;			/* number of brigtnesses per surface temperature */

int nstars;
star *stars;			/* array of stars */

double display_gamma;		/* gamma correction factor for display */
double v;			/* ship velocity (fraction of c) */
double g;			/* sqrt(1 - v^2) */
double a;			/* acceleration (1/s) */
double limit;			/* speed limit (fraction of c) */
double offset = 0.0;		/* time offset between start of simulation and start of acceleration (s) */
double elapsed = 0.0;		/* elapsed simulation time */
double origin = 0.0;		/* distance offset between start of simulation and start of acceleration (m) */
double dist = 0.0;		/* distance traveled (m) */
double timescale;		/* ratio of simulation time to real time */
double view_angle;		/* view angle of display window */
double view_plane;		/* view projection plane distance */
double spacing = 0.0;		/* spacing between stars in lattice */
int updating = 1;		/* are we doing periodic updates? */
int view = 0;			/* which view (forward, starboard, aft, port) */
int frames = 0;			/* frames drawn */
int paused = 0;			/* simulation paused? */
struct timeval lasttime;	/* last time display was updated */
XPoint *ptlist;			/* array of points to be drawn in batches */
XArc *arclist;			/* array of arcs to be drawn in batches */
int arcsize;			/* arcsize of stars (0 = draw points, >0 = draw circles) */

/* fill in X RGB color structure with color calculated from temperature T
   and brightness fraction f */
void color_from_temp(unsigned int T, double f, XColor *color)
{
	double red, green, blue, m;

	red = B(RED, T);
	green = B(GREEN, T);
	blue = B(BLUE, T);

	m = sqrt((red * red + green * green + blue * blue) / 3.0);
	red *= f / m;
	green *= f / m;
	blue *= f / m;
	
	if (red > 1.0) red = 1.0;
	if (green > 1.0) green = 1.0;
	if (blue > 1.0) blue = 1.0;

	color->red = (unsigned short) floor(65535.0 * pow(red, display_gamma));
	color->green = (unsigned short) floor(65535.0 * pow(green, display_gamma));
	color->blue = (unsigned short) floor(65535.0 * pow(blue, display_gamma));
}

/* set up colors and determine color allocation strategy based on visual type
   of display */
void init_colors(void)
{
	int ncolors = resdata.colors;
	int default_depth;
	XVisualInfo visual_info;
	int i, j, gotcolors, colors_allocated;
	XColor color;
	double temp;

	default_depth = DefaultDepth(display, screen);

	/* if the default depth is more than 8, try for TrueColor/DirectColor,
	   otherwise hope for PseudoColor */
	if (default_depth > 8) i = DirectColor;
	else i = PseudoColor;
	/* find the niftiest visual available at the default depth */
	while (i >= StaticGray && !XMatchVisualInfo(display, screen, default_depth, i, &visual_info)) i--;

	switch (visual_info.class) {
	case TrueColor:
	case DirectColor:
		/* system does not have a colormap, so initialize color array info
		   with zeros to indicate colors should be allocated on demand */
		ntemps = 0; nbrights = 0;
		colors = (unsigned long *) 0;
		/* allocate black and white explicitly (in case some
                   dork changes BlackPixel and WhitePixel) */
		color.red = color.green = color.blue = 0;
		XAllocColor(display, colormap, &color);
		black = color.pixel;
		color.red = color.green = color.blue = 65535;
		XAllocColor(display, colormap, &color);
		white = color.pixel;
		break;
	case PseudoColor:
	case StaticColor:
		/* system has a colormap, so try to allocate some colors for a
		   range of temperatures and brightnesses in the colormap */
		if (resdata.owncmap) ncolors = (1 << default_depth) - 2;
		if (ncolors > (1 << default_depth) - 2)
			/* limit ncolors to at most the number of colors in the colormap minus black and white */
			ncolors = (1 << default_depth) - 2;
		/* allocate black and white explicitly (in case some
                   dork changes BlackPixel and WhitePixel) */
		color.red = color.green = color.blue = 0;
		if (XAllocColor(display, colormap, &color)) black = color.pixel;
		else black = BlackPixel(display, screen);
		color.red = color.green = color.blue = 65535;
		if (XAllocColor(display, colormap, &color)) white = color.pixel;
		else white = WhitePixel(display, screen);
		colors_allocated = 0;
		/* attempt to allocate the requested number of colors until we get them
		   or fail completely */
		do {
			/* if there aren't enough colors, fail to monochrome mode */
			if (ncolors < 8) goto monochrome;
			/* try to find a number of brightnesses and temperatures that will fit
                           into the requested number of colors */
			nbrights = (int) sqrt(ncolors / 2.0);
			ntemps = ncolors / nbrights;
			ncolors = ntemps * nbrights;
			/* allocate color array */
			colors = calloc(ncolors, sizeof(*colors));
			gotcolors = 0;
			/* attempt to allocate the requested number of colors */
			for (i = 0; i < ntemps; i++) {
				temp = VIS_TEMP_MIN + ((VIS_TEMP_MAX - VIS_TEMP_MIN) * i) / (ntemps - 1);
				for (j = 0; j < nbrights; j++) {
					color_from_temp(temp, (double) (j + 1) / nbrights, &color);
					if (XAllocColor(display, colormap, &color)) {
						colors[i * nbrights + j] = color.pixel;
						gotcolors++;
					}
					else goto out_of_colors;
				}
			}
			out_of_colors:
			/* if we couldn't get all the colors requested */
			if (gotcolors < ncolors) {
				/* free all the colors allocated in this attempt */
				XFreeColors(display, colormap, colors, gotcolors, 0);
				/* free the color array */
				free(colors);
				printf("couldn't get %d colors, trying to make do with %d\n", ncolors, gotcolors);
				ncolors = gotcolors;
			}
			else colors_allocated = 1;
		} while (!colors_allocated);
		break;
	case GrayScale:
	case StaticGray:
		monochrome:
		ntemps = 1; nbrights = 2;
		ncolors = 2;
		black = BlackPixel(display, screen);
		white = WhitePixel(display, screen);
		break;
	}
}

unsigned long get_color(unsigned int T, int M)
{
	XColor color;

	/* clip M so it is less than VIS_MAG_RANGE */
	if (M < VIS_MAG_MAX) M = VIS_MAG_MAX;
	M = VIS_MAG_MIN - M;

	if (ntemps == 0) {
		/* DirectColor/TrueColor: compute color directly */
		if (M < 0) return black;
		color_from_temp(T, (double) M / VIS_MAG_RANGE, &color);
		XAllocColor(display, colormap, &color);
		return color.pixel;
	}
	else if (ntemps == 1) {
		/* monochrome: return black or white */
		if (M < 0) return black;
		else return white;
	}
	else {
		/* pseudocolor: find matching color in preallocated array */
		if (M < 0) return black;
		/* scale M to 0..nbrights-1 */
		M = (nbrights * M) / VIS_MAG_RANGE;
		if (M >= nbrights) M = nbrights - 1;
		/* clip T to (VIS_TEMP_MIN) .. (VIS_TEMP_MAX - 1) */
		if (T >= VIS_TEMP_MAX) T = VIS_TEMP_MAX - 1;
		else if (T < VIS_TEMP_MIN) T = VIS_TEMP_MIN;
		/* scale T to 0 .. (ntemps-1) */
		T = (ntemps * (T - VIS_TEMP_MIN)) / VIS_TEMP_RANGE;
		/* pick a color out of the preallocated array */
		return colors[T * nbrights + M];
	}
}

/* return a random real 0 <= x < 1 */
double random_unit(void)
{
	return random() / (double) (1U << 31);
}

/* create and fill in the array of star data */
void init_stars(void)
{
	/* try to allocate star data */
	nstars = resdata.stars;
	stars = (star *) calloc(nstars, sizeof(star));
	if (stars == 0) {
		printf("Can't make %d stars\n", nstars);
		exit(1);
	}
	if (spacing != 0.0) {
		/* create a cubic lattice of stars centered around the origin */
		int i, j, k, zsize, size, n;
		double x = 0.0, y = 0.0, z = 0.0;
		star s;

		/* compute the size of the cube edges */
		zsize = (int) rint(pow((double) nstars, 1 / 3.0));
		size = (int) floor(sqrt((double) nstars / zsize));
		/* compute nstars from the number of stars in the cube */
		nstars = zsize * size * size;
		s.T = 3000; s.M = 500;
		n = 0;
		for (k = zsize - 1; k >= 0; k -= 2) {
			s.z = z = k * C * YR * spacing / 2.0;
			/* blue (7000K) and red (3000K) stars alternate in layers */
			s.T = 10000 - s.T;
			for (j = size - 1; j >= 0; j -= 2) {
				s.y = y = j * C * YR * spacing / 2.0;
				for (i = size - 1; i >= 0; i -= 2) {
					s.x = x = i * C * YR * spacing / 2.0;
					stars[n++] = s;
					if (i > 0) {
						s.x = -x; stars[n++] = s; s.x = x;
					}
					if (j > 0) {
						s.y = -y; stars[n++] = s;
						if (i > 0) {
							s.x = -x; stars[n++] = s; s.x = x;
						}
						s.y = y;
					}
					if (k > 0) {
						s.z = -z; stars[n++] = s;
						if (i > 0) {
							s.x = -x; stars[n++] = s; s.x = x;
						}
						if (j > 0) {
							s.y = -y; stars[n++] = s;
							if (i > 0) {
								s.x = -x; stars[n++] = s; s.x = x;
							}
							s.y = y;
						}
						s.z = z;
					}
				}
			}
		}
	}
	else {
		/* create stars at random in a cylinder around the Z axis */
		int i;
		long seed = resdata.seed;

		if (seed == 0) seed = time(0);
		srandom(seed);
		printf("seed = %ld\n", seed);
		for (i = 0; i < nstars; i++) {
			if (i == 0) {
				/* make a Sun-like star as the first one */
				stars[0].T = 5800; stars[0].M = 477;
				stars[0].x = stars[0].y = stars[0].z = 0.0;
			}
			else {
				double t, m, r, theta;
				t = random_unit();
				m = random_unit();
				/* most of the time, magnitude correlates with temperature */
				if (random_unit() < 0.9) m = (3.0 * t + m) / 4.0;
				stars[i].T = (unsigned short) (TEMP_MIN + t * TEMP_RANGE);
				stars[i].M = (short) (ABS_MAG_MIN - m * ABS_MAG_RANGE);
				r = sqrt(random_unit()) * KNOWN_SPACE;
				theta = random_unit() * 2.0 * M_PI;
				stars[i].x = r * cos(theta); stars[i].y = r * sin(theta);
				stars[i].z = (2.0 * random_unit() - 1.0) * KNOWN_SPACE;
			}
		}
	}
}

/* draw a single star */
void draw_star(int x, int y)
{
	if (arcsize) XFillArc(display, window, gc, x - resdata.size, y - resdata.size, arcsize, arcsize, 0, 360 * 64);
	else XDrawPoint(display, window, gc, x, y);
}

/* plot n stars in the window */
void plot_stars(star *s, int n)
{
	for (; n--; s++) {
		double x = s->x, y = s->y, z = s->z, t, r, z_app, r_app, T_fact, M_diff, sx, sy;
		unsigned long color;

		r = sqrt(x * x + y * y + z * z);
		z_app = g * (v * r + z);
		/* handle rotations by swapping/negating coordinates */
		switch (view) {
		case 1:	t = x; x = -z_app; z_app = t; break; /* right window */
		case 2:	x = -x; z_app = -z_app; break; /* back window */
		case 3:	t = -x;	x = z_app; z_app = t; break; /* left window */
		}
		/* check if the star in front of us */
		if (z_app > 0.0) {
			sx = x * view_plane / z_app;
			sy = y * view_plane / z_app;
			/* check if the star is in window */
			if (fabs(sx) <= (width / 2) && fabs(sy) <= (height / 2)) {
				r_app = g * (r + v * z);
				T_fact = r_app / r;
				M_diff = -500.0 * (log10(T_fact / r) + 17.4893506303);
				color = get_color(T_fact * s->T, s->M + M_diff);
				/* only draw a star if it is not too dim to be seen */
				if (color != black) {
					XSetForeground(display, gc, color);
					s->sx = (short) sx + width / 2;
					s->sy = (short) sy + height / 2;
					if (arcsize) XFillArc(display, window, gc, s->sx - resdata.size, s->sy - resdata.size,
							      arcsize, arcsize, 0, 360 * 64);
					else XDrawPoint(display, window, gc, s->sx, s->sy);
					continue;
				}
			}
		}
		/* mark star as undrawn */
		s->sx = -1;
	}
}

/* erase n stars from the window */
void unplot_stars(star *s, int n)
{
	int ct = 0;

	XSetForeground(display, gc, black);
	if (arcsize) {
		for (; n--; s++) {
			if (s->sx >= 0) {
				/* if star was drawn in last update, add it to list of erasures */
				arclist[ct].x = s->sx - resdata.size;
				arclist[ct].y = s->sy - resdata.size;
				ct++;
			}
		}
		XFillArcs(display, window, gc, arclist, ct);
	}
	else {
		for (; n--; s++) {
			if (s->sx >= 0) {
				ptlist[ct].x = s->sx;
				ptlist[ct].y = s->sy;
				ct++;
			}
		}
		XDrawPoints(display, window, gc, ptlist, ct, CoordModeOrigin);
	}
}

/* show elapsed time, distance from origin, and frame velocity in window */
void show_stats(Bool draw)
{
	char msg[128];
	double t, d;
	char *ut, *ud;
	char *dirs[] = { "forward", "starboard", "aft", "port" };
	int y = fontinfo->ascent, row = fontinfo->ascent + fontinfo->descent;

	t = elapsed;
	if (fabs(t) < YR) {
		t /= DAY; ut = "days";
	}
	else {
		t /= YR; ut = "years";
	}
	d = dist;
	if (fabs(d) < C * YR) {
		d /= C * DAY; ud = "light-days";
	}
	else {
		d /= C * YR; ud = "light-years";
	}
	if (draw == True) XSetForeground(display, gc, white);
	else XSetForeground(display, gc, black);
	sprintf(msg, "%g %s", t, ut);
	XDrawString(display, window, gc, 0, y, msg, strlen(msg));
	y += row;
	sprintf(msg, "%g %s", d, ud);
	XDrawString(display, window, gc, 0, y, msg, strlen(msg));
	y += row;
	sprintf(msg, "%.9g c", v);
	XDrawString(display, window, gc, 0, y, msg, strlen(msg));
	y += row;
	sprintf(msg, "%s view", dirs[view]);
	XDrawString(display, window, gc, 0, y, msg, strlen(msg));
}

/* recalculate view parameters based on window size and view angle */
void update_view_plane(void)
{
	int min;

	if (width < height) min = width;
	else min = height;
	view_plane =  0.5 * min / tan(view_angle * M_PI / 360.0);
}

/* computes log((1 - x) / (1 + x)) */
double log1md1p(double x)
{
	return log1p(-v) - log1p(v);
}

/* computes cosh(x) - 1 */
double coshm1(double x)
{
	if (x < 0.25) {
		/* use series to compute cosh(x) - 1 accurately for very small values */
		int i = 0;
		double a = 1.0, f = 1.0, s = 0.0, s1;
		do {
			s1 = s;
			i += 2;
			a *= x * x;
			f *= i * (i - 1);
			s += a / f;
		} while (s != s1);
		return s;
	}
	else return (expm1(x) + expm1(-x)) / 2.0;
}

/* perform an incremental redraw of the display */
Bool update(XtPointer client_data)
{
	int i;
	double t, d;
	struct timeval thistime;

	/* recreate TimeOut if we are doing timed updates */
	if (resdata.fps > 0) {
		update_timeout_id = XtAppAddTimeOut(app_context, 1000 / resdata.fps, (XtTimerCallbackProc) &update, client_data);
	}
	/* erase stats from display */
	show_stats(False);
	/* calculate simulation time between last update and this one */
	gettimeofday(&thistime, 0);
	t = timescale * ((thistime.tv_sec - lasttime.tv_sec) * 1000000 + thistime.tv_usec - lasttime.tv_usec) / 1.0e6;
	elapsed += t;
	lasttime = thistime;
	/* if acceleration is non-zero, check velocity against limit */
	if ((a > 0.0 && v > limit) || (a < 0.0 && v < -limit)) {
		if (resdata.bounce == True) {
			/* reverse acceleration and recompute offset and origin */
			a = -a;
			offset = elapsed + 0.5 * log1md1p(v) / a;
			origin = C * coshm1(a * (elapsed - offset)) / a - dist;
		}
		else {
			/* stop acceleration */
			a = 0.0;
		}
	}
	/* if accelerating, compute new velocity and location and distance
	   moved since last update */
	if (a != 0.0) {
		v = tanh(a * (elapsed - offset));
		g = 1.0 / sqrt(1.0 - v * v);
		d = dist;
		dist = C * coshm1(a * (elapsed - offset)) / a - origin;
		d = dist - d;
	}
	/* else compute distance moved from current velocity */
	else {
		d = C * g * v * t;
		dist += d;
	}
	/* redraw new stats */
	show_stats(True);
	/* for all stars */
	for (i = 0; i < nstars; i += resdata.chunk) {
		int batch = resdata.chunk, j;
		
		if (nstars - i < batch)	batch = nstars - i;
		/* erase some stars */
		unplot_stars(&stars[i], batch);
		/* move stars by the distance traveled by the ship */
		for (j = 0; j < batch; j++) {
			stars[i+j].z -= d;
			if (spacing == 0.0) {
				/* wrap stars that fall off one end of the cylinder to the other end */
				if (stars[i+j].z < -KNOWN_SPACE) stars[i+j].z += 2.0 * KNOWN_SPACE;
				else if (stars[i+j].z > KNOWN_SPACE) stars[i+j].z -= 2.0 * KNOWN_SPACE;
			}
		}
		/* redraw those stars */
		plot_stars(&stars[i], batch);
	}
	frames++;
	/* return False to indicate this WorkProc shouldn't be removed */
	return False;
}

/* resume updating after a pause */
void start_updates(void)
{
	if (!updating) return;
	/* add update function as a WorkProc to be called when the Xt
	   event loop is idle */
	if (resdata.fps == 0) update_workproc_id = XtAppAddWorkProc(app_context, (XtWorkProc) &update, (XtPointer) drawing);
	/* add scheduled update as TimeOut to our update function */
	else update_timeout_id = XtAppAddTimeOut(app_context, 1000 / resdata.fps,
						 (XtTimerCallbackProc) &update, (XtPointer) drawing);
	/* reset time of last update */
	gettimeofday(&lasttime, 0);
}

/* stop updating to pause */
void stop_updates(void)
{
	if (!updating) return;
	if (resdata.fps == 0) XtRemoveWorkProc(update_workproc_id);
	else XtRemoveTimeOut(update_timeout_id);
}

/* load application resources, do some sanity checks, and calculate some initial values */
void load_resources(Widget w)
{
	/* load resources */
	XtVaGetApplicationResources(w, &resdata, resources, XtNumber(resources), 0);
	/* convert some string values to doubles */
	display_gamma = atof(resdata.gamma);
	view_angle = atof(resdata.viewangle);
	v = atof(resdata.velocity);
	limit = atof(resdata.limit);
	a = atof(resdata.accel) / C;
	timescale = atof(resdata.timescale);
	spacing = fabs(atof(resdata.lattice));
	if (display_gamma <= 0.0) {
		printf("Gamma really ought to be greater than 0.\n");
		exit(1);
	}
	if (view_angle <= 0.0 || view_angle >= 180.0) {
		printf("Please specify a view angle between 0 and 180 degrees.\n");
		exit(1);
	}
	resdata.fps = abs(resdata.fps);
	if (resdata.stars < 0) resdata.stars = -resdata.stars;
	if (resdata.stars == 0) resdata.stars = 1;
	if (resdata.chunk < 0) resdata.chunk = -resdata.chunk;
	if (resdata.chunk == 0) resdata.chunk = 1;
	if (resdata.size <= 1) arcsize = 0;
	else arcsize = 2 * resdata.size - 1;
	/* allocate array of XArcs/XPoints to do mass erasing of stars */
	if (arcsize) {
		int i;

		arclist = calloc(resdata.chunk, sizeof(XArc));
		if (arclist == 0) {
			printf("cannot allocate %d XArcs for chunking.\n", resdata.chunk);
			exit(1);
		}
		/* initialize XArcs with constants for width, height, start and end angles */
		for (i = 0; i < resdata.chunk; i++) {
			arclist[i].width = arclist[i].height = arcsize;
			arclist[i].angle1 = 0; arclist[i].angle2 = 360 * 64;
		}
	}
	else {
		ptlist = calloc(resdata.chunk, sizeof(XPoint));
		if (ptlist == 0) {
			printf("cannot allocate %d XPoints for chunking.\n", resdata.chunk);
			exit(1);
		}
	}
	if (fabs(v) >= 1.0 || fabs(limit) >= 1.0) {
		printf("299,792,458 meters per second:\nIt's not just a good idea, it's the law.\n");
		exit(1);
	}
	g = 1.0 / sqrt(1.0 - v * v);
	if (a != 0.0) {
		offset = 0.5 * log1md1p(v) / a;
		origin = C * coshm1(a * offset) / a;
	}
	if (resdata.colortest == True || timescale == 0.0) updating = 0;
	else updating = 1;
}

/* expose callback for expose() translation */
void expose(Widget widget, XEvent *event, String *parms, Cardinal *num_parms)
{
	int i, j;

	/* avoid multiple redraws for multiple contiguous expose events */
	if (event->type == Expose && event->xexpose.count != 0) return;
	/* clear window */
	XClearWindow(display, window);
	if (resdata.colortest == True) {
		/* draw color test pattern */
		for (i = 0; i < 64; i++) {
			int temp = TEMP_MIN + ((i - 1) * (TEMP_RANGE)) / 62;
			for (j = 0; j < 48; j++) {
				int mag = VIS_MAG_MIN - ((j - 1) * VIS_MAG_RANGE) / 46;
				XSetForeground(display, gc, get_color(temp, mag));
				draw_star((width * (2 * i + 1)) / 128, (height * (2 * j + 1)) / 96);
			}
		}
	}
	else {
		/* draw normal display */
		show_stats(True);
		plot_stars(stars, nstars);
	}
}

/* resize callback for resize() translation */
void resize(Widget widget, XEvent *event, String *parms, Cardinal *num_parms)
{
	int redraw;

	if (event->xconfigure.width > width || event->xconfigure.height > height) redraw = 0;
	/* window grew in size so expose event is generated after configure */
	else redraw = 1;
	/* window shrank and no expose event is generated; do it myself */

	width = event->xconfigure.width;
	height = event->xconfigure.height;
	update_view_plane();
	if (redraw) expose(widget, event, parms, num_parms);
}	  

/* callback for change_view() translation */
void change_view(Widget widget, XEvent *event, String *parms, Cardinal *num_parms)
{
	if (event->type == ButtonPress && event->xbutton.x < width / 2) view += 2;
	view = (view + 1) % 4;
	expose(widget, event, parms, num_parms);
}

/* callback for toggle_pause() translation */
void toggle_pause(Widget widget, XEvent *event, String *parms, Cardinal *num_parms)
{
	paused = !paused;
	if (paused) stop_updates();
	else start_updates();
}

/* callback for quit() translation */
void quit(Widget widget, XEvent *event, String *parms, Cardinal *num_parms)
{
	if (updating) printf("%.3g frames per second\n", timescale * frames / elapsed);
	exit(0);
}

/* action functions for translations */
XtActionsRec actions[] = {
	{ "expose", &expose },
	{ "resize", &resize },
	{ "toggle_pause", &toggle_pause },
	{ "change_view", &change_view },
	{ "quit", &quit },
};

/* translations for events */
char translation_table[] = 
"<Expose>: expose()\n"
"<ConfigureNotify>: resize()\n"
"<Btn1Down>: change_view()\n"
"<Btn2Down>: toggle_pause()\n"
"<Btn3Down>: quit()\n"
"<Key>r: change_view()\n"
"<Key>p: toggle_pause()\n"
"<Key>q: quit()\n";

int main(int argc, char *argv[])
{
	XWindowAttributes attributes;
	XGCValues gcvalues;

	/* create application context and Shell widget */
	top_level = XtVaAppInitialize(&app_context,"Xrel", options, XtNumber(options),
				      &argc, argv, fallback_resources,
				      XtNinput, True, 0);
	/* get display and screen we're using */
	display = XtDisplay(top_level);
	screen = XScreenNumberOfScreen(XtScreen(top_level));
	load_resources(top_level);
	if (argc > 1) {
		/* command line arguments were left over after Xt argument processing */
		printf("%s: simple relativistic space flight simulator\n\n", argv[0]);
		printf("These options are recognized in addition to the standard X toolkit options:\n");
		printf("-owncmap           use window-local colormap (%s)\n", resdata.owncmap == True ? "True" : "False");
		printf("-colortest         display color test pattern\n");
		printf("-max               resize to cover whole screen\n");
		printf("-colors [int]      maximum number of colors to allocate (%d)\n", resdata.colors);
		printf("-gamma [real]      gamma correction factor for display (%s)\n", resdata.gamma);
		printf("-size [int]        size (radius) of dots for stars (%d)\n", resdata.size);
		printf("-viewangle [real]  view angle of display (%s)\n", resdata.viewangle);
		printf("-fps [int]         frames per second to display [0 = maximum rate] (%d)\n", resdata.fps);
		printf("-stars [int]       number of stars in universe (%d)\n", resdata.stars);
		printf("-chunk [int]       number of stars to redraw in a chunk (%d)\n", resdata.chunk);
		printf("-seed [int]        seed for random star generation [0 = new seed] (%d)\n", resdata.seed);
		printf("-lattice [real]    lattice spacing for stars [0.0 = random generation] (%s)\n", resdata.lattice);
		printf("-velocity [real]   ship initial velocity as a fraction of c (%s)\n", resdata.velocity);
		printf("-accel [real]      ship acceleration in m/s^2 (%s)\n", resdata.accel);
		printf("-limit [real]      ship speed limit as fraction of c (%s)\n", resdata.limit);
		printf("-bounce            reverse acceleration at limit (%s)\n", resdata.bounce == True ? "True" : "False");
		printf("-timescale [real]  ratio of simulation time to real time (%s)\n", resdata.timescale);
		printf("\nWith the pointer in the window, press 'r' or left button to rotate view,\n");
		printf("'p' or middle button to pause/resume, 'q' or right button to quit.\n");
		exit(0);
	}
	/* if we want a private colormap, then create it and put it in
           the top_level colormap resource */
	if (resdata.owncmap == True) {
		colormap = XCreateColormap(display, RootWindow(display, screen), DefaultVisual(display, screen), AllocNone);
		XtVaSetValues(top_level, XtNcolormap, colormap, 0);
	}
	/* if the window should be maximum size for the screen, then
	   set its overrideRedirect resource */
	if (resdata.max == True) {
		XtVaSetValues(top_level, XtNoverrideRedirect, True, 0);
	}
	XtVaSetValues(top_level, XtNtranslations, XtParseTranslationTable(translation_table), 0);
	XtAppAddActions(app_context, actions, XtNumber(actions));
	/* create drawing window in shell widget */
	drawing = XtVaCreateManagedWidget("drawing", coreWidgetClass, top_level, 0);
	XtRealizeWidget(top_level);
	/* raise, resize, and install colormap for maximum size window */
	if (resdata.max == True) {
		XRaiseWindow(display, XtWindow(top_level));
		XMoveResizeWindow(display, XtWindow(top_level), 0, 0,
				  WidthOfScreen(XtScreen(top_level)), HeightOfScreen(XtScreen(top_level)));
		if (resdata.owncmap == True) XInstallColormap(display, colormap);
		XGrabKeyboard(display, XtWindow(top_level), True, GrabModeAsync, GrabModeAsync, CurrentTime);
	}
	window = XtWindow(drawing);
	/* get drawing window attributes */
	XGetWindowAttributes(display, window, &attributes);
	/* calculate initial view parameters */
	width = attributes.width;
	height = attributes.height;
	update_view_plane();
	/* allocate new colors in window's colormap (if necessary) */
	colormap = attributes.colormap;
	init_colors();
	/* set background to black allocated by init_colors */
	XtVaSetValues(top_level, XtNbackground, black, 0);
	XtVaSetValues(drawing, XtNbackground, black, 0);
	/* get info about font */
	fontinfo = XQueryFont(display, resdata.font);
	gcvalues.background = black;
	gcvalues.foreground = white;
	gcvalues.font = resdata.font;
	/* create drawing GC with our parameters */
	gc = XCreateGC(display, window, GCForeground|GCBackground|GCFont, &gcvalues);
	if (resdata.colortest == False) init_stars();
	start_updates();
	XtAppMainLoop(app_context);
	return 0;
}