Transit Method: Finding Exoplanets

// Exoplanet transit-method simulation.
// Top panel: edge-on view of a star with limb darkening and a planet
// orbiting in a circle; transits occur when the planet's silhouette
// crosses the stellar disk along the line of sight.
// Bottom panel: scrolling light curve of measured flux, showing
// U-shaped transit dips. mouseY scrubs the planet radius.

let W, H;
let topH, botH;
let starCx, starCy, starR;
let orbitA; // orbit semi-major axis in screen px

let phase; // orbit phase in [0, 2*pi)
let period; // seconds per orbit
let timeSec;

let lightCurve; // ring buffer of recent flux samples
let lcHead;
const LC_LEN = 600;
const LC_DT = 1 / 60; // each sample = one frame at 60fps

// Planet radius is scrubbed via mouseY.
// Rp/Rs ranges from ~0.005 (sub-Earth) up to ~0.5 (deep eclipse).
let rpOverRs;

// Photometric noise sigma (fractional flux).
const NOISE = 0.0015;

const LIMB_U1 = 0.5;
const LIMB_U2 = 0.2;

function computeFlux(planetX, planetY, planetZ, rp) {
  // The "observer" is at +infinity along +z (out of the screen).
  // A transit happens when planetZ > 0 (in front of the star) and
  // the projected (planetX, planetY) lies within the stellar disk.
  //
  // We integrate over the stellar disk with a quadratic limb-darkening
  // intensity profile and subtract the planet's covered area.

  // Distance from planet center to star center in the plane of sky:
  const dx = planetX - starCx;
  const dy = planetY - starCy;
  const d = Math.sqrt(dx * dx + dy * dy);

  // Normalize to stellar radius
  const dN = d / starR;
  const rpN = rp / starR;

  // If planet is behind the star or fully off-disk, no occultation.
  if (planetZ < 0) return 1.0;
  if (dN > 1 + rpN) return 1.0;

  // Total emitted flux of the star with quadratic limb darkening:
  //   I(mu) = 1 - u1(1-mu) - u2(1-mu)^2,  mu = cos(theta) = sqrt(1-r^2)
  // Total integrated flux Ftot = integral over disk of I dA, normalize so 1.
  // We'll Monte-Carlo / numerical-grid sample over the planet's covering
  // region for the lost flux. For real-time, use a small radial-bin grid
  // sampled over the planet disk.

  // Quick path: if planet fully on disk, sample over its area.
  // Sample using a polar grid centered on the planet, radius 0..rpN.
  const NR = 6;
  const NTH = 16;
  let lost = 0;
  let weight = 0;
  // We'll evaluate I(r_star) at each sample point inside the planet
  // that is also inside the star.
  for (let i = 0; i < NR; i++) {
    const rr = ((i + 0.5) / NR) * rpN;
    const dA = rr; // polar weight (constant times dr*dtheta)
    for (let j = 0; j < NTH; j++) {
      const th = (j / NTH) * 2 * Math.PI;
      const sx = dN + rr * Math.cos(th);
      const sy = rr * Math.sin(th);
      const r2 = sx * sx + sy * sy;
      if (r2 < 1) {
        const mu = Math.sqrt(1 - r2);
        const I = 1 - LIMB_U1 * (1 - mu) - LIMB_U2 * (1 - mu) * (1 - mu);
        lost += I * dA;
      }
      weight += dA;
    }
  }
  // Normalize the lost flux by the disk's mean intensity.
  // Mean intensity over the unit disk for quadratic LD:
  //   <I> = 1 - u1/3 - u2/6
  const meanI = 1 - LIMB_U1 / 3 - LIMB_U2 / 6;
  // Fraction of total stellar flux that is blocked.
  // Planet's projected area in units of stellar area = rpN^2.
  // Blocked area sampled at intensity proportional to lost/weight.
  const sampledI = weight > 0 ? lost / weight : 0;
  const blockedAreaFrac = rpN * rpN; // (Rp/Rs)^2
  const dF = blockedAreaFrac * (sampledI / meanI);
  return Math.max(0, 1 - dF);
}

