summaryrefslogtreecommitdiff
path: root/appl/wm/collide.b
diff options
context:
space:
mode:
Diffstat (limited to 'appl/wm/collide.b')
-rw-r--r--appl/wm/collide.b2180
1 files changed, 2180 insertions, 0 deletions
diff --git a/appl/wm/collide.b b/appl/wm/collide.b
new file mode 100644
index 00000000..1d8a2527
--- /dev/null
+++ b/appl/wm/collide.b
@@ -0,0 +1,2180 @@
+#
+# initially generated by c2l
+#
+
+implement Collide;
+
+include "draw.m";
+ draw: Draw;
+ Display, Image: import draw;
+
+Collide: module
+{
+ init: fn(nil: ref Draw->Context, argl: list of string);
+};
+
+include "sys.m";
+ sys: Sys;
+include "tk.m";
+ tk: Tk;
+ Toplevel: import tk;
+include "tkclient.m";
+ tkclient: Tkclient;
+include "math.m";
+ maths: Math;
+include "rand.m";
+ rand: Rand;
+include "daytime.m";
+ daytime: Daytime;
+include "bufio.m";
+include "arg.m";
+ arg: Arg;
+include "math/polyhedra.m";
+ polyhedra: Polyhedra;
+
+init(ctxt: ref Draw->Context, argl: list of string)
+{
+ sys = load Sys Sys->PATH;
+ draw = load Draw Draw->PATH;
+ tk = load Tk Tk->PATH;
+ tkclient = load Tkclient Tkclient->PATH;
+ maths = load Math Math->PATH;
+ rand = load Rand Rand->PATH;
+ arg = load Arg Arg->PATH;
+ daytime = load Daytime Daytime->PATH;
+ main(ctxt, argl);
+}
+
+π: con Math->Pi;
+∞: con real (1<<30);
+ε: con 0.001;
+√2: con 1.4142135623730950488016887242096980785696718753769486732;
+
+M1: con 1.0;
+M2: con 1.0;
+E: con 1.0; # coefficient of restitution/elasticity
+
+COLLIDE, REFLECT: con 1<<iota;
+
+MAXX, MAXY: con 512;
+
+RDisp: ref Draw->Image;
+black, white, red: ref Draw->Image;
+display: ref Draw->Display;
+toplev: ref Toplevel;
+
+Vector: adt{
+ x: real;
+ y: real;
+ z: real;
+};
+
+Line: adt{
+ a: Vector;
+ d: Vector; # normalized
+};
+
+Plane: adt{
+ id: int;
+ n: Vector; # normalized
+ d: real;
+ min: Vector;
+ max: Vector;
+ far: Vector;
+ v: array of Vector;
+};
+
+Object: adt{
+ id: int;
+ poly: ref Polyhedra->Polyhedron; # if any
+ c: ref Draw->Image; # colour
+ cb: ref Draw->Image; # border colour
+ l: ref Line; # initial point and direction
+ p: Vector; # current position
+ rp: Vector; # position after reflection
+ cp: Vector; # any collision point
+ rt: real; # time to reflection
+ ct: real; # time to collision
+ plane: ref Plane; # reflecting off
+ pmask: int; # plane mask
+ obj: cyclic ref Object; # colliding with
+ v: real; # speed
+ ω: real; # speed of rotation
+ roll: real; # roll
+ pitch: real; # pitch
+ yaw: real; # yaw
+ todo: int; # work to do
+};
+
+planes: list of ref Plane;
+
+V0: con Vector(real 0, real 0, real 0);
+VZ: con Vector(0.0, 0.0, 1.0);
+
+far: Vector;
+
+DOCIRCLE: con 1;
+POLY, FILLPOLY, CIRCLE, FILLCIRCLE, ELLIPSE, FILLELLIPSE: con iota;
+
+#
+# final object is centred on (0, 0, -objd)
+# viewer is at origin looking along (0 0 -1)
+#
+maxx, maxy: int;
+
+SCRW: con 320; # screen width
+SCRH: con 240; # screen height
+
+frac := 0.5; # % of screen for cube
+front := 0.5; # % of cube in front of screen
+hpar := 0.0; # horizontal parallax
+fov := -1.0; # field of view : 0 for parallel projection, -1 for unspecified
+objd := 500.0; # eye to middle of cube
+cubd := 100.0; # half side of cube
+icubd: real; # half side of inner cube
+icubd2: real; # square of above
+eyed := 32.0; # half eye to eye
+trkd := 5.0; # side/diameter of object
+trkd2: real; # square of above
+rpy := 0;
+roll := 0.0; # z
+pitch := 0.0; # y
+yaw := 0.0; # x
+
+scrd, objD, scrD: real; # screen distance
+left := 0; # left or right eye
+sx, sy, sz: real; # screen scale factors
+sf: real; # perspective scale factor
+fbpar: real; # -1 for front of cube, 1 for back
+vf := 1.0; # current velocity factor
+
+cmin, cmax: Vector; # cube extents
+
+# special transformation matrix without roll, pitch, yaw
+# this is needed so that spheres can be drawn as circles
+mod0 := array[4] of array of real;
+
+stereo := 0; # stereopsis
+
+SPHERE, ELLIPSOID, CUBE, POLYHEDRON: con iota;
+surr := CUBE; # surround
+
+poly := 0; # show polyhedra
+flat: int; # objects in one plane
+projx: int; # strange projection
+
+# ellipse parameters
+ef: Vector = (1.0, 0.8, 1.0);
+e2: Vector;
+
+# objects
+nobjs: int;
+objs: array of ref Object;
+me: ref Object;
+
+# circle drawing
+NC: con 72;
+cost, sint: array of real;
+
+# polyhedra
+polys: ref Polyhedra->Polyhedron;
+npolys: int;
+polyh: ref Polyhedra->Polyhedron;
+
+rgba(r: int, g: int, b: int, α: int): ref Image
+{
+ c := draw->setalpha((r<<24)|(g<<16)|(b<<8), α);
+ return display.newimage(((0, 0), (1, 1)), display.image.chans, 1, c);
+}
+
+random(a: int, b: int): int
+{
+ return a+rand->rand(b-a+1);
+}
+
+urand(): real
+{
+ M: con 1000;
+ return real random(0, M)/real M;
+}
+
+randomr(a: real, b: real): real
+{
+ return a+urand()*(b-a);
+}
+
+randomc(): ref Image
+{
+ r, g, b: int;
+
+ do{
+ r = random(0, 255);
+ g = random(0, 255);
+ b = random(0, 255);
+ }while(r+g+b < 384);
+ return rgba(r, g, b, 255);
+}
+
+randomv(a: real, b: real): Vector
+{
+ x := randomr(a, b);
+ y := randomr(a, b);
+ if(flat)
+ return (x, y, (a+b)/2.0);
+ return (x, y, randomr(a, b));
+}
+
+randomd(): Vector
+{
+ M: con 1000.0;
+ v := randomv(-M, M);
+ while(vlen(v) == 0.0)
+ v = randomv(-M, M);
+ return vnorm(v);
+}
+
+randomp(min: real, max: real): Vector
+{
+ case(surr){
+ SPHERE =>
+ return vmul(randomd(), max*maths->sqrt(urand()));
+ ELLIPSOID =>
+ return vmul(randomd(), max*vmin(ef)*maths->sqrt(urand()));
+ CUBE =>
+ return randomv(min, max);
+ * =>
+ v := randomv(min, max);
+ while(outside3(v, cmin, cmax))
+ v = randomv(min, max);
+ return v;
+ }
+}
+
+det(a: real, b: real, c: real, d: real): real
+{
+ return a*d-b*c;
+}
+
+simeq(a: real, b: real, c: real, d: real, e: real, f: real): (real, real)
+{
+ de := det(a, b, c, d);
+ return (det(e, b, f, d)/de, det(a, e, c, f)/de);
+}
+
+cksimeq(a: real, b: real, c: real, d: real, e: real, f: real): (int, real, real)
+{
+ ade := de := det(a, b, c, d);
+ if(ade < 0.0)
+ ade = -ade;
+ if(ade < ε)
+ return (0, 0.0, 0.0);
+ return (1, det(e, b, f, d)/de, det(a, e, c, f)/de);
+}
+
+ostring(o: ref Object): string
+{
+ return lstring(o.l) + "+" + vstring(o.p) + "+" + string o.v;
+}
+
+pstring(p: ref Plane): string
+{
+ return vstring(p.n) + "=" + string p.d;
+}
+
+lstring(l: ref Line): string
+{
+ return vstring(l.a) + "->" + vstring(l.d);
+}
+
+vstring(v: Vector): string
+{
+ return "(" + string v.x + " " + string v.y + " " + string v.z + ")";
+}
+
+vpt(x: real, y: real, z: real): Vector
+{
+ p: Vector;
+
+ p.x = x;
+ p.y = y;
+ p.z = z;
+ return p;
+}
+
+vdot(v1: Vector, v2: Vector): real
+{
+ return v1.x*v2.x+v1.y*v2.y+v1.z*v2.z;
+}
+
+vcross(v1: Vector, v2: Vector): Vector
+{
+ v: Vector;
+
+ v.x = v1.y*v2.z-v1.z*v2.y;
+ v.y = v1.z*v2.x-v1.x*v2.z;
+ v.z = v1.x*v2.y-v1.y*v2.x;
+ return v;
+}
+
+vadd(v1: Vector, v2: Vector): Vector
+{
+ v: Vector;
+
+ v.x = v1.x+v2.x;
+ v.y = v1.y+v2.y;
+ v.z = v1.z+v2.z;
+ return v;
+}
+
+vsub(v1: Vector, v2: Vector): Vector
+{
+ v: Vector;
+
+ v.x = v1.x-v2.x;
+ v.y = v1.y-v2.y;
+ v.z = v1.z-v2.z;
+ return v;
+}
+
+vmul(v1: Vector, s: real): Vector
+{
+ v: Vector;
+
+ v.x = s*v1.x;
+ v.y = s*v1.y;
+ v.z = s*v1.z;
+ return v;
+}
+
+vdiv(v1: Vector, s: real): Vector
+{
+ v: Vector;
+
+ v.x = v1.x/s;
+ v.y = v1.y/s;
+ v.z = v1.z/s;
+ return v;
+}
+
+vlen(v: Vector): real
+{
+ return maths->sqrt(vdot(v, v));
+}
+
+vlen2(v: Vector): (real, real)
+{
+ d2 := vdot(v, v);
+ d := maths->sqrt(d2);
+ return (d, d2);
+}
+
+vnorm(v: Vector): Vector
+{
+ d := maths->sqrt(vdot(v, v));
+ if(d == 0.0)
+ return v;
+ return vmul(v, real 1/d);
+}
+
+vnorm2(v: Vector): (real, Vector)
+{
+ d := maths->sqrt(vdot(v, v));
+ if(d == 0.0)
+ return (0.0, VZ);
+ return (d, vmul(v, real 1/d));
+}
+
+clip(x: real, d: real): real
+{
+ if(x < -d)
+ x = -d;
+ if(x > d)
+ x = d;
+ return x;
+}
+
+vclip(v: Vector, d: real): Vector
+{
+ c: Vector;
+
+ c.x = clip(v.x, d);
+ c.y = clip(v.y, d);
+ c.z = clip(v.z, d);
+ return c;
+}
+
+vout(v1: Vector, v2: Vector): int
+{
+ v := vsub(v2, v1);
+ return v.x < 0.0 || v.y < 0.0 || v.z < 0.0;
+}
+
+vmin(v: Vector): real
+{
+ m := v.x;
+ if(v.y < m)
+ m = v.y;
+ if(v.z < m)
+ m = v.z;
+ return m;
+}
+
+vvmul(v1: Vector, v2: Vector): Vector
+{
+ v: Vector;
+
+ v.x = v1.x*v2.x;
+ v.y = v1.y*v2.y;
+ v.z = v1.z*v2.z;
+ return v;
+}
+
+vvdiv(v1: Vector, v2: Vector): Vector
+{
+ v: Vector;
+
+ v.x = v1.x/v2.x;
+ v.y = v1.y/v2.y;
+ v.z = v1.z/v2.z;
+ return v;
+}
+
+vmuldiv(v1: Vector, v2: Vector, v3: Vector): real
+{
+ return vdot(vvdiv(v1, v3), v2);
+}
+
+newp(id: int, a: real, b: real, c: real, d: real, m: real, v: array of Vector)
+{
+ p := ref Plane;
+ p.id = id;
+ p.n = (a, b, c);
+ p.d = d;
+ m += ε;
+ p.min = (-m, -m, -m);
+ p.max = (+m, +m, +m);
+ p.v = v;
+ if(v != nil){
+ p.min = (∞, ∞, ∞);
+ p.max = (-∞, -∞, -∞);
+ for(i := 0; i < len v; i++){
+ vtx := v[i];
+ if(vtx.x < p.min.x)
+ p.min.x = vtx.x;
+ if(vtx.y < p.min.y)
+ p.min.y = vtx.y;
+ if(vtx.z < p.min.z)
+ p.min.z = vtx.z;
+ if(vtx.x > p.max.x)
+ p.max.x = vtx.x;
+ if(vtx.y > p.max.y)
+ p.max.y = vtx.y;
+ if(vtx.z > p.max.z)
+ p.max.z = vtx.z;
+ }
+ (x, y, z) := p.far = vmul(p.max, 2.0);
+ if(a != 0.0)
+ p.far.x = (d-b*y-c*z)/a;
+ else if(b != 0.0)
+ p.far.y = (d-c*z-a*x)/b;
+ else if(c != 0.0)
+ p.far.z = (d-a*x-b*y)/c;
+ else
+ fatal("null plane");
+ }
+ planes = p :: planes;
+}
+
+pinit()
+{
+ case(surr){
+ SPHERE or
+ ELLIPSOID =>
+ newp(0, 0.0, 0.0, 1.0, ∞, ∞, nil);
+ CUBE =>
+ newp(0, 1.0, 0.0, 0.0, -icubd, icubd, nil);
+ newp(1, 1.0, 0.0, 0.0, +icubd, icubd, nil);
+ newp(2, 0.0, 1.0, 0.0, -icubd, icubd, nil);
+ newp(3, 0.0, 1.0, 0.0, +icubd, icubd, nil);
+ newp(4, 0.0, 0.0, 1.0, -icubd, icubd, nil);
+ newp(5, 0.0, 0.0, 1.0, +icubd, icubd, nil);
+ * =>
+ p := polyh;
+ F := p.F;
+ v := p.v;
+ f := p.f;
+ fv := p.fv;
+ d := 0.0;
+ for(i := 0; i < F; i++){
+ n := vnorm(f[i]);
+ dn := vmul(n, cubd-icubd);
+ fvi := fv[i];
+ m := fvi[0];
+ av := array[m] of Vector;
+ for(j := 0; j < m; j++){
+ av[j] = vtx := vsub(vmul(v[fvi[j+1]], 2.0*cubd), dn);
+ d += vdot(n, vtx);
+ }
+ d /= real m;
+ newp(i, n.x, n.y, n.z, d, 0.0, av);
+ }
+ }
+}
+
+inside(v: Vector, vmin: Vector, vmax: Vector): int
+{
+ return !vout(vmin, v) && !vout(v, vmax);
+}
+
+inside2(v: Vector, p: ref Plane): int
+{
+ h := 0;
+ pt := p.far;
+ vs := p.v;
+ n := len p.v;
+ j := n-1;
+ for(i := 0; i < n; i++){
+ (ok, λ, μ) := cksimeq(vs[j].x-vs[i].x, v.x-pt.x, vs[j].y-vs[i].y, v.y-pt.y, v.x-vs[i].x, v.y-vs[i].y);
+ if(!ok)
+ (ok, λ, μ) = cksimeq(vs[j].y-vs[i].y, v.y-pt.y, vs[j].z-vs[i].z, v.z-pt.z, v.y-vs[i].y, v.z-vs[i].z);
+ if(!ok)
+ (ok, λ, μ) = cksimeq(vs[j].z-vs[i].z, v.z-pt.z, vs[j].x-vs[i].x, v.x-pt.x, v.z-vs[i].z, v.x-vs[i].x);
+ if(ok && μ >= 0.0 && λ >= 0.0 && λ < 1.0)
+ h++;
+ j = i;
+ }
+ return h&1;
+}
+
+inside3(v: Vector, lp: list of ref Plane): int
+{
+ h := 0;
+ l := ref Line;
+ l.a = v;
+ l.d = vnorm(vsub(far, v));
+ for( ; lp != nil; lp = tl lp){
+ (ok, nil, nil) := intersect(l, hd lp);
+ if(ok)
+ h++;
+ }
+ return h&1;
+}
+
+# outside of a face
+outside2(v: Vector, p: ref Plane): int
+{
+ if(surr == CUBE)
+ return vout(p.min, v) || vout(v, p.max);
+ else
+ return !inside2(v, p);
+}
+
+# outside of a polyhedron
+outside3(v: Vector, vmin: Vector, vmax: Vector): int
+{
+ case(surr){
+ SPHERE =>
+ return vout(vmin, v) || vout(v, vmax) || vdot(v, v) > icubd2 ;
+ ELLIPSOID =>
+ return vout(vmin, v) || vout(v, vmax) || vmuldiv(v, v, e2) > 1.0;
+ CUBE =>
+ return vout(vmin, v) || vout(v, vmax);
+ * =>
+ return !inside3(v, planes);
+ }
+}
+
+intersect(l: ref Line, p: ref Plane): (int, real, Vector)
+{
+ m := vdot(p.n, l.d);
+ if(m == real 0)
+ return (0, real 0, V0);
+ c := vdot(p.n, l.a);
+ λ := (p.d-c)/m;
+ if(λ < real 0)
+ return (0, λ, V0);
+ pt := vadd(l.a, vmul(l.d, λ));
+ if(outside2(pt, p))
+ return (0, λ, pt);
+ return (1, λ, pt);
+}
+
+reflection(tr: ref Object, lp: list of ref Plane)
+{
+ ok: int;
+ λ: real;
+
+ l := tr.l;
+ if(surr == SPHERE){
+ (ok, λ) = quadratic(1.0, 2.0*vdot(l.a, l.d), vdot(l.a, l.a)-icubd2);
+ if(!ok || λ < 0.0)
+ fatal("no sphere intersections");
+ tr.rp = vadd(l.a, vmul(l.d, λ));
+ tr.plane = hd lp; # anything
+ }
+ else if(surr == ELLIPSOID){
+ (ok, λ) = quadratic(vmuldiv(l.d, l.d, e2), 2.0*vmuldiv(l.a, l.d, e2), vmuldiv(l.a, l.a, e2)-1.0);
+ if(!ok || λ < 0.0)
+ fatal("no ellipsoid intersections");
+ tr.rp = vadd(l.a, vmul(l.d, λ));
+ tr.plane = hd lp; # anything
+ }
+ else{
+ p: ref Plane;
+ pt := V0;
+ λ = ∞;
+ for( ; lp != nil; lp = tl lp){
+ p0 := hd lp;
+ if((1<<p0.id)&tr.pmask)
+ continue;
+ (ok0, λ0, pt0) := intersect(l, p0);
+ if(ok0 && λ0 < λ){
+ λ = λ0;
+ p = p0;
+ pt = pt0;
+ }
+ }
+ if(λ == ∞)
+ fatal("no intersections");
+ tr.rp = pt;
+ tr.plane = p;
+ }
+ if(tr.v == 0.0)
+ tr.rt = ∞;
+ else
+ tr.rt = λ/tr.v;
+}
+
+reflect(tr: ref Object)
+{
+ l := tr.l;
+ if(surr == SPHERE)
+ n := vdiv(tr.rp, -icubd);
+ else if(surr == ELLIPSOID)
+ n = vnorm(vdiv(vvdiv(tr.rp, e2), -1.0));
+ else
+ n = tr.plane.n;
+ tr.l.a = tr.rp;
+ tr.l.d = vnorm(vsub(l.d, vmul(n, 2.0*vdot(n, l.d))));
+}
+
+impact(u2: real): (real, real)
+{
+ # u1 == 0
+ return simeq(M1, M2, -1.0, 1.0, M2*u2, -E*u2);
+}
+
+collision(t1: ref Object, t2: ref Object): (int, real, Vector, Vector)
+{
+ # stop t2
+ (v3, f) := vnorm2(vsub(vmul(t1.l.d, t1.v), vmul(t2.l.d, t2.v)));
+ if(v3 == 0.0)
+ return (0, 0.0, V0, V0);
+ ab := vsub(t2.p, t1.p);
+ (d, d2) := vlen2(ab);
+ cos := vdot(f, ab)/d;
+ cos2 := cos*cos;
+ if(cos < 0.0 || (disc := trkd2 - d2*(1.0-cos2)) < 0.0)
+ return (0, 0.0, V0, V0);
+ s := d*cos - maths->sqrt(disc);
+ t := s/v3;
+ s1 := t1.v*t;
+ s2 := t2.v*t;
+ cp1 := vadd(t1.p, vmul(t1.l.d, s1));
+ if(outside3(cp1, cmin, cmax))
+ return (0, 0.0, V0, V0);
+ cp2 := vadd(t2.p, vmul(t2.l.d, s2));
+ if(outside3(cp2, cmin, cmax))
+ return (0, 0.0, V0, V0);
+ return (1, t, cp1, cp2);
+}
+
+collisions(tr: ref Object)
+{
+ mincp1, mincp2: Vector;
+
+ n := nobjs;
+ t := objs;
+ tr0 := tr.obj;
+ mint := tr.ct;
+ tr1: ref Object;
+ for(i := 0; i < n; i++){
+ if((tr3 := t[i]) == tr || tr3 == tr0)
+ continue;
+ (c, tm, cp1, cp2) := collision(tr, tr3);
+ if(c && tm < mint && tm < tr3.ct){
+ mint = tm;
+ tr1 = tr3;
+ mincp1 = cp1;
+ mincp2 = cp2;
+ }
+ }
+ if(tr1 != nil){
+ tr.ct = mint;
+ tr1.ct = mint;
+ tr.obj = tr1;
+ tr2 := tr1.obj;
+ tr1.obj = tr;
+ tr.cp = mincp1;
+ tr1.cp = mincp2;
+ zerot(tr0, COLLIDE, 0);
+ zerot(tr2, COLLIDE, 0);
+ if(tr0 != nil && tr0 != tr2)
+ collisions(tr0);
+ if(tr2 != nil)
+ collisions(tr2);
+ }
+}
+
+collide(t1: ref Object, t2: ref Object)
+{
+ # stop t2
+ ov := vmul(t2.l.d, t2.v);
+ (v3, f) := vnorm2(vsub(vmul(t1.l.d, t1.v), ov));
+ ab := vsub(t2.cp, t1.cp);
+ α := vdot(f, ab)/vdot(ab, ab);
+ abr := vsub(f, vmul(ab, α));
+ (v2, v1) := impact(α*v3);
+ t1.l.a = t1.cp;
+ t2.l.a = t2.cp;
+ dir1 := vadd(vmul(ab, v1), vmul(abr, v3));
+ dir2 := vmul(ab, v2);
+ # start t2
+ (t1.v, t1.l.d) = vnorm2(vadd(dir1, ov));
+ (t2.v, t2.l.d) = vnorm2(vadd(dir2, ov));
+}
+
+deg2rad(d: real): real
+{
+ return π*d/180.0;
+}
+
+rad2deg(r: real): real
+{
+ return 180.0*r/π;
+}
+
+rp2d(r: real, p: real): Vector
+{
+ r = deg2rad(r);
+ cr := maths->cos(r);
+ sr := maths->sin(r);
+ p = deg2rad(p);
+ cp := maths->cos(p);
+ sp := maths->sin(p);
+ return (cr*cp, sr*cp, sp);
+}
+
+d2rp(v: Vector): (real, real)
+{
+ r := maths->atan2(v.y, v.x);
+ p := maths->asin(v.z);
+ return (rad2deg(r), rad2deg(p));
+}
+
+collideω(t1: ref Object, t2: ref Object)
+{
+ d1 := rp2d(t1.roll, t1.pitch);
+ d2 := rp2d(t2.roll, t2.pitch);
+ oω := vmul(d2, t2.ω);
+ (ω3, f) := vnorm2(vsub(vmul(d1, t1.ω), oω));
+ ab := vsub(t2.cp, t1.cp);
+ α := vdot(f, ab)/vdot(ab, ab);
+ abr := vsub(f, vmul(ab, α));
+ (ω2, ω1) := impact(α*ω3);
+ dir1 := vadd(vmul(ab, ω1), vmul(abr, ω3));
+ dir2 := vmul(ab, ω2);
+ (t1.ω, d1) = vnorm2(vadd(dir1, oω));
+ (t2.ω, d2) = vnorm2(vadd(dir2, oω));
+ (t1.roll, t1.pitch) = d2rp(d1);
+ (t2.roll, t2.pitch) = d2rp(d2);
+}
+
+plane(p1: Vector, p2: Vector, p3: Vector): (Vector, real)
+{
+ n := vnorm(vcross(vsub(p2, p1), vsub(p3, p1)));
+ d := vdot(n, p1);
+ return (n, d);
+}
+
+# angle subtended by the eyes at p in minutes
+angle(p: Vector): real
+{
+ l, r: Vector;
+
+ # left eye at (-eyed, 0, 0)
+ # right eye at (+eyed, 0, 0)
+ #
+ l = p;
+ l.x += eyed;
+ r = p;
+ r.x -= eyed;
+ return real 60*(real 180*maths->acos(vdot(l, r)/(maths->sqrt(vdot(l, l))*maths->sqrt(vdot(r, r))))/π);
+}
+
+# given coordinates relative to centre of cube
+disparity(p: Vector, b: Vector): real
+{
+ p.z -= objd;
+ b.z -= objd;
+ return angle(p)-angle(b);
+}
+
+# rotation about any of the axes
+# rotate(theta, &x, &y, &z) for x-axis
+# rotate(theta, &y, &z, &x) for y-axis
+# rotate(theta, &z, &x, &y) for z-axis
+#
+rotate(theta: int, x: real, y: real, z: real): (real, real, real)
+{
+ a := π*real theta/real 180;
+ c := maths->cos(a);
+ s := maths->sin(a);
+ oy := y;
+ oz := z;
+ y = c*oy-s*oz;
+ z = c*oz+s*oy;
+ return (x, y, z);
+}
+
+# solve the quadratic ax^2 + bx + c = 0, returning the larger root
+# * (a > 0)
+#
+quadratic(a: real, b: real, c: real): (int, real)
+{
+ d := b*b-real 4*a*c;
+ if(d < real 0)
+ return (0, 0.0); # no real roots
+ x := (maths->sqrt(d)-b)/(real 2*a);
+ return (1, x);
+}
+
+# calculate the values of objD, scrD given the required parallax
+dopar()
+{
+ a := real 1;
+ b, c: real;
+ f := real 2*front-real 1;
+ x: real;
+ s: int;
+ w := sx*real SCRW;
+ ok: int;
+
+ if(hpar == 0.0){ # natural parallax
+ objD = objd;
+ scrD = scrd;
+ return;
+ }
+ if(fbpar < f)
+ s = -1;
+ else
+ s = 1;
+ if(fbpar == f)
+ fatal("parallax value is zero at screen distance");
+ b = (fbpar+f)*cubd-(fbpar-f)*w*eyed*real s*frac/hpar;
+ c = fbpar*f*cubd*cubd;
+ (ok, x) = quadratic(a, b, c);
+ if(ok){
+ objD = x;
+ scrD = x+f*cubd;
+ if(objD > real 0 && scrD > real 0)
+ return;
+ }
+ fatal("unachievable parallax value");
+}
+
+# update graphics parameters
+update(init: int)
+{
+ if(fov != real 0){
+ if(objd == real 0)
+ fov = 180.0;
+ else
+ fov = real 2*(real 180*maths->atan(cubd/(frac*objd))/π);
+ }
+ scrd = objd+(real 2*front-real 1)*cubd;
+ if(init){
+ if(objd < ε)
+ objd = ε;
+ if(fov != real 0)
+ sf = real (1<<2)*cubd/(objd*frac);
+ else
+ sf = cubd/frac;
+ }
+ # dopar();
+ domats();
+}
+
+fovtodist()
+{
+ if(fov != real 0)
+ objd = cubd/(frac*maths->tan(π*(fov/real 2)/real 180));
+}
+
+getpolys()
+{
+ (n, p, b) := polyhedra->scanpolyhedra("/lib/polyhedra");
+ polyhedra->getpolyhedra(p, b);
+ polys = p;
+ npolys = n;
+ do{
+ for(i := 0; i < p.V; i++)
+ p.v[i] = vmul(p.v[i], 0.5);
+ for(i = 0; i < p.F; i++)
+ p.f[i] = vmul(p.f[i], 0.5);
+ p = p.nxt;
+ }while(p != polys);
+}
+
+randpoly(p: ref Polyhedra->Polyhedron, n: int): ref Polyhedra->Polyhedron
+{
+ r := random(0, n-1);
+ for( ; --r >= 0; p = p.nxt)
+ ;
+ return p;
+}
+
+drawpoly(p: ref Polyhedra->Polyhedron, typex: int)
+{
+ # V := p.V;
+ F := p.F;
+ v := p.v;
+ # f := p.f;
+ fv := p.fv;
+ for(i := 0; i < F; i++){
+ fvi := fv[i];
+ n := fvi[0];
+ m_begin(typex, n);
+ for(j := 0; j < n; j++){
+ vtx := v[fvi[j+1]];
+ m_vertex(vtx.x, vtx.y, vtx.z);
+ }
+ m_end();
+ }
+}
+
+# objects with unit sides/diameter
+H: con 0.5;
+
+square(typex: int)
+{
+ m_begin(typex, 4);
+ m_vertex(-H, -H, 0.0);
+ m_vertex(-H, +H, 0.0);
+ m_vertex(+H, +H, 0.0);
+ m_vertex(+H, -H, 0.0);
+ m_end();
+}
+
+diamond(typex: int)
+{
+ m_pushmatrix();
+ m_rotatez(45.0);
+ square(typex);
+ m_popmatrix();
+}
+
+circleinit()
+{
+ i: int;
+ a := 0.0;
+ cost = array[NC] of real;
+ sint = array[NC] of real;
+ for(i = 0; i < NC; i++){
+ cost[i] = H*maths->cos(a);
+ sint[i] = H*maths->sin(a);
+ a += (2.0*π)/real NC;
+ }
+}
+
+circle(typex: int)
+{
+ i: int;
+
+ if(DOCIRCLE){
+ m_begin(typex, 2);
+ m_circle(0.0, 0.0, 0.0, 0.5);
+ m_end();
+ return;
+ }
+ else{
+ m_begin(typex, NC);
+ for(i = 0; i < NC; i++)
+ m_vertex(cost[i], sint[i], 0.0);
+ m_end();
+ }
+}
+
+ellipse(typex: int)
+{
+ m_begin(typex, 4);
+ m_ellipse(0.0, 0.0, 0.0, vmul(ef, 0.5));
+ m_end();
+}
+
+hexahedron(typex: int)
+{
+ i, j, k: int;
+ V := array[8] of {
+ array[3] of {
+ -H, -H, -H,
+ },
+ array[3] of {
+ -H, -H, +H,
+ },
+ array[3] of {
+ -H, +H, -H,
+ },
+ array[3] of {
+ -H, +H, +H,
+ },
+ array[3] of {
+ +H, -H, -H,
+ },
+ array[3] of {
+ +H, -H, +H,
+ },
+ array[3] of {
+ +H, +H, -H,
+ },
+ array[3] of {
+ +H, +H, +H,
+ },
+ };
+ F := array[6] of {
+ array[4] of {
+ 0, 4, 6, 2,
+ },
+ array[4] of {
+ 0, 4, 5, 1,
+ },
+ array[4] of {
+ 0, 1, 3, 2,
+ },
+ array[4] of {
+ 1, 5, 7, 3,
+ },
+ array[4] of {
+ 2, 6, 7, 3,
+ },
+ array[4] of {
+ 4, 5, 7, 6,
+ },
+ };
+
+ for(i = 0; i < 6; i++){
+ m_begin(typex, 4);
+ for(j = 0; j < 4; j++){
+ k = F[i][j];
+ m_vertex(V[k][0], V[k][1], V[k][2]);
+ }
+ m_end();
+ }
+}
+
+zerot(tr: ref Object, zero: int, note: int)
+{
+ if(tr != nil){
+ if(zero&REFLECT){
+ tr.rt = ∞;
+ tr.plane = nil;
+ }
+ if(zero&COLLIDE){
+ tr.ct = ∞;
+ tr.obj = nil;
+ }
+ if(note)
+ tr.todo = zero;
+ }
+}
+
+newobj(t: array of ref Object, n: int, vel: int, velf: real): ref Object
+{
+ tr: ref Object;
+ p1: Vector;
+ again: int;
+
+ d := icubd-1.0;
+ cnt := 1024;
+ do{
+ p1 = randomp(-d, d);
+ again = 0;
+ for(i := 0; i < n; i++){
+ (nil, d2) := vlen2(vsub(t[i].p, p1));
+ if(d2 <= trkd2){
+ again = 1;
+ break;
+ }
+ }
+ cnt--;
+ }while(again && cnt > 0);
+ if(again)
+ return nil;
+ # p2 := randomp(-d, d);
+ p21 := randomd();
+ tr = ref Object;
+ tr.id = n;
+ tr.poly = nil;
+ if(poly){
+ if(n == 0)
+ tr.poly = randpoly(polys, npolys);
+ else
+ tr.poly = t[0].poly;
+ }
+ tr.c = randomc();
+ tr.cb = tr.c; # randomc();
+ if(vel)
+ tr.v = vf*velf*randomr(0.5, 2.0);
+ else
+ tr.v = 0.0;
+ tr.ω = vf*randomr(1.0, 10.0);
+ tr.roll = randomr(0.0, 360.0);
+ tr.pitch = randomr(0.0, 360.0);
+ tr.yaw = randomr(0.0, 360.0);
+ tr.l = ref Line(p1, vnorm(p21));
+ tr.p = p1;
+ tr.todo = 0;
+ zerot(tr, REFLECT|COLLIDE, 0);
+ tr.pmask = 0;
+ reflection(tr, planes);
+ return tr;
+}
+
+objinit(m: int, v: int)
+{
+ velf := real m/real v;
+ p := nobjs;
+ n := p+m;
+ v += p;
+ t := array[n] of ref Object;
+ for(i := 0; i < p; i++)
+ t[i] = objs[i];
+ for(i = p; i < n; i++){
+ t[i] = newobj(t, i, i < v, velf);
+ if(t[i] == nil)
+ return;
+ }
+ sort(t, n);
+ nobjs = n;
+ objs = t;
+ for(i = p; i < n; i++)
+ collisions(t[i]);
+}
+
+zobj: Object;
+
+newo(n: int, p: Vector, c: ref Draw->Image): ref Object
+{
+ o := ref Object;
+ *o = zobj;
+ o.id = n;
+ o.c = o.cb = c;
+ o.l = ref Line(p, VZ);
+ o.p = p;
+ zerot(o, REFLECT|COLLIDE, 0);
+ reflection(o, planes);
+ return o;
+}
+
+objinits(nil: int, nil: int)
+{
+ n := 16;
+ t := array[n] of ref Object;
+ r := trkd/2.0;
+ i := 0;
+ yc := 0.0;
+ for(y := 0; y < 5; y++){
+ xc := -real y*r;
+ for(x := 0; x <= y; x++){
+ t[i] = newo(i, (xc, yc, 0.0), red);
+ xc += trkd;
+ i++;
+ }
+ yc += trkd;
+ }
+ t[i] = newo(i, (0.0, -50.0, 0.0), white);
+ t[i].l.d = (0.0, 1.0, 0.0);
+ t[i].v = 1.0;
+ sort(t, n);
+ nobjs = n;
+ objs = t;
+ for(i = 0; i < n; i++)
+ collisions(t[i]);
+}
+
+initme(): ref Object
+{
+ t := newobj(nil, 0, 1, 1.0);
+ t.roll = t.pitch = t.yaw = 0.0;
+ t.v = t.ω = 0.0;
+ t.l.a = (0.0, 0.0, objd); # origin when translated
+ t.l.d = (0.0, 0.0, -1.0);
+ t.p = t.l.a;
+ zerot(t, REFLECT|COLLIDE, 0);
+ return t;
+}
+
+retime(s: real)
+{
+ r := 1.0/s;
+ n := nobjs;
+ t := objs;
+ for(i := 0; i < n; i++){
+ tr := t[i];
+ tr.v *= s;
+ tr.ω *= s;
+ tr.rt *= r;
+ tr.ct *= r;
+ }
+ me.v *= s;
+ me.ω *= s;
+ me.rt *= r;
+ me.ct *= r;
+}
+
+drawobjs()
+{
+ tr: ref Object;
+ p: Vector;
+
+ n := nobjs;
+ t := objs;
+
+ for(i := 0; i < n; i++){
+ tr = t[i];
+ tr.pmask = 0;
+ p = tr.p;
+ m_pushmatrix();
+ if(rpy && tr.poly == nil){
+ m_loadmatrix(mod0);
+ (p.x, p.y, p.z) = rotate(int yaw, p.x, p.y, p.z);
+ (p.y, p.z, p.x) = rotate(int pitch, p.y, p.z, p.x);
+ (p.z, p.x, p.y) = rotate(int roll, p.z, p.x, p.y);
+ }
+ m_translate(p.x, p.y, p.z);
+ m_scale(trkd, trkd, trkd);
+ if(tr.poly != nil){
+ m_rotatez(tr.roll);
+ m_rotatey(tr.pitch);
+ m_rotatex(tr.yaw);
+ tr.yaw += tr.ω;
+ }
+ m_matmul();
+ if(tr.cb != tr.c){
+ m_colour(tr.cb);
+ if(tr.poly != nil)
+ drawpoly(tr.poly, POLY);
+ else if(DOCIRCLE)
+ circle(CIRCLE);
+ else
+ circle(POLY);
+ }
+ m_colour(tr.c);
+ if(tr.poly != nil)
+ drawpoly(tr.poly, FILLPOLY);
+ else if(DOCIRCLE)
+ circle(FILLCIRCLE);
+ else
+ circle(FILLPOLY);
+ m_popmatrix();
+ }
+
+ tm := 1.0;
+ do{
+ δt := ∞;
+
+ for(i = 0; i < n; i++){
+ tr = t[i];
+ if(tr.rt < δt)
+ δt = tr.rt;
+ if(tr.ct < δt)
+ δt = tr.ct;
+ }
+
+ if(δt > tm)
+ δt = tm;
+
+ for(i = 0; i < n; i++){
+ tr = t[i];
+ if(tr.rt == δt){
+ tr1 := tr.obj;
+ reflect(tr);
+ tr.p = tr.rp;
+ if(δt > 0.0)
+ tr.pmask = (1<<tr.plane.id);
+ else
+ tr.pmask |= (1<<tr.plane.id);
+ zerot(tr, REFLECT|COLLIDE, 1);
+ zerot(tr1, COLLIDE, 1);
+ }
+ else if(tr.ct == δt){
+ tr1 := tr.obj ;
+ collide(tr, tr1);
+ if(0 && poly)
+ collideω(tr, tr1);
+ tr.p = tr.cp;
+ tr1.p = tr1.cp;
+ tr.pmask = tr1.pmask = 0;
+ zerot(tr, REFLECT|COLLIDE, 1);
+ zerot(tr1, REFLECT|COLLIDE, 1);
+ }
+ else if(tr.todo != (REFLECT|COLLIDE)){
+ tr.p = vclip(vadd(tr.p, vmul(tr.l.d, tr.v*δt)), icubd);
+ tr.rt -= δt;
+ tr.ct -= δt;
+ }
+ }
+
+ for(i = 0; i < n; i++){
+ tr = t[i];
+ if(tr.todo){
+ if(tr.todo&REFLECT)
+ reflection(tr, planes);
+ if(tr.todo&COLLIDE)
+ collisions(tr);
+ tr.todo = 0;
+ }
+ }
+
+ tm -= δt;
+
+ }while(tm > 0.0);
+
+ sort(t, n);
+}
+
+drawscene()
+{
+ m_pushmatrix();
+ m_scale(real 2*cubd, real 2*cubd, real 2*cubd);
+ m_colour(white);
+ m_matmul();
+ case(surr){
+ SPHERE =>
+ if(DOCIRCLE)
+ circle(CIRCLE);
+ else
+ circle(POLY);
+ ELLIPSOID =>
+ ellipse(ELLIPSE);
+ CUBE =>
+ if(flat)
+ square(POLY);
+ else
+ hexahedron(POLY);
+ * =>
+ drawpoly(polyh, POLY);
+ }
+ m_popmatrix();
+ drawobjs();
+}
+
+# ensure parallax doesn't alter between images
+adjpar(x: array of real, y: array of real, z: array of real)
+{
+ zed, incr: real;
+
+ y = nil;
+ if(z[0] < real 0)
+ zed = -z[0];
+ else
+ zed = z[0];
+ incr = eyed*zed*(real 1-scrD/(zed+objD-objd))/scrd;
+ if(!stereo || fov == real 0)
+ return;
+ if(left)
+ x[0] -= incr;
+ else
+ x[0] += incr;
+}
+
+view()
+{
+ m_mode(PROJ);
+ m_loadidentity();
+ m_scale(sx, sy, sz);
+ if(fov != real 0)
+ m_frustum(sf, real (1<<2), real (1<<20));
+ else
+ m_ortho(sf, real (1<<2), real (1<<20));
+ # m_print();
+ m_mode(MODEL);
+}
+
+model(rot: int)
+{
+ m_loadidentity();
+ m_translate(0.0, 0.0, -objd);
+ if(rpy && rot){
+ m_rotatez(roll);
+ m_rotatey(pitch);
+ m_rotatex(yaw);
+ }
+}
+
+# store projection and modelview matrices
+domats()
+{
+ model(0);
+ m_storematrix(mod0);
+ model(1);
+ view();
+}
+
+scale()
+{
+ if(maxx > maxy)
+ sx = real maxy/real maxx;
+ else
+ sx = 1.0;
+ if(maxy > maxx)
+ sy = real maxx/real maxy;
+ else
+ sy = 1.0;
+ sz = 1.0;
+}
+
+rescale(w: int, h: int)
+{
+ maxx = w;
+ maxy = h;
+ scale();
+ m_viewport(0, 0, maxx, maxy);
+}
+
+grinit()
+{
+ for(i := 0; i < 4; i++)
+ mod0[i] = array[4] of real;
+ far = (2.0*cubd, 2.0*cubd, 2.0*cubd);
+ icubd = cubd-trkd/2.0;
+ icubd2 = icubd*icubd;
+ trkd2 = trkd*trkd;
+ cmin = vpt(-icubd, -icubd, -icubd);
+ cmax = vpt(+icubd, +icubd, +icubd);
+ maxx = MAXX;
+ maxy = MAXY;
+ e2 = vmul(vvmul(ef, ef), icubd2);
+
+ m_init();
+ pinit();
+ circleinit();
+
+ m_viewport(0, 0, maxx, maxy);
+
+ scale();
+ if(fov > real 0)
+ fovtodist();
+ update(1);
+}
+
+newimage(win: ref Toplevel, init: int)
+{
+ maxx = int cmd(win, ".p cget -actwidth");
+ maxy = int cmd(win, ".p cget -actheight");
+ RDisp = display.newimage(((0, 0), (maxx, maxy)), win.image.chans, 0, Draw->Black);
+ tk->putimage(win, ".p", RDisp, nil);
+ RDisp.draw(RDisp.r, black, nil, (0, 0));
+ reveal();
+ rescale(maxx, maxy);
+ update(init);
+}
+
+reveal()
+{
+ cmd(toplev, ".p dirty; update");
+}
+
+usage()
+{
+ sys->fprint(sys->fildes(2), "usage: collide [-cse] [-f] [-op] [-b num] [-v num]\n");
+ exit;
+}
+
+main(ctxt: ref Draw->Context, args: list of string)
+{
+ rand->init(daytime->now());
+ daytime = nil;
+
+ b := v := random(16, 32);
+
+ arg->init(args);
+ while((o := arg->opt()) != 0){
+ case(o){
+ * =>
+ usage();
+ 's' =>
+ surr = SPHERE;
+ 'e' =>
+ surr = ELLIPSOID;
+ 'c' =>
+ surr = CUBE;
+ 'o' =>
+ fov = 0.0;
+ 'p' =>
+ fov = -1.0;
+ 'f' =>
+ flat = 1;
+ 'b' =>
+ b = v = int arg->arg();
+ 'v' =>
+ v = int arg->arg();
+ }
+ }
+ if(arg->argv() != nil)
+ usage();
+
+ if(b <= 0)
+ b = 1;
+ if(b > 100)
+ b = 100;
+ if(v <= 0)
+ v = 1;
+ if(v > b)
+ v = b;
+
+ if(poly || surr == POLYHEDRON){
+ polyhedra = load Polyhedra Polyhedra->PATH;
+ getpolys();
+ }
+ if(surr == POLYHEDRON)
+ polyh = randpoly(polys, npolys);
+
+ grinit();
+
+ tkclient->init();
+ (win, wmch) := tkclient->toplevel(ctxt, "", "Collide", Tkclient->Resize | Tkclient->Hide);
+ toplev = win;
+ sys->pctl(Sys->NEWPGRP, nil);
+ cmdch := chan of string;
+ tk->namechan(win, cmdch, "cmd");
+ for(i := 0; i < len winconfig; i++)
+ cmd(win, winconfig[i]);
+ cmd(win, "update");
+
+ tkclient->onscreen(win, nil);
+ tkclient->startinput(win, "kbd"::"ptr"::nil);
+
+ display = win.image.display;
+ newimage(win, 1);
+
+ black = display.color(Draw->Black);
+ white = display.color(Draw->White);
+ red = display.color(Draw->Red);
+
+ objinit(b, v);
+ me = initme();
+
+ pid := -1;
+ sync := chan of int;
+ cmdc := chan of int;
+ spawn animate(sync, cmdc);
+ pid = <- sync;
+
+ for(;;){
+ alt{
+ c := <-win.ctxt.kbd =>
+ tk->keyboard(win, c);
+ c := <-win.ctxt.ptr =>
+ tk->pointer(win, *c);
+ c := <-win.ctxt.ctl or
+ c = <-win.wreq =>
+ tkclient->wmctl(win, c);
+ c := <- wmch =>
+ case c{
+ "exit" =>
+ if(pid != -1)
+ kill(pid);
+ exit;
+ * =>
+ sync <-= 0;
+ tkclient->wmctl(win, c);
+ if(c[0] == '!')
+ newimage(win, 0);
+ sync <-= 1;
+ }
+ c := <- cmdch =>
+ case c{
+ "stop" =>
+ cmdc <-= 's';
+ "zoomin" =>
+ cmdc <-= 'z';
+ "zoomout" =>
+ cmdc <-= 'o';
+ "slow" =>
+ cmdc <-= '-';
+ "fast" =>
+ cmdc <-= '+';
+ "objs" =>
+ sync <-= 0;
+ b >>= 1;
+ if(b == 0)
+ b = 1;
+ objinit(b, b);
+ sync <-= 1;
+ }
+ }
+ }
+}
+
+sign(r: real): real
+{
+ if(r < 0.0)
+ return -1.0;
+ return 1.0;
+}
+
+abs(r: real): real
+{
+ if(r < 0.0)
+ return -r;
+ return r;
+}
+
+animate(sync: chan of int, cmd: chan of int)
+{
+ zd := objd/250.0;
+ δ := θ := 0.0;
+ f := 8;
+
+ sync <- = sys->pctl(0, nil);
+ for(;;){
+ σ := 1.0;
+ alt{
+ <- sync =>
+ <- sync;
+ c := <- cmd =>
+ case(c){
+ 's' =>
+ δ = θ = 0.0;
+ 'z' =>
+ δ = zd;
+ 'o' =>
+ δ = -zd;
+ 'r' =>
+ θ = 1.0;
+ '+' =>
+ σ = 1.25;
+ f++;
+ if(f > 16){
+ f--;
+ σ = 1.0;
+ }
+ else
+ vf *= σ;
+ '-' =>
+ σ = 0.8;
+ f--;
+ if(f < 0){
+ f++;
+ σ = 1.0;
+ }
+ else
+ vf *= σ;
+ }
+ * =>
+ sys->sleep(0);
+ }
+
+ RDisp.draw(RDisp.r, black, nil, (0, 0));
+ drawscene();
+ reveal();
+
+ if(δ != 0.0 || θ != 0.0){
+ objd -= δ;
+ me.l.a.z -= δ;
+ if(θ != 0.0){
+ roll += θ;
+ pitch += θ;
+ yaw += θ;
+ rpy = 1;
+ }
+ update(projx);
+ }
+ if(σ != 1.0)
+ retime(σ);
+ }
+}
+
+# usually almost sorted
+sort(ts: array of ref Object, n: int)
+{
+ done: int;
+ t: ref Object;
+ q, p: int;
+
+ q = n;
+ do{
+ done = 1;
+ q--;
+ for(p = 0; p < q; p++){
+ if(ts[p].p.z > ts[p+1].p.z){
+ t = ts[p];
+ ts[p] = ts[p+1];
+ ts[p+1] = t;
+ done = 0;
+ }
+ }
+ }while(!done);
+}
+
+kill(pid: int): int
+{
+ fd := sys->open("#p/"+string pid+"/ctl", Sys->OWRITE);
+ if(fd == nil)
+ return -1;
+ if(sys->write(fd, array of byte "kill", 4) != 4)
+ return -1;
+ return 0;
+}
+
+fatal(e: string)
+{
+ sys->fprint(sys->fildes(2), "%s\n", e);
+ raise "fatal";
+}
+
+cmd(top: ref Toplevel, s: string): string
+{
+ e := tk->cmd(top, s);
+ if (e != nil && e[0] == '!')
+ fatal(sys->sprint("tk error on '%s': %s", s, e));
+ return e;
+}
+
+winconfig := array[] of {
+ "frame .f",
+ "button .f.zoomin -text {zoomin} -command {send cmd zoomin}",
+ "button .f.zoomout -text {zoomout} -command {send cmd zoomout}",
+ "button .f.stop -text {stop} -command {send cmd stop}",
+ "pack .f.zoomin -side left",
+ "pack .f.zoomout -side right",
+ "pack .f.stop -side top",
+
+ "frame .f2",
+ "button .f2.slow -text {slow} -command {send cmd slow}",
+ "button .f2.fast -text {fast} -command {send cmd fast}",
+ "button .f2.objs -text {objects} -command {send cmd objs}",
+ "pack .f2.slow -side left",
+ "pack .f2.fast -side right",
+ "pack .f2.objs -side top",
+
+ "panel .p -width " + string MAXX + " -height " + string MAXY,
+
+ "pack .f -side top -fill x",
+ "pack .f2 -side top -fill x",
+ "pack .p -side bottom -fill both -expand 1",
+
+ "pack propagate . 0",
+ "update"
+};
+
+############################################################
+
+# gl.b
+
+#
+# initially generated by c2l
+#
+
+MODEL, PROJ: con iota;
+
+Matrix: type array of array of real;
+
+Mstate: adt{
+ matl: list of Matrix;
+ modl: list of Matrix;
+ prjl: list of Matrix;
+ mull: Matrix;
+ freel: list of Matrix;
+ vk: int;
+ vr: int;
+ vrr: int;
+ vc: ref Draw->Image;
+ ap: array of Draw->Point;
+ apn: int;
+ mx, cx, my, cy: real;
+ ignore: int;
+};
+
+ms: Mstate;
+
+m_new(): Matrix
+{
+ if(ms.freel != nil){
+ m := hd ms.freel;
+ ms.freel = tl ms.freel;
+ return m;
+ }
+ m := array[4] of array of real;
+ for(i := 0; i < 4; i++)
+ m[i] = array[4] of real;
+ return m;
+}
+
+m_free(m: Matrix)
+{
+ ms.freel = m :: ms.freel;
+}
+
+m_init()
+{
+ ms.modl = m_new() :: nil;
+ ms.prjl = m_new() :: nil;
+ ms.matl = ms.modl;
+ ms.mull = m_new();
+ ms.vk = 0;
+ ms.apn = 0;
+ ms.mx = ms.cx = ms.my = ms.cy = 0.0;
+ ms.ignore = 0;
+}
+
+m_print()
+{
+ m := hd ms.matl;
+
+ for(i := 0; i < 4; i++){
+ for(j := 0; j < 4; j++)
+ sys->print("%f ", m[i][j]);
+ sys->print("\n");
+ }
+ sys->print("\n");
+}
+
+m_mode(m: int)
+{
+ if(m == PROJ)
+ ms.matl = ms.prjl;
+ else
+ ms.matl = ms.modl;
+}
+
+m_pushmatrix()
+{
+ if(ms.matl == ms.modl){
+ ms.modl = m_new() :: ms.modl;
+ ms.matl = ms.modl;
+ }
+ else{
+ ms.prjl = m_new() :: ms.prjl;
+ ms.matl = ms.prjl;
+ }
+ s := hd tl ms.matl;
+ d := hd ms.matl;
+ for(i := 0; i < 4; i++)
+ for(j := 0; j < 4; j++)
+ d[i][j] = s[i][j];
+}
+
+m_popmatrix()
+{
+ if(ms.matl == ms.modl){
+ m_free(hd ms.modl);
+ ms.modl = tl ms.modl;
+ ms.matl = ms.modl;
+ }
+ else{
+ m_free(hd ms.prjl);
+ ms.prjl = tl ms.prjl;
+ ms.matl = ms.prjl;
+ }
+}
+
+m_loadidentity()
+{
+ i, j: int;
+ m := hd ms.matl;
+
+ for(i = 0; i < 4; i++){
+ for(j = 0; j < 4; j++)
+ m[i][j] = real 0;
+ m[i][i] = real 1;
+ }
+}
+
+m_translate(x: real, y: real, z: real)
+{
+ i: int;
+ m := hd ms.matl;
+
+ for(i = 0; i < 4; i++)
+ m[i][3] = x*m[i][0]+y*m[i][1]+z*m[i][2]+m[i][3];
+}
+
+m_scale(x: real, y: real, z: real)
+{
+ i: int;
+ m := hd ms.matl;
+
+ for(i = 0; i < 4; i++){
+ m[i][0] *= x;
+ m[i][1] *= y;
+ m[i][2] *= z;
+ }
+}
+
+# rotate about x, y or z axes
+rot(deg: real, j: int, k: int)
+{
+ i: int;
+ m := hd ms.matl;
+ rad := Math->Pi*deg/real 180;
+ s := maths->sin(rad);
+ c := maths->cos(rad);
+ a, b: real;
+
+ for(i = 0; i < 4; i++){
+ a = m[i][j];
+ b = m[i][k];
+ m[i][j] = c*a+s*b;
+ m[i][k] = c*b-s*a;
+ }
+}
+
+m_rotatex(a: real)
+{
+ rot(a, 1, 2);
+}
+
+m_rotatey(a: real)
+{
+ rot(a, 2, 0);
+}
+
+m_rotatez(a: real)
+{
+ rot(a, 0, 1);
+}
+
+# (l m n) normalized
+m_rotate(deg: real, l: real, m: real, n:real)
+{
+ i: int;
+ mx := hd ms.matl;
+ rad := Math->Pi*deg/real 180;
+ s := maths->sin(rad);
+ c := maths->cos(rad);
+ f := 1.0-c;
+ m0, m1, m2: real;
+
+ for(i = 0; i < 4; i++){
+ m0 = mx[i][0];
+ m1 = mx[i][1];
+ m2 = mx[i][2];
+ mx[i][0] = m0*(l*l*f+c)+m1*(l*m*f+n*s)+m2*(l*n*f-m*s);
+ mx[i][1] = m0*(l*m*f-n*s)+m1*(m*m*f+c)+m2*(m*n*f+l*s);
+ mx[i][2] = m0*(l*n*f+m*s)+m1*(m*n*f-l*s)+m2*(n*n*f+c);
+ }
+}
+
+# Frustum(-l, l, -l, l, n, f)
+m_frustum(l: real, n: real, f: real)
+{
+ i: int;
+ m := hd ms.matl;
+ r := n/l;
+ a, b: real;
+
+ f = ∞;
+ for(i = 0; i < 4; i++){
+ a = m[i][2];
+ b = m[i][3];
+ m[i][0] *= r;
+ m[i][1] *= r;
+ m[i][2] = a+b;
+ m[i][3] = 0.0;
+ # m[i][2] = -(a*(f+n)/(f-n)+b);
+ # m[i][3] = real -2*f*n*a/(f-n);
+ }
+}
+
+# Ortho(-l, l, -l, l, n, f)
+m_ortho(l: real, n: real, f: real)
+{
+ i: int;
+ m := hd ms.matl;
+ r := real 1/l;
+ # a: real;
+
+ n = 0.0;
+ f = ∞;
+ for(i = 0; i < 4; i++){
+ # a = m[i][2];
+ m[i][0] *= r;
+ m[i][1] *= r;
+ # m[i][2] *= real -2/(f-n);
+ # m[i][3] -= a*(f+n)/(f-n);
+ }
+}
+
+m_loadmatrix(u: array of array of real)
+{
+ m := hd ms.matl;
+
+ for(i := 0; i < 4; i++)
+ for(j := 0; j < 4; j++)
+ m[i][j] = u[i][j];
+}
+
+m_storematrix(u: array of array of real)
+{
+ m := hd ms.matl;
+
+ for(i := 0; i < 4; i++)
+ for(j := 0; j < 4; j++)
+ u[i][j] = m[i][j];
+}
+
+m_matmul()
+{
+ m, p, r: Matrix;
+
+ m = hd ms.modl;
+ p = hd ms.prjl;
+ r = ms.mull;
+ for(i := 0; i < 4; i++){
+ pr := p[i];
+ rr := r[i];
+ for(j := 0; j < 4; j++)
+ rr[j] = pr[0]*m[0][j]+pr[1]*m[1][j]+pr[2]*m[2][j]+pr[3]*m[3][j];
+ }
+}
+
+m_vertexo(x: real, y: real, z: real)
+{
+ m: Matrix;
+ mr: array of real;
+ w, x1, y1, z1: real;
+
+ m = hd ms.modl;
+ mr = m[0]; x1 = x*mr[0]+y*mr[1]+z*mr[2]+mr[3];
+ mr = m[1]; y1 = x*mr[0]+y*mr[1]+z*mr[2]+mr[3];
+ mr = m[2]; z1 = x*mr[0]+y*mr[1]+z*mr[2]+mr[3];
+ mr = m[3]; w = x*mr[0]+y*mr[1]+z*mr[2]+mr[3];
+ if(w != real 1){
+ x1 /= w;
+ y1 /= w;
+ z1 /= w;
+ }
+ if(z1 >= 0.0){
+ ms.ignore = 1;
+ return;
+ }
+ m = hd ms.prjl;
+ mr = m[0]; x = x1*mr[0]+y1*mr[1]+z1*mr[2]+mr[3];
+ mr = m[1]; y = x1*mr[0]+y1*mr[1]+z1*mr[2]+mr[3];
+ mr = m[2]; z = x1*mr[0]+y1*mr[1]+z1*mr[2]+mr[3];
+ mr = m[3]; w = x1*mr[0]+y1*mr[1]+z1*mr[2]+mr[3];
+ if(w == real 0){
+ ms.ignore = 1;
+ return;
+ }
+ if(w != real 1){
+ x /= w;
+ y /= w;
+ # z /= w;
+ }
+ ms.ap[ms.apn++] = (int (ms.mx*x+ms.cx), int (ms.my*y+ms.cy));
+}
+
+m_vertex(x: real, y: real, z: real): (real, real)
+{
+ m: Matrix;
+ mr: array of real;
+ w, x1, y1, z1: real;
+
+ m = ms.mull;
+ mr = m[0]; x1 = x*mr[0]+y*mr[1]+z*mr[2]+mr[3];
+ mr = m[1]; y1 = x*mr[0]+y*mr[1]+z*mr[2]+mr[3];
+ mr = m[2]; z1 = x*mr[0]+y*mr[1]+z*mr[2]+mr[3];
+ mr = m[3]; w = x*mr[0]+y*mr[1]+z*mr[2]+mr[3];
+ if(w == real 0){
+ ms.ignore = 1;
+ return (x1, y1);
+ }
+ if(w != real 1){
+ x1 /= w;
+ y1 /= w;
+ # z1 /= w;
+ }
+ if(z1 >= 0.0){
+ ms.ignore = 1;
+ return (x1, y1);
+ }
+ ms.ap[ms.apn++] = (int (ms.mx*x1+ms.cx), int (ms.my*y1+ms.cy));
+ return (x1, y1);
+}
+
+m_circle(x: real, y: real, z: real, r: real)
+{
+ (d, nil) := m_vertex(x, y, z);
+ (e, nil) := m_vertex(x+r, y, z);
+ d -= e;
+ if(d < 0.0)
+ d = -d;
+ ms.vr = int (ms.mx*d);
+}
+
+m_ellipse(x: real, y: real, z: real, v: Vector)
+{
+ m_circle(x, y, z, v.x);
+ (nil, d) := m_vertex(x, y, z);
+ (nil, e) := m_vertex(x, y+v.y, z);
+ d -= e;
+ if(d < 0.0)
+ d = -d;
+ ms.vrr = int (ms.my*d);
+}
+
+m_begin(k: int, n: int)
+{
+ ms.ignore = 0;
+ ms.vk = k;
+ ms.ap = array[n+1] of Draw->Point;
+ ms.apn = 0;
+}
+
+m_end()
+{
+ if(ms.ignore)
+ return;
+ case(ms.vk){
+ CIRCLE =>
+ RDisp.ellipse(ms.ap[0], ms.vr, ms.vr, 0, ms.vc, (0, 0));
+ FILLCIRCLE =>
+ RDisp.fillellipse(ms.ap[0], ms.vr, ms.vr, ms.vc, (0, 0));
+ ELLIPSE =>
+ RDisp.ellipse(ms.ap[0], ms.vr, ms.vrr, 0, ms.vc, (0, 0));
+ FILLELLIPSE =>
+ RDisp.fillellipse(ms.ap[0], ms.vr, ms.vrr, ms.vc, (0, 0));
+ POLY =>
+ ms.ap[len ms.ap-1] = ms.ap[0];
+ RDisp.poly(ms.ap, Draw->Endsquare, Draw->Endsquare, 0, ms.vc, (0, 0));
+ FILLPOLY =>
+ ms.ap[len ms.ap-1] = ms.ap[0];
+ RDisp.fillpoly(ms.ap, ~0, ms.vc, (0, 0));
+ }
+}
+
+m_colour(i: ref Draw->Image)
+{
+ ms.vc = i;
+}
+
+m_viewport(x1: int, y1: int, x2: int, y2: int)
+{
+ ms.mx = real (x2-x1)/2.0;
+ ms.cx = real (x2+x1)/2.0;
+ ms.my = real (y2-y1)/2.0;
+ ms.cy = real (y2+y1)/2.0;
+}
+
+# sys->print("%d %f (%f %f %f) %s\n", ok, λ, 1.0, 2.0*vdot(l.a, l.d), vdot(l.a, l.a)-icubd2, lstring(l));
+
+# sys->print("%d %f (%f %f %f) %s\n", ok, λ, vmuldiv(l.d, l.d, e2), 2.0*vmuldiv(l.a, l.d, e2), vmuldiv(l.a, l.a, e2)-1.0, lstring(l));
+
+# for(lp = lp0 ; lp != nil; lp = tl lp){
+# p := hd lp;
+# (ok, λ, pt) := intersect(l, p);
+# sys->print("%d %x %d %f %s %s %s\n", p.id, tr.pmask, ok, λ, vstring(pt), pstring(p), lstring(l));
+# } \ No newline at end of file
generated by cgit v1.2.3 (git 2.46.0) at 2025-12-10 22:07:06 +0000