/*
============================================================================
Cavitational Capacitive Drive (CCD) Process Handlers (mW/cm2)
Author: Mithun Padmakumar
Date: July 2026

Citation:
  Padmakumar, M., Rajan, D., & Steephen, J. E. (2026). Cavitational capacitive 
  drive: A computationally efficient model for ultrasonic neuromodulation. 
  Journal of Neural Engineering.

Description:
  Contains procedure definitions for inserting/updating the CCD mechanism, 
  managing simulation time steps, and pre-processing/saving simulation output 
  vectors (raw Vm, window-averaged Vm with metadata, full Cm waveform, AP times).
============================================================================
*/

// Procedures for CCD Mechanism Control

proc InsertCCD() {
    if ($1) {
        forall {
            insert ccd
            tbegin_ccd = tbegin
            tdur_ccd = tdurn
            usf_ccd = freq
            usi_ccd = usi
            for (x, 0) {
                setpointer c_ccd(x), cm(x)
            }
        }
        dt = 0.025 / freq
        bFUSS = 1
        printf("Activated CCD Focused Ultrasound Stimulation (FUSS).\n")
    } else {
        forall {
            uninsert ccd
        }
        dt = 0.025
        bFUSS = 0
        printf("Deactivated CCD FUSS (Baseline RS Neuron).\n")
    }
}

proc update_ccd() {
    if (bFUSS) {
        forall {
            tbegin_ccd = tbegin
            tdur_ccd = tdurn
            usf_ccd = freq
            usi_ccd = usi
        }
        dt = 0.025 / freq
        printf("Updated CCD parameters: Freq=%d kHz, Intensity=%g mW/cm2, tbegin=%g ms, tdurn=%g ms.\n", freq, usi, tbegin, tdurn)
    }
}

proc toggle_fuss() {
    InsertCCD($1)
}


// Procedures for Output Management

objref outFile
objref processed_v, processed_t, g_proc

proc saveRawVm() {
    strdef fname
    if (bFUSS) {
        sprint(fname, "Results/ccd_%g_%d.dat", usi, freq)
    } else {
        sprint(fname, "Results/baseline.dat")
    }

    outFile = new File()
    outFile.wopen(fname)
    for (i = 0; i < ts.size(); i += 1) {
        outFile.printf("%g\t%g\n", ts.x(i), vsoma.x(i))
    }
    outFile.close()

    print "Raw simulation result stored in ", fname
}

proc saveProcessedVm() {
    strdef fname, metafile
    processed_v = new Vector()
    processed_t = new Vector()

    if (bFUSS) {
        w = 40
        half_w = 20
        step_size = freq
        if (step_size < 1) step_size = 1

        // Perform 40-point centered window averaging and downsampling by step_size (freq)
        for (i = 0; i < ts.size(); i += step_size) {
            sum_v = 0
            count_v = 0
            start_idx = i - half_w
            if (start_idx < 0) start_idx = 0
            end_idx = i + half_w - 1
            if (end_idx >= vsoma.size()) end_idx = vsoma.size() - 1

            for (k = start_idx; k <= end_idx; k += 1) {
                sum_v += vsoma.x(k)
                count_v += 1
            }
            if (count_v > 0) {
                processed_v.append(sum_v / count_v)
                processed_t.append(ts.x(i))
            }
        }

        sprint(fname, "Results/ccd_processed_%g_%d.dat", usi, freq)
        sprint(metafile, "Results/ccd_processed_%g_%d_metadata.dat", usi, freq)
    } else {
        // When CCD is deactivated, dt = 0.025 ms already; raw and processed waveforms are identical.
        for (i = 0; i < ts.size(); i += 1) {
            processed_v.append(vsoma.x(i))
            processed_t.append(ts.x(i))
        }

        sprint(fname, "Results/baseline_processed.dat")
        sprint(metafile, "Results/baseline_processed_metadata.dat")
    }

    // Display processed trace in a new graph window
    g_proc = new Graph()
    g_proc.size(0, tstop, -80, 40)
    processed_v.plot(g_proc, processed_t, 2, 1)
    if (bFUSS) {
        g_proc.label(0.1, 0.9, "Processed Vm (Window Avg & Downsampled)")
    } else {
        g_proc.label(0.1, 0.9, "Processed Vm (Raw Vm at dt=0.025 ms)")
    }

    // Save processed waveform file
    outFile = new File()
    outFile.wopen(fname)
    for (i = 0; i < processed_t.size(); i += 1) {
        outFile.printf("%g\t%g\n", processed_t.x(i), processed_v.x(i))
    }
    outFile.close()

    // Save metadata file (stored ONLY when processed waveform is saved)
    outFile = new File()
    outFile.wopen(metafile)
    if (bFUSS) {
        outFile.printf("FUSS=CCD\n")
        outFile.printf("Frequency(kHz)=%d\n", freq)
        outFile.printf("Intensity(mW/cm2)=%g\n", usi)
        outFile.printf("Start(ms)=%g\n", tbegin)
        outFile.printf("Duration(ms)=%g\n", tdurn)
        outFile.printf("TSTOP(ms)=%g\n", tstop)
        outFile.printf("WindowAveraging=40-point centered moving average\n")
        outFile.printf("DownsamplingFactor=%d\n", freq)
        outFile.printf("EffectiveDt(ms)=0.025\n")
    } else {
        outFile.printf("FUSS=none\n")
        outFile.printf("TSTOP(ms)=%g\n", tstop)
        outFile.printf("WindowAveraging=none (baseline simulation at dt=0.025 ms)\n")
        outFile.printf("EffectiveDt(ms)=0.025\n")
    }
    outFile.close()

    print "Processed Vm stored in ", fname
    print "Metadata stored in ", metafile
}

proc saveCm() {
    strdef fname
    if (bFUSS) {
        sprint(fname, "Results/ccd_Cm_%g_%d.dat", usi, freq)
    } else {
        sprint(fname, "Results/baseline_Cm.dat")
    }

    outFile = new File()
    outFile.wopen(fname)
    for (i = 0; i < ts.size(); i += 1) {
        outFile.printf("%g\t%g\n", ts.x(i), cm_vec.x(i))
    }
    outFile.close()

    print "Full Cm vector stored in ", fname
}

proc saveAPC() {
    strdef fname
    if (bFUSS) {
        sprint(fname, "Results/APtimes_ccd_%g_%d.dat", usi, freq)
    } else {
        sprint(fname, "Results/APtimes_baseline.dat")
    }

    outFile = new File()
    outFile.wopen(fname)
    for (i = 0; i < ap_times.size(); i += 1) {
        outFile.printf("%g\n", ap_times.x(i))
    }
    outFile.close()

    print "AP times stored in ", fname
}