1: //
   2: // Heccer : a compartmental solver that implements efficient Crank-Nicolson
   3: // integration for neuronal models.
   4: //
   5: 
   6: //////////////////////////////////////////////////////////////////////////////
   7: //'
   8: //' Heccer : testbed C implementation
   9: //'
  10: //' Copyright (C) 2006-2008 Hugo Cornelis
  11: //'
  12: //' functional ideas .. Hugo Cornelis, hugo.cornelis@gmail.com
  13: //'
  14: //' coding ............ Hugo Cornelis, hugo.cornelis@gmail.com
  15: //'
  16: //////////////////////////////////////////////////////////////////////////////
  17: 
  18: 
  19: #include <math.h>
  20: #include <stdio.h>
  21: #include <stdlib.h>
  22: #include <string.h>
  23: 
  24: 
  25: #include "heccer/vclamp.h"
  26: 
  27: 
  28: /// **************************************************************************
  29: ///
  30: /// SHORT: VClampAddVariable()
  31: ///
  32: /// ARGS.:
  33: ///
  34: ///     pvc.........: voltage clamper.
  35: ///     pvVoltage...: pointer to the voltage variable, assumed is double *
  36: ///     pvInjector..: pointer to the variable for injected current.
  37: ///
  38: /// RTN..: int
  39: ///
  40: ///     success of operation.
  41: ///
  42: /// DESCR: Clamp the given voltage variable, using the given injector.
  43: ///
  44: /// **************************************************************************
  45: 
  46: int
  47: VClampAddVariable
  48: (struct VClamp * pvc, void *pvVoltage, void *pvInjector)
  49: {
  50:     //- set default result: ok
  51: 
  52:     int iResult = 1;
  53: 
  54:     if (pvc->iClampsActive >= 1)
  55:     {
  56:         return(0);
  57:     }
  58: 
  59:     //- set next variable
  60: 
  61:     pvc->pdVoltage = (double *)pvVoltage;
  62: 
  63:     pvc->pdInjector = (double *)pvInjector;
  64: 
  65:     pvc->iClampsActive++;
  66: 
  67:     //- return result
  68: 
  69:     return(iResult);
  70: }
  71: 
  72: 
  73: /// **************************************************************************
  74: ///
  75: /// SHORT: VClampFinish()
  76: ///
  77: /// ARGS.:
  78: ///
  79: ///     pvc...: voltage clamper.
  80: ///
  81: /// RTN..: int
  82: ///
  83: ///     success of operation.
  84: ///
  85: /// DESCR: Free the voltage clamper.
  86: ///
  87: /// **************************************************************************
  88: 
  89: int VClampFinish(struct VClamp * pvc)
  90: {
  91:     //- set default result: ok
  92: 
  93:     int iResult = 1;
  94: 
  95:     //- free all allocated memory
  96: 
  97:     free(pvc);
  98: 
  99:     //- return result
 100: 
 101:     return(iResult);
 102: }
 103: 
 104: 
 105: /// **************************************************************************
 106: ///
 107: /// SHORT: VClampInitiate()
 108: ///
 109: /// ARGS.:
 110: ///
 111: ///     pvc...: voltage clamper.
 112: ///
 113: /// RTN..: int
 114: ///
 115: ///     success of operation.
 116: ///
 117: /// DESCR: Initiate the voltage clamper.
 118: ///
 119: /// **************************************************************************
 120: 
 121: int VClampInitiate(struct VClamp * pvc)
 122: {
 123:     //- set default result: ok
 124: 
 125:     int iResult = 1;
 126: 
 127:     //- just clear out some varialbes...
 128: 
 129:     //t should initialize from a stream as heccer does.
 130: 
 131:     pvc->dEIntegral = 0;
 132:     pvc->dEPrevious = 0;
 133:     pvc->dOutput = 0;
 134: 
 135:     //- return result
 136: 
 137:     return(iResult);
 138: }
 139: 
 140: 
 141: /// **************************************************************************
 142: ///
 143: /// SHORT: VClampNew()
 144: ///
 145: /// ARGS.:
 146: ///
 147: ///     pcName..: name of this object.
 148: ///
 149: /// RTN..: struct VClamp *
 150: ///
 151: ///     voltage clamper, NULL for failure.
 152: ///
 153: /// DESCR: voltage clamper.
 154: ///
 155: /// **************************************************************************
 156: 
 157: struct VClamp * VClampNew(char *pcName)
 158: {
 159:     //- set default result: failure
 160: 
 161:     struct VClamp * pvcResult = NULL;
 162: 
 163:     //- allocate voltage clamper
 164: 
 165:     pvcResult = (struct VClamp *)calloc(1, sizeof(struct VClamp));
 166: 
 167:     if (!pvcResult)
 168:     {
 169:         return(NULL);
 170:     }
 171: 
 172:     //- set name
 173: 
 174:     pvcResult->pcName = calloc(1 + strlen(pcName), sizeof(char));
 175: 
 176:     strcpy(pvcResult->pcName, pcName);
 177: 
 178:     //- return result
 179: 
 180:     return(pvcResult);
 181: }
 182: 
 183: 
 184: /// **************************************************************************
 185: ///
 186: /// SHORT: VClampSetFields()
 187: ///
 188: /// ARGS.:
 189: ///
 190: ///     dInjected.....: injected current.
 191: ///     dC............: parallel initial capacitance.
 192: ///     dR............: initial resistance.
 193: ///     dCommand_init.: initial command voltage.
 194: ///     dGain.........: gain value.
 195: ///     dTau_i........: integrating time constant.
 196: ///     dTau_d........: derivative time constant.
 197: ///     dSaturation...: saturation value.
 198: ///
 199: /// RTN..: int
 200: ///
 201: ///     success of operation.
 202: ///
 203: /// DESCR: set operation fields of voltage clamper.
 204: ///
 205: /// NOTE.:
 206: ///
 207: ///     The initial command voltage should probably have a separate
 208: ///     setter method.
 209: ///
 210: /// **************************************************************************
 211: 
 212: int
 213: VClampSetFields
 214: (struct VClamp * pvc,
 215:  double dInjected,
 216:  double dC,
 217:  double dR,
 218:  double dCommand_init,
 219:  double dGain,
 220:  double dTau_i,
 221:  double dTau_d,
 222:  double dSaturation)
 223: {
 224:     //- set default result: ok
 225: 
 226:     int iResult = 1;
 227: 
 228:     //- set fields
 229: 
 230:     pvc->dInjected = dInjected;
 231: 
 232:     pvc->dC = dC;
 233:     pvc->dR = dR;
 234:     pvc->dCommand = dCommand_init;
 235: 
 236:     pvc->dGain = dGain;
 237:     pvc->dTau_i = dTau_i;
 238:     pvc->dTau_d = dTau_d;
 239:     pvc->dSaturation = dSaturation;
 240: 
 241:     //- return result
 242: 
 243:     return(iResult);
 244: }
 245: 
 246: 
 247: /// **************************************************************************
 248: ///
 249: /// SHORT: VClampSingleStep()
 250: ///
 251: /// ARGS.:
 252: ///
 253: ///     pvc...: voltage clamper.
 254: ///     dTime.: current simulation time.
 255: ///
 256: /// RTN..: int
 257: ///
 258: ///     success of operation.
 259: ///
 260: /// DESCR: Compute new currents to correct voltages.
 261: ///
 262: /// NOTE.:
 263: ///
 264: ///     Old current values are overwritten.
 265: ///
 266: /// **************************************************************************
 267: 
 268: int VClampSingleStep(struct VClamp * pvc, double dTime)
 269: {
 270:     //- set default result: ok
 271: 
 272:     int iResult = 1;
 273: 
 274:     //- calculate the time step
 275: 
 276:     double dStep = dTime - pvc->dPreviousTime;
 277: 
 278:     //- exponential euler integration for the RC circuit
 279: 
 280:     double dA = pvc->dInjected / pvc->dC;
 281: 
 282:     double dB = 1 / (pvc->dR * pvc->dC);
 283: 
 284:     double dD = exp( - dB * dStep);
 285: 
 286:     pvc->dCommand = (pvc->dCommand * dD + (dA / dB) * (1 - dD));
 287: 
 288:     //- first preserve the previous value
 289: 
 290:     pvc->dEPrevious = pvc->dE;
 291: 
 292:     pvc->dE = pvc->dCommand - *pvc->pdVoltage;
 293: 
 294:     pvc->dEDerivative = (pvc->dE - pvc->dEPrevious) / dStep;
 295: 
 296:     pvc->dEIntegral += 0.5 * (pvc->dE + pvc->dEPrevious) * dStep;
 297: 
 298:     //- compute output
 299: 
 300:     pvc->dOutput
 301:         = (pvc->dGain *
 302:            (pvc->dE
 303:             + pvc->dTau_d * pvc->dEDerivative
 304:             + 1 / pvc->dTau_i * pvc->dEIntegral));
 305: 
 306:     //- apply saturation
 307: 
 308:     if (pvc->dOutput > pvc->dSaturation)
 309:     {
 310:         pvc->dOutput = pvc->dSaturation;
 311: 
 312:         pvc->dEIntegral -= 0.5 * (pvc->dE + pvc->dEPrevious) * dStep;
 313:     }
 314:     else if (pvc->dOutput < - pvc->dSaturation)
 315:     {
 316:         pvc->dOutput = - pvc->dSaturation;
 317: 
 318:         pvc->dEIntegral -= 0.5 * (pvc->dE + pvc->dEPrevious) * dStep;
 319:     }
 320: 
 321:     //- set the output
 322: 
 323:     *pvc->pdInjector = pvc->dOutput;
 324: 
 325:     //- register the current simulation time for the next time we are scheduled
 326: 
 327:     pvc->dPreviousTime = dTime;
 328: 
 329:     //- return result
 330: 
 331:     return(iResult);
 332: }
 333: 
 334: 
 335: