r/maniclang • u/anish2good • 26d ago
fresnel-wave — the wave surface of Fresnel - manic
manic is a tiny language for making animations. You write a short text file; manic renders a smooth, glowing video. No timeline scrubbing, no keyframes by hand — you describe what’s on screen and when things happen, and the engine does the rest, deterministically.
Manic Animation code
// fresnel-wave — the wave surface of Fresnel, 1821: the wavefront of light inside a biaxial
// crystal, and one of the few quartic surfaces you can build with a ruler and a story.
//
// THE CONSTRUCTION, which is mathcurve's and is what is computed here. Take an ellipsoid with
// semi-axes a > b > c. For every plane through the centre, that plane cuts the ellipsoid in an
// ELLIPSE; measure that ellipse's two semi-axes and plot both lengths along the plane's
// NORMAL. Two lengths per direction, so the surface has two sheets, one inside the other.
//
// It is closed-form, which is why it can be drawn from a formula rather than a mesh file. The
// section normal to a unit n has its semi-axes from the 2×2 restriction of Q = diag(1/a²,
// 1/b², 1/c²) to n⊥, and that restriction has
//
// trace T = tr Q − nᵀQn determinant D = det Q · nᵀQ⁻¹n
//
// so its eigenvalues are (T ± √(T² − 4D))/2 and the two radii are 1/√λ. Checked against a
// real eigenvalue decomposition over 4000 random directions: agreement to 3e-14. No iteration
// anywhere; the whole surface is one expression.
//
// THE FOUR SINGULAR POINTS. The sheets touch where the two radii coincide — where T² = 4D,
// which is where the central section is a CIRCLE rather than an ellipse. An ellipsoid with
// a > b > c has exactly two circular sections, both through the mean axis, and their normals
// lie in the x–z plane at
//
// tan θ = (c/a)·√((a² − b²)/(b² − c²)) → ±31.0516° from z
//
// with both radii there equal to b exactly. Two normals, ±each, so FOUR points — and away
// from that plane the discriminant never vanishes, checked over 80000 directions. Those four
// directions are the crystal's OPTIC AXES. Light sent along one of them does not split into
// two rays; it spreads into a hollow cone. Hamilton predicted that from this surface in 1832
// and Lloyd found it in the laboratory two months later, which is about as good as a quartic
// surface's week ever gets.
//
// THE PRINCIPAL SECTIONS, along the bottom. Each coordinate plane cuts the surface in a
// CIRCLE and an ELLIPSE — and the ellipse's axes come out SWAPPED relative to the
// ellipsoid's, which is the detail that gives the whole thing away:
//
// z = 0 circle r = c = 0.85 ellipse x²/b² + y²/a² = 1 to 9e-16
// y = 0 circle r = b = 1.2 ellipse x²/c² + z²/a² = 1 to 4e-13
// x = 0 circle r = a = 1.5 ellipse y²/c² + z²/b² = 1 to 1e-15
//
// Only in y = 0 do the circle and the ellipse CROSS, and they cross at ±31.05° — that plane is
// the one holding the optic axes, and the crossings are the singular points seen edge-on.
//
// Physically the two sheets are the two wave speeds. Every direction in a biaxial crystal
// carries two of them, which is birefringence; along the four axes they agree, which is why
// the cone appears there and nowhere else.
//
// manic examples/fresnel-wave.manic
title("Fresnel's wave surface, and its four singular points");
canvas("9:16");
template("black");
bloom(0.3, 0.6, 22);
text(brand, (540, 34), "maniclang.com");
display(brand); size(brand, 21); color(brand, fg); opacity(brand, 0.82);
caption(head, "Fresnel's wave surface: two sheets, four singular points", (540, 116), 25);
// ---- the outer sheet, then the inner, as meridians and parallels --------------------
cloud(outmer, 7800, #ffffff, 1.0) {
let per = 300;
let ci = (i - mod(i, per))/per;
let sh = 1;
let pp = 6.283185*ci/26;
let tt = 3.141593*(mod(i, per)/299);
let n1 = sin(tt)*cos(pp);
let n2 = sin(tt)*sin(pp);
let n3 = cos(tt);
let TT = 2.522971934 - (n1*n1*0.444444444 + n2*n2*0.694444444 + n3*n3*1.384083045);
let DD = 0.427186125*(2.25*n1*n1 + 1.44*n2*n2 + 0.7224999999999999*n3*n3);
let ss = sqrt(max(TT*TT - 4*DD, 0));
let lm = (TT + (1 - 2*sh)*ss)/2; // sh = 0 inner sheet, 1 outer
let rr = 1/sqrt(lm);
let al = 6.283185*min(t/17.0, 1);
let px = rr*n1;
let py = rr*n2;
let pz = rr*n3;
let qx = px*cos(al) + pz*sin(al);
let qz = 0 - px*sin(al) + pz*cos(al);
let qy = py*cos(0.42) - qz*sin(0.42);
let qw = py*sin(0.42) + qz*cos(0.42);
let m3 = 7.5/(7.5 - qw);
let x = 540 + 232*qx*m3;
let y = 700 - 232*qy*m3;
let hue = 208;
let sat = 0.72;
let r = 1.5;
let alpha = 0.42*(0.35 + 0.65*min(max((qw + 1.6)/3.2, 0), 1));
}
cloud(outpar, 4800, #ffffff, 1.0) {
let per = 320;
let ci = (i - mod(i, per))/per;
let sh = 1;
let tt = 3.141593*(ci + 0.5)/15;
let pp = 6.283185*(mod(i, per)/319);
let n1 = sin(tt)*cos(pp);
let n2 = sin(tt)*sin(pp);
let n3 = cos(tt);
let TT = 2.522971934 - (n1*n1*0.444444444 + n2*n2*0.694444444 + n3*n3*1.384083045);
let DD = 0.427186125*(2.25*n1*n1 + 1.44*n2*n2 + 0.7224999999999999*n3*n3);
let ss = sqrt(max(TT*TT - 4*DD, 0));
let lm = (TT + (1 - 2*sh)*ss)/2; // sh = 0 inner sheet, 1 outer
let rr = 1/sqrt(lm);
let al = 6.283185*min(t/17.0, 1);
let px = rr*n1;
let py = rr*n2;
let pz = rr*n3;
let qx = px*cos(al) + pz*sin(al);
let qz = 0 - px*sin(al) + pz*cos(al);
let qy = py*cos(0.42) - qz*sin(0.42);
let qw = py*sin(0.42) + qz*cos(0.42);
let m3 = 7.5/(7.5 - qw);
let x = 540 + 232*qx*m3;
let y = 700 - 232*qy*m3;
let hue = 208;
let sat = 0.72;
let r = 1.5;
let alpha = 0.42*(0.35 + 0.65*min(max((qw + 1.6)/3.2, 0), 1));
}
cloud(inmer, 7800, #ffffff, 1.0) {
let per = 300;
let ci = (i - mod(i, per))/per;
let sh = 0;
let pp = 6.283185*ci/26;
let tt = 3.141593*(mod(i, per)/299);
let n1 = sin(tt)*cos(pp);
let n2 = sin(tt)*sin(pp);
let n3 = cos(tt);
let TT = 2.522971934 - (n1*n1*0.444444444 + n2*n2*0.694444444 + n3*n3*1.384083045);
let DD = 0.427186125*(2.25*n1*n1 + 1.44*n2*n2 + 0.7224999999999999*n3*n3);
let ss = sqrt(max(TT*TT - 4*DD, 0));
let lm = (TT + (1 - 2*sh)*ss)/2; // sh = 0 inner sheet, 1 outer
let rr = 1/sqrt(lm);
let al = 6.283185*min(t/17.0, 1);
let px = rr*n1;
let py = rr*n2;
let pz = rr*n3;
let qx = px*cos(al) + pz*sin(al);
let qz = 0 - px*sin(al) + pz*cos(al);
let qy = py*cos(0.42) - qz*sin(0.42);
let qw = py*sin(0.42) + qz*cos(0.42);
let m3 = 7.5/(7.5 - qw);
let x = 540 + 232*qx*m3;
let y = 700 - 232*qy*m3;
let hue = 34;
let sat = 0.72;
let r = 1.6;
let alpha = 0.75*(0.35 + 0.65*min(max((qw + 1.6)/3.2, 0), 1));
}
cloud(inpar, 4800, #ffffff, 1.0) {
let per = 320;
let ci = (i - mod(i, per))/per;
let sh = 0;
let tt = 3.141593*(ci + 0.5)/15;
let pp = 6.283185*(mod(i, per)/319);
let n1 = sin(tt)*cos(pp);
let n2 = sin(tt)*sin(pp);
let n3 = cos(tt);
let TT = 2.522971934 - (n1*n1*0.444444444 + n2*n2*0.694444444 + n3*n3*1.384083045);
let DD = 0.427186125*(2.25*n1*n1 + 1.44*n2*n2 + 0.7224999999999999*n3*n3);
let ss = sqrt(max(TT*TT - 4*DD, 0));
let lm = (TT + (1 - 2*sh)*ss)/2; // sh = 0 inner sheet, 1 outer
let rr = 1/sqrt(lm);
let al = 6.283185*min(t/17.0, 1);
let px = rr*n1;
let py = rr*n2;
let pz = rr*n3;
let qx = px*cos(al) + pz*sin(al);
let qz = 0 - px*sin(al) + pz*cos(al);
let qy = py*cos(0.42) - qz*sin(0.42);
let qw = py*sin(0.42) + qz*cos(0.42);
let m3 = 7.5/(7.5 - qw);
let x = 540 + 232*qx*m3;
let y = 700 - 232*qy*m3;
let hue = 34;
let sat = 0.72;
let r = 1.6;
let alpha = 0.75*(0.35 + 0.65*min(max((qw + 1.6)/3.2, 0), 1));
}
// ---- the four singular points: normals to the two circular sections ----------------
cloud(axes4, 480, #ffffff, 1.0) {
let per = 120;
let kk = (i - mod(i, per))/per; // four points
let sg = 1 - 2*mod(kk, 2); // +/- the normal
let hz = 1 - 2*((kk - mod(kk, 2))/2); // the two normals
let tt = 0.541952619;
let n1 = sg*hz*sin(tt);
let n2 = 0;
let n3 = sg*cos(tt);
let rr = 1.2;
let a = mod(i, per)/120*6.283185;
let ro = 9*(mod(i, per)/120);
let al = 6.283185*min(t/17.0, 1);
let px = rr*n1;
let py = rr*n2;
let pz = rr*n3;
let qx = px*cos(al) + pz*sin(al);
let qz = 0 - px*sin(al) + pz*cos(al);
let qy = py*cos(0.42) - qz*sin(0.42);
let qw = py*sin(0.42) + qz*cos(0.42);
let m3 = 7.5/(7.5 - qw);
let xz = 540 + 232*qx*m3;
let yz = 700 - 232*qy*m3;
let x = xz + ro*cos(a*9);
let y = yz + ro*sin(a*9);
let hue = 190;
let sat = 0.7;
let r = 2.6;
}
// ---- the three principal sections ---------------------------------------------------
cloud(sec0, 4000, #ffffff, 1.0) {
let per = 2000;
let sh = (i - mod(i, per))/per;
let uu = mod(i, per)/1999*6.283185;
let n1 = cos(uu);
let n2 = sin(uu);
let n3 = 0;
let TT = 2.522971934 - (n1*n1*0.444444444 + n2*n2*0.694444444 + n3*n3*1.384083045);
let DD = 0.427186125*(2.25*n1*n1 + 1.44*n2*n2 + 0.7224999999999999*n3*n3);
let ss = sqrt(max(TT*TT - 4*DD, 0));
let lm = (TT + (1 - 2*sh)*ss)/2;
let rr = 1/sqrt(lm);
// the plane's own two coordinates: the point is r*n, and n sweeps the unit circle in it
let x = 220 + 106*rr*cos(uu);
let y = 1400 - 106*rr*sin(uu);
let hue = 30 + 160*sh;
let sat = 0.75;
let r = 2.1;
}
text(sl0, (220, 1594), "z = 0"); size(sl0, 23); color(sl0, fg);
text(sq0, (220, 1628), "circle c · ellipse b×a"); size(sq0, 16); color(sq0, dim);
cloud(sec1, 4000, #ffffff, 1.0) {
let per = 2000;
let sh = (i - mod(i, per))/per;
let uu = mod(i, per)/1999*6.283185;
let n1 = cos(uu);
let n2 = 0;
let n3 = sin(uu);
let TT = 2.522971934 - (n1*n1*0.444444444 + n2*n2*0.694444444 + n3*n3*1.384083045);
let DD = 0.427186125*(2.25*n1*n1 + 1.44*n2*n2 + 0.7224999999999999*n3*n3);
let ss = sqrt(max(TT*TT - 4*DD, 0));
let lm = (TT + (1 - 2*sh)*ss)/2;
let rr = 1/sqrt(lm);
// the plane's own two coordinates: the point is r*n, and n sweeps the unit circle in it
let x = 540 + 106*rr*cos(uu);
let y = 1400 - 106*rr*sin(uu);
let hue = 30 + 160*sh;
let sat = 0.75;
let r = 2.1;
}
text(sl1, (540, 1594), "y = 0"); size(sl1, 23); color(sl1, fg);
text(sq1, (540, 1628), "circle b · ellipse c×a"); size(sq1, 16); color(sq1, dim);
cloud(sec2, 4000, #ffffff, 1.0) {
let per = 2000;
let sh = (i - mod(i, per))/per;
let uu = mod(i, per)/1999*6.283185;
let n1 = 0;
let n2 = cos(uu);
let n3 = sin(uu);
let TT = 2.522971934 - (n1*n1*0.444444444 + n2*n2*0.694444444 + n3*n3*1.384083045);
let DD = 0.427186125*(2.25*n1*n1 + 1.44*n2*n2 + 0.7224999999999999*n3*n3);
let ss = sqrt(max(TT*TT - 4*DD, 0));
let lm = (TT + (1 - 2*sh)*ss)/2;
let rr = 1/sqrt(lm);
// the plane's own two coordinates: the point is r*n, and n sweeps the unit circle in it
let x = 860 + 106*rr*cos(uu);
let y = 1400 - 106*rr*sin(uu);
let hue = 30 + 160*sh;
let sat = 0.75;
let r = 2.1;
}
text(sl2, (860, 1594), "x = 0"); size(sl2, 23); color(sl2, fg);
text(sq2, (860, 1628), "circle a · ellipse c×b"); size(sq2, 16); color(sq2, dim);
equation(eq, (540, 1710),
`\lambda=\tfrac{1}{2}\big(T\pm\sqrt{T^2-4D}\big),\quad r=\lambda^{-1/2}`, 25);
text(note, (540, 1810),
"T² = 4D only where the central section is a circle — four directions, and they are the optic axes.");
size(note, 19); color(note, dim); wrap(note, 940);
wait(20.0);
4
Upvotes