ref: 1018d7a7091fd5e0433f501a3474cdc3fec322cc
dir: /appl/wm/collide.b/
# # 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)); # }