function drawStar(ctx) {
  // Soft glow
  const glow = ctx.createRadialGradient(starCx, starCy, starR * 0.9, starCx, starCy, starR * 2.0);
  glow.addColorStop(0, 'rgba(255, 220, 140, 0.35)');
  glow.addColorStop(1, 'rgba(255, 220, 140, 0)');
  ctx.fillStyle = glow;
  ctx.fillRect(starCx - starR * 2.2, starCy - starR * 2.2, starR * 4.4, starR * 4.4);

  // Disk with limb darkening: build a small radial gradient.
  const grad = ctx.createRadialGradient(starCx, starCy, 0, starCx, starCy, starR);
  // Approximate quadratic LD as a few stops.
  for (let i = 0; i <= 6; i++) {
    const t = i / 6;
    // t = r/R, so mu = sqrt(1 - t^2)
    const mu = Math.sqrt(Math.max(0, 1 - t * t));
    const I = 1 - LIMB_U1 * (1 - mu) - LIMB_U2 * (1 - mu) * (1 - mu);
    const v = Math.max(0, Math.min(1, I));
    // Star color: warm yellow-white at center, deeper orange at limb.
    const r = Math.round(255 * v);
    const g = Math.round((220 + 35 * v) * v + (1 - v) * 140);
    const b = Math.round(180 * v);
    grad.addColorStop(t, `rgb(${r}, ${Math.min(255, g)}, ${b})`);
  }
  ctx.fillStyle = grad;
  ctx.beginPath();
  ctx.arc(starCx, starCy, starR, 0, Math.PI * 2);
  ctx.fill();
}

function drawPlanet(ctx, px, py, pz, rp) {
  // Only draw the planet silhouette when it's in front of the star (pz>0).
  // When pz < 0 (behind), we hide it.
  if (pz < 0) return;
  ctx.fillStyle = '#0b0c12';
  ctx.beginPath();
  ctx.arc(px, py, rp, 0, Math.PI * 2);
  ctx.fill();
  // Thin rim to make it readable on the bright disk
  ctx.strokeStyle = 'rgba(0,0,0,0.85)';
  ctx.lineWidth = 1;
  ctx.stroke();
}

function drawOrbit(ctx) {
  ctx.strokeStyle = 'rgba(120, 140, 180, 0.18)';
  ctx.setLineDash([4, 6]);
  ctx.lineWidth = 1;
  ctx.beginPath();
  // Edge-on: orbit is a horizontal line through the star center.
  ctx.moveTo(starCx - orbitA, starCy);
  ctx.lineTo(starCx + orbitA, starCy);
  ctx.stroke();
  ctx.setLineDash([]);
}

function drawLightCurve(ctx) {
  const x0 = 10;
  const y0 = topH + 8;
  const w = W - 20;
  const h = botH - 16;

  // Panel background
  ctx.fillStyle = '#08090d';
  ctx.fillRect(x0, y0, w, h);
  ctx.strokeStyle = 'rgba(120,140,180,0.25)';
  ctx.lineWidth = 1;
  ctx.strokeRect(x0 + 0.5, y0 + 0.5, w - 1, h - 1);

  // Determine vertical range. Show flux between (1 - depthMax*1.3) and 1.01.
  const expectedDepth = rpOverRs * rpOverRs;
  const depthMax = Math.max(0.0035, expectedDepth * 1.4);
  const fMin = 1 - depthMax;
  const fMax = 1 + 0.005;
  const yOf = (f) => y0 + h - ((f - fMin) / (fMax - fMin)) * h;

  // 1.00 reference line
  ctx.strokeStyle = 'rgba(180,200,230,0.35)';
  ctx.setLineDash([2, 4]);
  ctx.beginPath();
  ctx.moveTo(x0, yOf(1.0));
  ctx.lineTo(x0 + w, yOf(1.0));
  ctx.stroke();
  ctx.setLineDash([]);

  // Expected depth line
  if (expectedDepth > 0.0008) {
    ctx.strokeStyle = 'rgba(255,160,90,0.5)';
    ctx.setLineDash([3, 5]);
    ctx.beginPath();
    ctx.moveTo(x0, yOf(1 - expectedDepth));
    ctx.lineTo(x0 + w, yOf(1 - expectedDepth));
    ctx.stroke();
    ctx.setLineDash([]);
  }

  // Plot light curve. Newest sample at right edge.
  ctx.strokeStyle = '#9fd9ff';
  ctx.lineWidth = 1.25;
  ctx.beginPath();
  let started = false;
  for (let i = 0; i < LC_LEN; i++) {
    const idx = (lcHead + i) % LC_LEN;
    const f = lightCurve[idx];
    if (f <= 0) continue;
    const x = x0 + (i / (LC_LEN - 1)) * w;
    const y = yOf(f);
    if (!started) {
      ctx.moveTo(x, y);
      started = true;
    } else {
      ctx.lineTo(x, y);
    }
  }
  ctx.stroke();

  // Labels
  ctx.fillStyle = 'rgba(220,230,250,0.85)';
  ctx.font = '12px system-ui, sans-serif';
  ctx.textAlign = 'left';
  ctx.fillText('flux F', x0 + 6, y0 + 14);
  ctx.fillText('time →', x0 + w - 56, y0 + h - 6);

  // Depth label near the dashed orange line
  if (expectedDepth > 0.0008) {
    ctx.fillStyle = 'rgba(255,180,110,0.95)';
    const lblY = yOf(1 - expectedDepth);
    const lbl = `ΔF/F = ${(expectedDepth * 100).toFixed(3)}%`;
    ctx.fillText(lbl, x0 + 6, Math.max(y0 + 26, Math.min(y0 + h - 4, lblY - 4)));
  }

  // 1.00 label
  ctx.fillStyle = 'rgba(180,200,230,0.7)';
  ctx.fillText('1.00', x0 + w - 36, yOf(1.0) - 4);
}

