#pragma mode( separator(.,;) integer(h32) )

// ============================================================
// BINORM3D - beta 1a
// The bivariate normal surface, turned with a finger:
//   correlation and regression to the mean
// ============================================================
// Author: Roger Metcalf, 2026.  Written with AI assistance
//   (Anthropic Claude); concept, direction and all hardware
//   testing by the author.
//
// FREE FOR EDUCATIONAL AND PERSONAL USE.  Keep this header with
//   the file and credit the author if you pass it on.  Not for
//   sale or commercial redistribution.  Provided as is, with no
//   warranty: check any result you intend to rely on.
// DISCLAIMER: a teaching and demonstration aid. Not validated for
//   forensic, clinical or legal use; do not rely on its results
//   for legal purposes.
//
// Firmware: written for the 2026-09-09 firmware (2.4.2).  This
//   version has NOT yet been run on hardware or on the emulator,
//   only in the PC harness (ppl_run_v3e, 19 Oct 2026); beta 1
//   passed the emulator (HP Prime Virtual Calculator, build
//   2026-09-09, 17 Oct 2026).  The arithmetic was
//   checked against scipy (multivariate_normal, norm, chi2) and
//   the projection against a Python copy of it before the PPL was
//   written: verify_BiNorm3D.py.
//
// Run BiNorm3D() from the Program Catalog.
//
// beta 1a (19 Oct 2026): N1, N2 (fleet review 19 Oct 2026,
//   findings_B).  Line numbers are beta 1's.  No number changed.
//   - N1: help page 4 had lost a verb: "offenders chosen for a /
//     programme because their risk score was / extreme score closer
//     to the average next / time ..." (981-984).  It now reads
//     "... was / extreme will score closer to the average / next
//     time, programme or not.  That drop is / not the programme's
//     effect; only a comparison / group picked the same way shows
//     what it added." -- still six lines, none over 46 characters.
//   - N2: the menu line "Grid while dragging (now 15x15)" (872-873)
//     was 31 characters, and the emulator cut menu lines of 31 to 36
//     characters at the box's edge.  It now reads "Grid while
//     dragging: 15x15" (or 11x11), 26 characters; same item, same
//     action.
//
// beta 1 (16 Oct 2026): first version.
//   - WIREFRAME ONLY.  The surface is drawn with LINE_P alone, with
//     this program's own rotation and projection arithmetic.  The
//     3-D TRIANGLE_P mesh, FILLPOLY_P and alpha are not device-
//     verified yet (GfxProbe beta 1 tests them, T01-T07).  A shaded
//     TRIANGLE_P version can follow once GfxProbe T07 passes on the
//     G2 and the G1.
//   - ANGLE MODE.  SIN and COS follow the calculator's angle setting
//     (degrees, radians or grads).  This program never changes that
//     setting: it keeps every angle in RADIANS and converts at each
//     call with bnAU = ASIN(1) / (pi / 2), which is 57.30 in degrees
//     mode, 1 in radians and 63.66 in grads.
//   - COARSE GRID WHILE DRAGGING.  A full 3-D frame is 426 LINE_P
//     calls; with the slice drawn, 485 LINE_P and 49 PIXON_P (counted
//     in the PC harness).  That is past the ~450 segments the house
//     budget allows a live drag on a G1, so while a finger turns the
//     surface it is drawn on an 11 x 11 grid (226 LINE_P) and redrawn
//     at 15 x 15 when the finger lifts.  Menu > Grid while dragging
//     turns this off; Menu > Frame timing measures the 3-D frame at
//     both grids and the contour-with-slice frame on the calculator.
//     No frame time has been measured yet.
//   - GLYPHS.  rho, theta, phi and the degree sign (CHAR 961, 952,
//     966 and 176) are not on the verified list yet (GfxProbe T14
//     tests rho and theta); the corpus uses all four in TEXTOUT_P.
//     They are built in ONE place, bnInit: if any shows as a box,
//     change that line to plain text ("r", "th", "ph", " deg").
//   - MAKELIST keeps Home's K as its counter (the verified form), so
//     a run overwrites Home's K, as ChiTable and GaltonBoard do,
//     until GfxProbe T19 clears a LOCAL counter.
// ============================================================
//
// WHAT YOU SEE
//   3-D view: the joint density f(x, y) of two N(0, 1) variables X
//     and Y with correlation rho, as a wireframe over the square
//     -3..3 by -3..3.  Colour is height: blue low, magenta, red at
//     the top.  The z scale is fixed, so the peak visibly rises as
//     the bell narrows (the volume under it stays 1).
//   Contour view (Cont): the surface from above.  The ellipses hold
//     39 % and 86 % of all (x, y) pairs; the magenta line is the
//     regression line E[Y | X = x] = rho x; the dotted line is y = x.
//   Slice: the cut through the surface at x = x0.  From above it is
//     a red line you drag, with a side panel: the distribution of Y
//     among cases with X = x0, N(rho x0, (1 - rho^2)^0.5), against
//     N(0, 1) in grey.  Its centre rho x0 is always nearer 0 than
//     x0 itself: REGRESSION TO THE MEAN.  In 3-D the cut is drawn in
//     red on the surface, its top marked above the regression line.
//
// THE ARITHMETIC (checked against scipy, verify_BiNorm3D.py)
//   f(x, y) = exp(-Q / 2) / (2 pi (1 - rho^2)^0.5),
//   Q = (x^2 - 2 rho x y + y^2) / (1 - rho^2)
//   Y | X = x0  ~  N(rho x0, (1 - rho^2)^0.5)
//   share of pairs inside the ellipse Q = k^2:  1 - exp(-k^2 / 2)
//     (Q is chi-square with 2 df): 39.3 % for k = 1, 86.5 % for 2
//   ellipse Q = k^2 drawn as x = k (a cos t + b sin t) / 2^0.5,
//     y = k (a cos t - b sin t) / 2^0.5, a = (1 + rho)^0.5,
//     b = (1 - rho)^0.5
//   3-D view, orthographic, azimuth theta and elevation phi:
//     screen x = 160 + s (x cos theta - y sin theta)
//     screen y = cy - s ((x sin theta + y cos theta) sin phi
//                        + k z cos phi)
//     s = 22 px per unit, k = 24: the rho = 0 peak is 84 px tall,
//     the largest that keeps the rho = 0.9 peak (2.29 times taller)
//     on the screen at every tilt.
//
// KEYS AND TOUCH
//   3-D view: drag sideways = turn (0.6 degrees per pixel); drag
//     down = look from higher up, drag up = lower (phi 10..80
//     degrees); Left / Right turn and Up / Down raise or lower the
//     eye in 10-degree steps.
//   Contour view: tap or drag in the plot = slice at that x0;
//     Left / Right move x0 by 0.25.
//   Soft keys: rho- rho+ (also the - and + keys) | Cont (3-D) |
//     Slice | Reset (theta 35, phi 30, rho 0.60) | Menu (also Enter)
//   Menu: set rho exactly, the views, slice, reset, the grid while
//     dragging, frame timing, help, quit.  Esc = quit.
//
// HOME EXPORT (readable after the program ends)
//   BN_RHO  the last rho
//
// SECTIONS
//   state and geometry; small helpers; the cache (once per rho);
//   the projection; the 3-D view; the contour view and slice
//   panel; dragging; setting rho and x0; menu, timing, help;
//   init and main.
//
// LAYOUT NOTE
//   Font 1 is about 6.5 px per character.  Help pages start at
//   x = 12 and keep each line to 46 characters; every TEXTOUT_P
//   passes a width so nothing runs off the right-hand edge.  The
//   soft-key strip takes y 218-239, so the views stop at y 216 and
//   carry no footer line.
// ============================================================


EXPORT BN_RHO;

// ------------------------------------------------------------
// STATE SHARED BY EVERY ROUTINE (private to this program; all of
// it is set in bnInit at the start of every run)
// ------------------------------------------------------------
LOCAL bnRho, bnTh, bnPh, bnView, bnSlc, bnX0;
//  bnRho  the correlation, -0.9 .. 0.9
//  bnTh   azimuth theta in RADIANS: the turn about the vertical axis
//  bnPh   elevation phi in RADIANS (10..80 degrees): the eye's height
//  bnView 1 = the 3-D surface, 2 = the contour view from above
//  bnSlc  1 = the slice is on      bnX0  where it cuts: x = x0
LOCAL bnDrag, bnCrs, bnMsg, bnSlD, bnTX, bnTY;
//  bnDrag 1 while a finger is turning the surface
//  bnCrs  1 = draw the coarse 11 x 11 grid while dragging
//  bnMsg  a one-frame message, "" for none
//  bnSlD  1 = the 3-D slice curve needs recomputing
//  bnTX, bnTY  where the last tap landed (set by bnWaitAny)
LOCAL bnAU, bnPI, bnDR;
//  bnAU angle unit: radians times bnAU is what SIN and COS expect
//  bnPI pi;  bnDR one degree in radians

