MEMORY PALACE

/*
▖  ▖             ▄▖  ▜
▛▖▞▌█▌▛▛▌▛▌▛▘▌▌  ▙▌▀▌▐ ▀▌▛▘█▌
▌▝ ▌▙▖▌▌▌▙▌▌ ▙▌  ▌ █▌▐▖█▌▙▖▙▖
             ▄▌
multi-tap delay with chorus/reverb capabilities
by @ersatzben

- channels parameter = number of taps
- first tap occurs at time_ms, last tap at time_ms * span
- curve parameter controls (non)linearity of how taps are distributed across that period
- curve: 1 = linear, >1 bunched at start (early reflections), <1 bunched at end (reverse effect)
- all other controls as described below

SETUP
- two inputs
- two outputs

*/


History counter(0);

// ---------- TAP PARAMETERS --------------
Param channels(42, min=1, max=50);			// number of taps
Param time_ms(50, min=0.01, max=250);		// first tap time (ms)
Param span(7, min=2, max=10);				// last tap time (as multiple of first)
Param curve(3.2, min=0.1, max=5);			// distribution of taps
Param diffusion(0.25, min=0, max=1);       	// blur (jitter + allpass)
Param modRate(0.0001, min=0, max=10);		// modulation rate
Param modDepth(0.002, min=0, max=0.02);		// modulation depth

// ---------- FEEDBACK PARAMETERS ---------
Param fb(0.8, min=0, max=0.99999);			// feedback
Param damp(0.05, min=0, max=0.99);			// LPF damping
Param ap_coeff(0.6, min=-0.7, max=0.7);		// allpass coefficient
Param hp_hz(25, min=5, max=200);			// HPF bass control
Param width(1.0, min=0, max=1);				// stereo spread of taps
Param freeze(0.0, min=0, max=1);			// freeze (0 = off, 1 = on)

// ---------- MIXING PARAMETERS -----------
Param wet(1.0, min=0, max=1);				// dry/wet mix
Param inGain(4.0, min=0, max=8);			// input gain
Param drive(0.02, min=0, max=3);			// output soft saturation / drive
Param xmix(1.0, min=0, max=1);            	// L/R stereo crosstalk
Param xangle(0.25, min=0, max=1);         	// crosstalk rotation angle
Param xpasses(1, min=1, max=2);				// crosstalk passes (once/twice)

// ---------- RE-TIME XFADE ----------
Param xfade_ms(30, min=0, max=100);

// ---------- data ----------
Data buf(samplerate*5, 50);
Data dtime(50,1);
Data lp_z(50,1);
Data ap1_x_z(50,1);                       // AP1 x[n-1]
Data ap1_y_z(50,1);                       // AP1 y[n-1]
Data ap2_x_z(50,1);                       // AP2 x[n-1]
Data ap2_y_z(50,1);                       // AP2 y[n-1]
Data dry1_z(50,1);                        // 1-sample dry delay for orthonormal mix
Data seed(50,1);
Data tmp_fb(50,1);
Data mix_fb(50,1);
Data dc_lp(50,1);
Data prev_dt(50,1);
Data xf_p(50,1);
Data lfo_p(50,1);

// ---------- init ----------
if (counter == 0) {
    for (i=0; i<50; i+=1) {
        t = sin(i * 127.1) * 43758.5453;
        j = t - floor(t);
        poke(seed, j, i);
        poke(lp_z, 0, i);
        poke(ap1_x_z, 0, i);
        poke(ap1_y_z, 0, i);
        poke(ap2_x_z, 0, i);
        poke(ap2_y_z, 0, i);
        poke(dry1_z, 0, i);
        poke(tmp_fb, 0, i);
        poke(mix_fb, 0, i);
        poke(dc_lp, 0, i);
        poke(prev_dt, 1, i);
        poke(xf_p, 1, i);
        poke(lfo_p, fract(j + 0.123*i), i);
    }
}

// ---------- param shaping ----------
base    = slide(mstosamps(time_ms), samplerate/20, samplerate/20);
maxd    = base * slide(span, samplerate/8, samplerate/8);
ch      = max(1, floor(channels));
buflen  = dim(buf, 0);