function drawHUD(ctx) {
  // HUD lines in upper-right of top panel.
  const expectedDepth = rpOverRs * rpOverRs;
  ctx.fillStyle = 'rgba(220,230,250,0.92)';
  ctx.font = '13px ui-monospace, Menlo, monospace';
  ctx.textAlign = 'right';
  const xR = W - 12;
  let y = 18;
  const ln = (s) => { ctx.fillText(s, xR, y); y += 17; };
  ln(`R_p / R_s = ${rpOverRs.toFixed(3)}`);
  ln(`depth     = ${(expectedDepth * 100).toFixed(3)}%`);
  ln(`period T  = ${period.toFixed(2)} s`);
  ctx.textAlign = 'left';

  // Hint
  ctx.fillStyle = 'rgba(180,200,230,0.65)';
  ctx.font = '11px system-ui, sans-serif';
  ctx.fillText('drag Y to scrub planet size', 10, 16);

  // Mark transit period on the light curve with a faint horizontal bracket
  // (top of bottom panel)
  const x0 = 10;
  const y0 = topH + 8;
  const w = W - 20;
  // Pixels per second on the light curve: each sample is LC_DT; LC_LEN samples span w.
  const pxPerSec = w / (LC_LEN * LC_DT);
  const tBracket = period;
  if (tBracket * pxPerSec < w * 0.95) {
    const bx0 = x0 + 8;
    const by = y0 + 22;
    const bx1 = bx0 + tBracket * pxPerSec;
    ctx.strokeStyle = 'rgba(160,220,255,0.7)';
    ctx.lineWidth = 1;
    ctx.beginPath();
    ctx.moveTo(bx0, by - 4);
    ctx.lineTo(bx0, by + 4);
    ctx.moveTo(bx1, by - 4);
    ctx.lineTo(bx1, by + 4);
    ctx.moveTo(bx0, by);
    ctx.lineTo(bx1, by);
    ctx.stroke();
    ctx.fillStyle = 'rgba(160,220,255,0.95)';
    ctx.font = '12px system-ui, sans-serif';
    ctx.textAlign = 'left';
    ctx.fillText(`T = ${period.toFixed(2)}s`, bx0 + 4, by - 6);
  }
}

function layout() {
  topH = Math.round(H * 0.58);
  botH = H - topH;
  starR = Math.max(28, Math.min(W * 0.10, topH * 0.30));
  starCx = Math.round(W * 0.32);
  starCy = Math.round(topH * 0.5);
  orbitA = Math.min(W - starCx - 30, starR * 4.5);
}

function init({ canvas, ctx, width, height, input }) {
  W = width;
  H = height;
  layout();

  phase = -0.15; // start just before transit so a dip is visible quickly
  period = 6.0; // seconds per orbit
  timeSec = 0;

  lightCurve = new Float32Array(LC_LEN);
  for (let i = 0; i < LC_LEN; i++) lightCurve[i] = 1.0;
  lcHead = 0;

  rpOverRs = 0.10; // initial: roughly Jupiter/Sun ~ 0.1

  ctx.fillStyle = '#05060a';
  ctx.fillRect(0, 0, W, H);
}