// ---- the cache: everything that depends on rho alone ----
LOCAL bnGV, bnZ, bnSCX, bnSCY, bnPk, bnRamp;
//  bnGV  grid values: 15 for the fine grid (-3 .. 3 in steps of
//        3/7), then 11 for the coarse grid (steps of 0.6)
//  bnZ   heights f(x, y): 225 for the fine grid, row by row (y
//        outer, x inner), then 121 for the coarse grid
//  bnSCX colour of every segment along x: 210 fine, then 110 coarse
//  bnSCY colour of every segment along y: 210 fine, then 110 coarse
//  bnPk  the peak f(0, 0)    bnRamp 16 colours, blue..magenta..red
LOCAL bnEX, bnEY, bnCnd, bnSigC, bnPSc, bnMg;
//  bnEX, bnEY the two contour ellipses in pixels: 61 points for
//             k = 1, then 61 for k = 2
//  bnCnd  the slice density f(y | x0) at 0, 1, .. 120 pixels from
//         its centre (its SHAPE is the same for every x0)
//  bnSigC its SD, (1 - rho^2)^0.5   bnPSc panel pixels per density unit
//  bnMg   the N(0, 1) density at every pixel row of the contour plot
LOCAL bnSlZ, bnSlTop, bnCos, bnSin;
//  bnSlZ   the 3-D slice: f(x0, y) at 29 values of y
//  bnSlTop f(x0, rho x0), the top of that slice
//  bnCos, bnSin  cos t and sin t at 61 angles round an ellipse

// ---- one frame of the 3-D view (set by bnFrame and bnDraw3D) ----
LOCAL bnPX, bnPY, bnWA, bnWB, bnWC, bnWD;
//  bnPX, bnPY screen x and y of every grid point
//  bnWA .. bnWD  the parts of those that depend on x or on y alone
LOCAL bnFA, bnFB, bnFC, bnFD, bnFE, bnFY;
LOCAL bnFCt, bnFSt;
//  screen x = 160 + FA x - FB y,  screen y = FY - FC x - FD y - FE z
//  (bnFrame explains);  bnFCt, bnFSt = cos and sin of theta

// ---- geometry, shared by the drawing and the touch code ----
LOCAL bnS, bnK, bnCX3, bnR3, bnHM, bnBT;
LOCAL bnBB, bnCL, bnCR, bnCT, bnCB, bnCU;
LOCAL bnCXc, bnCYc, bnPL, bnPW;
//  3-D: bnS pixels per unit of x and y, bnK the height stretch; the
//       floor's centre is at x = bnCX3; bnR3 is half the floor's
//       diagonal in pixels and bnHM the tallest peak (rho = 0.9);
//       the picture fits between y = bnBT and y = bnBB
//  contour: plot box bnCL..bnCR by bnCT..bnCB, bnCU pixels per
//       unit, centre (bnCXc, bnCYc); the slice panel's density axis
//       is at x = bnPL and its tallest density is bnPW pixels wide

// ---- house colours and display glyphs ----
LOCAL bnBLUE, bnMAG, bnRED, bnGRN, bnORNG, bnGRAY;
LOCAL bnFOOT, bnGRID, bnFLR, bnLBLU, bnLMAG, bnLORN;
LOCAL bnLRED, bnLGRN;
LOCAL bnRS, bnTS, bnFS, bnDg, bnSq, bnMin;
LOCAL bnX0S;
//  bnRS rho, bnTS theta, bnFS phi, bnDg degree sign, bnSq a
//  superscript 2, bnMin a minus sign, bnX0S "x0" with a subscript
//  zero -- built with CHAR() in bnInit so the code stays ASCII.


// ------------------------------------------------------------
// SMALL HELPERS
// ------------------------------------------------------------

// A number as text with a FIXED number of decimals.  Mode 2 keeps
// trailing zeros; mode 1 would ignore the digits argument.
bnF(fmV, fmD)
BEGIN
  RETURN STRING(fmV, 2, fmD);
END;

// Wait for the finger to lift -- but never for more than a second.
// MOUSE(1) is -1 while nothing touches the screen.  Without the
// time limit a stuck or phantom touch would hold the program here
// for ever, and the only way out would be a reset.
bnDrain()
BEGIN
  LOCAL drT;
  drT := TICKS;                                  // milliseconds
  REPEAT UNTIL MOUSE(1) == -1 OR TICKS - drT > 1000;
END;

// Drain, wait for a key OR a tap, drain.  Returns the key code, or
// -1 for a tap (its x and y are left in bnTX and bnTY).
bnWaitAny()
BEGIN
  LOCAL wk, wm, wm1;
  REPEAT wk := GETKEY; UNTIL wk == -1;
  bnDrain();
  REPEAT
    wk := GETKEY;
    wm := MOUSE;
    wm1 := wm(1);
  UNTIL wk > -1 OR SIZE(wm1) > 0;
  bnTX := -1; bnTY := -1;
  IF SIZE(wm1) > 0 THEN bnTX := wm1(1); bnTY := wm1(2); wk := -1; END;
  bnDrain();
  REPEAT UNTIL GETKEY == -1;
  RETURN wk;
END;

// Keep theta between 0 and 360 degrees (2 pi) after a turn.
bnNormTh()
BEGIN
  WHILE bnTh < 0 DO bnTh := bnTh + 2 * bnPI; END;
  WHILE bnTh >= 2 * bnPI DO bnTh := bnTh - 2 * bnPI; END;
END;

// An angle held in radians, shown in whole degrees.
bnDeg(adV)
BEGIN
  LOCAL adD;
  adD := ROUND(adV / bnDR, 0);
  IF adD >= 360 THEN adD := adD - 360; END;
  RETURN STRING(adD, 2, 0) + bnDg;
END;

// The colour for a height given as a share of the peak (0 to 1):
// one of 16 steps from blue (low) through magenta to red (the top).
bnRampC(rcT)
BEGIN
  RETURN bnRamp(1 + IP(15 * MAX(0, MIN(1, rcT)) + 0.5));
END;


// ------------------------------------------------------------
// THE CACHE.  Everything that depends on rho alone is computed
// here, once each time rho changes, and never inside a drag:
//   the surface heights on both grids and the segment colours,
//   the two contour ellipses, the shape of the slice density.
//   f(x, y) = exp(-(x^2 - 2 rho x y + y^2) / (2 (1 - rho^2)))
//             / (2 pi (1 - rho^2)^0.5)
// ------------------------------------------------------------
bnCache()
BEGIN
  LOCAL chI, chJ, chN, chOf, chG, chSo;
  LOCAL chX, chY, chQ, chC, chP, chM;
  LOCAL chA, chB, chU, chV;
  chQ := 1 / (2 * (1 - bnRho^2));           // divides the quadratic form
  chC := 1 / (2 * bnPI * (1 - bnRho^2)^0.5);  // the height at (0, 0)
  bnPk := chC;
  FOR chM FROM 1 TO 2 DO
    // chM = 1: the fine 15 x 15 grid; chM = 2: the coarse 11 x 11
    // grid, stored after it in the same lists
    IF chM == 1 THEN
      chN := 15; chOf := 0; chG := 0; chSo := 0;
    ELSE
      chN := 11; chOf := 225; chG := 15; chSo := 210;
    END;
    // the heights, row by row: point (i, j) is x = value i, y = value j
    FOR chJ FROM 1 TO chN DO
      chY := bnGV(chG + chJ);
      FOR chI FROM 1 TO chN DO
        chX := bnGV(chG + chI);
        bnZ(chOf + (chJ - 1) * chN + chI) := chC * EXP(-(chX * chX - 2 * bnRho * chX * chY + chY * chY) * chQ);
      END;
    END;
    // segment colours, from the mean height of the two ends as a
    // share of the peak.  Along x: point p to point p + 1.
    FOR chJ FROM 1 TO chN DO
      FOR chI FROM 1 TO chN - 1 DO
        chP := chOf + (chJ - 1) * chN + chI;
        bnSCX(chSo + (chJ - 1) * (chN - 1) + chI) := bnRampC((bnZ(chP) + bnZ(chP + 1)) / (2 * chC));
      END;
    END;
    // along y: point p to point p + n (the next row)
    FOR chI FROM 1 TO chN DO
      FOR chJ FROM 1 TO chN - 1 DO
        chP := chOf + (chJ - 1) * chN + chI;
        bnSCY(chSo + (chI - 1) * (chN - 1) + chJ) := bnRampC((bnZ(chP) + bnZ(chP + chN)) / (2 * chC));
      END;
    END;
  END;
  // the contour ellipses Q = k^2 for k = 1 and 2, parametric form:
  //   x = k (a cos t + b sin t) / 2^0.5,  y = k (a cos t - b sin t) / 2^0.5
  // with a = (1 + rho)^0.5 and b = (1 - rho)^0.5; 61 points each
  chA := (1 + bnRho)^0.5 / 2^0.5;
  chB := (1 - bnRho)^0.5 / 2^0.5;
  FOR chM FROM 1 TO 2 DO
    FOR chI FROM 1 TO 61 DO
      chU := chM * (chA * bnCos(chI) + chB * bnSin(chI));
      chV := chM * (chA * bnCos(chI) - chB * bnSin(chI));
      bnEX((chM - 1) * 61 + chI) := IP(bnCXc + bnCU * chU + 0.5);
      bnEY((chM - 1) * 61 + chI) := IP(bnCYc - bnCU * chV + 0.5);
    END;
  END;
  // the slice panel: Y given X = x0 is N(rho x0, (1 - rho^2)^0.5).
  // Only its centre moves with x0, so its shape is cached by the
  // distance in pixels from the centre: 0 .. 120 px (4.4 SD at most).
  bnSigC := (1 - bnRho^2)^0.5;
  FOR chI FROM 0 TO 120 DO
    chY := chI / bnCU;
    bnCnd(chI + 1) := EXP(-chY * chY / (2 * bnSigC * bnSigC)) / (bnSigC * (2 * bnPI)^0.5);
  END;
  // scale the panel so the slice density's peak is bnPW px wide
  bnPSc := bnPW * bnSigC * (2 * bnPI)^0.5;
  bnSlD := 1;
