new file mode 100644
--- /dev/null
+++ b/packages/services/gfx/mw/current/src/mwin/winlib/graph3d.c
@@ -0,0 +1,540 @@
+/*
+ * Copyright (c) 1999 Greg Haerr <greg@censoft.com>
+ *
+ * 3D Graphics Library for Micro-Windows
+ */
+#define MWINCLUDECOLORS
+#include "windows.h"
+#include "device.h"
+#include "graph3d.h"
+#define USEBLIT		1	/* =1 to use memDC's*/
+
+static int	nxpix;
+static int	nypix;
+static vec1	xscale;
+static vec1	yscale;
+static vec3	eye;
+static vec3	direct;
+static double 	Q[5][5];
+static HDC	hdc;
+static HDC	hdcMem;
+static HBITMAP	hbmp, hbmpOrg;
+
+/* setup eye, direction, calc observation matrix Q*/
+void
+look3(vec1 x, vec1 y, vec1 z)
+{
+	eye.x = x;
+	eye.y = y;
+	eye.z = z;
+	direct.x = -eye.x;
+	direct.y = -eye.y;
+	direct.z = -eye.z;
+	findQ();
+}
+
+void
+init3(HDC hDC, HWND memhwnd)
+{
+	HBRUSH	hbr;
+
+	hdc = hDC;
+	if(hdc) {
+		nxpix = hdc->hwnd->clirect.right - hdc->hwnd->clirect.left;
+		nypix = hdc->hwnd->clirect.bottom - hdc->hwnd->clirect.top;
+		xscale = (vec1)(nxpix-1) / nxpix * nxpix/2;
+		yscale = (vec1)(nypix-1) / nypix * nypix/2;
+
+		if(memhwnd) {
+			hdcMem = CreateCompatibleDC(NULL);
+			if(hdcMem) {
+				hbmp = CreateCompatibleBitmap(hdcMem,
+					nxpix, nypix);
+				hbmpOrg = SelectObject(hdcMem, hbmp);
+				hdc = hdcMem;
+			}
+			hbr = (HBRUSH)GetClassLong(memhwnd, GCL_HBRBACKGROUND);
+			FillRect(hdc, NULL, hbr);
+		}
+		/* create pen for setcolor3() color override*/
+		SelectObject(hdc, CreatePen(PS_SOLID, 1, BLACK));
+	}
+}
+
+void
+paint3(HDC hDC)
+{
+	if(hdcMem) {
+		BitBlt(hDC, 0, 0, nxpix, nypix, hdcMem, 0, 0, SRCCOPY);
+		DeleteObject(SelectObject(hdcMem, hbmpOrg));
+		DeleteDC(hdcMem);
+	}
+	hdcMem = NULL;
+}
+
+int
+fx(vec1 x)
+{
+	return (int)(x * xscale + nxpix*0.5 - 0.5);
+}
+
+int
+fy(vec1 y)
+{
+	return (int)(y * yscale + nypix*0.5 - 0.5);
+}
+
+void
+moveto3(vec2 pt)
+{
+	MoveToEx(hdc, fx(pt.x), fy(pt.y), NULL);
+}
+
+void
+setcolor3(MWCOLORVAL c)
+{
+	if(hdc)
+		hdc->pen->color = c;
+}
+
+void
+lineto3(vec2 pt)
+{
+	LineTo(hdc, fx(pt.x), fy(pt.y));
+}
+
+void
+polyfill(int n, vec2 points[])
+{
+	int	i;
+	int	xoff, yoff;
+	MWPOINT	pv[MAXPOLY];
+
+	if(!hdc)
+		return;
+
+	/* calc window offset*/
+	xoff = hdc->hwnd->clirect.left;
+	yoff = hdc->hwnd->clirect.top;
+
+	/* only plot non-trivial polygons*/
+	if(n > 2) {
+		for(i=0; i<n; ++i) {
+			pv[i].x = fx(points[i].x) + xoff;
+			pv[i].y = fy(points[i].y) + yoff;
+			/* fix: floating round error, y intercept difference
+			 * with GdLine
+			 */
+			/*pv[i].x = fx(points[i].x + xoff);*/
+			/*pv[i].y = fy(points[i].y + yoff);*/
+		}
+		GdSetForeground(GdFindColor(hdc->pen->color));
+		GdFillPoly(hdc->psd, n, pv);
+	}
+}
+
+void
+square(void)
+{
+	vec2	pt0, pt1, pt2, pt3;
+
+	pt0.x = -1; pt0.y = 1;
+	pt1.x = -1; pt1.y = -1;
+	pt2.x = 1; pt2.y = -1;
+	pt3.x = 1; pt3.y = 1;
+	moveto3(pt0);
+	lineto3(pt1);
+	lineto3(pt2);
+	lineto3(pt3);
+	lineto3(pt0);
+}
+
+void
+circle3(vec1 r)
+{
+	vec1 	theta = 0;
+	vec1	thinc = 2*pi/100;
+	int 	i;
+	vec2	pt;
+
+	pt.x = r;
+	pt.y = 0.0;
+	moveto3(pt);
+
+	for(i=0; i<100; ++i) {
+		theta = theta + thinc;
+		pt.x = r*cos(theta);
+		pt.y = r*sin(theta);
+		lineto3(pt);
+	}
+}
+
+void
+daisy(vec1 r,int points)
+{
+	int	i, j;
+	vec1	theta = 0;
+	vec1	thinc;
+	vec2	pt[100];
+
+	/* calculate n points on a circle*/
+	thinc = 2*pi/points;
+	for(i=0; i<points; ++i) {
+		pt[i].x = r*cos(theta);
+		pt[i].y = r*sin(theta);
+		theta += thinc;
+	}
+
+	/* join point i to point j for all 0 <= i < j < n */
+	for(i=0; i<points-1; ++i) {
+		for(j=i+1; j<points; ++j) {
+			moveto3(pt[i]);
+			lineto3(pt[j]);
+		}
+	}
+}
+
+void
+rose(vec1 r,int levels,int points)
+{
+	int	i, j, m, n;
+	vec1	r1, theta, thinc;
+	vec2	inner[100];
+	vec2	outer[100];
+	vec2	triangle[3];
+
+	m = levels;
+	n = points;
+	thinc = 2*pi/n;
+
+	/* initial inner circle*/
+	for(i=0; i<n; ++i) {
+		inner[i].x = 0.0;
+		inner[i].y = 0.0;
+	}
+
+	/* loop thru m levels*/
+	for(j=1; j<=m; ++j) {
+		theta = -j*pi/n;
+		r1 = r * (vec1)j/m;
+
+		/* calc n points on outer circle*/
+		for(i=0; i<n; ++i) {
+			theta += thinc;
+			outer[i].x = r1*cos(theta);
+			outer[i].y = r1*sin(theta);
+		}
+
+		/* construct/draw triangles with vertices on
+		 * inner and outer circles
+		 */
+		for(i=0; i<n; ++i) {
+			triangle[0] = outer[i];
+			triangle[1] = outer[(i+1) % n];
+			triangle[2] = inner[i];
+
+			/* fill triangle in red*/
+			setcolor3(RED);
+			polyfill(3, triangle);
+
+#if 1
+			/* outline triangle in white*/
+			setcolor3(WHITE);
+			moveto3(triangle[0]);
+			lineto3(triangle[1]);
+			lineto3(triangle[2]);
+			lineto3(triangle[0]);
+#endif
+		}
+
+		/* copy points on outer circle to inner arrays*/
+		for(i=0; i<n; ++i)
+			inner[i] = outer[i];
+	}
+}
+
+/* draw a triangle with cordners v0, v1, v2*/
+void
+triangle(vec2 v0, vec2 v1, vec2 v2)
+{
+	vec2	poly[3];
+
+	poly[0] = v0;
+	poly[1] = v1;
+	poly[2] = v2;
+
+	setcolor3(GREEN);
+	polyfill(3, poly);
+	setcolor3(BLACK);
+	moveto3(poly[2]);
+	lineto3(poly[0]);
+	lineto3(poly[1]);
+	lineto3(poly[2]);
+}
+
+/* draw a quadrilateral with corners v0, v1, v2, v3*/
+void
+quadrilateral(vec2 v0, vec2 v1, vec2 v2, vec2 v3)
+{
+	vec2	poly[4];
+
+	poly[0] = v0;
+	poly[1] = v1;
+	poly[2] = v2;
+	poly[3] = v3;
+	setcolor3(GREEN);
+	polyfill(4, poly);
+	setcolor3(BLACK);
+	moveto3(poly[3]);
+	lineto3(poly[0]);
+	lineto3(poly[1]);
+	lineto3(poly[2]);
+	lineto3(poly[3]);
+}
+
+/* find intersection of lines v0 to v1 and v2 to v3*/
+static int
+patch(vec2 v0, vec2 v1, vec2 v2, vec2 v3)
+{
+	vec1	denom;
+	vec1	mu;
+	vec2	v4;
+
+	denom = (v1.x-v0.x)*(v3.y-v2.y) - (v1.y-v0.y)*(v3.x-v2.x);
+	if(fabs(denom) > epsilon) {
+		mu = ((v2.x-v0.x)*(v3.y-v2.y) - (v2.y-v0.y)*(v3.x-v2.x))/denom;
+
+		/* if intersection between lines v0 to v1 and v2 to v3,
+		 * call it v4 and form triangles v0,v2,v4 and v1,v3,v4
+		 */
+		if(mu >= 0 && mu <= 1) {
+			v4.x = (1-mu)*v0.x + mu*v1.x;
+			v4.y = (1-mu)*v0.y + mu*v1.y;
+			triangle(v0, v2, v4);
+			triangle(v1, v3, v4);
+			return 0;
+		}
+	}
+
+	/* else find intersection of lines v0 to v2 and v1 to v3*/
+	denom = (v2.x-v0.x)*(v3.y-v1.y) - (v2.y-v0.y)*(v3.x-v1.x);
+	if(fabs(denom) > epsilon) {
+		mu = ((v1.x-v0.x)*(v3.y-v1.y) - (v1.y-v0.y)*(v3.x-v1.x))/denom;
+
+		/* if intersection between v0 and v1, call it v4
+		 * and form triangles v0,v1,v4 and v2,v3,v4
+		 */
+		if(mu >= 0 && mu <= 1) {
+			v4.x = (1-mu)*v0.x + mu*v2.x;
+			v4.y = (1-mu)*v0.y + mu*v2.y;
+			triangle(v0, v1, v4);
+			triangle(v2, v3, v4);
+			return 0;
+		}
+	}
+
+	/* there are no proper intersections so form quadrilateral v0,v1,v3,v2*/
+	quadrilateral(v0, v1, v3, v2);
+	return 1;
+}
+
+/* plotted function*/
+static vec1
+plotfn(vec1 x, vec1 z)
+{
+	vec1	t;
+
+	/* y = 4sin(sqrt(x*x+z*z))/sqrt(x*x+z*z) */
+	t = sqrt(x*x + z*z);
+	if(fabs(t) < epsilon)
+		return 4.0;
+	return 4.0 * sin(t) / t;
+}
+
+/* draw mathematical function plotfn*/
+void
+drawgrid(vec1 xmin, vec1 xmax, int nx, vec1 zmin, vec1 zmax, int nz)
+{
+	int	i, j;
+	vec1	xi, xstep, yij;
+	vec1	zj, zstep;
+	vec2	v[2][100];
+	double	S[5][5];
+
+	/* scale it down*/
+	scale3(1.0/(xmax-xmin)*2, 1.0/(xmax-xmin)*2, 1.0/(zmax-zmin), S);
+	mult3(Q, S, Q);
+
+	/* grid from xmin to xmax in nx steps and zmin to xmax in nz steps*/
+	xstep = (xmax-xmin)/nx;
+	zstep = (zmax-zmin)/nz;
+	xi = xmin;
+	zj = zmin;
+
+	/* calc grid points on first fixed-z line, fine the y-height
+	 * and transfrorm the points (xi,yij,zj) into observed
+	 * position.  Observed first set stored in v[0,1..nx]
+	 */
+	for(i=0; i<=nx; ++i) {
+		yij = plotfn(xi, zj);
+		v[0][i].x = Q[1][1]*xi + Q[1][2]*yij + Q[1][3]*zj;
+		v[0][i].y = Q[2][1]*xi + Q[2][2]*yij + Q[2][3]*zj;
+		xi += xstep;
+	}
+
+	/* run thru consecutive fixed-z lines (the second set)*/
+	for(j=0; j<nz; ++j) {
+		xi = xmin;
+		zj += zstep;
+
+		/* calc grid points on this second set, find the
+		 * y-height and transform the points (xi,yij,zj)
+		 * into observed position.  Observed second set
+		 * stored in v[1,0..nx]
+		 */
+		for(i=0; i<=nx; ++i) {
+			yij = plotfn(xi, zj);
+			v[1][i].x = Q[1][1]*xi + Q[1][2]*yij + Q[1][3]*zj;
+			v[1][i].y = Q[2][1]*xi + Q[2][2]*yij + Q[2][3]*zj;
+			xi += xstep;
+		}
+
+		/* run thru the nx patches formed by these two sets*/
+		for(i=0; i<nx; ++i)
+			patch(v[0][i], v[0][i+1], v[1][i], v[1][i+1]);
+
+		/* copy second set into first set*/
+		for(i=0; i<=nx; ++i)
+			v[0][i] = v[1][i];
+	}
+}
+
+/* returns the angle whose tangent is y/x.
+ * all anomalies such as x=0 are also checked
+ */
+vec1
+angle(vec1 x, vec1 y)
+{
+	if(fabs(x) < epsilon)
+		if(fabs(y) < epsilon)
+			return 0.0;
+		else 
+			if(y > 0.0)
+				return pi*0.5;
+			else return pi*1.5;
+	else
+		if(x < 0.0)
+			return atan(y/x) + pi;
+		else return atan(y/x);
+}
+
+/* calc 3d scaling matrix A giving scaling vector sx,sy,sz.
+ * one unit on the x axis becomes sx units, one unit on y, sy,
+ * and one unit on the z axis becomes sz units
+ */
+void
+scale3(vec1 sx, vec1 sy, vec1 sz, double A[][5])
+{
+	int	i, j;
+
+	for(i=1; i<5; ++i)
+		for(j=1; j<5; ++j)
+			A[i][j] = 0.0;
+	A[1][1] = sx;
+	A[2][2] = sy;
+	A[3][3] = sz;
+	A[4][4] = 1.0;
+}
+
+/* calc 3d axes translation matrix A
+ * origin translated by vectdor tx,ty,tz
+ */
+void
+tran3(vec1 tx, vec1 ty, vec1 tz, double A[][5])
+{
+	int	i, j;
+
+	for(i=1; i<5; ++i) {
+		for(j=1; j<5; ++j)
+			A[i][j] = 0.0;
+		A[i][i] = 1.0;
+	}
+	A[1][4] = -tx;
+	A[2][4] = -ty;
+	A[3][4] = -tz;
+}
+
+/* calc 3d axes rotation matrix A.  The axes are
+ * rotated anti-clockwise through an angle theta radians
+ * about an axis specified by m: m=1 means x, m=2 y, m=3 z axis
+ */
+void
+rot3(int m, vec1 theta, double A[][5])
+{
+	int	i, j, m1, m2;
+	vec1	c, s;
+
+	for(i=1; i<5; ++i) 
+		for(j=1; j<5; ++j)
+			A[i][j] = 0.0;
+	A[m][m] = 1.0;
+	A[4][4] = 1.0;
+	m1 = (m % 3) + 1;
+	m2 = (m1 % 3) + 1;
+	c = cos(theta);
+	s = sin(theta);
+	A[m1][m1] = c;
+	A[m2][m2] = c;
+	A[m1][m2] = s;
+	A[m2][m1] = s;
+}
+
+/* calc the matrix product C of two matrices A and B*/
+void
+mult3(double A[][5], double B[][5], double C[][5])
+{
+	int	i, j, k;
+	vec1	ab;
+
+	for(i=1; i<5; ++i) 
+		for(j=1; j<5; ++j) {
+			ab = 0;
+			for(k=1; k<5; ++k)
+				ab += A[i][k] * B[k][j];
+			C[i][j] = ab;
+		}
+}
+
+/* calc observation matrix Q for given observer*/
+void
+findQ(void)
+{
+	vec1	alpha, beta, gamma, v, w;
+	double	E[5][5];
+	double	F[5][5];
+	double	G[5][5];
+	double	H[5][5];
+	double	U[5][5];
+
+	/* calc translation matrix F*/
+	tran3(eye.x, eye.y, eye.z, F);
+
+	/* calc rotation matrix G*/
+	alpha = angle(-direct.x, -direct.y);
+	rot3(3, alpha, G);
+
+	/* calc rotation matrix H*/
+	v = sqrt(direct.x*direct.x + direct.y*direct.y);
+	beta = angle(-direct.z, v);
+	rot3(2, beta, H);
+
+	/* calc rotation matrix U*/
+	w = sqrt(v*v + direct.z*direct.z);
+	gamma = angle(-direct.x*w, direct.y*direct.z);
+	rot3(3, -gamma, U);
+
+	/* combine the transformations to find Q*/
+	mult3(G, F, Q);
+	mult3(H, Q, E);
+	mult3(U, E, Q);
+}