function tick({ ctx, dt, width, height, input }) {
  if (width !== W || height !== H) {
    W = width;
    H = height;
    layout();
  }

  // Update Rp/Rs from mouseY. Top of canvas = small planet, bottom = huge.
  const my = input.mouseY;
  if (Number.isFinite(my) && my >= 0 && my <= H) {
    const t = Math.max(0, Math.min(1, my / H));
    // Log scale so the interesting Earth-Jupiter range (0.01-0.1) gets
    // most of the slider, with deep eclipses at the far end.
    const lo = Math.log(0.005);
    const hi = Math.log(0.5);
    rpOverRs = Math.exp(lo + t * (hi - lo));
  }

  // Advance orbit.
  timeSec += dt;
  phase += (2 * Math.PI / period) * dt;
  if (phase > Math.PI) phase -= 2 * Math.PI;

  // Planet position in 3D (edge-on orbit in the x-z plane).
  // Define x along the screen horizontal, z toward observer (+z = front).
  // When phase = 0, the planet is at +x (right of the star) behind no transit.
  // Pick orbit so transit happens when planet's projected x crosses the star
  // and z > 0. Use:
  //   xWorld = orbitA * sin(phase) (left-right)
  //   zWorld = orbitA * cos(phase) (toward/away from observer)
  // Then a transit occurs when cos(phase) > 0 (planet in front)
  // and |sin(phase)| * orbitA < starR + planetR_px.
  const xWorld = orbitA * Math.sin(phase);
  const zWorld = orbitA * Math.cos(phase);
  const planetX = starCx + xWorld;
  const planetY = starCy; // edge-on: stays on the same horizontal line
  // For drawing, scale planet radius in screen px from Rp/Rs.
  const planetRpx = rpOverRs * starR;

  // Compute flux (true model).
  const fluxTrue = computeFlux(planetX, planetY, zWorld, planetRpx);
  // Add measurement noise.
  // Box-Muller would be cleaner; cheap uniform sum is fine here.
  let noise = 0;
  for (let i = 0; i < 4; i++) noise += Math.random() - 0.5;
  noise *= NOISE; // sigma ~ NOISE * sqrt(4/12) = NOISE/sqrt(3)
  const fluxMeasured = fluxTrue + noise;

  // Push to ring buffer (newest at the "right edge" = (lcHead + LC_LEN - 1)).
  lightCurve[lcHead] = fluxMeasured;
  lcHead = (lcHead + 1) % LC_LEN;

  // --- Render ---
  // Top panel background
  ctx.fillStyle = '#03050b';
  ctx.fillRect(0, 0, W, topH);
  // Subtle starfield
  // (deterministic, low cost: skip per-frame redraw of stars by using a hash)
  ctx.fillStyle = 'rgba(255,255,255,0.6)';
  // Use a few pseudo-random stars
  for (let i = 0; i < 60; i++) {
    const sx = ((i * 9301 + 49297) % 233280) / 233280 * W;
    const sy = ((i * 67129 + 1009) % 233280) / 233280 * topH;
    const sr = ((i * 4111 + 7) % 100) / 100 * 0.9 + 0.2;
    ctx.fillRect(sx, sy, sr, sr);
  }

  drawOrbit(ctx);
  // Draw planet first if behind the star, then star, then planet if in front,
  // so silhouette occludes the disk properly.
  if (zWorld < 0) drawPlanet(ctx, planetX, planetY, zWorld, planetRpx);
  drawStar(ctx);
  if (zWorld >= 0) drawPlanet(ctx, planetX, planetY, zWorld, planetRpx);

  drawHUD(ctx);

  // Bottom panel
  ctx.fillStyle = '#05070c';
  ctx.fillRect(0, topH, W, botH);
  drawLightCurve(ctx);

  // Divider
  ctx.strokeStyle = 'rgba(120,140,180,0.3)';
  ctx.lineWidth = 1;
  ctx.beginPath();
  ctx.moveTo(0, topH + 0.5);
  ctx.lineTo(W, topH + 0.5);
  ctx.stroke();
}