END;

// The 3-D slice: the surface along x = x0 at 29 values of y.  Its
// top is at y = rho x0 -- above the regression line, not the ridge.
// Recomputed only when x0 or rho has changed (bnSlD = 1).
bnSliceC()
BEGIN
  LOCAL scM, scY, scQ;
  scQ := 1 / (2 * (1 - bnRho^2));
  FOR scM FROM 1 TO 29 DO
    scY := -3 + (scM - 1) * 6 / 28;
    bnSlZ(scM) := bnPk * EXP(-(bnX0 * bnX0 - 2 * bnRho * bnX0 * scY + scY * scY) * scQ);
  END;
  scY := bnRho * bnX0;
  bnSlTop := bnPk * EXP(-(bnX0 * bnX0 - 2 * bnRho * bnX0 * scY + scY * scY) * scQ);
  bnSlD := 0;
END;


// ------------------------------------------------------------
// THE PROJECTION (orthographic: no perspective).
// Turn the surface by theta about the vertical axis, then look at
// it from elevation phi:
//   screen x = 160 + s (x cos theta - y sin theta)
//   screen y = cy  - s ((x sin theta + y cos theta) sin phi
//                       + k z cos phi)
// x cos theta - y sin theta is the left-right position after the
// turn; x sin theta + y cos theta is the depth (far points sit
// higher on the screen when we look down at angle phi); the height
// z shows by cos phi.  bnFrame works out the constants once per
// frame; its four SIN / COS calls convert radians with bnAU.
// ------------------------------------------------------------
bnFrame()
BEGIN
  LOCAL frCt, frSt, frSp, frCp, frT, frB;
  frCt := COS(bnTh * bnAU); frSt := SIN(bnTh * bnAU);
  frSp := SIN(bnPh * bnAU); frCp := COS(bnPh * bnAU);
  bnFA := bnS * frCt;
  bnFB := bnS * frSt;
  bnFC := bnS * frSt * frSp;
  bnFD := bnS * frCt * frSp;
  bnFE := bnS * bnK * frCp;
  bnFCt := frCt; bnFSt := frSt;
  // Where the floor's centre goes.  It depends on phi only, so the
  // picture does not jump while theta or rho change: room above for
  // the tallest peak (rho = 0.9) or the far corner, room below for
  // the near corner; centred between bnBT and bnBB when it all fits.
  frT := MAX(bnHM * frCp, bnR3 * frSp);
  frB := bnR3 * frSp;
  IF frT + frB <= bnBB - bnBT THEN
    bnFY := bnBT + frT + (bnBB - bnBT - frT - frB) / 2;
  ELSE
    bnFY := bnBB - frB;
  END;
END;

// One point (x, y, z) of the 3-D view on the screen: the floor,
// the axes and the slice use these; the grid uses bnDraw3D's lists.
bnPrjX(pjX, pjY)
BEGIN
  RETURN IP(bnCX3 + 0.5 + bnFA * pjX - bnFB * pjY);
END;

bnPrjY(pjX, pjY, pjZ)
BEGIN
  RETURN IP(bnFY + 0.5 - bnFC * pjX - bnFD * pjY - bnFE * pjZ);
END;

// The soft keys, painted on the screen after every blit.  The
// labels never produce a key code; a tap with y >= 218 is read as
// slot MIN(6, IP(x / 53.4) + 1) by the main loop.
bnMenuBar()
BEGIN
  IF bnView == 1 THEN
    DRAWMENU({bnRS + bnMin, bnRS + "+", "Cont", "Slice", "Reset", "Menu"});
  ELSE
    DRAWMENU({bnRS + bnMin, bnRS + "+", "3-D", "Slice", "Reset", "Menu"});
  END;
END;