fb_smooth        = slide(fb, samplerate/10, samplerate/10);
damp_smooth      = slide(damp, samplerate/10, samplerate/10);
wet_smooth       = slide(wet, samplerate/20, samplerate/20);
width_smooth     = slide(width, samplerate/20, samplerate/20);
inGain_smooth    = slide(inGain, samplerate/10, samplerate/10);
drive_smooth     = slide(drive, samplerate/15, samplerate/15);
modRate_smooth   = slide(modRate, samplerate/5, samplerate/5);
modDepth_smooth  = slide(modDepth, samplerate/10, samplerate/10);
freeze_smooth    = slide(freeze, samplerate/50, samplerate/50);
xmix_smooth      = slide(xmix, samplerate/15, samplerate/15);
xangle_smooth    = slide(xangle, samplerate/15, samplerate/15);
curve_smooth     = slide(curve, samplerate/8, samplerate/8);
diffusion_smooth = slide(diffusion, samplerate/12, samplerate/12);

// AP coeff: tighter clamp + smoothing
ap_coeff_smooth  = slide(clamp(ap_coeff, -0.75, 0.75), samplerate/30, samplerate/30);

// freeze as continuous mix
freeze_mix = clamp(freeze_smooth, 0, 1);
in_atten   = 1 - freeze_mix;
freeze_fb  = mix(fb_smooth, 0.9995, freeze_mix);

// ---------- delay table with diffusion-driven jitter (stronger) ----------
for (i=0; i<ch; i+=1) {
    t01    = (ch==1) ? 0 : (i / (ch-1));
    shaped = mix(base, maxd, pow(t01, curve_smooth));
    j      = peek(seed, i);
    jofs   = (j - 0.5) * (0.2 * diffusion_smooth) * shaped;   // was 0.1
    dt     = clamp(shaped + jofs, 1, buflen-2);
    poke(dtime, dt, i);
}

// pointer
w = wrap(counter, 0, buflen);

// ---------- read + per-channel filtering ----------
accL = 0; accR = 0;

hp_fc = max(0, hp_hz) * (1 - freeze_mix);
hp_a  = (hp_fc <= 0) ? 0 : (1 - exp(-twopi * hp_fc / samplerate));

xf_inc = (xfade_ms <= 0) ? 1 : (1000 / (xfade_ms * samplerate));

xinL = inGain_smooth * in1 * in_atten;
xinR = inGain_smooth * in2 * in_atten;

for (i=0; i<ch; i+=1) {
    dt_now  = peek(dtime, i);
    dt_prev = peek(prev_dt, i);

    ph  = peek(lfo_p, i) + modRate_smooth / samplerate;
    ph  = ph - floor(ph);
    poke(lfo_p, ph, i);
    lfo = sin(twopi * (peek(seed,i) + ph));
    md  = modDepth_smooth * (1 - freeze_mix);

    xf   = peek(xf_p, i);
    dt_thresh = 8;
    trig = (abs(dt_now - dt_prev) > dt_thresh) && (xf >= 0.999);
    xf   = trig ? 0 : min(1, xf + xf_inc);
    poke(xf_p, xf, i);

    dmod_old = dt_prev * md * lfo;
    dmod_new = dt_now  * md * lfo;

    rpos_old = wrap(w - (dt_prev + dmod_old), 0, buflen);
    rpos_new = wrap(w - (dt_now  + dmod_new), 0, buflen);

    y_old = peek(buf, rpos_old, i, 1);
    y_new = peek(buf, rpos_new, i, 1);

    y = mix(y_old, y_new, xf);
    y = (abs(y) < 1e-20) ? 0 : y;

    if (xf >= 0.999) { poke(prev_dt, dt_now, i); }

    // damping with time tilt
    tsec   = dt_now / samplerate;
    tf     = clamp(1.0 / (1.0 + 2.0 * tsec), 0.1, 1.0);
    a_base = 1 - (1 - damp_smooth) * tf;
    a      = a_base * (1 - freeze_mix);

    z0   = peek(lp_z, i);
    y_lp = (1 - a) * y + a * z0;
    poke(lp_z, y_lp, i);

    // -------- STRONGER ALLPASS DIFFUSION (two cascaded stages) --------
    ap_depth = diffusion_smooth * (1 - freeze_mix);          // 0..1

    // nonlinear depth mapping: hits harder near 1
    depth_sh  = sqrt(ap_depth);                               // 0..1, faster rise
    g_base    = ap_coeff_smooth * (0.3 + 0.7 * depth_sh);     // scale with depth
    g1        =  g_base;                                      // stage 1
    g2        = -0.9 * g_base;                                // stage 2, opposite sign

    // tiny contraction per stage
    ap_loss1  = 1 - (0.0015 * ap_depth);
    ap_loss2  = 1 - (0.0015 * ap_depth);

    // Stage 1: y1 = -g1*x + x[n-1] + g1*y[n-1]
    x1_n   = y_lp;
    x1_nm1 = peek(ap1_x_z, i);  x1_nm1 = (abs(x1_nm1) < 1e-20) ? 0 : x1_nm1;
    y1_nm1 = peek(ap1_y_z, i);  y1_nm1 = (abs(y1_nm1) < 1e-20) ? 0 : y1_nm1;
    y1_raw = ap_loss1 * (-g1 * x1_n + x1_nm1 + g1 * y1_nm1);
    poke(ap1_x_z, x1_n, i);
    poke(ap1_y_z, y1_raw, i);

    // Stage 2: y2 = -g2*y1 + y1[n-1] + g2*y2[n-1]
    x2_n   = y1_raw;
    x2_nm1 = peek(ap2_x_z, i);  x2_nm1 = (abs(x2_nm1) < 1e-20) ? 0 : x2_nm1;
    y2_nm1 = peek(ap2_y_z, i);  y2_nm1 = (abs(y2_nm1) < 1e-20) ? 0 : y2_nm1;
    y2_raw = ap_loss2 * (-g2 * x2_n + x2_nm1 + g2 * y2_nm1);
    poke(ap2_x_z, x2_n, i);
    poke(ap2_y_z, y2_raw, i);

    // dry leg phase-align then orthonormal mix with cascaded AP
    y_lp_d = peek(dry1_z, i);
    poke(dry1_z, y_lp, i);

    phi      = ap_depth * (pi * 0.5);                      // 0..90°
    c_ap     = cos(phi);
    s_ap     = sin(phi);
    y_mix_ap = c_ap * y_lp_d + s_ap * y2_raw;              // at 0 → dry; at 1 → AP cascade

    // DC blocker on feedback candidate
    dc0  = peek(dc_lp, i);
    dc1  = dc0 + hp_a * (y_mix_ap - dc0);
    y_fb = (hp_fc <= 0) ? y_mix_ap : (y_mix_ap - dc1);
    poke(dc_lp, dc1, i);

    // store for rotation mixing
    poke(tmp_fb, y_fb, i);

    // stereo sum via unity-gain M/S width law
    side = (i % 2);
    Lc = (side==0) ? y : 0;
    Rc = (side==1) ? y : 0;
    M  = 0.5 * (Lc + Rc);
    S  = 0.5 * (Lc - Rc);
    S *= width_smooth;
    accL += M + S;
    accR += M - S;
}

