ref: ec35f468e0eba87c9f09cbbe5fa8af2591e6f914
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));
# }