// === CORONAL MASS EJECTION ===
// Physically-grounded CME simulation:
// Parker spiral field · Magnetic flux rope · Filament eruption
// Shock front · Solar wind · SDO/SOHO false-color palette

const scale     = addControl("scale",    "Scale",           8,  60,  28);
const eruptSpeed= addControl("espeed",   "Eruption Speed",  0,  3.0, 1.0);
const fieldTwist= addControl("twist",    "Flux Rope Twist", 0,  5.0, 2.2);
const windSpeed = addControl("wind",     "Solar Wind",      0,  3.0, 1.1);
const flareInt  = addControl("flare",    "Flare Intensity", 0,  2.0, 1.0);

// ── Physical constants (normalized) ─────────────────────
const PHI = 1.6180339887;
const SR2 = 1.4142135623;
const SR3 = 1.7320508075;
const SR5 = 2.2360679774;
const SR7 = 2.6457513110;

const T = time;
const es = eruptSpeed;

// ── Eruption lifecycle ───────────────────────────────────
// Slow build → filament suspension → catastrophic release → propagation
// Driven by irrational period so it never identically repeats
const cycleT    = T * es * 0.08 * PHI;
const cycleMod  = (cycleT % (2*Math.PI));

// Phase weights
const buildPhase   = Math.max(0, Math.sin(cycleMod * 0.5));            // slow energy build
const filamentPhase= Math.max(0, Math.sin(cycleMod * 0.5 - 0.8));     // filament rises
const eruptPhase   = Math.max(0, Math.sin(cycleMod * 0.5 - 1.4));     // CME launch
const propagPhase  = Math.max(0, Math.sin(cycleMod * 0.5 - 2.1));     // outward travel
const decayPhase   = Math.max(0, Math.sin(cycleMod * 0.5 - 2.8));     // dissipation

// Eruption front distance from sun (grows with time during eruption)
const eruptFront = eruptPhase * 2.2 + propagPhase * 3.8 + decayPhase * 1.2;

// ── Halton seed ──────────────────────────────────────────
const idx = i + 1;
let hx=0,hb2=0.5,hi2=idx;
while(hi2>0){hx+=(hi2&1)*hb2;hi2>>=1;hb2*=0.5;}
let hy=0,hf3=1/3,hi3=idx;
while(hi3>0){const r=hi3%3;hy+=r*hf3;hi3=Math.floor(hi3/3);hf3/=3;}
let hz=0,hf5=1/5,hi5=idx;
while(hi5>0){const r=hi5%5;hz+=r*hf5;hi5=Math.floor(hi5/5);hf5/=5;}

const n  = i / count;
const sx = hx*2-1, sy = hy*2-1, sz = hz*2-1;

// ── Particle layer assignment ────────────────────────────
// 5 populations: solar surface · corona · filament · CME rope · solar wind
// Each occupies a band of index space, softly blended

const surfaceW  = Math.max(0, 1.0 - n*8.0);                          // ~0-12%
const coronaW   = Math.exp(-Math.pow((n-0.18)*9.0, 2));              // ~10-26%
const filamentW = Math.exp(-Math.pow((n-0.38)*7.0, 2)) * filamentPhase;
const ropeW     = Math.exp(-Math.pow((n-0.60)*5.0, 2)) * eruptPhase;
const windW     = Math.max(0, n-0.55) * 1.8;                         // ~55-100%

const layerSum  = surfaceW+coronaW+filamentW+ropeW+windW+0.001;
const wSurf  = surfaceW /layerSum;
const wCor   = coronaW  /layerSum;
const wFil   = filamentW/layerSum;
const wRope  = ropeW    /layerSum;
const wWind  = windW    /layerSum;

// ── LAYER 0: Solar surface (photosphere/chromosphere) ────
// Granulation cells + spicules rising from surface
const sTheta = hx * Math.PI * 2.0;
const sPhi2  = (hy - 0.5) * Math.PI;
const sR     = 1.0 + 0.04*Math.sin(sTheta*7.0 + T*es*0.4*SR2)
                   + 0.02*Math.sin(sTheta*13.0 - T*es*0.6*PHI)
                   + 0.03*buildPhase*Math.sin(sTheta*3.0+sPhi2*4.0+T*es*0.3); // pre-eruption bulge
const sX = sR*Math.cos(sPhi2)*Math.cos(sTheta);
const sY = sR*Math.sin(sPhi2);
const sZ = sR*Math.cos(sPhi2)*Math.sin(sTheta);

// ── LAYER 1: Corona ──────────────────────────────────────
// Follows coronal streamer belt structure
// Streamers concentrate near solar equator
const cTheta = hx * Math.PI * 2.0;
const cPhi   = (hy - 0.5) * Math.PI * 0.7; // equatorial bias
const cR     = 1.3 + hz*0.9
             + 0.15*Math.sin(cTheta*2.0 + T*es*0.15*SR3)  // helmet streamers
             + 0.08*buildPhase*Math.sin(cTheta*1.0+T*es*0.1); // energy buildup
