// Modified from calcrxc.hoc and stim.hoc of Ted Carnevale
// $Id: stim.hoc,v 1.5 2009/02/24 00:55:27 ted Exp ted $
//
// Calculates the transfer resistances between extracellular stimulating|recording
// electrode(s) and a model neuron.  Relies on the principle of reciprocity, which
// assumes that the intervening bath and tissue can be treated as linear.  Suppose a
// stimulus current of amplitude Is, applied to a particular configuration of
// extracellular electrode(s), produces a potential Vext(x,y,z) at location (x,y,z).
// Then the transfer resistance between the electrode(s) and (x,y,z) is 
//   rx(x,y,z) = Vext(x,y,z)/Is
//
// Insert the extracellular and xtra mechanisms in all sections that are subject to 
// the extracellular field. Compute the transfer resistance rx for every section that
// contains xtra, as illustrated below. Construct a stimulus waveform template and 
// copy it to a Vector. For each internal node along the axon, use this Vector to 
// drive is_xtra(x). The xtra mechanism uses the rx values to convert the stimulus
// current waveform into the proper amplitude and sign of the local extracellular 
// field.
// 
// If rho, b, or c is changed, new_elec() must be invoked to make the changes take
// effect.
//
// To a first approximation, a monopolar stimulating electrode that delivers current 
// I produces a field in which potential V is given by 
//   V = I rho / 4 PI r
// where r is the distance from the center of the electrode.
// 
// The principle of superposition may be applied to deal with an arbitrary number of
// monopolar electrodes, or even surface electrodes with different shapes and areas,
// which are located at arbitrary positions, and deliver arbitrary stimulus currents.
//
// The stimulus is constructed from a waveform template that is copied to a Vector. 
// For each section that has the xtra mechanism, this Vector is used to drive
// is_xtra. The transfer resistance rx_xtra takes care of the amplitude and sign of 
// the local extracellular field. This works with fixed dt and adaptive integration.


// resistivity of medium (ohm cm)
rho = 110  // 60 in Greenberg99; 35.4 for sea water

// monopolar electrode position (um)
stimX = 295
stimY = 300    //295
stimZ = 35  // >= 25 for electrode radius + soma radius

// rectangular pulse 
stimDel = 1200    // ms
stimDur = 0.5     // ms
stimAmp = -0.020  // mA

objref stim_amp, stim_time
stim_amp = new Vector()
stim_time = new Vector()


/////////////////////////////////////////////////////////////////////

proc setrx() {
    // xyc coords as arguments
    forall {
        if (ismembrane("xtra")) {
            // avoid nodes at 0 and 1 ends: not override values at internal nodes
            for (x,0) {
                r = sqrt( \
                    (x_xtra(x) - xe)^2 + (y_xtra(x) - ye)^2 + (z_xtra(x) - ze)^2 )
                rx_xtra(x) = (rho / 4 / PI) * (1/r) * 0.01
            }
        }
    }
}

create sElec  // bogus section to show extracell stim/rec electrode location
objref pElec  // bogus PointProcess just to show stim location
objref gElec  // will be a Shape that shows extracellular electrode location
gElec = new Shape(0)  // create it but don't map it to the screen yet
gElec.view(-100, -100, 900, 900, 200, 450, 201.6, 201.28)

proc drawelec() {
    sElec {
        pt3dclear()
        pt3dadd($1-0.5, $2, $3, 1)
        pt3dadd($1+0.5, $2, $3, 1)
        pElec = new IClamp(0.5)
    }
    gElec.point_mark(pElec, 2)  // red
}

proc setelec() {
    xe = $1
    ye = $2
    ze = $3
    setrx(xe, ye, ze)
    drawelec(xe, ye, ze)
}

proc stim_waveform() {
    // this uses interpolated play
    // index    0  1    2    3        4        5
    // stim vec 0, 0,   1,   1,       0        0
    // time vec 0, stimDel, stimDel, stimDel+stimDur, stimDel+stimDur, stimDel+stimDur+1
    //  really  0, $1,  $1,  $1+$2,   $1+$2,   $1+$2+1
    // first the stim vector
    stim_amp.resize(6)
    stim_amp.fill(0)
    stim_amp.x[2]=1
    stim_amp.x[3]=1
    stim_amp.mul($3)
    // now the time vector
    stim_time.resize(6)
    stim_time.x[1]=$1
    stim_time.x[2]=$1
    stim_time.x[3]=$1+$2
    stim_time.x[4]=$1+$2
    stim_time.x[5]=$1+$2+1
}

ATTACHED__ = 0
proc attach_stim() {
    // since is_xtra is GLOBAL, we only need to specify Vector.play()
    // for one instance of xtra, i.e. at just one internal node
    // of only one section that contains xtra
    forall {
        if (ATTACHED__ == 0) {
            if (ismembrane("xtra")) {
                stim_amp.play(&is_xtra, stim_time, 1) // "interpolated" play
                ATTACHED__ = 1
            }
        }
    }
}

proc setstim() {
    del = $1
    dur = $2
    amp = $3
    stim_waveform(del, dur, amp)
    attach_stim()
}


/////////////////////////////////////////////////////////////////////

setelec(stimX, stimY, stimZ)  // put stim electrode at (x, y, z)

xpanel("Extracellular Electrode Location", 0)
    xlabel("xyz coords in um")
    xvalue("x", "stimX", 1, "setelec(stimX,stimY,stimZ)", 0, 1)
    xvalue("y", "stimY", 1, "setelec(stimX,stimY,stimZ)", 0, 1)
    xvalue("z", "stimZ", 1, "setelec(stimX,stimY,stimZ)", 0, 1)
xpanel(0,450)

setstim(stimDel, stimDur, stimAmp)

xpanel("Extracellular Stimulus Current", 0)
    xvalue("del (ms)", "stimDel", 1, "setstim(stimDel,stimDur,stimAmp)", 0, 1)
    xvalue("dur (ms)", "stimDur", 1, "setstim(stimDel,stimDur,stimAmp)", 0, 1)
    xvalue("amp (mA)", "stimAmp", 1, "setstim(stimDel,stimDur,stimAmp)", 0, 1)
xpanel(0,580)