PROCEDURE ComputePolynomialCoeffs( READONLY p, q, r: LR3.T; vp, vq, vr: LONGREAL; v0: LONGREAL; VAR c0: SPPolyPWFunction.Coeffs; (* 3 coefficients *) VAR c1: SPPolyPWFunction.Coeffs; (* 1 coefficient *) ) = (* Given points "p", "q", "r" on the sphere, and values "vp", "vq", "vr", compute the coefficients of the affine (first-degree non-homogeneous) polynomial "P((x,y,z)) = A*x + B*y + C*z + D" such that "P(p) = vp", "P(q) = vq", "P(r) = vr". *) BEGIN WITH m = LR3x3.T{ LR3.T{p[0], p[1], p[2]}, LR3.T{q[0], q[1], q[2]}, LR3.T{r[0], r[1], r[2]} }, b = LR3.T{vp - v0, vq - v0, vr - v0}, mInv = LR3x3.Inv(m), coef = LR3x3.MapCol(mInv,b) DO c0[0]:= coef[0]; c0[1]:= coef[1]; c0[2]:= coef[2]; c1[0]:= v0; Wr.PutText(stderr, "coeffs = "); FOR k := 0 TO 2 DO Wr.PutText(stderr, FLR(c0[k], 15, 8)); END; Wr.PutText(stderr, FLR(c1[0], 15, 8)); Wr.PutText(stderr, "\n"); END END ComputePolynomialCoeffs; PROCEDURE ChoosePolynomial( tri: SPTriang.T; READONLY vv: ARRAY OF LONGREAL; v0: LONGREAL; ): T = VAR SP: SPPolyPWFunction.T; BEGIN WITH rtri = NEW(REF SPTriang.T), side = tri.side^ DO rtri^ := tri; SP := NEW(T, tri := rtri, supp := LR4.T{1.0d0, 0.0d0, 0.0d0, 0.0d0}, (* positive everywhere *) data := NEW(REF SPPWFunction.Triangles, NUMBER(side)) ); FOR i:= 0 TO LAST(side) DO WITH a = side[i], p = Dest(Onext(a)), q = Org(a), r = Dest(a), cp = p.c, vp = vv[p.num], cq = q.c, vq = vv[q.num], cr = r.c, vr = vv[r.num], t = NEW(SPPolyPWFunction.TriangleData) DO SP.data[i] := t; t.face := i; t.deg := 1; t.c0 := NEW(REF SPPolyPWFunction.Coeffs, 3); t.c1 := NEW(REF SPPolyPWFunction.Coeffs, 1); ComputePolynomialCoeffs(cp, cq, cr, vp, vq, vr, v0, t.c0^, t.c1^) END END END; RETURN SP END ChoosePolynomial; PROCEDURE ChooseVertexValues(READONLY out: ARRAY OF Arc): REF ARRAY OF LONGREAL = BEGIN WITH t = NEW(REF ARRAY OF LONGREAL, NUMBER(out)), tt = t^ DO FOR i := 0 TO LAST(tt) DO tt[i]:= 3.0d0 (* WITH p = Org(out[i]).c DO tt[i]:= p[2] END *) END; RETURN t END; END ChooseVertexValues; PROCEDURE ChooseConstantValue(): LONGREAL = BEGIN RETURN 3.0d0 END ChooseConstantValue; pp.getKeyword("-name"); o.name := pp.getNext(); IF pp.keywordPresent("-icosahedron") THEN o.nSites := 12; o.icosahedron := TRUE ELSE pp.getKeyword("-nSites"); o.nSites := pp.getNextInt(4, 100000); o.icosahedron := FALSE END; IF pp.keywordPresent("-obs") THEN o.obs := ParsePointOption(pp) ELSE o.obs := HLR3.Point{c := LR4.T{0.0d0, 2.0d0, 5.0d0, 2.0d0}} END; IF pp.keywordPresent("-up") THEN o.upp := ParsePointOption(pp) ELSE o.upp := HLR3.Point{c := LR4.T{0.0d0, 0.0d0, 0.0d0, 1.0d0}} END; PROCEDURE ParsePointOption(pp: ParseParams.T): HLR3.Point RAISES {ParseParams.Error} = VAR p: HLR3.Point; BEGIN FOR i := 0 TO 3 DO p.c[i] := pp.getNextLongReal(-1.0d20, +1.0d20) END; RETURN p END ParsePointOption; <*UNUSED*> PROCEDURE Beep() = BEGIN Wr.PutChar(stderr, '.') END Beep; PROCEDURE OutputDelaunay( e: Arc; name: TEXT; obs, upp: HLR3.Point; ) = BEGIN Wr.PutText(stderr, "Writing " & name & "...\n"); WITH tri = SPTriang.BuildTables(e), arc = tri.arc^, side = tri.side^ DO SPTriang.PerspDraw(arc, name, obs, upp); SPTriang.OutputX3D(side, name); SPTriang.Print(tri, name); SPTriang.Check(arc); END END OutputDelaunay; PROCEDURE MakeSites(n: CARDINAL; icosahedron: BOOLEAN): REF ARRAY OF REF Site = (* Creates an array of sites *) CONST IcoA = 1.6180339887498948482d0; IcoB = 1.0d0; BEGIN WITH rand = NEW(Random.Default).init(fixed := TRUE), rs = NEW(REF ARRAY OF REF Site, n), s = rs^ DO FOR i := 0 TO LAST(s) DO s[i] := NEW(REF Site); s[i].num := i; WITH c = s[i].c DO IF icosahedron THEN CASE i OF | 00 => c[0] := 0.0d0; c[1] := +IcoA; c[2] := +IcoB; | 01 => c[0] := 0.0d0; c[1] := +IcoA; c[2] := -IcoB; | 02 => c[0] := +IcoA; c[1] := -IcoB; c[2] := 0.0d0; | 03 => c[0] := -IcoA; c[1] := -IcoB; c[2] := 0.0d0; | 04 => c[0] := +IcoB; c[1] := 0.0d0; c[2] := +IcoA; | 05 => c[0] := +IcoB; c[1] := 0.0d0; c[2] := -IcoA; | 06 => c[0] := 0.0d0; c[1] := -IcoA; c[2] := +IcoB; | 07 => c[0] := 0.0d0; c[1] := -IcoA; c[2] := -IcoB; | 08 => c[0] := -IcoB; c[1] := 0.0d0; c[2] := +IcoA; | 09 => c[0] := -IcoB; c[1] := 0.0d0; c[2] := -IcoA; | 10 => c[0] := +IcoA; c[1] := +IcoB; c[2] := 0.0d0; | 11 => c[0] := -IcoA; c[1] := +IcoB; c[2] := 0.0d0; ELSE <* ASSERT FALSE *> END ELSE REPEAT c[0] := 2.0d0*rand.longreal() - 1.0d0; c[1] := 2.0d0*rand.longreal() - 1.0d0; c[2] := 2.0d0*rand.longreal() - 1.0d0; UNTIL c[0]*c[0] + c[1]*c[1] + c[2]*c[2] <= 1.0d0; END; c := LR3.Dir(c); PrintCoords(stderr, c); PrintCoords(stderr, c); END; Wr.PutText(stderr, "\n"); END; RETURN rs END END MakeSites; s = MakeSites(o.nSites, o.icosahedron), e = SPTriang.BuildInc(s^), tri = SPTriang.BuildTables(e), PROCEDURE Small(READONLY p, q, r: Coords) : BOOLEAN = VAR dpq, dpr, dqr: LONGREAL; BEGIN dpq:= Math.sqrt((p[0] - q[0])*(p[0] - q[0]) + (p[1] - q[1])*(p[1] - q[1]) + (p[2] - q[2])*(p[2] - q[2])); dpr:= Math.sqrt((p[0] - r[0])*(p[0] - r[0]) + (p[1] - r[1])*(p[1] - r[1]) + (p[2] - r[2])*(p[2] - r[2])); dqr:= Math.sqrt((q[0] - r[0])*(q[0] - r[0]) + (q[1] - r[1])*(q[1] - r[1]) + (q[2] - r[2])*(q[2] - r[2])); IF (dpq <= o.tol) AND (dpr <= o.tol) AND (dqr <= o.tol) THEN RETURN TRUE ELSE RETURN FALSE; END; END Small; PROCEDURE BuildInc( READONLY site: ARRAY OF REF Site; name: TEXT; obs, upp: HLR3.Point; ): Arc = VAR e: Arc; ss: ARRAY [0..3] OF REF Site := SUBARRAY(site, 0, 4); BEGIN <* ASSERT NUMBER(site) >= 4 *> e := BuildTetrahedron(ss); FOR i := 4 TO LAST(site) DO WITH pname = name & "-partial-" & Fmt.Pad(Fmt.Int(i), 3, '0') DO OutputDelaunay(e, pname, obs, upp); END; e := InsertSite(site[i], e) END; RETURN e; END BuildInc;