const cX = cR*Math.cos(cPhi)*Math.cos(cTheta);
const cY = cR*Math.sin(cPhi)*0.6; // flattened toward equator
const cZ = cR*Math.cos(cPhi)*Math.sin(cTheta);

// ── LAYER 2: Filament ────────────────────────────────────
// Cool dense plasma suspended in magnetic dip
// Rises as eruption approaches, then launches
const fAngle = hx * Math.PI * 2.0;
const fHeight= 1.15 + filamentPhase * 0.6 * (1.0 + 0.3*Math.sin(T*es*SR5*0.2));
const fWidth = 0.35 + 0.1*Math.sin(T*es*0.12*PHI);
// Filament has a twisted spine — flux rope embryo
const fSpine = fAngle * fieldTwist * 0.4;
const fX = (fHeight + hz*fWidth*0.3)*Math.cos(fAngle)
          + hy*fWidth*Math.sin(fSpine)*0.4;
const fY = hy*fWidth*Math.cos(fSpine) + filamentPhase*0.2*Math.sin(fAngle*2.0+T*es*0.3);
const fZ = (fHeight + hz*fWidth*0.3)*Math.sin(fAngle)
          - hy*fWidth*Math.sin(fSpine)*0.4;

// ── LAYER 3: CME Flux Rope ───────────────────────────────
// The core of the CME — twisted magnetic tube propagating outward
// Parker spiral bends trajectory as it travels
const rAngle  = hx * Math.PI * 2.0;  // angle around rope axis
const rAlong  = sy;                   // position along rope length (-1..1)
const rTube   = 0.28 + 0.08*Math.sin(rAngle*3.0 + T*es*SR2*0.3); // tube radius

// Rope center follows eruption front outward
const ropeDist= 1.4 + eruptFront * 0.8;
// Parker spiral: azimuthal lag increases with distance
// Ω_sun * r / v_wind — we normalize this
const parkerLag = ropeDist * 0.35 * windSpeed;
const ropeAzimuth = T*es*0.05 + parkerLag; // base azimuth + parker wind

// Rope axis direction (tilted slightly from ecliptic)
const axisX = Math.cos(ropeAzimuth);
const axisZ = Math.sin(ropeAzimuth);
const axisY = 0.18; // slight north-south tilt

// Flux rope: helical winding around propagating axis
const helixPhase = rAngle + rAlong*fieldTwist*PHI + T*es*SR3*0.15;
const perpX = -Math.sin(ropeAzimuth);
const perpZ =  Math.cos(ropeAzimuth);

const rX = axisX*ropeDist
          + axisX*rAlong*0.9
          + perpX*rTube*Math.cos(helixPhase)
          + 0.0*Math.sin(helixPhase); // keep in XZ plane mostly
const rY = axisY*ropeDist
          + axisY*rAlong*0.9
          + rTube*Math.sin(helixPhase)*0.8;
const rZ = axisZ*ropeDist
          + axisZ*rAlong*0.9
          + perpZ*rTube*Math.cos(helixPhase);

// ── LAYER 4: Solar wind + shock front ───────────────────
// Ambient solar wind follows Parker spiral outward
// Shock front piles up just ahead of CME
const wTheta   = hx * Math.PI * 2.0;
const wR_base  = 1.0 + hb2 * 4.5;  // radial distance 1..5.5
// Parker spiral angle: φ = Ω*r/v
const parkerAngle = wTheta + wR_base * 0.28 * windSpeed
                  + T*es*0.06*SR2;
// Shock compression: density piles up at eruption front distance
const shockDist  = 1.35 + eruptFront * 0.82;
const shockProx  = Math.exp(-Math.pow((wR_base - shockDist)*2.5, 2));
const wR         = wR_base + shockProx * 0.3 * eruptPhase; // pile-up

const wX = wR*Math.cos(hz*0.6)*Math.cos(parkerAngle);
const wY = wR*Math.sin(hz*0.5)*0.4;
const wZ = wR*Math.cos(hz*0.6)*Math.sin(parkerAngle);

// ── BLEND ALL LAYERS ────────────────────────────────────
let fx = wSurf*sX + wCor*cX + wFil*fX + wRope*rX + wWind*wX;
let fy = wSurf*sY + wCor*cY + wFil*fY + wRope*rY + wWind*wY;
let fz = wSurf*sZ + wCor*cZ + wFil*fZ + wRope*rZ + wWind*wZ;

// ── Magnetic turbulence ──────────────────────────────────
// Small-scale field perturbations — Alfvén waves
const turbScale = 0.04 * (1.0 + eruptPhase*2.0);
const tx2 = Math.sin(fy*4.1+T*es*SR5*0.4)*Math.cos(fz*3.7-T*es*SR2*0.3);
const ty2 = Math.sin(fz*3.9-T*es*PHI*0.35)*Math.cos(fx*4.3+T*es*SR3*0.25);
const tz2 = Math.sin(fx*4.5+T*es*SR7*0.3)*Math.cos(fy*3.3-T*es*PHI*0.28);
fx += tx2*turbScale; fy += ty2*turbScale; fz += tz2*turbScale;