// ------------------------------------------------------------
// THE 3-D VIEW: floor, wireframe, slice, readouts -- one frame in
// G1, one blit.  While a finger drags (bnDrag = 1) and the coarse
// option is on, the 11 x 11 grid is drawn instead of 15 x 15.
// ------------------------------------------------------------
bnDraw3D()
BEGIN
  LOCAL d3N, d3Of, d3G, d3So, d3I, d3J;
  LOCAL d3Q, d3P, d3X, d3Y, d3Lx, d3Ly;
  LOCAL d3Cx, d3Cy, d3M, d3V;
  bnFrame();
  IF bnSlc == 1 AND bnSlD == 1 THEN bnSliceC(); END;
  DIMGROB_P(G1, 320, 240, #FFFFFFh);
  // ---- which grid ----
  d3N := 15; d3Of := 0; d3G := 0; d3So := 0;
  IF bnDrag == 1 AND bnCrs == 1 THEN
    d3N := 11; d3Of := 225; d3G := 15; d3So := 210;
  END;
  // ---- project every grid point (x, y, f(x, y)) ----
  // The x-part and the y-part of each coordinate are worked out
  // once per grid value; a point then costs two lookups and a sum.
  FOR d3I FROM 1 TO d3N DO
    d3V := bnGV(d3G + d3I);
    bnWA(d3I) := bnFA * d3V;
    bnWB(d3I) := bnFB * d3V;
    bnWC(d3I) := bnFC * d3V;
    bnWD(d3I) := bnFD * d3V;
  END;
  d3Cx := bnCX3 + 0.5; d3Cy := bnFY + 0.5;
  FOR d3J FROM 1 TO d3N DO
    FOR d3I FROM 1 TO d3N DO
      d3P := (d3J - 1) * d3N + d3I;
      bnPX(d3P) := IP(d3Cx + bnWA(d3I) - bnWB(d3J));
      bnPY(d3P) := IP(d3Cy - bnWC(d3I) - bnWD(d3J) - bnFE * bnZ(d3Of + d3P));
    END;
  END;
  // ---- the floor (z = 0): the square -3..3 by -3..3 and the axes ----
  LINE_P(G1, bnPrjX(-3, -3), bnPrjY(-3, -3, 0), bnPrjX(3, -3), bnPrjY(3, -3, 0), bnFLR);
  LINE_P(G1, bnPrjX(3, -3), bnPrjY(3, -3, 0), bnPrjX(3, 3), bnPrjY(3, 3, 0), bnFLR);
  LINE_P(G1, bnPrjX(3, 3), bnPrjY(3, 3, 0), bnPrjX(-3, 3), bnPrjY(-3, 3, 0), bnFLR);
  LINE_P(G1, bnPrjX(-3, 3), bnPrjY(-3, 3, 0), bnPrjX(-3, -3), bnPrjY(-3, -3, 0), bnFLR);
  LINE_P(G1, bnPrjX(-3, 0), bnPrjY(-3, 0, 0), bnPrjX(3.5, 0), bnPrjY(3.5, 0, 0), bnGRAY);
  LINE_P(G1, bnPrjX(0, -3), bnPrjY(0, -3, 0), bnPrjX(0, 3.5), bnPrjY(0, 3.5, 0), bnGRAY);
  TEXTOUT_P("x", G1, bnPrjX(3.85, 0) - 3, bnPrjY(3.85, 0, 0) - 6, 1, bnGRAY, 10);
  TEXTOUT_P("y", G1, bnPrjX(0, 3.85) - 3, bnPrjY(0, 3.85, 0) - 6, 1, bnGRAY, 10);
  // ---- with the slice on, two lines on the floor: the ridge
  //      y = x (grey dots) and the regression line y = rho x (magenta)
  IF bnSlc == 1 THEN
    FOR d3M FROM 0 TO 48 DO
      d3V := -3 + d3M * 0.125;
      PIXON_P(G1, bnPrjX(d3V, d3V), bnPrjY(d3V, d3V, 0), bnGRAY);
    END;
    LINE_P(G1, bnPrjX(-3, -3 * bnRho), bnPrjY(-3, -3 * bnRho, 0), bnPrjX(3, 3 * bnRho), bnPrjY(3, 3 * bnRho, 0), bnMAG);
  END;
  // ---- the wireframe, far lines first so near ones are drawn over
  // them.  Lines along x (one per y value): the far end of the
  // floor is large y when cos theta > 0.
  FOR d3Q FROM 1 TO d3N DO
    IF bnFCt > 0 THEN d3J := d3N + 1 - d3Q; ELSE d3J := d3Q; END;
    FOR d3I FROM 1 TO d3N - 1 DO
      d3P := (d3J - 1) * d3N + d3I;
      LINE_P(G1, bnPX(d3P), bnPY(d3P), bnPX(d3P + 1), bnPY(d3P + 1), bnSCX(d3So + (d3J - 1) * (d3N - 1) + d3I));
    END;
  END;
  // lines along y (one per x value): far is large x when sin theta > 0
  FOR d3Q FROM 1 TO d3N DO
    IF bnFSt > 0 THEN d3I := d3N + 1 - d3Q; ELSE d3I := d3Q; END;
    FOR d3J FROM 1 TO d3N - 1 DO
      d3P := (d3J - 1) * d3N + d3I;
      LINE_P(G1, bnPX(d3P), bnPY(d3P), bnPX(d3P + d3N), bnPY(d3P + d3N), bnSCY(d3So + (d3I - 1) * (d3N - 1) + d3J));
    END;
  END;
  // ---- the slice: the surface cut along x = x0, in red ----
  IF bnSlc == 1 THEN
    // where the cut meets the floor
    LINE_P(G1, bnPrjX(bnX0, -3), bnPrjY(bnX0, -3, 0), bnPrjX(bnX0, 3), bnPrjY(bnX0, 3, 0), bnRED);
    // the cut itself, drawn twice (one pixel apart) to stand out
    d3Lx := bnPrjX(bnX0, -3); d3Ly := bnPrjY(bnX0, -3, bnSlZ(1));
    FOR d3M FROM 2 TO 29 DO
      d3V := -3 + (d3M - 1) * 6 / 28;
      d3X := bnPrjX(bnX0, d3V); d3Y := bnPrjY(bnX0, d3V, bnSlZ(d3M));
      LINE_P(G1, d3Lx, d3Ly, d3X, d3Y, bnRED);
      LINE_P(G1, d3Lx, d3Ly + 1, d3X, d3Y + 1, bnRED);
      d3Lx := d3X; d3Ly := d3Y;
    END;
    // its top stands on the regression line: a magenta post from
    // the floor at (x0, rho x0) up to the surface
    d3X := bnPrjX(bnX0, bnRho * bnX0);
    d3Y := bnPrjY(bnX0, bnRho * bnX0, bnSlTop);
    LINE_P(G1, d3X, bnPrjY(bnX0, bnRho * bnX0, 0), d3X, d3Y, bnMAG);
    RECT_P(G1, d3X - 2, d3Y - 2, d3X + 2, d3Y + 2, bnMAG, bnMAG);
  END;
  // ---- readouts in the corners, clear of the floor ----
  TEXTOUT_P(bnRS + " = " + bnF(bnRho, 2), G1, 4, 3, 2, bnBLUE, 84);
  TEXTOUT_P("peak " + bnF(bnPk, 3), G1, 4, 21, 1, bnGRAY, 80);
  TEXTOUT_P(bnTS + " = " + bnDeg(bnTh), G1, 262, 4, 1, bnGRAY, 56);
  TEXTOUT_P(bnFS + " = " + bnDeg(bnPh), G1, 262, 17, 1, bnGRAY, 56);
  IF DIM(bnMsg) > 0 THEN
    TEXTOUT_P(bnMsg, G1, 4, 34, 1, bnRED, 120);
  ELSE
    IF bnDrag == 1 AND bnCrs == 1 THEN
      TEXTOUT_P("11" + CHAR(215) + "11 while dragging", G1, 4, 34, 1, bnGRAY, 120);
    END;
  END;
  IF bnSlc == 1 THEN
    TEXTOUT_P("cut at " + bnX0S + " = " + bnF(bnX0, 2), G1, 4, 190, 1, bnRED, 100);
    TEXTOUT_P("top at y = " + bnF(bnRho * bnX0, 2), G1, 4, 203, 1, bnMAG, 100);
  END;
  BLIT_P(G0, 0, 0, G1);
  bnMenuBar();
END;


// ------------------------------------------------------------
// THE CONTOUR VIEW: the surface from above, the two ellipses, the
// regression line, y = x, and with the slice on the red line at x0
// and the side panel f(y | x0).  The panel's rows are the plot's
// rows, so y reads straight across.  All curves come from the
// cache; nothing here evaluates a density.
// ------------------------------------------------------------
bnDrawC()
BEGIN
  LOCAL dcI, dcX, dcY, dcS, dcR, dcO;
  LOCAL dcW, dcLm, dcLc, dcRm, dcRx, dcM;
  DIMGROB_P(G1, 320, 240, #FFFFFFh);
  TEXTOUT_P("From above, " + bnRS + " = " + bnF(bnRho, 2), G1, 22, 2, 1, bnBLUE, 170);
  // ---- grid lines and tick labels at -3 .. 3 ----
  FOR dcI FROM 0 TO 6 DO
    dcX := bnCL + dcI * bnCU;
    dcY := bnCT + dcI * bnCU;
    LINE_P(G1, dcX, bnCT, dcX, bnCB, bnGRID);
    LINE_P(G1, bnCL, dcY, bnCR, dcY, bnGRID);
    dcS := STRING(dcI - 3, 2, 0);
    TEXTOUT_P(dcS, G1, dcX - IP(3 * DIM(dcS)), bnCB + 2, 1, bnGRAY, 16);
    dcS := STRING(3 - dcI, 2, 0);
    TEXTOUT_P(dcS, G1, 2, dcY - 5, 1, bnGRAY, 18);
  END;
  // ---- the axes and the frame ----
  LINE_P(G1, bnCL, bnCYc, bnCR, bnCYc, bnFLR);
  LINE_P(G1, bnCXc, bnCT, bnCXc, bnCB, bnFLR);
  LINE_P(G1, bnCL, bnCT, bnCR, bnCT, bnGRAY);
  LINE_P(G1, bnCR, bnCT, bnCR, bnCB, bnGRAY);
  LINE_P(G1, bnCR, bnCB, bnCL, bnCB, bnGRAY);
  LINE_P(G1, bnCL, bnCB, bnCL, bnCT, bnGRAY);
  TEXTOUT_P("x", G1, bnCR - 8, bnCYc + 2, 1, bnGRAY, 8);
  TEXTOUT_P("y", G1, bnCXc + 4, bnCT + 1, 1, bnGRAY, 8);
  // ---- y = x, dotted (LINE_P has no dash style) ----
  FOR dcI FROM 0 TO 6 * bnCU STEP 3 DO
    PIXON_P(G1, bnCL + dcI, bnCB - dcI, bnGRAY);
  END;
  // ---- the regression line E[Y | X = x] = rho x, magenta ----
  LINE_P(G1, bnCL, IP(bnCYc + 3 * bnCU * bnRho + 0.5), bnCR, IP(bnCYc - 3 * bnCU * bnRho + 0.5), bnMAG);
  // ---- the ellipses Q = 1 and Q = 4, from the cache ----
  FOR dcM FROM 1 TO 2 DO
    FOR dcI FROM 1 TO 60 DO
      dcO := (dcM - 1) * 61 + dcI;
      LINE_P(G1, bnEX(dcO), bnEY(dcO), bnEX(dcO + 1), bnEY(dcO + 1), bnBLUE);
    END;
    // the share inside, above the ellipse's top point (rho k, k)
    dcX := IP(bnCXc + bnCU * bnRho * dcM + 0.5) - 13;
    dcY := IP(bnCYc - bnCU * dcM + 0.5) - 12;
    // a steep bell (|rho| > 0.75) leaves no room above the inner
    // ellipse, so its label goes on the long axis, 1.65 units out,
    // between the two ellipses (checked for every rho in 0.01 steps)
    IF dcM == 1 AND ABS(bnRho) > 0.75 THEN
      dcX := IP(bnCXc + bnCU * 1.16673 + 0.5) - 13;
      IF bnRho > 0 THEN
        dcY := IP(bnCYc - bnCU * 1.16673 + 0.5) - 6;
      ELSE
        dcY := IP(bnCYc + bnCU * 1.16673 + 0.5) - 6;
      END;
    END;
    IF dcM == 1 THEN dcS := "39 %"; ELSE dcS := "86 %"; END;
    TEXTOUT_P(dcS, G1, dcX, dcY, 1, bnBLUE, 30);
  END;
  IF bnSlc == 1 THEN
    dcX := IP(bnCXc + bnCU * bnX0 + 0.5);                // the slice, x = x0
    dcRm := IP(bnCYc - bnCU * bnRho * bnX0 + 0.5);       // row of y = rho x0
    dcRx := IP(bnCYc - bnCU * bnX0 + 0.5);               // row of y = x0
    // ---- the panel: f(y | x0) shaded, row by row (every row, so
    //      the shading has no gaps), from the cached shape ----
    FOR dcR FROM bnCT TO bnCB DO
      dcO := ABS(dcR - dcRm);
      IF dcO <= 120 THEN
        dcW := IP(bnCnd(dcO + 1) * bnPSc + 0.5);
        IF dcW >= 1 THEN LINE_P(G1, bnPL + 1, dcR, bnPL + dcW, dcR, bnLMAG); END;
      END;
    END;
    LINE_P(G1, bnPL, bnCT, bnPL, bnCB, bnGRAY);
    // the outlines, every second row: N(0, 1) grey, f(y | x0) magenta
    dcLm := bnPL; dcLc := bnPL;
    FOR dcR FROM bnCT TO bnCB STEP 2 DO
      dcW := IP(bnMg(dcR - bnCT + 1) * bnPSc + 0.5);
      IF dcR > bnCT THEN LINE_P(G1, dcLm, dcR - 2, bnPL + dcW, dcR, bnGRAY); END;
      dcLm := bnPL + dcW;
      dcO := ABS(dcR - dcRm);
      dcW := 0;
      IF dcO <= 120 THEN dcW := IP(bnCnd(dcO + 1) * bnPSc + 0.5); END;
      IF dcR > bnCT THEN LINE_P(G1, dcLc, dcR - 2, bnPL + dcW, dcR, bnMAG); END;
      dcLc := bnPL + dcW;
    END;
    // ---- across to the panel: y = x0 (red dots), where Y would sit
    //      if it kept pace with x0; y = rho x0 (magenta), the centre
    //      of the slice, on the regression line ----
    FOR dcI FROM dcX TO 316 STEP 3 DO PIXON_P(G1, dcI, dcRx, bnRED); END;
    LINE_P(G1, dcX, dcRm, 316, dcRm, bnMAG);
    // ---- the slice line, two pixels wide, and its two points ----
    LINE_P(G1, dcX, bnCT, dcX, bnCB, bnRED);
    LINE_P(G1, dcX + 1, bnCT, dcX + 1, bnCB, bnRED);
    RECT_P(G1, dcX - 2, dcRx - 2, dcX + 2, dcRx + 2, bnRED, bnRED);
    RECT_P(G1, dcX - 2, dcRm - 2, dcX + 2, dcRm + 2, bnMAG, bnMAG);
    // ---- labels: the panel title, and each guide labelled on the
    //      side away from the other one ----
    TEXTOUT_P("f(y | " + bnX0S + ")", G1, 198, 2, 1, bnMAG, 66);
    TEXTOUT_P("N(0, 1)", G1, 268, 2, 1, bnGRAY, 48);
    IF dcRx <= dcRm THEN
      TEXTOUT_P(bnX0S, G1, 300, MAX(bnCT, dcRx - 11), 1, bnRED, 16);
      TEXTOUT_P(bnRS + bnX0S, G1, 296, MIN(bnCB - 10, dcRm + 2), 1, bnMAG, 22);
    ELSE
      TEXTOUT_P(bnX0S, G1, 300, MIN(bnCB - 10, dcRx + 2), 1, bnRED, 16);
      TEXTOUT_P(bnRS + bnX0S, G1, 296, MAX(bnCT, dcRm - 11), 1, bnMAG, 22);
    END;
    // ---- the numbers and the verdict ----
    TEXTOUT_P(bnX0S + " = " + bnF(bnX0, 2) + ":  E[Y | " + bnX0S + "] = " + bnRS + bnX0S + " = " + bnF(bnRho * bnX0, 2) + ",  SD " + bnF(bnSigC, 2), G1, 4, 194, 1, #000000h, 312);
    CASE
      IF DIM(bnMsg) > 0 THEN TEXTOUT_P(bnMsg, G1, 4, 206, 1, bnRED, 312); END;
      IF ABS(bnX0) < 0.001 THEN
        TEXTOUT_P(bnX0S + " = 0 is the mean itself: nothing to regress", G1, 4, 206, 1, bnGRAY, 312);
      END;
      IF ABS(bnRho) < 0.001 THEN
        TEXTOUT_P(bnRS + " = 0: " + bnX0S + " tells nothing, E[Y | " + bnX0S + "] = 0", G1, 4, 206, 1, bnMAG, 312);
      END;
      DEFAULT
        TEXTOUT_P("REGRESSION TO THE MEAN: " + bnF(bnRho * bnX0, 2) + " is nearer 0", G1, 4, 206, 1, bnRED, 312);
    END;
  ELSE
    // ---- no slice: a legend in the panel's place ----
    TEXTOUT_P("Share of all (x, y)", G1, 198, 22, 1, bnBLUE, 118);
    TEXTOUT_P("pairs inside each:", G1, 198, 34, 1, bnBLUE, 118);
    TEXTOUT_P("inner 39 %,", G1, 198, 46, 1, bnBLUE, 118);
    TEXTOUT_P("outer 86 %", G1, 198, 58, 1, bnBLUE, 118);
    TEXTOUT_P("Magenta: the", G1, 198, 76, 1, bnMAG, 118);
    TEXTOUT_P("regression line", G1, 198, 88, 1, bnMAG, 118);
    TEXTOUT_P("E[Y | x] = " + bnRS + "x", G1, 198, 100, 1, bnMAG, 118);
    TEXTOUT_P("Dotted: y = x", G1, 198, 118, 1, bnGRAY, 118);
    TEXTOUT_P("Tap the plot to", G1, 198, 138, 1, bnRED, 118);
    TEXTOUT_P("slice it at x = " + bnX0S, G1, 198, 150, 1, bnRED, 118);
    CASE
      IF bnRho > 0.001 THEN
        TEXTOUT_P("Tilted along y = x: high x goes with high y", G1, 4, 194, 1, #000000h, 312);
      END;
      IF bnRho < -0.001 THEN
        TEXTOUT_P("Tilted along y = " + bnMin + "x: high x goes with low y", G1, 4, 194, 1, #000000h, 312);
      END;
      DEFAULT
        TEXTOUT_P("Circles: x tells nothing about y", G1, 4, 194, 1, #000000h, 312);
    END;
    IF DIM(bnMsg) > 0 THEN
      TEXTOUT_P(bnMsg, G1, 4, 206, 1, bnRED, 312);
    ELSE
      TEXTOUT_P("The magenta line is flatter than the ridge", G1, 4, 206, 1, bnGRAY, 312);
    END;
  END;
  BLIT_P(G0, 0, 0, G1);
  bnMenuBar();
END;

// Draw whichever view is current.
bnShow()
BEGIN
  IF bnView == 1 THEN bnDraw3D(); ELSE bnDrawC(); END;
END;


// ------------------------------------------------------------
// SETTING RHO AND x0
// ------------------------------------------------------------

// A new rho: rounded to two decimals, kept within -0.9 .. 0.9 (with
// a word on screen when it stops), then the cache is rebuilt.
bnSetRho(srV)
BEGIN
  LOCAL srN;
  srN := ROUND(srV, 2);
  IF srN > 0.9 THEN
    srN := 0.9;
    bnMsg := bnRS + " stops at 0.9 here";
  END;
  IF srN < -0.9 THEN
    srN := -0.9;
    bnMsg := bnRS + " stops at " + bnMin + "0.9 here";
  END;
  IF srN <> bnRho THEN
    bnRho := srN;
    bnCache();
  END;
  BN_RHO := bnRho;
END;

// A new x0 for the slice: kept within -2.9 .. 2.9 and snapped to
// steps of 0.05, so the readouts are round numbers.
bnSetX0(sxV)
BEGIN
  LOCAL sxN;
  sxN := ROUND(MAX(-2.9, MIN(2.9, sxV)) * 20, 0) / 20;
  IF sxN <> bnX0 THEN
    bnX0 := sxN;
    bnSlD := 1;
  END;
END;

// Reset: theta 35 degrees, phi 30 degrees, rho 0.60 (x0 back to 1.5).
bnReset()
BEGIN
  bnTh := 35 * bnDR;
  bnPh := 30 * bnDR;
  bnSetRho(0.6);
  bnSetX0(1.5);
END;


// ------------------------------------------------------------
// DRAGGING
// ------------------------------------------------------------

// Turn the surface with a finger.  dgPx, dgPy is the first MOUSE
// sample (the press).  While the finger stays down, every new sample
// sets theta and phi from the distance moved since the press, 0.6
// degrees per pixel; movement under 8 px is a tap, not a drag.  The
// view is redrawn per sample from the cache (coarse if bnCrs = 1)
// and once more at full detail by the main loop after the lift.
// Any key ends the drag too, so a stuck touch cannot hold it.
bnDrag3D(dgPx, dgPy)
BEGIN
  LOCAL dgM, dgM1, dgK, dgT0, dgP0, dgMv;
  LOCAL dgDx, dgDy, dgNt, dgNp;
  dgT0 := bnTh; dgP0 := bnPh; dgMv := 0;
  REPEAT
    dgM := MOUSE;
    dgM1 := dgM(1);
    dgK := GETKEY;
    IF SIZE(dgM1) > 0 THEN
      dgDx := dgM1(1) - dgPx;
      dgDy := dgM1(2) - dgPy;
      IF ABS(dgDx) > 8 OR ABS(dgDy) > 8 THEN dgMv := 1; END;
      IF dgMv == 1 THEN
        // sideways = turn; down = look from higher up (phi grows)
        dgNt := dgT0 + 0.6 * dgDx * bnDR;
        dgNp := MAX(10 * bnDR, MIN(80 * bnDR, dgP0 + 0.6 * dgDy * bnDR));
        IF dgNt <> bnTh OR dgNp <> bnPh THEN
          bnTh := dgNt; bnPh := dgNp;
          bnDrag := 1;
          bnDraw3D();
        END;
      END;
    END;
  UNTIL SIZE(dgM1) == 0 OR dgK > -1;
  bnDrag := 0;
  bnNormTh();
  bnDrain();                         // capped: at most one second
  RETURN dgMv;
END;

// Move the slice with a finger (contour view).  A press inside the
// plot box puts the slice there and switches it on; it then follows
// the finger, redrawn only when x0 really changes (0.05 steps).
bnDragC(dsPx, dsPy)
BEGIN
  LOCAL dsM, dsM1, dsK, dsOld;
  IF dsPx < bnCL OR dsPx > bnCR OR dsPy < bnCT OR dsPy > bnCB THEN
    bnDrain();
    RETURN 0;
  END;
  bnSlc := 1;
  bnSetX0((dsPx - bnCXc) / bnCU);
  bnDrawC();
  dsOld := bnX0;
  REPEAT
    dsM := MOUSE;
    dsM1 := dsM(1);
    dsK := GETKEY;
    IF SIZE(dsM1) > 0 THEN
      bnSetX0((dsM1(1) - bnCXc) / bnCU);
      IF bnX0 <> dsOld THEN
        bnDrawC();
        dsOld := bnX0;
      END;
    END;
  UNTIL SIZE(dsM1) == 0 OR dsK > -1;
  bnDrain();
  RETURN 1;
END;


// ------------------------------------------------------------
// MENU, RHO BY HAND, FRAME TIMING
// ------------------------------------------------------------

// Type rho exactly.  INPUT needs a variable NAME, so the value goes
// through the local arV.
bnAskRho()
BEGIN
  LOCAL arV, arOK, arL;
  arV := bnRho;
  arL := bnRS + " =";
  arOK := INPUT({arV}, "Correlation", {arL}, {"Between -0.9 and 0.9; two decimals are kept"});
  IF arOK == 1 THEN
    IF arV < -0.9 OR arV > 0.9 THEN
      MSGBOX(bnRS + " must be between -0.9 and 0.9 here.  At 1 or -1 the bell collapses onto a line: no spread across it at all.");
    ELSE
      bnSetRho(arV);
    END;
  END;
END;

// Frame timing, for the G1 test: five 3-D frames at 15 x 15, five
// at 11 x 11 (as while dragging), then five contour frames with the
// slice on (as while dragging x0), in milliseconds (TICKS).
bnTiming()
BEGIN
  LOCAL tmI, tmT, tmF, tmC, tmV, tmK;
  LOCAL tmS, tmP, tmL;
  tmV := bnView; tmK := bnCrs; tmL := bnSlc;
  bnView := 1;
  bnDrag := 0;
  tmT := TICKS;
  FOR tmI FROM 1 TO 5 DO
    bnTh := bnTh + 2 * bnDR;
    bnDraw3D();
  END;
  tmF := (TICKS - tmT) / 5;
  bnDrag := 1; bnCrs := 1;
  tmT := TICKS;
  FOR tmI FROM 1 TO 5 DO
    bnTh := bnTh + 2 * bnDR;
    bnDraw3D();
  END;
  tmC := (TICKS - tmT) / 5;
  bnDrag := 0;
  // the contour view with the slice: the frame an x0 drag redraws
  bnView := 2; bnSlc := 1;
  tmT := TICKS;
  FOR tmI FROM 1 TO 5 DO bnDrawC(); END;
  tmP := (TICKS - tmT) / 5;
  bnCrs := tmK; bnView := tmV; bnSlc := tmL;
  bnNormTh();
  tmS := "off";
  IF bnSlc == 1 THEN tmS := "on"; END;
  tmF := ROUND(tmF, 0); tmC := ROUND(tmC, 0); tmP := ROUND(tmP, 0);
  MSGBOX("3-D frame (slice " + tmS + "): " + STRING(tmF, 2, 0) + " ms at 15" + CHAR(215) + "15, " + STRING(tmC, 2, 0) + " ms at 11" + CHAR(215) + "11." + CHAR(10) + "From above with the slice: " + STRING(tmP, 2, 0) + " ms." + CHAR(10) + "Over about 250 ms a drag feels sticky.");
END;

// The Enter / Menu choices.  Returns 0 to quit, 1 to carry on.
bnMenu()
BEGIN
  LOCAL muI, muS1, muS2, muS3;
  muS1 := "Contour view (from above)";
  IF bnView == 2 THEN muS1 := "3-D view"; END;
  muS2 := "Slice on";
  IF bnSlc == 1 THEN muS2 := "Slice off"; END;
  // beta 1a (N2): 26 characters, not 31 -- the emulator cut menu lines
  // of 31 to 36 characters at the box's edge
  muS3 := "Grid while dragging: 15" + CHAR(215) + "15";
  IF bnCrs == 1 THEN muS3 := "Grid while dragging: 11" + CHAR(215) + "11"; END;
  muI := 0;
  CHOOSE(muI, "BiNorm3D", {
    "Set " + bnRS + " exactly (now " + bnF(bnRho, 2) + ")",
    muS1,
    muS2,
    "Reset the view",
    muS3,
    "Frame timing (G1 test)",
    "Help",
    "Quit"});
  IF muI == 1 THEN bnAskRho(); END;
  IF muI == 2 THEN bnView := 3 - bnView; END;
  IF muI == 3 THEN bnSlc := 1 - bnSlc; END;
  IF muI == 4 THEN bnReset(); END;
  IF muI == 5 THEN bnCrs := 1 - bnCrs; END;
  IF muI == 6 THEN bnTiming(); END;
  IF muI == 7 THEN bnHelp(); END;
  IF muI == 8 THEN RETURN 0; END;
  RETURN 1;
END;


// ------------------------------------------------------------
// HELP: five drawn pages.  Any key or tap turns the page; Esc
// goes back to the surface.
// ------------------------------------------------------------

// A help page's frame and title.
bnPage(pgT, pgC, pgF)
BEGIN
  DIMGROB_P(G1, 320, 240, #FFFFFFh);
  RECT_P(G1, 4, 4, 315, 235, pgC, pgF);
  TEXTOUT_P(pgT, G1, 12, 8, 2, pgC, 296);
END;

// One line of help text (46 characters at most).
bnLn(lnS, lnY, lnC)
BEGIN
  TEXTOUT_P(lnS, G1, 12, lnY, 1, lnC, 296);
END;

// The footer, the blit and the wait; returns the key (-1 for a tap).
bnPgEnd(peS)
BEGIN
  TEXTOUT_P(peS, G1, 12, 220, 1, bnFOOT, 296);
  BLIT_P(G0, 0, 0, G1);
  RETURN bnWaitAny();
END;

bnHelp()
BEGIN
  LOCAL hpK;
  // ---- 1: what rho does ----
  bnPage("What " + bnRS + " does to the surface", bnBLUE, bnLBLU);
  bnLn("X and Y each follow N(0, 1); " + bnRS + " is their", 32, #000000h);
  bnLn("correlation.  The surface is their joint", 45, #000000h);
  bnLn("density f(x, y): the volume under it is 1.", 58, #000000h);
  bnLn(bnRS + " = 0: a round bell.  X says nothing about Y.", 75, bnBLUE);
  bnLn(bnRS + " > 0: it stretches along y = x (high x with", 88, bnBLUE);
  bnLn("high y); " + bnRS + " < 0: along y = " + bnMin + "x.", 101, bnBLUE);
  bnLn("The volume stays 1, so a narrower bell is a", 118, bnMAG);
  bnLn("taller one: the peak f(0, 0) is 0.159 at", 131, bnMAG);
  bnLn(bnRS + " = 0, 0.199 at 0.6 and 0.365 at 0.9.", 144, bnMAG);
  bnLn("Colour = height: blue low, then magenta, red", 161, bnGRAY);
  bnLn("at the top.  The height scale never changes,", 174, bnGRAY);
  bnLn("so the peaks can be compared across " + bnRS + ".", 187, bnGRAY);
  hpK := bnPgEnd("1/5, tap or key = next page, Esc = back");
  IF hpK == 4 THEN RETURN; END;
  // ---- 2: why the contours tilt ----
  bnPage("Why the contours tilt", bnGRN, bnLGRN);
  bnLn("Cont looks straight down.  Each ellipse is a", 32, #000000h);
  bnLn("level of the surface, Q = k" + bnSq + ", where", 45, #000000h);
  bnLn("Q = (x" + bnSq + " " + bnMin + " 2" + bnRS + "xy + y" + bnSq + ") / (1 " + bnMin + " " + bnRS + bnSq + ").", 58, #000000h);
  bnLn("The cross term " + bnMin + "2" + bnRS + "xy tilts them: along y = x", 75, bnBLUE);
  bnLn("the variance is 1 + " + bnRS + ", across it 1 " + bnMin + " " + bnRS + ".", 88, bnBLUE);
  bnLn("Q is chi-square with 2 df, so the share inside", 105, bnMAG);
  bnLn("is 1 " + bnMin + " e^(" + bnMin + "k" + bnSq + "/2): 39 % inside k = 1 and", 118, bnMAG);
  bnLn("86 % inside k = 2, whatever " + bnRS + " is.", 131, bnMAG);
  bnLn("The magenta line E[Y | x] = " + bnRS + "x is the", 148, bnGRAY);
  bnLn("regression line.  It runs through the left", 161, bnGRAY);
  bnLn("and right ends of every ellipse, not along", 174, bnGRAY);
  bnLn("the long axis.", 187, bnGRAY);
  hpK := bnPgEnd("2/5, tap or key = next page, Esc = back");
  IF hpK == 4 THEN RETURN; END;
  // ---- 3: what a slice is ----
  bnPage("What a slice is", bnMAG, bnLMAG);
  bnLn("Slice cuts the surface along the line x = " + bnX0S + ".", 32, #000000h);
  bnLn("Scaled to area 1, the cut is the distribution", 45, #000000h);
  bnLn("of Y among the cases that have X = " + bnX0S + ":", 58, #000000h);
  bnLn("Y | X = " + bnX0S + " is N(" + bnRS + bnX0S + ", (1 " + bnMin + " " + bnRS + bnSq + ")^0.5).", 75, bnMAG);
  bnLn("Its centre " + bnRS + bnX0S + " lies on the regression line", 88, bnMAG);
  bnLn("y = " + bnRS + "x; the SD is the same for every " + bnX0S + ".", 101, bnMAG);
  bnLn("The ridge of the bell runs along y = x, but", 118, bnRED);
  bnLn("the slices peak on the flatter magenta line.", 131, bnRED);
  bnLn("The panel: f(y | " + bnX0S + ") shaded, N(0, 1) grey for", 148, bnGRAY);
  bnLn("comparison -- narrower, and centred at " + bnRS + bnX0S + ".", 161, bnGRAY);
  bnLn("Drag the red line, or use Left and Right.", 174, bnGRAY);
  hpK := bnPgEnd("3/5, tap or key = next page, Esc = back");
  IF hpK == 4 THEN RETURN; END;
  // ---- 4: regression to the mean ----
  bnPage("Regression to the mean", bnRED, bnLRED);
  bnLn("Pick the cases with an extreme " + bnX0S + ", say 2 SD", 32, #000000h);
  bnLn("above the mean.  Their average y is " + bnRS + bnX0S + ":", 45, #000000h);
  bnLn("1.2 SD at " + bnRS + " = 0.6 -- nearer the mean, though", 58, #000000h);
  bnLn("nothing was done to them.  With |" + bnRS + "| < 1, " + bnRS + bnX0S, 71, #000000h);
  bnLn("is always pulled toward 0; the weaker " + bnRS + ",", 84, #000000h);
  bnLn("the stronger the pull.", 97, #000000h);
  // beta 1a (N1): the sentence had lost its verb ("... risk score
  // was extreme score closer ..."); "will" restored, the rest reflowed,
  // still six lines of at most 46 characters
  bnLn("Criminal justice: offenders chosen for a", 114, bnBLUE);
  bnLn("programme because their risk score was", 127, bnBLUE);
  bnLn("extreme will score closer to the average", 140, bnBLUE);
  bnLn("next time, programme or not.  That drop is", 153, bnBLUE);
  bnLn("not the programme's effect; only a comparison", 166, bnBLUE);
  bnLn("group picked the same way shows what it added.", 179, bnBLUE);
  hpK := bnPgEnd("4/5, tap or key = next page, Esc = back");
  IF hpK == 4 THEN RETURN; END;
  // ---- 5: keys and touch ----
  bnPage("Keys and touch", bnORNG, bnLORN);
  bnLn("Drag: sideways turns, downward raises the eye.", 32, #000000h);
  bnLn("Arrows: Left, Right turn; Up, Down raise and", 45, #000000h);
  bnLn("lower the eye; 10" + bnDg + " a step.", 58, #000000h);
  bnLn(bnRS + bnMin + " " + bnRS + "+ (or the " + bnMin + " and + keys): " + bnRS + " by 0.1.", 71, #000000h);
  bnLn("Cont: the view from above; 3-D goes back.", 84, #000000h);
  bnLn("Slice: on or off.  From above, tap or drag", 97, #000000h);
  bnLn("the plot (or Left, Right) to move " + bnX0S + ".", 110, #000000h);
  bnLn("Reset: " + bnTS + " = 35" + bnDg + ", " + bnFS + " = 30" + bnDg + ", " + bnRS + " = 0.60.", 123, #000000h);
  bnLn("Menu (Enter): set " + bnRS + " exactly, the views, the", 140, bnGRAY);
  bnLn("slice, reset, the grid while dragging, frame", 153, bnGRAY);
  bnLn("timing, help, quit.  Esc = quit.", 166, bnGRAY);
  bnLn("A drag draws an 11" + CHAR(215) + "11 grid, 15" + CHAR(215) + "15 on lifting.", 183, bnGRAY);
  bnLn("Home after a run: BN_RHO holds " + bnRS + ".", 200, bnGRN);
  bnPgEnd("5/5, tap or key = back to the surface");
END;


// ------------------------------------------------------------
// INIT AND MAIN
// ------------------------------------------------------------
bnInit()
BEGIN
  LOCAL ivI, ivT, ivU, ivY;
  // ---- house colours ----
  bnBLUE := #003CC8h; bnMAG := #C800C8h; bnRED := #FF0000h;
  bnGRN := #008C3Ch; bnORNG := #FF4500h; bnGRAY := #5A5A5Ah;
  bnFOOT := #1E90FFh; bnGRID := #E1E4EBh; bnLBLU := #EBF2FFh;
  bnLMAG := #FAEBFAh; bnLORN := #FFF0DCh; bnLRED := #FFEBEBh;
  bnLGRN := #E8F6EBh;
  bnFLR := #A0A0A0h;                     // floor and axes: a mid grey
  // ---- display glyphs.  If one shows as a box on the calculator,
  //      change ONLY this line: "r", "th", "ph", " deg". ----
  bnRS := CHAR(961); bnTS := CHAR(952); bnFS := CHAR(966); bnDg := CHAR(176);
  // verified glyphs: superscript 2, minus sign, subscript 0
  bnSq := CHAR(178); bnMin := CHAR(8722);
  bnX0S := "x" + CHAR(8320);
  // ---- constants and the angle unit ----
  bnPI := 3.14159265358979;
  bnDR := bnPI / 180;
  // ASIN(1) is a right angle in the user's own unit: 90 (degrees),
  // 1.5707963267949 (radians) or 100 (grads).  Radians times bnAU
  // is therefore what SIN and COS expect, whatever the setting.
  bnAU := ASIN(1) / 1.5707963267949;
  // ---- geometry ----
  bnS := 22; bnK := 24; bnCX3 := 160;
  bnR3 := 3 * 2^0.5 * bnS;                       // half the floor's diagonal
  bnHM := bnS * bnK / (2 * bnPI * 0.19^0.5);     // the rho = 0.9 peak, px
  bnBT := 3; bnBB := 215;
  bnCL := 22; bnCT := 18; bnCU := 27;
  bnCR := bnCL + 6 * bnCU; bnCB := bnCT + 6 * bnCU;
  bnCXc := bnCL + 3 * bnCU; bnCYc := bnCT + 3 * bnCU;
  bnPL := 196; bnPW := 104;
  // ---- the grid values: 15 fine (steps of 3/7), 11 coarse (0.6) ----
  bnGV := MAKELIST(0, K, 1, 26, 1);
  FOR ivI FROM 1 TO 15 DO bnGV(ivI) := -3 + (ivI - 1) * 6 / 14; END;
  FOR ivI FROM 1 TO 11 DO bnGV(15 + ivI) := -3 + (ivI - 1) * 0.6; END;
  // ---- room for the cache and the frame lists ----
  bnZ := MAKELIST(0, K, 1, 346, 1);
  bnSCX := MAKELIST(0, K, 1, 320, 1); bnSCY := MAKELIST(0, K, 1, 320, 1);
  bnPX := MAKELIST(0, K, 1, 225, 1); bnPY := MAKELIST(0, K, 1, 225, 1);
  bnWA := MAKELIST(0, K, 1, 15, 1); bnWB := MAKELIST(0, K, 1, 15, 1);
  bnWC := MAKELIST(0, K, 1, 15, 1); bnWD := MAKELIST(0, K, 1, 15, 1);
  bnEX := MAKELIST(0, K, 1, 122, 1); bnEY := MAKELIST(0, K, 1, 122, 1);
  bnCnd := MAKELIST(0, K, 1, 121, 1); bnSlZ := MAKELIST(0, K, 1, 29, 1);
  // ---- cos t and sin t round an ellipse, 6 degrees apart ----
  bnCos := MAKELIST(0, K, 1, 61, 1); bnSin := MAKELIST(0, K, 1, 61, 1);
  FOR ivI FROM 1 TO 61 DO
    ivT := 2 * bnPI * (ivI - 1) / 60;
    bnCos(ivI) := COS(ivT * bnAU);
    bnSin(ivI) := SIN(ivT * bnAU);
  END;
  // ---- the height colours: blue (0,60,200) to magenta (200,0,200)
  //      to red (255,0,0), in 16 steps ----
  bnRamp := MAKELIST(0, K, 1, 16, 1);
  FOR ivI FROM 1 TO 16 DO
    ivT := (ivI - 1) / 15;
    IF ivT <= 0.5 THEN
      ivU := 2 * ivT;
      bnRamp(ivI) := RGB(IP(200 * ivU + 0.5), IP(60 * (1 - ivU) + 0.5), 200);
    ELSE
      ivU := 2 * ivT - 1;
      bnRamp(ivI) := RGB(IP(200 + 55 * ivU + 0.5), 0, IP(200 * (1 - ivU) + 0.5));
    END;
  END;
  // ---- N(0, 1) at every pixel row of the contour plot (fixed) ----
  bnMg := MAKELIST(0, K, 1, bnCB - bnCT + 1, 1);
  FOR ivI FROM bnCT TO bnCB DO
    ivY := (bnCYc - ivI) / bnCU;
    bnMg(ivI - bnCT + 1) := EXP(-ivY * ivY / 2) / (2 * bnPI)^0.5;
  END;
  // ---- the starting state ----
  bnView := 1; bnSlc := 0; bnDrag := 0; bnCrs := 1; bnMsg := "";
  bnTX := -1; bnTY := -1;
  bnTh := 35 * bnDR; bnPh := 30 * bnDR;
  bnRho := 0.6; bnX0 := 1.5;
  bnCache();
  bnSliceC();
  BN_RHO := bnRho;
END;

// The opening page: what the picture is and how to touch it.
bnWelcome()
BEGIN
  DIMGROB_P(G1, 320, 240, #FFFFFFh);
  RECT_P(G1, 4, 4, 315, 235, bnORNG, bnLORN);
  TEXTOUT_P("Bivariate normal", G1, 14, 10, 3, bnBLUE, 292);
  TEXTOUT_P("Two variables X and Y, each N(0, 1), with", G1, 14, 42, 1, #000000h, 292);
  TEXTOUT_P("correlation " + bnRS + ".  The surface shows how often", G1, 14, 55, 1, #000000h, 292);
  TEXTOUT_P("each pair (x, y) turns up.", G1, 14, 68, 1, #000000h, 292);
  TEXTOUT_P("Drag a finger to turn it; the arrows work too.", G1, 14, 88, 1, bnGRN, 292);
  TEXTOUT_P(bnRS + bnMin + " and " + bnRS + "+ change the correlation.", G1, 14, 101, 1, bnGRN, 292);
  TEXTOUT_P("Cont: the view from above, with the ellipses", G1, 14, 121, 1, bnMAG, 292);
  TEXTOUT_P("and the regression line.  Slice: cut the bell", G1, 14, 134, 1, bnMAG, 292);
  TEXTOUT_P("at x = " + bnX0S + " and watch REGRESSION TO THE MEAN.", G1, 14, 147, 1, bnMAG, 292);
  TEXTOUT_P("Enter: menu, with five pages of help.", G1, 14, 167, 1, bnGRAY, 292);
  TEXTOUT_P("Any key or tap to begin", G1, 14, 214, 2, bnFOOT, 292);
  BLIT_P(G0, 0, 0, G1);
  bnWaitAny();
END;


// ------------------------------------------------------------
// MAIN: the live surface
// ------------------------------------------------------------
EXPORT BiNorm3D()
BEGIN
  LOCAL mnK, mnM, mnM1, mnGo, mnAct, mnX;
  LOCAL mnY;
  bnInit();
  bnWelcome();
  mnGo := 1;
  WHILE mnGo == 1 DO
    // a key or a finger left over from a dialog or a drag must not
    // count as a new press: drain both before the next frame
    REPEAT mnK := GETKEY; UNTIL mnK == -1;
    bnDrain();
    bnShow();
    bnMsg := "";                       // a message lasts one frame
    REPEAT
      mnK := GETKEY;
      mnM := MOUSE;
      mnM1 := mnM(1);                  // the first finger, if any
    UNTIL mnK > -1 OR SIZE(mnM1) > 0;
    // mnAct: 1 rho down, 2 rho up, 3 Cont / 3-D, 4 Slice, 5 Reset,
    // 6 Menu; 0 nothing
    mnAct := 0;
    IF SIZE(mnM1) > 0 THEN
      mnX := mnM1(1); mnY := mnM1(2);
      IF mnY >= 218 THEN
        // the soft-key strip: six 53-px slots, left to right
        mnAct := MIN(6, IP(mnX / 53.4) + 1);
        bnDrain();
      ELSE
        IF bnView == 1 THEN bnDrag3D(mnX, mnY); ELSE bnDragC(mnX, mnY); END;
      END;
    ELSE
      // the keys remain as shortcuts for the soft keys
      IF mnK == 4 THEN mnGo := 0; END;                 // Esc   = quit
      IF mnK == 45 THEN mnAct := 1; END;               // -     = rho down
      IF mnK == 50 THEN mnAct := 2; END;               // +     = rho up
      IF mnK == 30 THEN mnAct := 6; END;               // Enter = Menu
      IF bnView == 1 THEN
        // arrows turn (theta) and tilt (phi) in 10-degree steps
        IF mnK == 7 THEN bnTh := bnTh - 10 * bnDR; bnNormTh(); END;
        IF mnK == 8 THEN bnTh := bnTh + 10 * bnDR; bnNormTh(); END;
        IF mnK == 2 THEN bnPh := MIN(80 * bnDR, bnPh + 10 * bnDR); END;
        IF mnK == 12 THEN bnPh := MAX(10 * bnDR, bnPh - 10 * bnDR); END;
      ELSE
        // from above, Left and Right move the slice
        IF mnK == 7 THEN bnSlc := 1; bnSetX0(bnX0 - 0.25); END;
        IF mnK == 8 THEN bnSlc := 1; bnSetX0(bnX0 + 0.25); END;
      END;
    END;
    IF mnAct == 1 THEN bnSetRho(bnRho - 0.1); END;
    IF mnAct == 2 THEN bnSetRho(bnRho + 0.1); END;
    IF mnAct == 3 THEN bnView := 3 - bnView; END;
    IF mnAct == 4 THEN bnSlc := 1 - bnSlc; END;
    IF mnAct == 5 THEN bnReset(); END;
    IF mnAct == 6 THEN mnGo := bnMenu(); END;
  END;
  BN_RHO := bnRho;
  // Hand the summary back as the program's value: Home shows it in
  // the history and no Terminal screen opens.  BN_RHO holds rho.
  RETURN "BiNorm3D: " + bnRS + " = " + bnF(bnRho, 2) + ", E[Y | X = x] = " + bnF(bnRho, 2) + " x; " + bnRS + " in BN_RHO";
END;