// ---------- pure rotation diffusion (gain-invariant) ----------
theta = (xmix_smooth * xangle_smooth) * (pi * 0.5);
c = cos(theta);
s = sin(theta);

// pass 1
for (i=0; i<ch; i+=2) {
    j  = (i + 1) % ch;
    u  = peek(tmp_fb, i);
    v  = peek(tmp_fb, j);
    u2 = c*u - s*v;
    v2 = s*u + c*v;
    poke(mix_fb, u2, i);
    poke(mix_fb, v2, j);
}

// pass 2
if (xpasses > 1) {
    for (i=1; i<ch; i+=2) {
        j  = (i + 1) % ch;
        u  = peek(mix_fb, i);
        v  = peek(mix_fb, j);
        u2 = c*u - s*v;
        v2 = s*u + c*v;
        poke(mix_fb, u2, i);
        poke(mix_fb, v2, j);
    }
}

// optional: slight fb trim with diffusion for stability
freeze_fb *= (1 - 0.0008 * diffusion_smooth);

// ---------- write ----------
for (i=0; i<ch; i+=1) {
    y_mix   = peek(mix_fb, i);
    x_in_ch = ((i % 2)==0) ? xinL : xinR;
    poke(buf, x_in_ch + freeze_fb * y_mix, w, i);
}

// ---------- outputs ----------
procL = accL / ch;
procR = accR / ch;

k = max(1e-6, drive_smooth);
procL_driven = tanh(k * procL) / tanh(k);
procR_driven = tanh(k * procR) / tanh(k);
lim = 0.99;
procL_driven = clamp(procL_driven, -lim, lim);
procR_driven = clamp(procR_driven, -lim, lim);

out1  = mix(in1, procL_driven, wet_smooth);
out2  = mix(in2, procR_driven, wet_smooth);

// ---------- advance ----------
counter += 1;

/*
Copyright (c) 2025 Ben Johnson

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and
associated documentation files (the "Software"), to deal in the Software without restriction,
including without limitation the rights to use, copy, modify, merge, publish, distribute,
sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is
furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or
substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT
NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT
OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
*/