// ── Breath / pulsation (solar oscillations — p-modes) ───
const fDist = Math.sqrt(fx*fx+fy*fy+fz*fz)+0.001;
const pMode = 1.0
  + 0.025*Math.sin(fDist*5.5 - T*es*SR2*3.1)
  + 0.015*Math.sin(fDist*9.3 - T*es*PHI*2.4)
  + 0.010*Math.sin(fDist*14.7- T*es*SR3*1.9);

target.set(fx*scale*pMode, fy*scale*pMode, fz*scale*pMode);

// ══ SDO / SOHO FALSE-COLOR PALETTE ═══════════════════════
// Each AIA channel highlights different temperature plasma:
// 304Å  = chromosphere  = deep red/brick    (~50,000K)
// 171Å  = corona        = golden yellow     (~600,000K)
// 193Å  = hot corona    = electric teal     (~1.2MK)
// 131Å  = flare plasma  = cyan-white        (~10MK)
// Shock front = blue-white compression

const surfAngle2 = Math.atan2(fy, fx);
const elevAngle  = Math.atan2(fz, Math.sqrt(fx*fx+fy*fy));

// Temperature proxy: increases with layer (surface→corona→rope→shock)
const tempProxy = wSurf*0.05 + wCor*0.25 + wFil*0.15
                + wRope*(0.5+eruptPhase*0.4) + wWind*(0.3+shockProx*0.6);

// SDO channel hue mapping
// 0.0  = deep red (304Å chromosphere)
// 0.11 = orange-gold (171Å quiet corona)
// 0.48 = electric teal (193Å active corona)
// 0.55 = cyan (131Å flare)
// 0.62 = blue-white (shock)
const hue304  = 0.02;
const hue171  = 0.10;
const hue193  = 0.48;
const hue131  = 0.54;
const hueShock= 0.61;

// Blend hues by temperature
const h1 = hue304  + (hue171-hue304) *Math.min(1,tempProxy*4.0);
const h2 = hue171  + (hue193-hue171) *Math.min(1,Math.max(0,tempProxy*3.0-0.5));
const h3 = hue193  + (hue131-hue193) *Math.min(1,Math.max(0,tempProxy*2.5-1.0));
const h4 = hue131  + (hueShock-hue131)*Math.min(1,Math.max(0,shockProx*eruptPhase*3.0));
const hue = ((h1+h2+h3+h4)*0.25
  + 0.03*Math.sin(surfAngle2*2.0+T*es*SR3*0.2)
  + 0.02*Math.sin(elevAngle*3.0+T*es*PHI*0.15)) % 1.0;

// Saturation: flare core and shock front most vivid
const boundary3 = Math.abs(Math.sin(fDist*4.8+T*es*SR2*0.25));
const flareSat  = flareInt*(wRope*eruptPhase + shockProx*0.6);
const sat = 0.55 + 0.35*boundary3 + 0.3*flareSat;

// Luminance: surface granulation dark, corona glows, flare blindingly bright
const coronaGlow  = wCor*0.4*(1.0+buildPhase*0.5);
const flareBloom  = flareInt*eruptPhase*(wRope*0.7+shockProx*0.5);
const filamentDark= wFil*0.3*(1.0-filamentPhase*0.4); // filaments are cool/dark
const lum = 0.06
  + wSurf*0.35*(0.8+0.2*Math.sin(surfAngle2*8.0+T*es*0.5)) // granulation
  + coronaGlow
  + flareBloom*0.8
  + filamentDark
  + 0.15*Math.abs(Math.sin(fDist*6.1+T*es*SR5*0.3))*boundary3;

color.setHSL((hue+1)%1.0, Math.min(1.0,sat), Math.min(0.98,lum));

// ── HUD ─────────────────────────────────────────────────
if(i===0){
  const phaseLabel =
    eruptPhase > 0.5  ? "⚡ ERUPTION" :
    filamentPhase>0.4 ? "↑ FILAMENT RISING" :
    propagPhase >0.4  ? "→ PROPAGATING" :
    decayPhase  >0.4  ? "~ DISSIPATING" :
                        "● ENERGY BUILD";
  setInfo(
    "Coronal Mass Ejection  —  " + phaseLabel,
    "Physically-grounded CME: Parker spiral solar wind · Magnetic flux rope (helical field lines) · " +
    "Filament eruption · Bow shock compression. " +
    "False-color SDO/AIA palette: red=chromosphere · gold=quiet corona · teal=active corona · cyan=flare · blue=shock. " +
    "Eruption Speed · Flux Rope Twist · Solar Wind · Flare Intensity"
  );
  annotate("sun",     new THREE.Vector3(0,0,0),                             "☀ Photosphere");
  annotate("corona",  new THREE.Vector3(scale*0.25, scale*0.25, 0),         "◌ Corona");
  annotate("rope",    new THREE.Vector3(scale*0.5*Math.cos(T*es*0.05), 0, scale*0.5*Math.sin(T*es*0.05)), "⟳ Flux Rope");
  annotate("shock",   new THREE.Vector3(-scale*0.6, scale*0.1, scale*0.3),  "▷ Shock Front");
  annotate("wind",    new THREE.Vector3(scale*0.7, -scale*0.4, -scale*0.3), "~ Parker Spiral");
}