diff --git a/moose-core/biophysics/CMakeLists.txt b/moose-core/biophysics/CMakeLists.txt
index 27c34e12eebb6095c6a4b2ff17e5710ce703a973..c36def8e3799492b1dd42eeed476f476515719a8 100644
--- a/moose-core/biophysics/CMakeLists.txt
+++ b/moose-core/biophysics/CMakeLists.txt
@@ -33,8 +33,12 @@ set(BIOPHYSICS_SRCS
SynChan.cpp
NMDAChan.cpp
testBiophysics.cpp
- IzhikevichNrn.cpp
- DifShell.cpp
+ IzhikevichNrn.cpp
+ DifShellBase.cpp
+ DifShell.cpp
+ DifBufferBase.cpp
+ DifBuffer.cpp
+ MMPump.cpp
Leakage.cpp
VectorTable.cpp
MarkovRateTable.cpp
diff --git a/moose-core/biophysics/CaConcBase.cpp b/moose-core/biophysics/CaConcBase.cpp
index 2d66b70940c9820f6b958d32594b64f220056b77..c484b24827292cb1d2f8b031c066eb2c6562de87 100644
--- a/moose-core/biophysics/CaConcBase.cpp
+++ b/moose-core/biophysics/CaConcBase.cpp
@@ -1,11 +1,11 @@
/**********************************************************************
-** This program is part of 'MOOSE', the
-** Messaging Object Oriented Simulation Environment.
-** Copyright (C) 2003-2014 Upinder S. Bhalla. and NCBS
-** It is made available under the terms of the
-** GNU Lesser General Public License version 2.1
-** See the file COPYING.LIB for the full notice.
-**********************************************************************/
+ ** This program is part of 'MOOSE', the
+ ** Messaging Object Oriented Simulation Environment.
+ ** Copyright (C) 2003-2014 Upinder S. Bhalla. and NCBS
+ ** It is made available under the terms of the
+ ** GNU Lesser General Public License version 2.1
+ ** See the file COPYING.LIB for the full notice.
+ **********************************************************************/
// #include <cfloat>
#include "header.h"
@@ -23,178 +23,178 @@
*/
// Static function.
SrcFinfo1< double >* CaConcBase::concOut() {
- static SrcFinfo1< double > concOut( "concOut",
- "Concentration of Ca in pool" );
- return &concOut;
+ static SrcFinfo1< double > concOut( "concOut",
+ "Concentration of Ca in pool" );
+ return &concOut;
}
const Cinfo* CaConcBase::initCinfo()
{
- ///////////////////////////////////////////////////////
- // Shared message definitions
- ///////////////////////////////////////////////////////
- static DestFinfo process( "process",
- "Handles process call",
- new ProcOpFunc< CaConcBase >( &CaConcBase::process ) );
- static DestFinfo reinit( "reinit",
- "Handles reinit call",
- new ProcOpFunc< CaConcBase >( &CaConcBase::reinit ) );
-
- static Finfo* processShared[] =
- {
- &process, &reinit
- };
-
- static SharedFinfo proc( "proc",
- "Shared message to receive Process message from scheduler",
- processShared, sizeof( processShared ) / sizeof( Finfo* ) );
+ ///////////////////////////////////////////////////////
+ // Shared message definitions
+ ///////////////////////////////////////////////////////
+ static DestFinfo process( "process",
+ "Handles process call",
+ new ProcOpFunc< CaConcBase >( &CaConcBase::process ) );
+ static DestFinfo reinit( "reinit",
+ "Handles reinit call",
+ new ProcOpFunc< CaConcBase >( &CaConcBase::reinit ) );
+
+ static Finfo* processShared[] =
+ {
+ &process, &reinit
+ };
+
+ static SharedFinfo proc( "proc",
+ "Shared message to receive Process message from scheduler",
+ processShared, sizeof( processShared ) / sizeof( Finfo* ) );
-///////////////////////////////////////////////////////
-// Field definitions
-///////////////////////////////////////////////////////
- static ElementValueFinfo< CaConcBase, double > Ca( "Ca",
- "Calcium concentration.",
- &CaConcBase::setCa,
- &CaConcBase::getCa
- );
- static ElementValueFinfo< CaConcBase, double > CaBasal( "CaBasal",
- "Basal Calcium concentration.",
- &CaConcBase::setCaBasal,
- &CaConcBase::getCaBasal
- );
- static ElementValueFinfo< CaConcBase, double > Ca_base( "Ca_base",
- "Basal Calcium concentration, synonym for CaBasal",
- &CaConcBase::setCaBasal,
- &CaConcBase::getCaBasal
- );
- static ElementValueFinfo< CaConcBase, double > tau( "tau",
- "Settling time for Ca concentration",
- &CaConcBase::setTau,
- &CaConcBase::getTau
- );
- static ElementValueFinfo< CaConcBase, double > B( "B",
- "Volume scaling factor. "
- "Deprecated. This is a legacy field from GENESIS and exposes "
- "internal calculations. Please do not use. \n"
- "B = 1/(vol * F) so that it obeys:\n"
- "dC/dt = B*I_Ca - C/tau",
- &CaConcBase::setB,
- &CaConcBase::getB
- );
- static ElementValueFinfo< CaConcBase, double > thick( "thick",
- "Thickness of Ca shell, assumed cylindrical. Legal range is "
- "between zero and the radius. If outside this range it is "
- "taken as the radius. Default zero, ie, the shell is the entire "
- "thickness of the cylinder",
- &CaConcBase::setThickness,
- &CaConcBase::getThickness
- );
- static ElementValueFinfo< CaConcBase, double > length( "length",
- "Length of Ca shell, assumed cylindrical",
- &CaConcBase::setLength,
- &CaConcBase::getLength
- );
- static ElementValueFinfo< CaConcBase, double > diameter( "diameter",
- "Diameter of Ca shell, assumed cylindrical",
- &CaConcBase::setDiameter,
- &CaConcBase::getDiameter
- );
- static ElementValueFinfo< CaConcBase, double > ceiling( "ceiling",
- "Ceiling value for Ca concentration. If Ca > ceiling, Ca = ceiling. If ceiling <= 0.0, there is no upper limit on Ca concentration value.",
- &CaConcBase::setCeiling,
- &CaConcBase::getCeiling
- );
- static ElementValueFinfo< CaConcBase, double > floor( "floor",
- "Floor value for Ca concentration. If Ca < floor, Ca = floor",
- &CaConcBase::setFloor,
- &CaConcBase::getFloor
- );
-
-///////////////////////////////////////////////////////
-// MsgDest definitions
-///////////////////////////////////////////////////////
-
- static DestFinfo current( "current",
- "Calcium Ion current, due to be converted to conc.",
- new EpFunc1< CaConcBase, double >( &CaConcBase::current )
- );
-
- static DestFinfo currentFraction( "currentFraction",
- "Fraction of total Ion current, that is carried by Ca2+.",
- new EpFunc2< CaConcBase, double, double >( &CaConcBase::currentFraction )
- );
-
- static DestFinfo increase( "increase",
- "Any input current that increases the concentration.",
- new EpFunc1< CaConcBase, double >( &CaConcBase::increase )
- );
-
- static DestFinfo decrease( "decrease",
- "Any input current that decreases the concentration.",
- new EpFunc1< CaConcBase, double >( &CaConcBase::decrease )
- );
-
- static DestFinfo basal( "basal",
- "Synonym for assignment of basal conc.",
- new EpFunc1< CaConcBase, double >( &CaConcBase::setCaBasal )
- );
-
- static Finfo* CaConcBaseFinfos[] =
- {
- &proc, // Shared
- concOut(), // Src
- &Ca, // Value
- &CaBasal, // Value
- &Ca_base, // Value
- &tau, // Value
- &B, // Value
- &thick, // Value
- &diameter, // Value
- &length, // Value
- &ceiling, // Value
- &floor, // Value
- ¤t, // Dest
- ¤tFraction, // Dest
- &increase, // Dest
- &decrease, // Dest
- &basal, // Dest
- };
-
- // We want the Ca updates before channel updates, so along with compts.
- // static SchedInfo schedInfo[] = { { process, 0, 0 } };
-
- static string doc[] =
- {
- "Name", "CaConcBase",
- "Author", "Upinder S. Bhalla, 2014, NCBS",
- "Description",
- "CaConcBase: Base class for Calcium concentration "
- "pool. Takes current from a channel and keeps track of "
- "calcium buildup and depletion by a "
- "single exponential process. ",
- };
+ ///////////////////////////////////////////////////////
+ // Field definitions
+ ///////////////////////////////////////////////////////
+ static ElementValueFinfo< CaConcBase, double > Ca( "Ca",
+ "Calcium concentration.",
+ &CaConcBase::setCa,
+ &CaConcBase::getCa
+ );
+ static ElementValueFinfo< CaConcBase, double > CaBasal( "CaBasal",
+ "Basal Calcium concentration.",
+ &CaConcBase::setCaBasal,
+ &CaConcBase::getCaBasal
+ );
+ static ElementValueFinfo< CaConcBase, double > Ca_base( "Ca_base",
+ "Basal Calcium concentration, synonym for CaBasal",
+ &CaConcBase::setCaBasal,
+ &CaConcBase::getCaBasal
+ );
+ static ElementValueFinfo< CaConcBase, double > tau( "tau",
+ "Settling time for Ca concentration",
+ &CaConcBase::setTau,
+ &CaConcBase::getTau
+ );
+ static ElementValueFinfo< CaConcBase, double > B( "B",
+ "Volume scaling factor. "
+ "Deprecated. This is a legacy field from GENESIS and exposes "
+ "internal calculations. Please do not use. \n"
+ "B = 1/(vol * F) so that it obeys:\n"
+ "dC/dt = B*I_Ca - C/tau",
+ &CaConcBase::setB,
+ &CaConcBase::getB
+ );
+ static ElementValueFinfo< CaConcBase, double > thick( "thick",
+ "Thickness of Ca shell, assumed cylindrical. Legal range is "
+ "between zero and the radius. If outside this range it is "
+ "taken as the radius. Default zero, ie, the shell is the entire "
+ "thickness of the cylinder",
+ &CaConcBase::setThickness,
+ &CaConcBase::getThickness
+ );
+ static ElementValueFinfo< CaConcBase, double > length( "length",
+ "Length of Ca shell, assumed cylindrical",
+ &CaConcBase::setLength,
+ &CaConcBase::getLength
+ );
+ static ElementValueFinfo< CaConcBase, double > diameter( "diameter",
+ "Diameter of Ca shell, assumed cylindrical",
+ &CaConcBase::setDiameter,
+ &CaConcBase::getDiameter
+ );
+ static ElementValueFinfo< CaConcBase, double > ceiling( "ceiling",
+ "Ceiling value for Ca concentration. If Ca > ceiling, Ca = ceiling. If ceiling <= 0.0, there is no upper limit on Ca concentration value.",
+ &CaConcBase::setCeiling,
+ &CaConcBase::getCeiling
+ );
+ static ElementValueFinfo< CaConcBase, double > floor( "floor",
+ "Floor value for Ca concentration. If Ca < floor, Ca = floor",
+ &CaConcBase::setFloor,
+ &CaConcBase::getFloor
+ );
+
+ ///////////////////////////////////////////////////////
+ // MsgDest definitions
+ ///////////////////////////////////////////////////////
+
+ static DestFinfo current( "current",
+ "Calcium Ion current, due to be converted to conc.",
+ new EpFunc1< CaConcBase, double >( &CaConcBase::current )
+ );
+
+ static DestFinfo currentFraction( "currentFraction",
+ "Fraction of total Ion current, that is carried by Ca2+.",
+ new EpFunc2< CaConcBase, double, double >( &CaConcBase::currentFraction )
+ );
+
+ static DestFinfo increase( "increase",
+ "Any input current that increases the concentration.",
+ new EpFunc1< CaConcBase, double >( &CaConcBase::increase )
+ );
+
+ static DestFinfo decrease( "decrease",
+ "Any input current that decreases the concentration.",
+ new EpFunc1< CaConcBase, double >( &CaConcBase::decrease )
+ );
+
+ static DestFinfo basal( "basal",
+ "Synonym for assignment of basal conc.",
+ new EpFunc1< CaConcBase, double >( &CaConcBase::setCaBasal )
+ );
+
+ static Finfo* CaConcBaseFinfos[] =
+ {
+ &proc, // Shared
+ concOut(), // Src
+ &Ca, // Value
+ &CaBasal, // Value
+ &Ca_base, // Value
+ &tau, // Value
+ &B, // Value
+ &thick, // Value
+ &diameter, // Value
+ &length, // Value
+ &ceiling, // Value
+ &floor, // Value
+ ¤t, // Dest
+ ¤tFraction, // Dest
+ &increase, // Dest
+ &decrease, // Dest
+ &basal, // Dest
+ };
+
+ // We want the Ca updates before channel updates, so along with compts.
+ // static SchedInfo schedInfo[] = { { process, 0, 0 } };
+
+ static string doc[] =
+ {
+ "Name", "CaConcBase",
+ "Author", "Upinder S. Bhalla, 2014, NCBS",
+ "Description",
+ "CaConcBase: Base class for Calcium concentration "
+ "pool. Takes current from a channel and keeps track of "
+ "calcium buildup and depletion by a "
+ "single exponential process. ",
+ };
- static ZeroSizeDinfo< int > dinfo;
-
- static Cinfo CaConcBaseCinfo(
- "CaConcBase",
- Neutral::initCinfo(),
- CaConcBaseFinfos,
- sizeof( CaConcBaseFinfos )/sizeof(Finfo *),
- &dinfo,
- doc,
- sizeof(doc)/sizeof(string)
- );
-
- return &CaConcBaseCinfo;
+ static ZeroSizeDinfo< int > dinfo;
+
+ static Cinfo CaConcBaseCinfo(
+ "CaConcBase",
+ Neutral::initCinfo(),
+ CaConcBaseFinfos,
+ sizeof( CaConcBaseFinfos )/sizeof(Finfo *),
+ &dinfo,
+ doc,
+ sizeof(doc)/sizeof(string)
+ );
+
+ return &CaConcBaseCinfo;
}
///////////////////////////////////////////////////
static const Cinfo* caConcCinfo = CaConcBase::initCinfo();
CaConcBase::CaConcBase()
- :
- thickness_( 0.0 )
+ :
+ thickness_( 0.0 )
{;}
///////////////////////////////////////////////////
@@ -203,98 +203,98 @@ CaConcBase::CaConcBase()
void CaConcBase::setCa( const Eref& e, double Ca )
{
- vSetCa( e, Ca );
+ vSetCa( e, Ca );
}
double CaConcBase::getCa( const Eref& e ) const
{
- return vGetCa( e );
+ return vGetCa( e );
}
void CaConcBase::setCaBasal( const Eref& e, double CaBasal )
{
- vSetCaBasal( e, CaBasal );
+ vSetCaBasal( e, CaBasal );
}
double CaConcBase::getCaBasal( const Eref& e ) const
{
- return vGetCaBasal( e );
+ return vGetCaBasal( e );
}
void CaConcBase::setTau( const Eref& e, double tau )
{
- vSetTau( e, tau );
+ vSetTau( e, tau );
}
double CaConcBase::getTau( const Eref& e ) const
{
- return vGetTau( e );
+ return vGetTau( e );
}
void CaConcBase::setB( const Eref& e, double B )
{
- vSetB( e, B );
+ vSetB( e, B );
}
double CaConcBase::getB( const Eref& e ) const
{
- return vGetB( e );
+ return vGetB( e );
}
void CaConcBase::setCeiling( const Eref& e, double ceiling )
{
- vSetCeiling( e, ceiling );
+ vSetCeiling( e, ceiling );
}
double CaConcBase::getCeiling( const Eref& e ) const
{
- return vGetCeiling( e );
+ return vGetCeiling( e );
}
void CaConcBase::setFloor( const Eref& e, double floor )
{
- vSetFloor( e, floor );
+ vSetFloor( e, floor );
}
double CaConcBase::getFloor( const Eref& e ) const
{
- return vGetFloor( e );
+ return vGetFloor( e );
}
void CaConcBase::updateDimensions( const Eref& e )
{
- double vol = PI * diameter_ * diameter_ * length_ * 0.25;
- if ( thickness_ > 0 && thickness_ < diameter_/2.0 ) {
- double coreRadius = diameter_ / 2.0 - thickness_;
- vol -= PI * coreRadius * coreRadius * length_;
- }
- double B = 1.0 / ( FaradayConst * vol );
- vSetB( e, B );
+ double vol = PI * diameter_ * diameter_ * length_ * 0.25;
+ if ( thickness_ > 0 && thickness_ < diameter_/2.0 ) {
+ double coreRadius = diameter_ / 2.0 - thickness_;
+ vol -= PI * coreRadius * coreRadius * length_;
+ }
+ double B = 1.0 / ( FaradayConst * vol );
+ vSetB( e, B );
}
void CaConcBase::setThickness( const Eref& e, double thickness )
{
- thickness_ = thickness;
- updateDimensions( e );
+ thickness_ = thickness;
+ updateDimensions( e );
}
double CaConcBase::getThickness( const Eref& e ) const
{
- return thickness_;
+ return thickness_;
}
void CaConcBase::setDiameter( const Eref& e, double diameter )
{
- diameter_ = diameter;
- updateDimensions( e );
+ diameter_ = diameter;
+ updateDimensions( e );
}
double CaConcBase::getDiameter( const Eref& e ) const
{
- return diameter_;
+ return diameter_;
}
void CaConcBase::setLength( const Eref& e, double length )
{
- length_ = length;
- updateDimensions( e );
+ length_ = length;
+ updateDimensions( e );
}
double CaConcBase::getLength( const Eref& e ) const
{
- return length_;
+ return length_;
}
///////////////////////////////////////////////////
@@ -303,32 +303,32 @@ double CaConcBase::getLength( const Eref& e ) const
void CaConcBase::reinit( const Eref& e, ProcPtr p )
{
- vReinit( e, p );
+ vReinit( e, p );
}
void CaConcBase::process( const Eref& e, ProcPtr p )
{
- vProcess( e, p );
+ vProcess( e, p );
}
void CaConcBase::current( const Eref& e, double I )
{
- vCurrent( e, I );
+ vCurrent( e, I );
}
void CaConcBase::currentFraction( const Eref& e, double I, double fraction )
{
- vCurrentFraction( e, I, fraction );
+ vCurrentFraction( e, I, fraction );
}
void CaConcBase::increase( const Eref& e, double I )
{
- vIncrease( e, I );
+ vIncrease( e, I );
}
void CaConcBase::decrease( const Eref& e, double I )
{
- vDecrease( e, I );
+ vDecrease( e, I );
}
/////////////////////////////////////////////////////
@@ -342,48 +342,48 @@ void CaConcBase::vSetSolver( const Eref& e, Id hsolve )
// static func
void CaConcBase::zombify( Element* orig, const Cinfo* zClass,
- Id hsolve )
+ Id hsolve )
{
- if ( orig->cinfo() == zClass )
- return;
- unsigned int start = orig->localDataStart();
- unsigned int num = orig->numLocalData();
- if ( num == 0 )
- return;
- vector< double > data( num * 9 );
-
- unsigned int j = 0;
- for ( unsigned int i = 0; i < num; ++i ) {
- Eref er( orig, i + start );
- const CaConcBase* cb =
- reinterpret_cast< const CaConcBase* >( er.data() );
- data[j + 0] = cb->getCa( er );
- data[j + 1] = cb->getCaBasal( er );
- data[j + 2] = cb->getTau( er );
- data[j + 3] = cb->getB( er );
- data[j + 4] = cb->getCeiling( er );
- data[j + 5] = cb->getFloor( er );
- data[j + 6] = cb->getThickness( er );
- data[j + 7] = cb->getLength( er );
- data[j + 8] = cb->getDiameter( er );
- j += 9;
- }
- orig->zombieSwap( zClass );
- j = 0;
- for ( unsigned int i = 0; i < num; ++i ) {
- Eref er( orig, i + start );
- CaConcBase* cb = reinterpret_cast< CaConcBase* >( er.data() );
- cb->vSetSolver( er, hsolve );
- cb->setCa( er, data[j + 0] );
- cb->setCaBasal( er, data[j + 1] );
- cb->setTau( er, data[j + 2] );
- cb->setB( er, data[j + 3] );
- cb->setCeiling( er, data[j + 4] );
- cb->setFloor( er, data[j + 5] );
- cb->setThickness( er, data[j + 6] );
- cb->setLength( er, data[j + 7] );
- cb->setDiameter( er, data[j + 8] );
- j += 7;
- }
+ if ( orig->cinfo() == zClass )
+ return;
+ unsigned int start = orig->localDataStart();
+ unsigned int num = orig->numLocalData();
+ if ( num == 0 )
+ return;
+ vector< double > data( num * 9 );
+
+ unsigned int j = 0;
+ for ( unsigned int i = 0; i < num; ++i ) {
+ Eref er( orig, i + start );
+ const CaConcBase* cb =
+ reinterpret_cast< const CaConcBase* >( er.data() );
+ data[j + 0] = cb->getCa( er );
+ data[j + 1] = cb->getCaBasal( er );
+ data[j + 2] = cb->getTau( er );
+ data[j + 3] = cb->getB( er );
+ data[j + 4] = cb->getCeiling( er );
+ data[j + 5] = cb->getFloor( er );
+ data[j + 6] = cb->getThickness( er );
+ data[j + 7] = cb->getLength( er );
+ data[j + 8] = cb->getDiameter( er );
+ j += 9;
+ }
+ orig->zombieSwap( zClass );
+ j = 0;
+ for ( unsigned int i = 0; i < num; ++i ) {
+ Eref er( orig, i + start );
+ CaConcBase* cb = reinterpret_cast< CaConcBase* >( er.data() );
+ cb->vSetSolver( er, hsolve );
+ cb->setCa( er, data[j + 0] );
+ cb->setCaBasal( er, data[j + 1] );
+ cb->setTau( er, data[j + 2] );
+ cb->setB( er, data[j + 3] );
+ cb->setCeiling( er, data[j + 4] );
+ cb->setFloor( er, data[j + 5] );
+ cb->setThickness( er, data[j + 6] );
+ cb->setLength( er, data[j + 7] );
+ cb->setDiameter( er, data[j + 8] );
+ j += 9; //was 7?
+ }
}
diff --git a/moose-core/biophysics/DifBuffer.cpp b/moose-core/biophysics/DifBuffer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..60763e82b86ce5794edf170948758d40d31f214c
--- /dev/null
+++ b/moose-core/biophysics/DifBuffer.cpp
@@ -0,0 +1,467 @@
+// DifBuffer.cpp ---
+//
+// Filename: DifBuffer.cpp
+// Description:
+// Author: Subhasis Ray
+// Maintainer:
+// Created: Mon Feb 16 12:02:11 2015 (-0500)
+// Version:
+// Package-Requires: ()
+// Last-Updated: Mon Feb 23 13:07:56 2015 (-0500)
+// By: Subhasis Ray
+// Update #: 130
+// URL:
+// Doc URL:
+// Keywords:
+// Compatibility:
+//
+//
+
+// Commentary:
+//
+//
+//
+//
+
+// Change Log:
+// 5/25/16 Completing DifBuffer -- Asia J-Szmek (GMU)
+// 9/21/16 rewrote DifBuffer to account for DifBufferBase (AJS)
+//
+//
+// This program is free software: you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation, either version 3 of the License, or (at
+// your option) any later version.
+//
+// This program is distributed in the hope that it will be useful, but
+// WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+// General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with GNU Emacs. If not, see <http://www.gnu.org/licenses/>.
+//
+//
+
+// Code:
+
+#include "header.h"
+#include "DifBufferBase.h"
+#include "DifBuffer.h"
+#include "ElementValueFinfo.h"
+#include "../utility/numutil.h"
+#include <cmath>
+const double DifBuffer::EPSILON = 1.0e-10;
+
+
+const Cinfo * DifBuffer::initCinfo()
+{
+
+ static string doc[] = {
+ "Name", "DifBuffer",
+ "Author", "Subhasis Ray (ported from GENESIS2)",
+ "Description", "Models diffusible buffer where total concentration is constant. It is"
+ " coupled with a DifShell.",
+ };
+ static Dinfo<DifBuffer> dinfo;
+ static Cinfo difBufferCinfo(
+ "DifBuffer",
+ DifBufferBase::initCinfo(),
+ 0,
+ 0,
+ &dinfo,
+ doc,
+ sizeof(doc)/sizeof(string));
+
+ return &difBufferCinfo;
+}
+
+static const Cinfo * difBufferCinfo = DifBuffer::initCinfo();
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Class functions
+////////////////////////////////////////////////////////////////////////////////
+
+DifBuffer::DifBuffer() :
+ activation_(0),
+ Af_(0),
+ Bf_(0),
+ bFree_(0),
+ bBound_(0),
+ prevFree_(0),
+ prevBound_(0),
+ bTot_(0),
+ kf_(0),
+ kb_(0),
+ D_(0),
+ shapeMode_(0),
+ length_(0),
+ diameter_(0),
+ thickness_(0),
+ volume_(0),
+ outerArea_(0),
+ innerArea_(0)
+
+{}
+
+////////////////////////////////////////////////////////////////////////////////
+// Field access functions
+////////////////////////////////////////////////////////////////////////////////
+double DifBuffer::vGetActivation(const Eref& e) const
+{
+ return activation_;
+}
+
+void DifBuffer::vSetActivation(const Eref& e,double value)
+{
+ if ( value < 0.0 ) {
+ cerr << "Error: DifBuffer: Activation cannot be negative!\n";
+ return;
+ }
+ activation_ = value;
+}
+
+
+double DifBuffer::vGetBFree(const Eref& e) const
+{
+ return bFree_;
+}
+void DifBuffer::vSetBFree(const Eref& e,double value)
+{
+ if ( value < 0.0 ) {
+ cerr << "Error: DifBuffer: Free Buffer cannot be negative!\n";
+ return;
+ }
+ if (value > bTot_){
+ cerr << "Error: DifBuffer: Free Buffer cannot exceed total buffer!\n";
+ return;
+ }
+ bFree_ = value;
+ bBound_ = bTot_ - bFree_;
+ prevFree_= bFree_;
+ prevBound_ = bBound_;
+
+}
+
+double DifBuffer::vGetBBound(const Eref& e) const
+{
+ return bBound_;
+}
+
+void DifBuffer::vSetBBound(const Eref& e,double value)
+{
+ if ( value < 0.0 ) {
+ cerr << "Error: DifBuffer: Bound Buffer cannot be negative!\n";
+ return;
+ }
+ if (value > bTot_){
+ cerr << "Error: DifBuffer: Bound Buffer cannot exceed total buffer!\n";
+ return;
+ }
+ bBound_ = value;
+ bFree_ = bTot_ - bBound_;
+ prevFree_= bFree_;
+ prevBound_ = bBound_;
+}
+
+
+double DifBuffer::vGetBTot(const Eref& e) const
+{
+ return bTot_;
+}
+
+void DifBuffer::vSetBTot(const Eref& e,double value)
+{
+ if ( value < 0.0 ) {
+ cerr << "Error: DifBuffer: Total buffer concentration cannot be negative!\n";
+ return;
+ }
+ bTot_ = value;
+ bFree_ = bTot_;
+ bBound_ = 0;
+}
+
+
+double DifBuffer::vGetKf(const Eref& e) const
+{
+ return kf_;
+}
+
+void DifBuffer::vSetKf(const Eref& e,double value)
+{
+ if ( value < 0.0 ) {
+ cerr << "Error: DifBuffer: Kf cannot be negative!\n";
+ return;
+ }
+ kf_ = value;
+}
+
+
+double DifBuffer::vGetKb(const Eref& e) const
+{
+ return kb_;
+}
+
+void DifBuffer::vSetKb(const Eref& e,double value)
+{
+ if ( value < 0.0 ) {
+ cerr << "Error: DifBuffer: Kb cannot be negative!\n";
+ return;
+ }
+ kb_ = value;
+}
+
+double DifBuffer::vGetD(const Eref& e) const
+{
+ return D_;
+}
+
+void DifBuffer::vSetD(const Eref& e,double value)
+{
+
+ if ( value < 0.0 ) {
+ cerr << " Error: DifBuffer: Diffusion constant, D, cannot be negative!\n";
+ return;
+ }
+ D_ = value;
+}
+
+
+void DifBuffer::vSetShapeMode(const Eref& e, unsigned int shapeMode )
+{
+ if ( shapeMode != 0 && shapeMode != 1 && shapeMode != 3 ) {
+ cerr << "Error: DifBuffer: I only understand shapeModes 0, 1 and 3.\n";
+ return;
+ }
+ shapeMode_ = shapeMode;
+}
+
+unsigned int DifBuffer::vGetShapeMode(const Eref& e) const
+{
+ return shapeMode_;
+}
+
+void DifBuffer::vSetLength(const Eref& e, double length )
+{
+ if ( length < 0.0 ) {
+ cerr << "Error: DifBuffer: length cannot be negative!\n";
+ return;
+ }
+
+ length_ = length;
+}
+
+double DifBuffer::vGetLength(const Eref& e ) const
+{
+ return length_;
+}
+
+void DifBuffer::vSetDiameter(const Eref& e, double diameter )
+{
+ if ( diameter < 0.0 ) {
+ cerr << "Error: DifBuffer: diameter cannot be negative!\n";
+ return;
+ }
+
+ diameter_ = diameter;
+}
+
+double DifBuffer::vGetDiameter(const Eref& e ) const
+{
+ return diameter_;
+}
+
+void DifBuffer::vSetThickness( const Eref& e, double thickness )
+{
+ if ( thickness < 0.0 ) {
+ cerr << "Error: DifBuffer: thickness cannot be negative!\n";
+ return;
+ }
+
+ thickness_ = thickness;
+}
+
+double DifBuffer::vGetThickness(const Eref& e) const
+{
+ return thickness_;
+}
+
+void DifBuffer::vSetVolume(const Eref& e, double volume )
+{
+ if ( shapeMode_ != 3 )
+ cerr << "Warning: DifBuffer: Trying to set volume, when shapeMode is not USER-DEFINED\n";
+
+ if ( volume < 0.0 ) {
+ cerr << "Error: DifBuffer: volume cannot be negative!\n";
+ return;
+ }
+
+ volume_ = volume;
+}
+
+double DifBuffer::vGetVolume(const Eref& e ) const
+{
+ return volume_;
+}
+
+void DifBuffer::vSetOuterArea(const Eref& e, double outerArea )
+{
+ if (shapeMode_ != 3 )
+ cerr << "Warning: DifBuffer: Trying to set outerArea, when shapeMode is not USER-DEFINED\n";
+
+ if ( outerArea < 0.0 ) {
+ cerr << "Error: DifBuffer: outerArea cannot be negative!\n";
+ return;
+ }
+
+ outerArea_ = outerArea;
+}
+
+double DifBuffer::vGetOuterArea(const Eref& e ) const
+{
+ return outerArea_;
+}
+
+void DifBuffer::vSetInnerArea(const Eref& e, double innerArea )
+{
+ if ( shapeMode_ != 3 )
+ cerr << "Warning: DifBuffer: Trying to set innerArea, when shapeMode is not USER-DEFINED\n";
+
+ if ( innerArea < 0.0 ) {
+ cerr << "Error: DifBuffer: innerArea cannot be negative!\n";
+ return;
+ }
+
+ innerArea_ = innerArea;
+}
+
+double DifBuffer::vGetInnerArea(const Eref& e) const
+{
+ return innerArea_;
+}
+
+
+
+
+void DifBuffer::vBuffer(const Eref& e,
+ double C )
+{
+ activation_ = C;
+ //cout<<"Buffer"<< C<<" ";
+}
+
+
+double DifBuffer::integrate( double state, double dt, double A, double B )
+{
+ if ( B > EPSILON ) {
+ double x = exp( -B * dt );
+ return state * x + ( A / B ) * ( 1 - x );
+ }
+ return state + A * dt ;
+}
+
+
+void DifBuffer::vProcess( const Eref & e, ProcPtr p )
+{
+ /**
+ * Send ion concentration and thickness to adjacent DifShells. They will
+ * then compute their incoming fluxes.
+ */
+
+
+
+ Af_ += kb_ * bBound_;
+ Bf_ += kf_ * activation_;
+
+ bFree_ = integrate(bFree_,p->dt,Af_,Bf_);
+ bBound_ = bTot_ - bFree_;
+ prevFree_ = bFree_;
+ prevBound_ = bBound_;
+
+ /**
+ * Send ion concentration to ion buffers. They will send back information on
+ * the reaction (forward / backward rates ; free / bound buffer concentration)
+ * immediately, which this DifShell will use to find amount of ion captured
+ * or released in the current time-step.
+ */
+ innerDifSourceOut()->send( e, prevFree_, thickness_ );
+ outerDifSourceOut()->send( e, prevFree_, thickness_ );
+ reactionOut()->send(e,kf_,kb_,bFree_,bBound_);
+ Af_ = 0;
+ Bf_= 0;
+
+}
+
+void DifBuffer::vReinit( const Eref& e, ProcPtr p )
+{
+
+ Af_ = 0;
+ Bf_= 0;
+ double rOut = diameter_/2.;
+
+ double rIn = rOut - thickness_;
+
+ if (rIn<0)
+ rIn = 0.;
+ switch ( shapeMode_ )
+ {
+ /*
+ * Onion Shell
+ */
+ case 0:
+ if ( length_ == 0.0 ) { // Spherical shell
+ volume_ = 4./3.* M_PI * ( rOut * rOut * rOut - rIn * rIn * rIn );
+ outerArea_ = 4*M_PI * rOut * rOut;
+ innerArea_ = 4*M_PI * rIn * rIn;
+ } else { // Cylindrical shell
+ volume_ = ( M_PI * length_ ) * ( rOut * rOut - rIn * rIn );
+ outerArea_ = 2*M_PI * rOut * length_;
+ innerArea_ = 2*M_PI * rIn * length_;
+ }
+
+ break;
+
+ /*
+ * Cylindrical Slice
+ */
+ case 1:
+ volume_ = M_PI * diameter_ * diameter_ * thickness_ / 4.0;
+ outerArea_ = M_PI * diameter_ * diameter_ / 4.0;
+ innerArea_ = outerArea_;
+ break;
+
+ /*
+ * User defined
+ */
+ case 3:
+ // Nothing to be done here. Volume and inner-, outer areas specified by
+ // user.
+ break;
+
+ default:
+ assert( 0 );
+ }
+
+ bFree_ = bTot_/(1+activation_*kf_/kb_);
+ prevFree_ = bFree_;
+ bBound_ = bTot_ - bFree_;
+ prevBound_ = bBound_;
+ innerDifSourceOut()->send( e, prevFree_, thickness_ );
+ outerDifSourceOut()->send( e, prevFree_, thickness_ );
+
+}
+
+void DifBuffer::vFluxFromIn(const Eref& e,double innerC, double innerThickness)
+{
+ double dif = 2 * D_ / volume_* innerArea_ / (thickness_ + innerThickness);
+ Af_ += dif * innerC;
+ Bf_ += dif;
+}
+
+void DifBuffer::vFluxFromOut(const Eref& e,double outerC, double outerThickness)
+{
+ double dif = 2 * D_ / volume_* outerArea_ / (thickness_ + outerThickness);
+ Af_ += dif * outerC;
+ Bf_ += dif;
+}
diff --git a/moose-core/biophysics/DifBuffer.h b/moose-core/biophysics/DifBuffer.h
new file mode 100644
index 0000000000000000000000000000000000000000..60f574458385487a732b3653d568e99fbba1a420
--- /dev/null
+++ b/moose-core/biophysics/DifBuffer.h
@@ -0,0 +1,96 @@
+/**********************************************************************
+ ** This program is part of 'MOOSE', the
+ ** Multiscale Object Oriented Simulation Environment.
+ ** copyright (C) 2003-2008
+ ** Upinder S. Bhalla, Niraj Dudani and NCBS
+ ** It is made available under the terms of the
+ ** GNU Lesser General Public License version 2.1
+ ** See the file COPYING.LIB for the full notice.
+ *****/
+#ifndef _DifBuffer_h
+#define _DifBuffer_h
+
+class DifBuffer: public DifBufferBase{
+ public:
+ DifBuffer();
+
+ void vBuffer(const Eref& e,double C);
+ void vReinit( const Eref & e, ProcPtr p );
+ void vProcess(const Eref & e, ProcPtr p );
+ void vFluxFromOut(const Eref& e,double outerC, double outerThickness );
+ void vFluxFromIn( const Eref& e,double innerC, double innerThickness );
+ //Field access functions
+
+ double vGetActivation(const Eref& e) const;
+ void vSetActivation(const Eref& e,double value);
+
+ double vGetBFree(const Eref& e) const;
+ void vSetBFree(const Eref& e,double value);
+
+ double vGetBBound(const Eref& e) const;
+ void vSetBBound(const Eref& e,double value);
+
+ double vGetBTot(const Eref& e) const; // total buffer concentration in mM (free + bound)
+ void vSetBTot(const Eref& e,double value);
+
+ double vGetKf(const Eref& e) const; // forward rate constant in 1/(mM*sec)
+ void vSetKf(const Eref& e,double value);
+
+ double vGetKb(const Eref& e) const; // backward rate constant in 1/sec
+ void vSetKb(const Eref& e,double value);
+
+ double vGetD(const Eref& e) const; // diffusion constant of buffer molecules, m^2/sec
+ void vSetD(const Eref& e,double value);
+ void vSetShapeMode(const Eref& e, unsigned int shapeMode );
+ unsigned int vGetShapeMode(const Eref& e) const;
+
+ void vSetLength(const Eref& e, double length );
+ double vGetLength(const Eref& e) const;
+
+ void vSetDiameter(const Eref& e, double diameter );
+ double vGetDiameter(const Eref& e) const;
+
+ void vSetThickness(const Eref& e, double thickness );
+ double vGetThickness(const Eref& e) const;
+
+ void vSetVolume(const Eref& e, double volume );
+ double vGetVolume(const Eref& e) const;
+
+ void vSetOuterArea(const Eref& e, double outerArea );
+ double vGetOuterArea(const Eref& e) const;
+
+ void vSetInnerArea(const Eref& e, double innerArea );
+ double vGetInnerArea(const Eref& e) const;
+
+ static const Cinfo * initCinfo();
+
+ private:
+
+
+ double integrate(double state, double dt, double A, double B );
+
+ double activation_; //ion concentration from incoming CONCEN message
+ double Af_;
+ double Bf_;
+ double bFree_; //free buffer concentration
+ double bBound_; //bound buffer concentration
+ double prevFree_;
+ double prevBound_;
+ //double prevFree_; //bFree at previous time
+ //double prevBound_; //bBound at previous time
+ double bTot_; //bFree+bBound
+ double kf_; //forward rate constant
+ double kb_; //backward rate constant
+ double D_; //diffusion constant
+ unsigned int shapeMode_;
+ double length_;
+ double diameter_;
+ double thickness_;
+ double volume_;
+ double outerArea_;
+ double innerArea_;
+ static const double EPSILON;
+};
+
+#endif // _DifBuffer_h
+/* DifBuffer.h ends here */
diff --git a/moose-core/biophysics/DifBufferBase.cpp b/moose-core/biophysics/DifBufferBase.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..40628a21e44d95612293230d45824ccb52c1b638
--- /dev/null
+++ b/moose-core/biophysics/DifBufferBase.cpp
@@ -0,0 +1,461 @@
+#include "header.h"
+#include "DifBufferBase.h"
+#include "ElementValueFinfo.h"
+#include "../utility/numutil.h"
+#include <cmath>
+
+SrcFinfo4< double, double, double, double >* DifBufferBase::reactionOut()
+{
+ static SrcFinfo4< double, double, double, double > reactionOut(
+ "reactionOut",
+ "Sends out reaction rates (forward and backward), and concentrations"
+ " (free-buffer and bound-buffer molecules).");
+ return &reactionOut;
+}
+
+
+SrcFinfo2< double, double >* DifBufferBase::innerDifSourceOut(){
+ static SrcFinfo2< double, double > sourceOut("innerDifSourceOut",
+ "Sends out source information.");
+ return &sourceOut;
+}
+
+SrcFinfo2< double, double >* DifBufferBase::outerDifSourceOut(){
+ static SrcFinfo2< double, double > sourceOut("outerDifSourceOut",
+ "Sends out source information.");
+ return &sourceOut;
+}
+
+const Cinfo * DifBufferBase::initCinfo()
+{
+ static DestFinfo process( "process",
+ "Handles process call",
+ new ProcOpFunc< DifBufferBase >( &DifBufferBase::process) );
+ static DestFinfo reinit( "reinit",
+ "Reinit happens only in stage 0",
+ new ProcOpFunc< DifBufferBase >( &DifBufferBase::reinit));
+
+ static Finfo* processShared[] = {
+ &process,
+ &reinit
+ };
+
+ static SharedFinfo proc(
+ "proc",
+ "Here we create 2 shared finfos to attach with the Ticks. This is because we want to perform DifBufferBase "
+ "computations in 2 stages, much as in the Compartment object. "
+ "In the first stage we send out the concentration value to other DifBufferBases and Buffer elements. We also",
+ processShared,
+ sizeof( processShared ) / sizeof( Finfo* ));
+
+ static DestFinfo concentration("concentration",
+ "Receives concentration (from DifShell).",
+ new EpFunc1<DifBufferBase, double>(&DifBufferBase::buffer));
+ static Finfo* bufferShared[] = {
+ &concentration, DifBufferBase::reactionOut()
+ };
+ static SharedFinfo buffer( "buffer",
+ "This is a shared message with DifShell. "
+ "During stage 0:\n "
+ " - DifBufferBase sends ion concentration\n"
+ " - Buffer updates buffer concentration and sends it back immediately using a call-back.\n"
+ " - DifShell updates the time-derivative ( dC / dt ) \n"
+ "\n"
+ "During stage 1: \n"
+ " - DifShell advances concentration C \n\n"
+ "This scheme means that the Buffer does not need to be scheduled, and it does its computations when "
+ "it receives a cue from the DifShell. May not be the best idea, but it saves us from doing the above "
+ "computations in 3 stages instead of 2." ,
+ bufferShared,
+ sizeof( bufferShared ) / sizeof( Finfo* ));
+
+ static DestFinfo fluxFromOut( "fluxFromOut",
+ "Destination message",
+ new EpFunc2< DifBufferBase, double, double > ( &DifBufferBase::fluxFromOut ));
+
+ static Finfo* innerDifShared[] = {
+ &fluxFromOut,
+ DifBufferBase::innerDifSourceOut()
+
+ };
+
+ static SharedFinfo innerDif( "innerDif",
+ "This shared message (and the next) is between DifBufferBases: adjoining shells exchange information to "
+ "find out the flux between them. "
+ "Using this message, an inner shell sends to, and receives from its outer shell." ,
+ innerDifShared,
+ sizeof( innerDifShared ) / sizeof( Finfo* ));
+
+ static DestFinfo fluxFromIn( "fluxFromIn", "",
+ new EpFunc2< DifBufferBase, double, double> ( &DifBufferBase::fluxFromIn) );
+
+ static Finfo* outerDifShared[] = {
+ &fluxFromIn,
+ DifBufferBase::outerDifSourceOut(),
+
+ };
+
+ static SharedFinfo outerDif( "outerDif",
+ "Using this message, an outer shell sends to, and receives from its inner shell." ,
+ outerDifShared,
+ sizeof( outerDifShared ) / sizeof( Finfo* ));
+
+ ////////////////////////////
+ // Field defs
+ ////////////////////////////
+ static ElementValueFinfo<DifBufferBase, double> activation("activation",
+ "Ion concentration from incoming conc message.",
+ &DifBufferBase::setActivation,
+ &DifBufferBase::getActivation);
+ static ElementValueFinfo<DifBufferBase, double> kf("kf",
+ "Forward rate constant of buffer molecules 1/mM/s (?)",
+ &DifBufferBase::setKf,
+ &DifBufferBase::getKf);
+ static ElementValueFinfo<DifBufferBase, double> kb("kb",
+ "Backward rate constant of buffer molecules. 1/s",
+ &DifBufferBase::setKb,
+ &DifBufferBase::getKb);
+ static ElementValueFinfo<DifBufferBase, double> D("D",
+ "Diffusion constant of buffer molecules. m^2/s",
+ &DifBufferBase::setD,
+ &DifBufferBase::getD);
+ static ElementValueFinfo<DifBufferBase, double> bFree("bFree",
+ "Free buffer concentration",
+ &DifBufferBase::setBFree,
+ &DifBufferBase::getBFree);
+ static ElementValueFinfo<DifBufferBase, double> bBound("bBound",
+ "Bound buffer concentration",
+ &DifBufferBase::setBBound,
+ &DifBufferBase::getBBound);
+ static ElementValueFinfo<DifBufferBase, double> bTot("bTot",
+ "Total buffer concentration.",
+ &DifBufferBase::setBTot,
+ &DifBufferBase::getBTot);
+ static ElementValueFinfo<DifBufferBase, double> length("length",
+ "Length of shell",
+ &DifBufferBase::setLength,
+ &DifBufferBase::getLength);
+ static ElementValueFinfo<DifBufferBase, double> diameter("diameter",
+ "Diameter of shell",
+ &DifBufferBase::setDiameter,
+ &DifBufferBase::getDiameter);
+ static ElementValueFinfo<DifBufferBase, unsigned int> shapeMode("shapeMode",
+ "shape of the shell: SHELL=0, SLICE=SLAB=1, USERDEF=3",
+ &DifBufferBase::setShapeMode,
+ &DifBufferBase::getShapeMode);
+
+ static ElementValueFinfo<DifBufferBase, double> thickness("thickness",
+ "Thickness of shell",
+ &DifBufferBase::setThickness,
+ &DifBufferBase::getThickness);
+
+ static ElementValueFinfo<DifBufferBase, double> innerArea("innerArea",
+ "Inner area of shell",
+ &DifBufferBase::setInnerArea,
+ &DifBufferBase::getInnerArea);
+ static ElementValueFinfo<DifBufferBase, double> outerArea("outerArea",
+ "Outer area of shell",
+ &DifBufferBase::setOuterArea,
+ &DifBufferBase::getOuterArea);
+ static ElementValueFinfo< DifBufferBase, double> volume( "volume", "",
+ &DifBufferBase::setVolume,
+ &DifBufferBase::getVolume );
+
+ ////
+ // DestFinfo
+ ////
+ static Finfo * difBufferFinfos[] = {
+ //////////////////////////////////////////////////////////////////
+ // Field definitions
+ //////////////////////////////////////////////////////////////////
+
+ &activation,
+ &D,
+ &bFree,
+ &bBound,
+ &bTot,
+ &kf,
+ &kb,
+ //&prevFree,
+ //&prevBound,
+ &length,
+ &diameter,
+ &shapeMode,
+ &thickness,
+ &innerArea,
+ &outerArea,
+ &volume,
+ //////////////////////////////////////////////////////////////////
+ // SharedFinfo definitions
+ /////////////////////////////////////////////////////////////////
+ &proc,
+ &buffer,
+ &innerDif,
+ &outerDif,
+ //
+ reactionOut(),
+ innerDifSourceOut(),
+ outerDifSourceOut(),
+ //////////////////////////////////////////////////////////////////
+ // DestFinfo definitions
+ //////////////////////////////////////////////////////////////////
+ &concentration,
+ };
+
+ static string doc[] = {
+ "Name", "DifBufferBase",
+ "Author", "Subhasis Ray (ported from GENESIS2)",
+ "Description", "Models diffusible buffer where total concentration is constant. It is"
+ " coupled with a DifShell.",
+ };
+ static ZeroSizeDinfo<int> dinfo;
+ static Cinfo difBufferCinfo(
+ "DifBufferBase",
+ Neutral::initCinfo(),
+ difBufferFinfos,
+ sizeof(difBufferFinfos)/sizeof(Finfo*),
+ &dinfo,
+ doc,
+ sizeof(doc)/sizeof(string));
+
+ return &difBufferCinfo;
+}
+
+static const Cinfo * difBufferCinfo = DifBufferBase::initCinfo();
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Class functions
+////////////////////////////////////////////////////////////////////////////////
+
+DifBufferBase::DifBufferBase()
+{ ; }
+
+
+double DifBufferBase::getActivation(const Eref& e) const
+{
+ return vGetActivation(e);
+}
+
+void DifBufferBase::setActivation(const Eref& e,double value)
+{
+ vSetActivation(e,value);
+}
+
+
+double DifBufferBase::getBFree(const Eref& e) const
+{
+ return vGetBFree(e);
+}
+
+void DifBufferBase::setBFree(const Eref& e,double value)
+{
+ vSetBFree(e,value);
+}
+
+double DifBufferBase::getBBound(const Eref& e) const
+{
+ return vGetBBound(e);
+}
+void DifBufferBase::setBBound(const Eref& e,double value)
+{
+ vSetBBound(e,value);
+}
+
+double DifBufferBase::getBTot(const Eref& e) const
+{
+ return vGetBTot(e);
+}
+
+void DifBufferBase::setBTot(const Eref& e,double value)
+{
+ vSetBTot(e,value);
+}
+
+
+double DifBufferBase::getKf(const Eref& e) const
+{
+ return vGetKf(e);
+}
+
+void DifBufferBase::setKf(const Eref& e,double value)
+{
+ vSetKf(e,value);
+}
+
+double DifBufferBase::getKb(const Eref& e) const
+{
+ return vGetKb(e);
+}
+
+void DifBufferBase::setKb(const Eref& e,double value)
+{
+ vSetKb(e,value);
+}
+
+double DifBufferBase::getD(const Eref& e) const
+{
+ return vGetD(e);
+}
+
+void DifBufferBase::setD(const Eref& e,double value)
+{
+ vSetD(e,value);
+}
+
+void DifBufferBase::setShapeMode(const Eref& e, unsigned int shapeMode )
+{
+ vSetShapeMode(e,shapeMode);
+}
+
+unsigned int DifBufferBase::getShapeMode(const Eref& e) const
+{
+ return vGetShapeMode(e);
+}
+
+void DifBufferBase::setLength(const Eref& e, double length )
+{
+ vSetLength(e,length);
+}
+
+double DifBufferBase::getLength(const Eref& e ) const
+{
+ return vGetLength(e);
+}
+
+void DifBufferBase::setDiameter(const Eref& e, double diameter )
+{
+ vSetDiameter(e,diameter);
+}
+
+double DifBufferBase::getDiameter(const Eref& e ) const
+{
+ return vGetDiameter(e);
+}
+
+void DifBufferBase::setThickness( const Eref& e, double thickness )
+{
+ vSetThickness(e,thickness);
+}
+
+double DifBufferBase::getThickness(const Eref& e) const
+{
+ return vGetThickness(e);
+}
+
+void DifBufferBase::setVolume(const Eref& e, double volume )
+{
+ vSetVolume(e,volume);
+}
+
+double DifBufferBase::getVolume(const Eref& e ) const
+{
+ return vGetVolume(e);
+}
+
+void DifBufferBase::setOuterArea(const Eref& e, double outerArea )
+{
+ vSetOuterArea(e,outerArea);
+}
+
+double DifBufferBase::getOuterArea(const Eref& e ) const
+{
+ return vGetOuterArea(e);
+}
+
+void DifBufferBase::setInnerArea(const Eref& e, double innerArea )
+{
+ vSetInnerArea(e,innerArea);
+}
+
+double DifBufferBase::getInnerArea(const Eref& e) const
+{
+ return vGetInnerArea(e);
+}
+
+
+
+void DifBufferBase::buffer(const Eref& e,double C)
+{
+ vBuffer(e,C);
+}
+
+void DifBufferBase::reinit( const Eref& e, ProcPtr p )
+{
+ vReinit( e, p );
+}
+
+void DifBufferBase::process( const Eref& e, ProcPtr p )
+{
+ vProcess( e, p );
+}
+void DifBufferBase:: fluxFromOut(const Eref& e,double outerC, double outerThickness )
+{
+ vFluxFromOut(e,outerC,outerThickness);
+}
+void DifBufferBase:: fluxFromIn(const Eref& e,double innerC, double innerThickness )
+{
+ vFluxFromIn(e,innerC,innerThickness);
+}
+void DifBufferBase::vSetSolver( const Eref& e, Id hsolve )
+{;}
+
+void DifBufferBase::zombify( Element* orig, const Cinfo* zClass,
+ Id hsolve )
+{
+ if ( orig->cinfo() == zClass )
+ return;
+ unsigned int start = orig->localDataStart();
+ unsigned int num = orig->numLocalData();
+ if ( num == 0 )
+ return;
+ unsigned int len = 14;
+ vector< double > data( num * len );
+
+ unsigned int j = 0;
+
+ for ( unsigned int i = 0; i < num; ++i ) {
+ Eref er( orig, i + start );
+ const DifBufferBase* ds =
+ reinterpret_cast< const DifBufferBase* >( er.data() );
+ data[j + 0] = ds->getActivation( er );
+ data[j + 1] = ds->getBFree( er );
+ data[j + 2] = ds->getBBound( er );
+ data[j + 3] = ds->getBTot( er );
+ data[j + 4] = ds->getKf( er );
+ data[j + 5] = ds->getKb( er );
+ data[j + 6] = ds->getD( er );
+ data[j + 7] = ds->getShapeMode( er );
+ data[j + 8] = ds->getLength( er );
+ data[j + 9] = ds->getDiameter( er );
+ data[j + 10] = ds->getThickness( er );
+ data[j + 11] = ds->getVolume( er );
+ data[j + 12] = ds->getOuterArea( er );
+ data[j + 13] = ds->getInnerArea( er );
+ j += len;
+ }
+ orig->zombieSwap( zClass );
+ j = 0;
+ for ( unsigned int i = 0; i < num; ++i ) {
+ Eref er( orig, i + start );
+ DifBufferBase* ds =
+ reinterpret_cast< DifBufferBase* >( er.data() );
+ ds->vSetSolver(er,hsolve);
+ ds->setActivation(er, data[j+0]);
+ ds->setBFree(er, data[j + 1]);
+ ds->setBBound(er, data[j + 2]);
+ ds->setBTot(er, data[j + 3]);
+ ds->setKf(er, data[j + 4]);
+ ds->setKb(er, data[j + 5]);
+ ds->setD(er, data[j + 6]);
+ ds->setShapeMode(er, data[j + 7]);
+ ds->setLength(er, data[j + 8]);
+ ds->setDiameter(er, data[j + 9]);
+ ds->setThickness(er, data[j + 10]);
+ ds->setVolume(er, data[j + 11]);
+ ds->setOuterArea(er, data[j + 12]);
+ ds->setInnerArea(er, data[j + 13]);
+ j += len; //??
+ }
+
+}
diff --git a/moose-core/biophysics/DifBufferBase.h b/moose-core/biophysics/DifBufferBase.h
new file mode 100644
index 0000000000000000000000000000000000000000..e1bc23b616caa153d9cde8affe5b8808fb0fa1d4
--- /dev/null
+++ b/moose-core/biophysics/DifBufferBase.h
@@ -0,0 +1,134 @@
+/**********************************************************************
+ ** This program is part of 'MOOSE', the
+ ** Multiscale Object Oriented Simulation Environment.
+ ** copyright (C) 2003-2008
+ ** Upinder S. Bhalla, Niraj Dudani and NCBS
+ ** It is made available under the terms of the
+ ** GNU Lesser General Public License version 2.1
+ ** See the file COPYING.LIB for the full notice.
+ **********************************************************************/
+
+#ifndef _DIFBUFFER_BASE_H
+#define _DIFBUFFER_BASE_H
+/*This is base class for DifBuffer*/
+class DifBufferBase
+{
+public:
+ DifBufferBase();
+
+ void buffer(const Eref& e,double C);
+ void reinit( const Eref & e, ProcPtr p );
+ void process(const Eref & e, ProcPtr p );
+ void fluxFromOut(const Eref& e,double outerC, double outerThickness );
+ void fluxFromIn( const Eref& e,double innerC, double innerThickness );
+ virtual void vBuffer(const Eref& e,double C) = 0;
+ virtual void vReinit( const Eref & e, ProcPtr p ) = 0;
+ virtual void vProcess(const Eref & e, ProcPtr p ) = 0;
+ virtual void vFluxFromOut(const Eref& e,double outerC, double outerThickness ) = 0;
+ virtual void vFluxFromIn( const Eref& e,double innerC, double innerThickness ) = 0;
+
+
+ double getActivation(const Eref& e) const;
+ void setActivation(const Eref& e,double value);
+
+ double getBFree(const Eref& e) const;
+ void setBFree(const Eref& e,double value);
+
+ double getBBound(const Eref& e) const;
+ void setBBound(const Eref& e,double value);
+
+ double getBTot(const Eref& e) const; // total buffer concentration in mM (free + bound)
+ void setBTot(const Eref& e,double value);
+
+ double getKf(const Eref& e) const; // forward rate constant in 1/(mM*sec)
+ void setKf(const Eref& e,double value);
+
+ double getKb(const Eref& e) const; // backward rate constant in 1/sec
+ void setKb(const Eref& e,double value);
+
+ double getD(const Eref& e) const; // diffusion constant of buffer molecules, m^2/sec
+ void setD(const Eref& e,double value);
+
+
+ unsigned int getShapeMode(const Eref& e) const;
+ void setShapeMode(const Eref& e,unsigned int value); // variables SHELL=0, SLICE=SLAB=1, USERDEF=3.
+
+ double getLength(const Eref& e) const; // shell length
+ void setLength(const Eref& e,double value);
+
+ double getDiameter(const Eref& e) const; // shell diameter
+ void setDiameter(const Eref& e,double value);
+
+ double getThickness(const Eref& e) const; // shell thickness
+ void setThickness(const Eref& e,double value);
+
+ void setOuterArea( const Eref& e,double outerArea );
+ double getOuterArea(const Eref& e) const; // area of upper (outer) shell surface
+
+ void setInnerArea( const Eref& e,double innerArea );
+ double getInnerArea(const Eref& e) const; // area of lower (inner) shell surface
+
+ double getVolume(const Eref& e) const; // shell volume
+ void setVolume(const Eref& e,double volume); //
+
+
+ virtual double vGetActivation(const Eref& e) const = 0;
+ virtual void vSetActivation(const Eref& e,double value) = 0;
+
+ virtual double vGetBFree(const Eref& e) const = 0;
+ virtual void vSetBFree(const Eref& e,double value) = 0;
+
+ virtual double vGetBBound(const Eref& e) const = 0;
+ virtual void vSetBBound(const Eref& e,double value) = 0;
+
+
+ virtual double vGetBTot(const Eref& e) const = 0; // total buffer concentration in mM (free + bound)
+ virtual void vSetBTot(const Eref& e,double value) = 0;
+
+ virtual double vGetKf(const Eref& e) const = 0; // forward rate constant in 1/(mM*sec)
+ virtual void vSetKf(const Eref& e,double value) = 0;
+
+ virtual double vGetKb(const Eref& e) const = 0; // backward rate constant in 1/sec
+ virtual void vSetKb(const Eref& e,double value) = 0;
+
+ virtual double vGetD(const Eref& e) const = 0; // diffusion constant of buffer molecules, m^2/sec
+ virtual void vSetD(const Eref& e,double value) = 0;
+
+ virtual void vSetShapeMode(const Eref& e, unsigned int shapeMode ) = 0;
+ virtual unsigned int vGetShapeMode(const Eref& e) const = 0;
+
+ virtual void vSetLength(const Eref& e, double length ) = 0;
+ virtual double vGetLength(const Eref& e) const = 0;
+
+ virtual void vSetDiameter(const Eref& e, double diameter ) = 0;
+ virtual double vGetDiameter(const Eref& e) const = 0;
+
+ virtual void vSetThickness(const Eref& e, double thickness ) = 0;
+ virtual double vGetThickness(const Eref& e) const = 0;
+
+ virtual void vSetVolume(const Eref& e, double volume ) = 0;
+ virtual double vGetVolume(const Eref& e) const = 0;
+
+ virtual void vSetOuterArea(const Eref& e, double outerArea ) = 0;
+ virtual double vGetOuterArea(const Eref& e) const = 0;
+
+ virtual void vSetInnerArea(const Eref& e, double innerArea ) = 0;
+ virtual double vGetInnerArea(const Eref& e) const = 0;
+
+ static SrcFinfo4< double, double, double, double >* reactionOut();
+ static SrcFinfo2< double, double >* innerDifSourceOut();
+ static SrcFinfo2< double, double >* outerDifSourceOut();
+ static const Cinfo * initCinfo();
+ virtual void vSetSolver( const Eref& e, Id hsolve );
+ static void zombify( Element* orig, const Cinfo* zClass,
+ Id hsolve );
+private:
+
+
+
+};
+
+
+
+
+#endif //_DIFBUFFER_BASE_H
diff --git a/moose-core/biophysics/DifShell.cpp b/moose-core/biophysics/DifShell.cpp
index 98fe9a2b1c9330390f9656288a0263de7207f9de..353683251aad3853eb9a4cc84e5768db6326e971 100644
--- a/moose-core/biophysics/DifShell.cpp
+++ b/moose-core/biophysics/DifShell.cpp
@@ -1,253 +1,49 @@
/**********************************************************************
-** This program is part of 'MOOSE', the
-** Multiscale Object Oriented Simulation Environment.
-** copyright (C) 2003-2008
-** Upinder S. Bhalla, Niraj Dudani and NCBS
-** It is made available under the terms of the
-** GNU Lesser General Public License version 2.1
-** See the file COPYING.LIB for the full notice.
-**********************************************************************/
+ ** This program is part of 'MOOSE', the
+ ** Multiscale Object Oriented Simulation Environment.
+ ** copyright (C) 2003-2008
+ ** Upinder S. Bhalla, Niraj Dudani and NCBS
+ ** It is made available under the terms of the
+ ** GNU Lesser General Public License version 2.1
+ ** See the file COPYING.LIB for the full notice.
+ **********************************************************************/
#include "header.h"
+#include "DifShellBase.h"
#include "DifShell.h"
-#include "../utility/numutil.h"
-#include <cmath>
-static SrcFinfo1< double >* concentrationOut()
-{
- static SrcFinfo1< double > concentrationOut("concentrationOut",
- "Sends out concentration");
- return &concentrationOut;
-}
-static SrcFinfo2< double, double >* innerDifSourceOut(){
- static SrcFinfo2< double, double > sourceOut("innerDifSourceOut",
- "Sends out source information.");
- return & sourceOut;
-}
-
-static SrcFinfo2< double, double >* outerDifSourceOut(){
- static SrcFinfo2< double, double > sourceOut("outerDifSourceOut",
- "Sends out source information.");
- return & sourceOut;
-}
+const double DifShell::EPSILON = 1.0e-10;
+const double DifShell::F = 96485.3415; /* C / mol like in genesis */
const Cinfo* DifShell::initCinfo()
{
- static DestFinfo process( "process",
- "Handles process call",
- new ProcOpFunc< DifShell >( &DifShell::process_0 ) );
- static DestFinfo reinit( "reinit",
- "Reinit happens only in stage 0",
- new ProcOpFunc< DifShell >( &DifShell::reinit_0 ));
-
- static Finfo* processShared_0[] = {
- &process,
- &reinit
- };
-
- static SharedFinfo process_0(
- "process_0",
- "Here we create 2 shared finfos to attach with the Ticks. This is because we want to perform DifShell "
- "computations in 2 stages, much as in the Compartment object. "
- "In the first stage we send out the concentration value to other DifShells and Buffer elements. We also "
- "receive fluxes and currents and sum them up to compute ( dC / dt ). "
- "In the second stage we find the new C value using an explicit integration method. "
- "This 2-stage procedure eliminates the need to store and send prev_C values, as was common in GENESIS.",
- processShared_0,
- sizeof( processShared_0 ) / sizeof( Finfo* ));
-
- static DestFinfo process1( "process",
- "Handle process call",
- new ProcOpFunc< DifShell >( &DifShell::process_1 ) );
- static DestFinfo reinit1( "reinit",
- "Reinit happens only in stage 0",
- new ProcOpFunc< DifShell >( &DifShell::reinit_1)
- );
- static Finfo* processShared_1[] = {
- &process1, &reinit1
- };
-
- static SharedFinfo process_1( "process_1",
- "Second process call",
- processShared_1,
- sizeof( processShared_1 ) / sizeof( Finfo* ) );
-
-
- static DestFinfo reaction( "reaction",
- "Here the DifShell receives reaction rates (forward and backward), and concentrations for the "
- "free-buffer and bound-buffer molecules.",
- new OpFunc4< DifShell, double, double, double, double >( &DifShell::buffer ));
- static Finfo* bufferShared[] = {
- concentrationOut(), &reaction
+
+ static string doc[] =
+ {
+ "Name", "DifShell",
+ "Author", "Niraj Dudani. Ported to async13 by Subhasis Ray. Rewritten by Asia Jedrzejewska-Szmek",
+ "Description", "DifShell object: Models diffusion of an ion (typically calcium) within an "
+ "electric compartment. A DifShell is an iso-concentration region with respect to "
+ "the ion. Adjoining DifShells exchange flux of this ion, and also keep track of "
+ "changes in concentration due to pumping, buffering and channel currents, by "
+ "talking to the appropriate objects.",
};
- static SharedFinfo buffer( "buffer",
- "This is a shared message from a DifShell to a Buffer (FixBuffer or DifBuffer). "
- "During stage 0:\n "
- " - DifShell sends ion concentration\n"
- " - Buffer updates buffer concentration and sends it back immediately using a call-back.\n"
- " - DifShell updates the time-derivative ( dC / dt ) \n"
- "\n"
- "During stage 1: \n"
- " - DifShell advances concentration C \n\n"
- "This scheme means that the Buffer does not need to be scheduled, and it does its computations when "
- "it receives a cue from the DifShell. May not be the best idea, but it saves us from doing the above "
- "computations in 3 stages instead of 2." ,
- bufferShared,
- sizeof( bufferShared ) / sizeof( Finfo* ));
-
- static DestFinfo fluxFromOut( "fluxFromOut",
- "Destination message",
- new OpFunc2< DifShell, double, double > ( &DifShell::fluxFromOut ));
-
- static Finfo* innerDifShared[] = {
- innerDifSourceOut(),
- &fluxFromOut
- };
- static SharedFinfo innerDif( "innerDif",
- "This shared message (and the next) is between DifShells: adjoining shells exchange information to "
- "find out the flux between them. "
- "Using this message, an inner shell sends to, and receives from its outer shell." ,
- innerDifShared,
- sizeof( innerDifShared ) / sizeof( Finfo* ));
-
- static DestFinfo fluxFromIn( "fluxFromIn", "",
- new OpFunc2< DifShell, double, double> ( &DifShell::fluxFromIn ) );
- static Finfo* outerDifShared[] = {
- &fluxFromIn,
- outerDifSourceOut(),
- };
-
- static SharedFinfo outerDif( "outerDif",
- "Using this message, an outer shell sends to, and receives from its inner shell." ,
- outerDifShared,
- sizeof( outerDifShared ) / sizeof( Finfo* ));
-
- static ReadOnlyValueFinfo< DifShell, double> C( "C",
- "Concentration C is computed by the DifShell and is read-only",
- &DifShell::getC);
- static ValueFinfo< DifShell, double> Ceq( "Ceq", "",
- &DifShell::setCeq,
- &DifShell::getCeq);
- static ValueFinfo< DifShell, double> D( "D", "",
- &DifShell::setD,
- &DifShell::getD);
- static ValueFinfo< DifShell, double> valence( "valence", "",
- &DifShell::setValence,
- &DifShell::getValence);
- static ValueFinfo< DifShell, double> leak( "leak", "",
- &DifShell::setLeak,
- &DifShell::getLeak);
- static ValueFinfo< DifShell, unsigned int> shapeMode( "shapeMode", "",
- &DifShell::setShapeMode,
- &DifShell::getShapeMode);
- static ValueFinfo< DifShell, double> length( "length", "",
- &DifShell::setLength,
- &DifShell::getLength);
- static ValueFinfo< DifShell, double> diameter( "diameter", "",
- &DifShell::setDiameter,
- &DifShell::getDiameter );
- static ValueFinfo< DifShell, double> thickness( "thickness", "",
- &DifShell::setThickness,
- &DifShell::getThickness );
- static ValueFinfo< DifShell, double> volume( "volume", "",
- &DifShell::setVolume,
- &DifShell::getVolume );
- static ValueFinfo< DifShell, double> outerArea( "outerArea", "",
- &DifShell::setOuterArea,
- &DifShell::getOuterArea);
- static ValueFinfo< DifShell, double> innerArea( "innerArea", "",
- &DifShell::setInnerArea,
- &DifShell::getInnerArea );
+ static Dinfo< DifShell > dinfo;
+ static Cinfo difShellCinfo(
+ "DifShell",
+ DifShellBase::initCinfo(),
+ 0,
+ 0,
+ &dinfo,
+ doc,
+ sizeof( doc ) / sizeof( string ));
-
- static DestFinfo influx( "influx", "",
- new OpFunc1< DifShell, double > (&DifShell::influx ));
- static DestFinfo outflux( "outflux", "",
- new OpFunc1< DifShell, double >( &DifShell::influx ));
- static DestFinfo fInflux( "fInflux", "",
- new OpFunc2< DifShell, double, double >( &DifShell::fInflux ));
- static DestFinfo fOutflux( "fOutflux", "",
- new OpFunc2< DifShell, double, double> (&DifShell::fOutflux ));
- static DestFinfo storeInflux( "storeInflux", "",
- new OpFunc1< DifShell, double >( &DifShell::storeInflux ) );
- static DestFinfo storeOutflux( "storeOutflux", "",
- new OpFunc1< DifShell, double > ( &DifShell::storeOutflux ) );
- static DestFinfo tauPump( "tauPump","",
- new OpFunc2< DifShell, double, double >( &DifShell::tauPump ) );
- static DestFinfo eqTauPump( "eqTauPump", "",
- new OpFunc1< DifShell, double >( &DifShell::eqTauPump ) );
- static DestFinfo mmPump( "mmPump", "",
- new OpFunc2< DifShell, double, double >( &DifShell::mmPump ) );
- static DestFinfo hillPump( "hillPump", "",
- new OpFunc3< DifShell, double, double, unsigned int >( &DifShell::hillPump ) );
- static Finfo* difShellFinfos[] = {
- //////////////////////////////////////////////////////////////////
- // Field definitions
- //////////////////////////////////////////////////////////////////
- &C,
- &Ceq,
- &D,
- &valence,
- &leak,
- &shapeMode,
- &length,
- &diameter,
- &thickness,
- &volume,
- &outerArea,
- &innerArea,
- //////////////////////////////////////////////////////////////////
- // MsgSrc definitions
- //////////////////////////////////////////////////////////////////
-
- //////////////////////////////////////////////////////////////////
- // SharedFinfo definitions
- //////////////////////////////////////////////////////////////////
- &process_0,
- &process_1,
- &buffer,
- &innerDif,
- &outerDif,
- //////////////////////////////////////////////////////////////////
- // DestFinfo definitions
- //////////////////////////////////////////////////////////////////
- &influx,
- &outflux,
- &fInflux,
- &fOutflux,
- &storeInflux,
- &storeOutflux,
- &tauPump,
- &eqTauPump,
- &mmPump,
- &hillPump,
- };
-
- static string doc[] =
- {
- "Name", "DifShell",
- "Author", "Niraj Dudani. Ported to async13 by Subhasis Ray.",
- "Description", "DifShell object: Models diffusion of an ion (typically calcium) within an "
- "electric compartment. A DifShell is an iso-concentration region with respect to "
- "the ion. Adjoining DifShells exchange flux of this ion, and also keep track of "
- "changes in concentration due to pumping, buffering and channel currents, by "
- "talking to the appropriate objects.",
- };
- static Dinfo< DifShell > dinfo;
- static Cinfo difShellCinfo(
- "DifShell",
- Neutral::initCinfo(),
- difShellFinfos,
- sizeof( difShellFinfos ) / sizeof( Finfo* ),
- &dinfo,
- doc,
- sizeof( doc ) / sizeof( string ));
-
- return &difShellCinfo;
+ return &difShellCinfo;
}
+//Cinfo *object* corresponding to the class.
static const Cinfo* difShellCinfo = DifShell::initCinfo();
////////////////////////////////////////////////////////////////////////////////
@@ -255,499 +51,432 @@ static const Cinfo* difShellCinfo = DifShell::initCinfo();
////////////////////////////////////////////////////////////////////////////////
/// Faraday's constant (Coulomb / Mole)
-const double DifShell::F_ = 0.0;
+
DifShell::DifShell() :
- dCbyDt_( 0.0 ),
- C_( 0.0 ),
- Ceq_( 0.0 ),
- D_( 0.0 ),
- valence_( 0.0 ),
- leak_( 0.0 ),
- shapeMode_( 0 ),
- length_( 0.0 ),
- diameter_( 0.0 ),
- thickness_( 0.0 ),
- volume_( 0.0 ),
- outerArea_( 0.0 ),
- innerArea_( 0.0 )
+ dCbyDt_( 0.0 ),
+ Cmultiplier_(0.0),
+ C_( 0.0 ),
+ Ceq_( 0.0 ),
+ prevC_(0.0),
+ D_( 0.0 ),
+ valence_( 0.0 ),
+ leak_( 0.0 ),
+ shapeMode_(0),
+ length_(0),
+ diameter_(0),
+ thickness_(0),
+ volume_(0),
+ outerArea_(0),
+ innerArea_(0)
{ ; }
////////////////////////////////////////////////////////////////////////////////
// Field access functions
////////////////////////////////////////////////////////////////////////////////
/// C is a read-only field
-double DifShell::getC() const
+
+
+void DifShell::vSetC( const Eref& e, double C)
+{
+ if ( C < 0.0 ) {
+ cerr << "Error: DifShell: C cannot be negative!\n";
+ return;
+ }
+
+ C_ = C;
+ prevC_ = C_;
+}
+double DifShell::vGetC(const Eref& e) const
{
- return C_;
+ return C_;
}
-void DifShell::setCeq( double Ceq )
+void DifShell::vSetCeq( const Eref& e, double Ceq )
{
- if ( Ceq < 0.0 ) {
- cerr << "Error: DifShell: Ceq cannot be negative!\n";
- return;
- }
+ if ( Ceq < 0.0 ) {
+ cerr << "Error: DifShell: Ceq cannot be negative!\n";
+ return;
+ }
- Ceq_ = Ceq;
+ Ceq_ = Ceq;
}
-double DifShell::getCeq( ) const
+double DifShell::vGetCeq(const Eref& e ) const
{
- return Ceq_;
+ return Ceq_;
}
-void DifShell::setD( double D )
+void DifShell::vSetD(const Eref& e, double D )
{
- if ( D < 0.0 ) {
- cerr << "Error: DifShell: D cannot be negative!\n";
- return;
- }
+ if ( D < 0.0 ) {
+ cerr << "Error: DifShell: D cannot be negative!\n";
+ return;
+ }
- D_ = D;
+ D_ = D;
}
-double DifShell::getD( ) const
+double DifShell::vGetD(const Eref& e ) const
{
- return D_;
+ return D_;
}
-void DifShell::setValence( double valence )
+void DifShell::vSetValence(const Eref& e, double valence )
{
- if ( valence < 0.0 ) {
- cerr << "Error: DifShell: valence cannot be negative!\n";
- return;
- }
+ if ( valence < 0.0 ) {
+ cerr << "Error: DifShell: valence cannot be negative!\n";
+ return;
+ }
- valence_ = valence;
+ valence_ = valence;
}
-double DifShell::getValence( ) const
+double DifShell::vGetValence(const Eref& e ) const
{
- return valence_;
+ return valence_;
}
-void DifShell::setLeak( double leak )
+void DifShell::vSetLeak(const Eref& e, double leak )
{
- leak_ = leak;
+ leak_ = leak;
}
-double DifShell::getLeak( ) const
+double DifShell::vGetLeak(const Eref& e ) const
{
- return leak_;
+ return leak_;
}
-void DifShell::setShapeMode( unsigned int shapeMode )
+
+void DifShell::vSetShapeMode(const Eref& e, unsigned int shapeMode )
{
- if ( shapeMode != 0 && shapeMode != 1 && shapeMode != 3 ) {
- cerr << "Error: DifShell: I only understand shapeModes 0, 1 and 3.\n";
- return;
- }
- shapeMode_ = shapeMode;
+ if ( shapeMode != 0 && shapeMode != 1 && shapeMode != 3 ) {
+ cerr << "Error: DifShell: I only understand shapeModes 0, 1 and 3.\n";
+ return;
+ }
+ shapeMode_ = shapeMode;
}
-unsigned int DifShell::getShapeMode() const
+unsigned int DifShell::vGetShapeMode(const Eref& e) const
{
- return shapeMode_;
+ return shapeMode_;
}
-void DifShell::setLength( double length )
+void DifShell::vSetLength(const Eref& e, double length )
{
- if ( length < 0.0 ) {
- cerr << "Error: DifShell: length cannot be negative!\n";
- return;
- }
+ if ( length < 0.0 ) {
+ cerr << "Error: DifShell: length cannot be negative!\n";
+ return;
+ }
- length_ = length;
+ length_ = length;
}
-double DifShell::getLength( ) const
+double DifShell::vGetLength(const Eref& e ) const
{
- return length_;
+ return length_;
}
-void DifShell::setDiameter( double diameter )
+void DifShell::vSetDiameter(const Eref& e, double diameter )
{
- if ( diameter < 0.0 ) {
- cerr << "Error: DifShell: diameter cannot be negative!\n";
- return;
- }
+ if ( diameter < 0.0 ) {
+ cerr << "Error: DifShell: diameter cannot be negative!\n";
+ return;
+ }
- diameter_ = diameter;
+ diameter_ = diameter;
}
-double DifShell::getDiameter( ) const
+double DifShell::vGetDiameter(const Eref& e ) const
{
- return diameter_;
+ return diameter_;
}
-void DifShell::setThickness( double thickness )
+void DifShell::vSetThickness( const Eref& e, double thickness )
{
- if ( thickness < 0.0 ) {
- cerr << "Error: DifShell: thickness cannot be negative!\n";
- return;
- }
+ if ( thickness < 0.0 ) {
+ cerr << "Error: DifShell: thickness cannot be negative!\n";
+ return;
+ }
- thickness_ = thickness;
+ thickness_ = thickness;
}
-double DifShell::getThickness() const
+double DifShell::vGetThickness(const Eref& e) const
{
- return thickness_;
+ return thickness_;
}
-void DifShell::setVolume( double volume )
+void DifShell::vSetVolume(const Eref& e, double volume )
{
- if ( shapeMode_ != 3 )
- cerr << "Warning: DifShell: Trying to set volume, when shapeMode is not USER-DEFINED\n";
+ if ( shapeMode_ != 3 )
+ cerr << "Warning: DifShell: Trying to set volume, when shapeMode is not USER-DEFINED\n";
- if ( volume < 0.0 ) {
- cerr << "Error: DifShell: volume cannot be negative!\n";
- return;
- }
+ if ( volume < 0.0 ) {
+ cerr << "Error: DifShell: volume cannot be negative!\n";
+ return;
+ }
- volume_ = volume;
+ volume_ = volume;
}
-double DifShell::getVolume( ) const
+double DifShell::vGetVolume(const Eref& e ) const
{
- return volume_;
+ return volume_;
}
-void DifShell::setOuterArea( double outerArea )
+void DifShell::vSetOuterArea(const Eref& e, double outerArea )
{
- if (shapeMode_ != 3 )
- cerr << "Warning: DifShell: Trying to set outerArea, when shapeMode is not USER-DEFINED\n";
+ if (shapeMode_ != 3 )
+ cerr << "Warning: DifShell: Trying to set outerArea, when shapeMode is not USER-DEFINED\n";
- if ( outerArea < 0.0 ) {
- cerr << "Error: DifShell: outerArea cannot be negative!\n";
- return;
- }
+ if ( outerArea < 0.0 ) {
+ cerr << "Error: DifShell: outerArea cannot be negative!\n";
+ return;
+ }
- outerArea_ = outerArea;
+ outerArea_ = outerArea;
}
-double DifShell::getOuterArea( ) const
+double DifShell::vGetOuterArea(const Eref& e ) const
{
- return outerArea_;
+ return outerArea_;
}
-void DifShell::setInnerArea( double innerArea )
+void DifShell::vSetInnerArea(const Eref& e, double innerArea )
{
- if ( shapeMode_ != 3 )
- cerr << "Warning: DifShell: Trying to set innerArea, when shapeMode is not USER-DEFINED\n";
+ if ( shapeMode_ != 3 )
+ cerr << "Warning: DifShell: Trying to set innerArea, when shapeMode is not USER-DEFINED\n";
- if ( innerArea < 0.0 ) {
- cerr << "Error: DifShell: innerArea cannot be negative!\n";
- return;
- }
+ if ( innerArea < 0.0 ) {
+ cerr << "Error: DifShell: innerArea cannot be negative!\n";
+ return;
+ }
- innerArea_ = innerArea;
-}
-
-double DifShell::getInnerArea() const
-{
- return innerArea_;
-}
-
-////////////////////////////////////////////////////////////////////////////////
-// Dest functions
-////////////////////////////////////////////////////////////////////////////////
-
-void DifShell::reinit_0( const Eref& e, ProcPtr p )
-{
- localReinit_0( e, p );
-}
-
-void DifShell::process_0( const Eref& e, ProcPtr p )
-{
- localProcess_0( e, p );
-}
-
-void DifShell::process_1(const Eref& e, ProcPtr p )
-{
- localProcess_1( e, p );
-}
-
-void DifShell::reinit_1(const Eref& e, ProcPtr p )
-{
- ;
-}
-
-void DifShell::buffer(
- double kf,
- double kb,
- double bFree,
- double bBound )
-{
- localBuffer( kf, kb, bFree, bBound );
-}
-
-void DifShell::fluxFromOut(
- double outerC,
- double outerThickness )
-{
- localFluxFromOut( outerC, outerThickness );
-}
-
-void DifShell::fluxFromIn(
- double innerC,
- double innerThickness )
-{
- localFluxFromIn( innerC, innerThickness );
-}
-
-void DifShell::influx(
- double I )
-{
- localInflux( I );
+ innerArea_ = innerArea;
}
-void DifShell::outflux(
- double I )
+double DifShell::vGetInnerArea(const Eref& e) const
{
- localOutflux( I );
+ return innerArea_;
}
-void DifShell::fInflux(
- double I,
- double fraction )
-{
-localFInflux( I, fraction );
-}
-void DifShell::fOutflux(
- double I,
- double fraction )
-{
-localFOutflux( I, fraction );
-}
-
-void DifShell::storeInflux(
- double flux )
-{
-localStoreInflux( flux );
-}
-
-void DifShell::storeOutflux(
- double flux )
-{
-localStoreOutflux( flux );
-}
-
-void DifShell::tauPump(
- double kP,
- double Ceq )
-{
-localTauPump( kP, Ceq );
-}
-
-void DifShell::eqTauPump(
- double kP )
-{
-localEqTauPump( kP );
-}
-
-void DifShell::mmPump(
- double vMax,
- double Kd )
-{
-localMMPump( vMax, Kd );
-}
-
-void DifShell::hillPump(
- double vMax,
- double Kd,
- unsigned int hill )
-{
-localHillPump( vMax, Kd, hill );
-}
////////////////////////////////////////////////////////////////////////////////
// Local dest functions
////////////////////////////////////////////////////////////////////////////////
-
-void DifShell::localReinit_0( const Eref& e, ProcPtr p )
+double DifShell::integrate( double state, double dt, double A, double B )
{
- dCbyDt_ = leak_;
-
- double Pi = M_PI;
- double dOut = diameter_;
- double dIn = diameter_ - thickness_;
-
- switch ( shapeMode_ )
- {
- /*
- * Onion Shell
- */
- case 0:
- if ( length_ == 0.0 ) { // Spherical shell
- volume_ = ( Pi / 6.0 ) * ( dOut * dOut * dOut - dIn * dIn * dIn );
- outerArea_ = Pi * dOut * dOut;
- innerArea_ = Pi * dIn * dIn;
- } else { // Cylindrical shell
- volume_ = ( Pi * length_ / 4.0 ) * ( dOut * dOut - dIn * dIn );
- outerArea_ = Pi * dOut * length_;
- innerArea_ = Pi * dIn * length_;
- }
+ if ( B > EPSILON ) {
+ double x = exp( -B * dt );
+ return state * x + ( A / B ) * ( 1 - x );
+ }
+
+ return state + A * dt ;
+}
+
+void DifShell::vReinit( const Eref& e, ProcPtr p )
+{
+ dCbyDt_ = leak_;
+ Cmultiplier_ = 0;
+
+ double rOut = diameter_/2.;
+
+ double rIn = rOut - thickness_;
+
+ if (rIn <0)
+ rIn = 0.;
+
+ switch ( shapeMode_ )
+ {
+ /*
+ * Onion Shell
+ */
+ case 0:
+ if ( length_ == 0.0 ) { // Spherical shell
+ volume_ = 4./3.* M_PI * ( rOut * rOut * rOut - rIn * rIn * rIn );
+ outerArea_ = 4*M_PI * rOut * rOut;
+ innerArea_ = 4*M_PI * rIn * rIn;
+ } else { // Cylindrical shell
+ volume_ = ( M_PI * length_ ) * ( rOut * rOut - rIn * rIn );
+ outerArea_ = 2*M_PI * rOut * length_;
+ innerArea_ = 2*M_PI * rIn * length_;
+ }
- break;
+ break;
- /*
- * Cylindrical Slice
- */
- case 1:
- volume_ = Pi * diameter_ * diameter_ * thickness_ / 4.0;
- outerArea_ = Pi * diameter_ * diameter_ / 4.0;
- innerArea_ = outerArea_;
- break;
+ /*
+ * Cylindrical Slice
+ */
+ case 1:
+ volume_ = M_PI * diameter_ * diameter_ * thickness_ / 4.0;
+ outerArea_ = M_PI * diameter_ * diameter_ / 4.0;
+ innerArea_ = outerArea_;
+ break;
- /*
- * User defined
- */
- case 3:
- // Nothing to be done here. Volume and inner-, outer areas specified by
- // user.
- break;
+ /*
+ * User defined
+ */
+ case 3:
+ // Nothing to be done here. Volume and inner-, outer areas specified by
+ // user.
+ break;
- default:
- assert( 0 );
- }
+ default:
+ assert( 0 );
+ }
+ C_= Ceq_;
+ prevC_ = Ceq_;
+ concentrationOut()->send( e, C_ );
+ innerDifSourceOut()->send( e, prevC_, thickness_ );
+ outerDifSourceOut()->send( e, prevC_, thickness_ );
}
-void DifShell::localProcess_0( const Eref & e, ProcPtr p )
+void DifShell::vProcess( const Eref & e, ProcPtr p )
{
- /**
- * Send ion concentration and thickness to adjacent DifShells. They will
- * then compute their incoming fluxes.
- */
- innerDifSourceOut()->send( e, C_, thickness_ );
- outerDifSourceOut()->send( e, C_, thickness_ );
-
- /**
- * Send ion concentration to ion buffers. They will send back information on
- * the reaction (forward / backward rates ; free / bound buffer concentration)
- * immediately, which this DifShell will use to find amount of ion captured
- * or released in the current time-step.
- */
- concentrationOut()->send( e, C_ );
-}
+ /**
+ * Send ion concentration and thickness to adjacent DifShells. They will
+ * then compute their incoming fluxes.
+ */
+
+
+ C_ = integrate(C_,p->dt,dCbyDt_,Cmultiplier_);
+
+ /**
+ * Send ion concentration to ion buffers. They will send back information on
+ * the reaction (forward / backward rates ; free / bound buffer concentration)
+ * immediately, which this DifShell will use to find amount of ion captured
+ * or released in the current time-step.
+ */
+ prevC_ = C_;
-void DifShell::localProcess_1( const Eref& e, ProcPtr p )
-{
- C_ += dCbyDt_ * p->dt;
- dCbyDt_ = leak_;
-}
+ //cout<<"Shell "<< C_<<" ";
+ dCbyDt_ = leak_;
+ Cmultiplier_ = 0;
+ innerDifSourceOut()->send( e, prevC_, thickness_ );
+ outerDifSourceOut()->send( e, prevC_, thickness_ );
+ concentrationOut()->send( e, C_ );
-void DifShell::localBuffer(
- double kf,
- double kb,
- double bFree,
- double bBound )
+}
+void DifShell::vBuffer(const Eref& e,
+ double kf,
+ double kb,
+ double bFree,
+ double bBound )
{
- dCbyDt_ += -kf * bFree * C_ + kb * bBound;
+ dCbyDt_ += kb * bBound;
+ Cmultiplier_ += kf * bFree;
}
-void DifShell::localFluxFromOut( double outerC, double outerThickness )
+void DifShell::vFluxFromOut(const Eref& e, double outerC, double outerThickness )
{
- double dx = ( outerThickness + thickness_ ) / 2.0;
-
- /**
- * We could pre-compute ( D / Volume ), but let us leave the optimizations
- * for the solver.
- */
- dCbyDt_ += ( D_ / volume_ ) * ( outerArea_ / dx ) * ( outerC - C_ );
+ double diff =2.* D_ /volume_ * outerArea_ / (outerThickness + thickness_) ;
+ //influx from outer shell
+ /**
+ * We could pre-compute ( D / Volume ), but let us leave the optimizations
+ * for the solver.
+ */
+
+ dCbyDt_ += diff * outerC;
+ Cmultiplier_ += diff ;
}
-void DifShell::localFluxFromIn( double innerC, double innerThickness )
+void DifShell::vFluxFromIn(const Eref& e, double innerC, double innerThickness )
{
- double dx = ( innerThickness + thickness_ ) / 2.0;
-
- dCbyDt_ += ( D_ / volume_ ) * ( innerArea_ / dx ) * ( innerC - C_ );
+ //influx from inner shell
+ //double dx = ( innerThickness + thickness_ ) / 2.0;
+ double diff = 2.* D_/volume_ * innerArea_ / (innerThickness + thickness_);
+ dCbyDt_ += diff * innerC ;
+ Cmultiplier_ += diff ;
}
-void DifShell::localInflux( double I )
+void DifShell::vInflux(const Eref& e, double I )
{
- /**
- * I: Amperes
- * F_: Faraday's constant: Coulomb / mole
- * valence_: charge on ion: dimensionless
- */
- dCbyDt_ += I / ( F_ * valence_ * volume_ );
+ /**
+ * I: Amperes
+ * F_: Faraday's constant: Coulomb / mole
+ * valence_: charge on ion: dimensionless
+ */
+ dCbyDt_ += I / ( F * valence_ * volume_ );
}
/**
* Same as influx, except subtracting.
*/
-void DifShell::localOutflux( double I )
+void DifShell::vOutflux(const Eref& e, double I )
{
- dCbyDt_ -= I / ( F_ * valence_ * volume_ );
+ dCbyDt_ -= I / ( F * valence_ * volume_ );
}
-void DifShell::localFInflux( double I, double fraction )
+void DifShell::vFInflux(const Eref& e, double I, double fraction )
{
- dCbyDt_ += fraction * I / ( F_ * valence_ * volume_ );
+ dCbyDt_ += fraction * I / ( F * valence_ * volume_ );
}
-void DifShell::localFOutflux( double I, double fraction )
+void DifShell::vFOutflux(const Eref& e, double I, double fraction )
{
- dCbyDt_ -= fraction * I / ( F_ * valence_ * volume_ );
+ dCbyDt_ -= fraction * I / ( F * valence_ * volume_ );
}
-void DifShell::localStoreInflux( double flux )
+void DifShell::vStoreInflux(const Eref& e, double flux )
{
- dCbyDt_ += flux / volume_;
+ dCbyDt_ += flux / volume_;
}
-void DifShell::localStoreOutflux( double flux )
+void DifShell::vStoreOutflux(const Eref& e, double flux )
{
- dCbyDt_ -= flux / volume_;
+ dCbyDt_ -= flux / volume_;
}
-void DifShell::localTauPump( double kP, double Ceq )
+void DifShell::vTauPump(const Eref& e, double kP, double Ceq )
{
- dCbyDt_ += -kP * ( C_ - Ceq );
+ //dCbyDt_ += -kP * ( C_ - Ceq );
+ dCbyDt_ += kP*Ceq;
+ Cmultiplier_ -= kP;
}
/**
- * Same as tauPump, except that we use the local value of Ceq.
+ * Same as tauPump, except that we use the v value of Ceq.
*/
-void DifShell::localEqTauPump( double kP )
-{
- dCbyDt_ += -kP * ( C_ - Ceq_ );
-}
-
-void DifShell::localMMPump( double vMax, double Kd )
-{
- dCbyDt_ += -( vMax / volume_ ) * ( C_ / ( C_ + Kd ) );
-}
-
-void DifShell::localHillPump( double vMax, double Kd, unsigned int hill )
-{
- double ch;
- switch( hill )
- {
- case 0:
- ch = 1.0;
- break;
- case 1:
- ch = C_;
- break;
- case 2:
- ch = C_ * C_;
- break;
- case 3:
- ch = C_ * C_ * C_;
- break;
- case 4:
- ch = C_ * C_;
- ch = ch * ch;
- break;
- default:
- ch = pow( C_, static_cast< double >( hill ) );
- };
+void DifShell::vEqTauPump(const Eref& e, double kP )
+{
+ dCbyDt_ += kP*Ceq_;
+ Cmultiplier_ -= kP;
+}
+
+void DifShell::vMMPump(const Eref& e, double vMax, double Kd )
+{
+ Cmultiplier_ += ( vMax / volume_ ) / ( C_ + Kd ) ;
+}
+
+void DifShell::vHillPump(const Eref& e, double vMax, double Kd, unsigned int hill )
+{
+ //This needs fixing
+ double ch;
+ switch( hill )
+ {
+ case 0:
+ ch = 1.0;
+ break;
+ case 1:
+ ch = C_;
+ break;
+ case 2:
+ ch = C_ * C_;
+ break;
+ case 3:
+ ch = C_ * C_ * C_;
+ break;
+ case 4:
+ ch = C_ * C_;
+ ch = ch * ch;
+ break;
+ default:
+ ch = pow( C_, static_cast< double >( hill ) );
+ };
- dCbyDt_ += -( vMax / volume_ ) * ( ch / ( ch + Kd ) );
+ dCbyDt_ += -( vMax / volume_ ) * ( ch / ( ch + Kd ) );
}
////////////////////////////////////////////////////////////////////////////////
diff --git a/moose-core/biophysics/DifShell.h b/moose-core/biophysics/DifShell.h
index 42f70600742dd163dff44138423428df38b73ed7..42cf643d69ca04e9af18f9cd5da8ea1729dcdf83 100644
--- a/moose-core/biophysics/DifShell.h
+++ b/moose-core/biophysics/DifShell.h
@@ -1,156 +1,105 @@
/**********************************************************************
-** This program is part of 'MOOSE', the
-** Multiscale Object Oriented Simulation Environment.
-** copyright (C) 2003-2008
-** Upinder S. Bhalla, Niraj Dudani and NCBS
-** It is made available under the terms of the
-** GNU Lesser General Public License version 2.1
-** See the file COPYING.LIB for the full notice.
-**********************************************************************/
-
-#ifndef _DifShell_h
-#define _DifShell_h
-
-class DifShell
-{
- public:
- DifShell();
-
- /////////////////////////////////////////////////////////////
- // Field access functions
- /////////////////////////////////////////////////////////////
- double getC( ) const;
-
- void setCeq(double Ceq );
- double getCeq() const;
-
- void setD( double D );
- double getD() const;
-
- void setValence( double valence );
- double getValence() const;
-
- void setLeak( double leak );
- double getLeak() const;
-
- void setShapeMode( unsigned int shapeMode );
- unsigned int getShapeMode() const;
-
- void setLength( double length );
- double getLength() const;
-
- void setDiameter( double diameter );
- double getDiameter() const;
-
- void setThickness( double thickness );
- double getThickness() const;
-
- void setVolume( double volume );
- double getVolume() const;
-
- void setOuterArea( double outerArea );
- double getOuterArea() const;
-
- void setInnerArea( double innerArea );
- double getInnerArea() const;
-
- /////////////////////////////////////////////////////////////
- // Dest functions
- /////////////////////////////////////////////////////////////
- void reinit_0( const Eref & e, ProcPtr p );
-
- void process_0(const Eref & e, ProcPtr p );
-
- void process_1(const Eref & e, ProcPtr p );
-
- void reinit_1(const Eref & e, ProcPtr p ); // dummyFunc
-
- void buffer(
- double kf,
- double kb,
- double bFree,
- double bBound );
-
- void fluxFromOut(
- double outerC,
- double outerThickness );
-
- void fluxFromIn(
- double innerC,
- double innerThickness );
-
- void influx(
- double I );
-
- void outflux(
- double I );
-
- void fInflux(
- double I,
- double fraction );
-
- void fOutflux(
- double I,
- double fraction );
-
- void storeInflux(
- double flux );
-
- void storeOutflux(
- double flux );
-
- void tauPump(
- double kP,
- double Ceq );
-
- void eqTauPump(
- double kP );
-
- void mmPump(
- double vMax,
- double Kd );
-
- void hillPump(
- double vMax,
- double Kd,
- unsigned int hill );
-
- static const Cinfo * initCinfo();
+ ** This program is part of 'MOOSE', the
+ ** Multiscale Object Oriented Simulation Environment.
+ ** copyright (C) 2003-2008
+ ** Upinder S. Bhalla, Niraj Dudani and NCBS
+ ** It is made available under the terms of the
+ ** GNU Lesser General Public License version 2.1
+ ** See the file COPYING.LIB for the full notice.
+ **********************************************************************/
+
+#ifndef _DIFSHELL_H
+#define _DIFSHELL_H
+
+class DifShell: public DifShellBase{
+ public:
+ DifShell();
+
+ /////////////////////////////////////////////////////////////
+ // Dest functions
+ /////////////////////////////////////////////////////////////
+ void vReinit( const Eref & e, ProcPtr p );
+ void vProcess(const Eref & e, ProcPtr p );
+ void vBuffer(const Eref& e, double kf, double kb, double bFree, double bBound );
+ void vFluxFromOut(const Eref& e, double outerC, double outerThickness );
+ void vFluxFromIn(const Eref& e, double innerC, double innerThickness );
+ void vInflux(const Eref& e, double I );
+ void vOutflux(const Eref& e, double I );
+ void vFInflux(const Eref& e, double I, double fraction );
+ void vFOutflux(const Eref& e, double I, double fraction );
+ void vStoreInflux(const Eref& e, double flux );
+ void vStoreOutflux(const Eref& e, double flux );
+ void vTauPump(const Eref& e, double kP, double Ceq );
+ void vEqTauPump(const Eref& e, double kP );
+ void vMMPump(const Eref& e, double vMax, double Kd );
+ void vHillPump(const Eref& e, double vMax, double Kd, unsigned int hill );
+
+ /////////////////////////////////////////////////////////////
+ // Field access functions
+ /////////////////////////////////////////////////////////////
+
+ void vSetC(const Eref& e,double C);
+ double vGetC( const Eref& e) const;
+
+ void vSetCeq(const Eref& e,double Ceq );
+ double vGetCeq(const Eref& e) const;
+
+ void vSetD(const Eref& e, double D );
+ double vGetD(const Eref& e) const;
+
+ void vSetValence(const Eref& e, double valence );
+ double vGetValence(const Eref& e) const;
+
+ void vSetLeak(const Eref& e, double leak );
+ double vGetLeak(const Eref& e) const;
+
+ void vSetShapeMode(const Eref& e, unsigned int shapeMode );
+ unsigned int vGetShapeMode(const Eref& e) const;
+
+ void vSetLength(const Eref& e, double length );
+ double vGetLength(const Eref& e) const;
+
+ void vSetDiameter(const Eref& e, double diameter );
+ double vGetDiameter(const Eref& e) const;
+
+ void vSetThickness(const Eref& e, double thickness );
+ double vGetThickness(const Eref& e) const;
+
+ void vSetVolume(const Eref& e, double volume );
+ double vGetVolume(const Eref& e) const;
+
+ void vSetOuterArea(const Eref& e, double outerArea );
+ double vGetOuterArea(const Eref& e) const;
+
+ void vSetInnerArea(const Eref& e, double innerArea );
+ double vGetInnerArea(const Eref& e) const;
+ static const Cinfo * initCinfo();
+
- private:
- void localReinit_0( const Eref & e, ProcPtr p );
- void localProcess_0( const Eref & e, ProcPtr p );
- void localProcess_1( const Eref & e, ProcPtr p );
- void localBuffer( double kf, double kb, double bFree, double bBound );
- void localFluxFromOut( double outerC, double outerThickness );
- void localFluxFromIn( double innerC, double innerThickness );
- void localInflux( double I );
- void localOutflux( double I );
- void localFInflux( double I, double fraction );
- void localFOutflux( double I, double fraction );
- void localStoreInflux( double flux );
- void localStoreOutflux( double flux );
- void localTauPump( double kP, double Ceq );
- void localEqTauPump( double kP );
- void localMMPump( double vMax, double Kd );
- void localHillPump( double vMax, double Kd, unsigned int hill );
-
- double dCbyDt_;
- double C_;
- double Ceq_;
- double D_;
- double valence_;
- double leak_;
- unsigned int shapeMode_;
- double length_;
- double diameter_;
- double thickness_;
- double volume_;
- double outerArea_;
- double innerArea_;
-
- /// Faraday's constant (Coulomb / Mole)
- static const double F_;
+ private:
+
+ double integrate( double state, double dt, double A, double B );
+
+ double dCbyDt_;
+ double Cmultiplier_;
+ double C_;
+ double Ceq_;
+ double prevC_;
+ double D_;
+ double valence_;
+ double leak_;
+ unsigned int shapeMode_;
+ double length_;
+ double diameter_;
+ double thickness_;
+ double volume_;
+ double outerArea_;
+ double innerArea_;
+
+ static const double EPSILON;
+ /// Faraday's constant (Coulomb / Mole)
+ static const double F;
+
};
-#endif // _DifShell_h
+#endif // _DIFSHELL_H
diff --git a/moose-core/biophysics/DifShellBase.cpp b/moose-core/biophysics/DifShellBase.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e6d1e5338b9242c03fcf6a7c49bc475faeabbea9
--- /dev/null
+++ b/moose-core/biophysics/DifShellBase.cpp
@@ -0,0 +1,513 @@
+/**********************************************************************
+ ** This program is part of 'MOOSE', the
+ ** Multiscale Object Oriented Simulation Environment.
+ ** copyright (C) 2003-2008
+ ** Upinder S. Bhalla, Niraj Dudani and NCBS
+ ** It is made available under the terms of the
+ ** GNU Lesser General Public License version 2.1
+ ** See the file COPYING.LIB for the full notice.
+ ****/
+
+#include "header.h"
+#include "DifShellBase.h"
+#include "ElementValueFinfo.h"
+
+SrcFinfo1< double >* DifShellBase::concentrationOut()
+{
+ static SrcFinfo1< double > concentrationOut("concentrationOut",
+ "Sends out concentration");
+ return &concentrationOut;
+}
+
+SrcFinfo2< double, double >* DifShellBase::innerDifSourceOut(){
+ static SrcFinfo2< double, double > sourceOut("innerDifSourceOut",
+ "Sends out source information.");
+ return &sourceOut;
+}
+
+SrcFinfo2< double, double >* DifShellBase::outerDifSourceOut(){
+ static SrcFinfo2< double, double > sourceOut("outerDifSourceOut",
+ "Sends out source information.");
+ return &sourceOut;
+}
+
+const Cinfo* DifShellBase::initCinfo()
+{
+
+ static DestFinfo process( "process",
+ "Handles process call",
+ new ProcOpFunc< DifShellBase>(&DifShellBase::process ) );
+ static DestFinfo reinit( "reinit",
+ "Reinit happens only in stage 0",
+ new ProcOpFunc< DifShellBase>( &DifShellBase::reinit ));
+
+ static Finfo* processShared[] = {
+ &process, &reinit
+ };
+
+ static SharedFinfo proc(
+ "proc",
+ "Shared message to receive Process message from scheduler",
+ processShared, sizeof( processShared ) / sizeof( Finfo* ));
+
+
+
+
+ static DestFinfo reaction( "reaction",
+ "Here the DifShell receives reaction rates (forward and backward), and concentrations for the "
+ "free-buffer and bound-buffer molecules.",
+ new EpFunc4< DifShellBase, double, double, double, double >( &DifShellBase::buffer ));
+
+ static Finfo* bufferShared[] = {
+ DifShellBase::concentrationOut(), &reaction
+ };
+
+ static SharedFinfo buffer( "buffer",
+ "This is a shared message from a DifShell to a Buffer (FixBuffer or DifBuffer). " ,
+ bufferShared,
+ sizeof( bufferShared ) / sizeof( Finfo* ));
+ /////
+
+
+
+
+ static DestFinfo fluxFromOut( "fluxFromOut",
+ "Destination message",
+ new EpFunc2< DifShellBase, double, double > ( &DifShellBase::fluxFromOut ));
+
+ static Finfo* innerDifShared[] = {
+ &fluxFromOut, DifShellBase::innerDifSourceOut(),
+ };
+ static SharedFinfo innerDif( "innerDif",
+ "This shared message (and the next) is between DifShells: adjoining shells exchange information to "
+ "find out the flux between them. "
+ "Using this message, an inner shell sends to, and receives from its outer shell." ,
+ innerDifShared,
+ sizeof( innerDifShared ) / sizeof( Finfo* ));
+
+ static DestFinfo fluxFromIn( "fluxFromIn", "",
+ new EpFunc2< DifShellBase, double, double> ( &DifShellBase::fluxFromIn ) );
+
+ static Finfo* outerDifShared[] = {
+ &fluxFromIn, DifShellBase::outerDifSourceOut(),
+ };
+
+ static SharedFinfo outerDif( "outerDif",
+ "Using this message, an outer shell sends to, and receives from its inner shell." ,
+ outerDifShared,
+ sizeof( outerDifShared ) / sizeof( Finfo* ));
+
+ static ElementValueFinfo< DifShellBase, double> C( "C",
+ "Concentration C",// is computed by the DifShell",
+ &DifShellBase::setC,
+ &DifShellBase::getC);
+ static ElementValueFinfo< DifShellBase, double> Ceq( "Ceq", "",
+ &DifShellBase::setCeq,
+ &DifShellBase::getCeq);
+ static ElementValueFinfo< DifShellBase, double> D( "D", "",
+ &DifShellBase::setD,
+ &DifShellBase::getD);
+ static ElementValueFinfo< DifShellBase, double> valence( "valence", "",
+ &DifShellBase::setValence,
+ &DifShellBase::getValence);
+ static ElementValueFinfo< DifShellBase, double> leak( "leak", "",
+ &DifShellBase::setLeak,
+ &DifShellBase::getLeak);
+ static ElementValueFinfo< DifShellBase, unsigned int> shapeMode( "shapeMode", "",
+ &DifShellBase::setShapeMode,
+ &DifShellBase::getShapeMode);
+ static ElementValueFinfo< DifShellBase, double> length( "length", "",
+ &DifShellBase::setLength,
+ &DifShellBase::getLength);
+ static ElementValueFinfo< DifShellBase, double> diameter( "diameter", "",
+ &DifShellBase::setDiameter,
+ &DifShellBase::getDiameter );
+ static ElementValueFinfo< DifShellBase, double> thickness( "thickness", "",
+ &DifShellBase::setThickness,
+ &DifShellBase::getThickness );
+ static ElementValueFinfo< DifShellBase, double> volume( "volume", "",
+ &DifShellBase::setVolume,
+ &DifShellBase::getVolume );
+ static ElementValueFinfo< DifShellBase, double> outerArea( "outerArea", "",
+ &DifShellBase::setOuterArea,
+ &DifShellBase::getOuterArea);
+ static ElementValueFinfo< DifShellBase, double> innerArea( "innerArea", "",
+ &DifShellBase::setInnerArea,
+ &DifShellBase::getInnerArea );
+
+ static DestFinfo mmPump( "mmPump", "Here DifShell receives pump outflux",
+ new EpFunc2< DifShellBase, double, double >( &DifShellBase::mmPump ) );
+ static DestFinfo influx( "influx", "",
+ new EpFunc1< DifShellBase, double > (&DifShellBase::influx ));
+ static DestFinfo outflux( "outflux", "",
+ new EpFunc1< DifShellBase, double >( &DifShellBase::influx ));
+ static DestFinfo fInflux( "fInflux", "",
+ new EpFunc2< DifShellBase, double, double >( &DifShellBase::fInflux ));
+ static DestFinfo fOutflux( "fOutflux", "",
+ new EpFunc2< DifShellBase, double, double> (&DifShellBase::fOutflux ));
+ static DestFinfo storeInflux( "storeInflux", "",
+ new EpFunc1< DifShellBase, double >( &DifShellBase::storeInflux ) );
+ static DestFinfo storeOutflux( "storeOutflux", "",
+ new EpFunc1< DifShellBase, double > ( &DifShellBase::storeOutflux ) );
+ static DestFinfo tauPump( "tauPump","",
+ new EpFunc2< DifShellBase, double, double >( &DifShellBase::tauPump ) );
+ static DestFinfo eqTauPump( "eqTauPump", "",
+ new EpFunc1< DifShellBase, double >( &DifShellBase::eqTauPump ) );
+ static DestFinfo hillPump( "hillPump", "",
+ new EpFunc3< DifShellBase, double, double, unsigned int >( &DifShellBase::hillPump ) );
+ static Finfo* difShellBaseFinfos[] = {
+ //////////////////////////////////////////////////////////////////
+ // Field definitions
+ //////////////////////////////////////////////////////////////////
+ &C,
+ &Ceq,
+ &D,
+ &valence,
+ &leak,
+ &shapeMode,
+ &length,
+ &diameter,
+ &thickness,
+ &volume,
+ &outerArea,
+ &innerArea,
+ //////////////////////////////////////////////////////////////////
+ // MsgSrc definitions
+ //////////////////////////////////////////////////////////////////
+
+ //////////////////////////////////////////////////////////////////
+ // SharedFinfo definitions
+ //////////////////////////////////////////////////////////////////
+ &proc,
+ &buffer,
+ concentrationOut(),
+ innerDifSourceOut(),
+ outerDifSourceOut(),
+ &innerDif,
+ &outerDif,
+ //////////////////////////////////////////////////////////////////
+ // DestFinfo definitions
+ //////////////////////////////////////////////////////////////////
+ &influx,
+ &outflux,
+ &fInflux,
+ &fOutflux,
+ &storeInflux,
+ &storeOutflux,
+ &tauPump,
+ &eqTauPump,
+ &mmPump,
+ &hillPump,
+ };
+
+ static string doc[] =
+ {
+ "Name", "DifShellBase",
+ "Author", "Niraj Dudani. Ported to async13 by Subhasis Ray/Asia Jedrzejewska-Szmek",
+ "Description", "DifShell object: Models diffusion of an ion (typically calcium) within an "
+ "electric compartment. A DifShell is an iso-concentration region with respect to "
+ "the ion. Adjoining DifShells exchange flux of this ion, and also keep track of "
+ "changes in concentration due to pumping, buffering and channel currents, by "
+ "talking to the appropriate objects.",
+ };
+ static ZeroSizeDinfo< int > dinfo;
+ //static Dinfo< DifShellBase> dinfo;
+ static Cinfo difShellBaseCinfo(
+ "DifShellBase",
+ Neutral::initCinfo(),
+ difShellBaseFinfos,
+ sizeof( difShellBaseFinfos ) / sizeof( Finfo* ),
+ &dinfo,
+ doc,
+ sizeof( doc ) / sizeof( string ));
+
+ return &difShellBaseCinfo;
+}
+
+//Cinfo *object* corresponding to the class.
+static const Cinfo* difShellBaseCinfo = DifShellBase::initCinfo();
+
+DifShellBase::DifShellBase()
+{ ; }
+
+void DifShellBase::setC( const Eref& e, double C)
+{
+ vSetC(e,C);
+}
+double DifShellBase::getC(const Eref& e) const
+{
+ return vGetC(e);
+}
+
+void DifShellBase::setCeq( const Eref& e, double Ceq )
+{
+ vSetCeq(e,Ceq);
+}
+
+double DifShellBase::getCeq(const Eref& e ) const
+{
+ return vGetCeq(e);
+}
+
+void DifShellBase::setD(const Eref& e, double D )
+{
+ vSetD(e,D);
+}
+
+double DifShellBase::getD(const Eref& e ) const
+{
+ return vGetD(e);
+}
+
+void DifShellBase::setValence(const Eref& e, double valence )
+{
+ vSetValence(e,valence);
+}
+
+double DifShellBase::getValence(const Eref& e ) const
+{
+ return vGetValence(e);
+}
+
+void DifShellBase::setLeak(const Eref& e, double leak )
+{
+ vSetLeak(e,leak);
+}
+
+double DifShellBase::getLeak(const Eref& e ) const
+{
+ return vGetLeak(e);
+}
+
+void DifShellBase::setShapeMode(const Eref& e, unsigned int shapeMode )
+{
+ vSetShapeMode(e,shapeMode);
+}
+
+unsigned int DifShellBase::getShapeMode(const Eref& e) const
+{
+ return vGetShapeMode(e);
+}
+
+void DifShellBase::setLength(const Eref& e, double length )
+{
+ vSetLength(e,length);
+}
+
+double DifShellBase::getLength(const Eref& e ) const
+{
+ return vGetLength(e);
+}
+
+void DifShellBase::setDiameter(const Eref& e, double diameter )
+{
+ vSetDiameter(e,diameter);
+}
+
+double DifShellBase::getDiameter(const Eref& e ) const
+{
+ return vGetDiameter(e);
+}
+
+void DifShellBase::setThickness( const Eref& e, double thickness )
+{
+ vSetThickness(e,thickness);
+}
+
+double DifShellBase::getThickness(const Eref& e) const
+{
+ return vGetThickness(e);
+}
+
+void DifShellBase::setVolume(const Eref& e, double volume )
+{
+ vSetVolume(e,volume);
+}
+
+double DifShellBase::getVolume(const Eref& e ) const
+{
+ return vGetVolume(e);
+}
+
+void DifShellBase::setOuterArea(const Eref& e, double outerArea )
+{
+ vSetOuterArea(e,outerArea);
+}
+
+double DifShellBase::getOuterArea(const Eref& e ) const
+{
+ return vGetOuterArea(e);
+}
+
+void DifShellBase::setInnerArea(const Eref& e, double innerArea )
+{
+ vSetInnerArea(e,innerArea);
+}
+
+double DifShellBase::getInnerArea(const Eref& e) const
+{
+ return vGetInnerArea(e);
+}
+
+////////////////////////////////////////////////////////////////////////////////
+// Dest functions
+////////////////////////////////////////////////////////////////////////////////
+
+void DifShellBase::reinit( const Eref& e, ProcPtr p )
+{
+ vReinit( e, p );
+}
+
+void DifShellBase::process( const Eref& e, ProcPtr p )
+{
+ vProcess( e, p );
+}
+
+void DifShellBase::buffer(
+ const Eref& e,
+ double kf,
+ double kb,
+ double bFree,
+ double bBound )
+{
+ vBuffer( e, kf, kb, bFree, bBound );
+}
+
+void DifShellBase::fluxFromOut(const Eref& e,
+ double outerC,
+ double outerThickness )
+{
+ vFluxFromOut(e, outerC, outerThickness );
+}
+
+void DifShellBase::fluxFromIn(
+ const Eref& e,
+ double innerC,
+ double innerThickness )
+{
+ vFluxFromIn( e, innerC, innerThickness );
+}
+
+void DifShellBase::influx(const Eref& e,
+ double I )
+{
+ vInflux( e, I );
+}
+
+void DifShellBase::outflux(const Eref& e,
+ double I )
+{
+ vOutflux(e, I );
+}
+
+void DifShellBase::fInflux(const Eref& e,
+ double I,
+ double fraction )
+{
+ vFInflux(e, I, fraction );
+}
+
+void DifShellBase::fOutflux(const Eref& e,
+ double I,
+ double fraction )
+{
+ vFOutflux(e, I, fraction );
+}
+
+void DifShellBase::storeInflux(const Eref& e,
+ double flux )
+{
+ vStoreInflux( e, flux );
+}
+
+void DifShellBase::storeOutflux(const Eref& e,
+ double flux )
+{
+ vStoreOutflux(e, flux );
+}
+
+void DifShellBase::tauPump(const Eref& e,
+ double kP,
+ double Ceq )
+{
+ vTauPump( e, kP, Ceq );
+}
+
+void DifShellBase::eqTauPump(const Eref& e,
+ double kP )
+{
+ vEqTauPump(e, kP );
+}
+
+void DifShellBase::mmPump(const Eref& e,
+ double vMax,
+ double Kd )
+{
+ vMMPump(e, vMax, Kd );
+}
+
+void DifShellBase::hillPump(const Eref& e,
+ double vMax,
+ double Kd,
+ unsigned int hill )
+{
+ vHillPump(e, vMax, Kd, hill );
+}
+void DifShellBase::vSetSolver( const Eref& e, Id hsolve )
+{;}
+
+void DifShellBase::zombify( Element* orig, const Cinfo* zClass,
+ Id hsolve )
+{
+ if ( orig->cinfo() == zClass )
+ return;
+ unsigned int start = orig->localDataStart();
+ unsigned int num = orig->numLocalData();
+ if ( num == 0 )
+ return;
+ unsigned int len = 12;
+ vector< double > data( num * len );
+
+ unsigned int j = 0;
+
+ for ( unsigned int i = 0; i < num; ++i ) {
+ Eref er( orig, i + start );
+ const DifShellBase* ds =
+ reinterpret_cast< const DifShellBase* >( er.data() );
+ data[j + 0] = ds->getC( er );
+ data[j + 1] = ds->getCeq( er );
+ data[j + 2] = ds->getD( er );
+ data[j + 3] = ds->getValence( er );
+ data[j + 4] = ds->getLeak( er );
+ data[j + 5] = ds->getShapeMode( er );
+ data[j + 6] = ds->getLength( er );
+ data[j + 7] = ds->getDiameter( er );
+ data[j + 8] = ds->getThickness( er );
+ data[j + 9] = ds->getVolume( er );
+ data[j + 10] = ds->getOuterArea( er );
+ data[j + 11] = ds->getInnerArea( er );
+ j += len;
+ }
+ orig->zombieSwap( zClass );
+ j = 0;
+ for ( unsigned int i = 0; i < num; ++i ) {
+ Eref er( orig, i + start );
+ DifShellBase* ds =
+ reinterpret_cast< DifShellBase* >( er.data() );
+ ds->vSetSolver(er,hsolve);
+ ds->setC(er, data[j+0]);
+ ds->setCeq(er, data[j + 1]);
+ ds->setD(er, data[j + 2]);
+ ds->setValence(er, data[j + 3]);
+ ds->setLeak(er, data[j + 4]);
+ ds->setShapeMode(er, data[j + 5]);
+ ds->setLength(er, data[j + 6]);
+ ds->setDiameter(er, data[j + 7]);
+ ds->setThickness(er, data[j + 8]);
+ ds->setVolume(er, data[j + 9]);
+ ds->setOuterArea(er, data[j + 10]);
+ ds->setInnerArea(er, data[j + 11]);
+ j += len; //??
+ }
+
+}
diff --git a/moose-core/biophysics/DifShellBase.h b/moose-core/biophysics/DifShellBase.h
new file mode 100644
index 0000000000000000000000000000000000000000..db74c53696c1d133d4bfbf96432da1e6a46b4ccb
--- /dev/null
+++ b/moose-core/biophysics/DifShellBase.h
@@ -0,0 +1,144 @@
+/**********************************************************************
+ ** This program is part of 'MOOSE', the
+ ** Multiscale Object Oriented Simulation Environment.
+ ** copyright (C) 2003-2008
+ ** Upinder S. Bhalla, Niraj Dudani and NCBS
+ ** It is made available under the terms of the
+ ** GNU Lesser General Public License version 2.1
+ ** See the file COPYING.LIB for the full notice.
+ **********************************************************************/
+
+#ifndef _DIFSHELL_BASE_H
+#define _DIFSHELL_BASE_H
+/*This is base class for DifShell*/
+
+class DifShellBase
+{
+ public:
+ DifShellBase();
+ /////////////////////////////////////////////////////////////
+ // Dest functions
+ /////////////////////////////////////////////////////////////
+ void reinit( const Eref & e, ProcPtr p );
+ void process(const Eref & e, ProcPtr p );
+ void buffer( const Eref& e, double kf, double kb, double bFree, double bBound );
+ void fluxFromOut(const Eref& e, double outerC, double outerThickness );
+ void fluxFromIn(const Eref& e, double innerC, double innerThickness );
+ void influx(const Eref& e, double I );
+ void outflux(const Eref& e, double I );
+ void fInflux(const Eref& e, double I, double fraction );
+ void fOutflux(const Eref& e, double I, double fraction );
+ void storeInflux(const Eref& e, double flux );
+ void storeOutflux(const Eref& e, double flux );
+ void tauPump(const Eref& e, double kP, double Ceq );
+ void eqTauPump(const Eref& e, double kP );
+ void mmPump(const Eref& e, double vMax, double Kd );
+ void hillPump(const Eref& e, double vMax, double Kd, unsigned int hill);
+
+ virtual void vReinit(const Eref & e, ProcPtr p ) = 0;
+ virtual void vProcess(const Eref & e, ProcPtr p ) = 0;
+ virtual void vBuffer(const Eref& e, double kf, double kb, double bFree, double bBound ) = 0;
+ virtual void vFluxFromOut(const Eref& e, double outerC, double outerThickness ) = 0;
+ virtual void vFluxFromIn(const Eref& e, double innerC, double innerThickness ) = 0;
+ virtual void vInflux(const Eref& e, double I ) = 0;
+ virtual void vOutflux(const Eref& e, double I ) = 0;
+ virtual void vFInflux(const Eref& e, double I, double fraction ) = 0;
+ virtual void vFOutflux(const Eref& e, double I, double fraction ) = 0;
+ virtual void vStoreInflux(const Eref& e, double flux ) = 0;
+ virtual void vStoreOutflux(const Eref& e, double flux ) = 0;
+ virtual void vTauPump(const Eref& e, double kP, double Ceq ) = 0;
+ virtual void vEqTauPump(const Eref& e, double kP ) = 0;
+ virtual void vMMPump(const Eref& e, double vMax, double Kd ) = 0;
+ virtual void vHillPump(const Eref& e, double vMax, double Kd, unsigned int hill ) = 0;
+ /////////////////////////////////////////////////////////////
+ // Field access functions
+ /////////////////////////////////////////////////////////////
+
+ void setC(const Eref& e,double C);
+ double getC( const Eref& e) const;
+
+ void setCeq(const Eref& e,double Ceq );
+ double getCeq(const Eref& e) const;
+
+ void setD(const Eref& e, double D );
+ double getD(const Eref& e) const;
+
+ void setValence(const Eref& e, double valence );
+ double getValence(const Eref& e) const;
+
+ void setLeak(const Eref& e, double leak );
+ double getLeak(const Eref& e) const;
+
+ void setShapeMode(const Eref& e, unsigned int shapeMode );
+ unsigned int getShapeMode(const Eref& e) const;
+
+ void setLength(const Eref& e, double length );
+ double getLength(const Eref& e) const;
+
+ void setDiameter(const Eref& e, double diameter );
+ double getDiameter(const Eref& e) const;
+
+ void setThickness(const Eref& e, double thickness );
+ double getThickness(const Eref& e) const;
+
+ void setVolume(const Eref& e, double volume );
+ double getVolume(const Eref& e) const;
+
+ void setOuterArea(const Eref& e, double outerArea );
+ double getOuterArea(const Eref& e) const;
+
+ void setInnerArea(const Eref& e, double innerArea );
+ double getInnerArea(const Eref& e) const;
+
+ virtual void vSetC(const Eref& e,double C) = 0;
+ virtual double vGetC( const Eref& e) const = 0;
+
+ virtual void vSetCeq(const Eref& e,double Ceq ) = 0;
+ virtual double vGetCeq(const Eref& e) const = 0;
+
+ virtual void vSetD(const Eref& e, double D ) = 0;
+ virtual double vGetD(const Eref& e) const = 0;
+
+ virtual void vSetValence(const Eref& e, double valence ) = 0;
+ virtual double vGetValence(const Eref& e) const = 0;
+
+ virtual void vSetLeak(const Eref& e, double leak ) = 0;
+ virtual double vGetLeak(const Eref& e) const = 0;
+
+ virtual void vSetShapeMode(const Eref& e, unsigned int shapeMode ) = 0;
+ virtual unsigned int vGetShapeMode(const Eref& e) const = 0;
+
+ virtual void vSetLength(const Eref& e, double length ) = 0;
+ virtual double vGetLength(const Eref& e) const = 0;
+
+ virtual void vSetDiameter(const Eref& e, double diameter ) = 0;
+ virtual double vGetDiameter(const Eref& e) const = 0;
+
+ virtual void vSetThickness(const Eref& e, double thickness ) = 0;
+ virtual double vGetThickness(const Eref& e) const = 0;
+
+ virtual void vSetVolume(const Eref& e, double volume ) = 0;
+ virtual double vGetVolume(const Eref& e) const = 0;
+
+ virtual void vSetOuterArea(const Eref& e, double outerArea ) = 0;
+ virtual double vGetOuterArea(const Eref& e) const = 0;
+
+ virtual void vSetInnerArea(const Eref& e, double innerArea ) = 0;
+ virtual double vGetInnerArea(const Eref& e) const = 0;
+
+
+ virtual void vSetSolver( const Eref& e, Id hsolve );
+ static void zombify( Element* orig, const Cinfo* zClass,
+ Id hsolve );
+ static SrcFinfo1< double >* concentrationOut();
+ static SrcFinfo2< double, double >* innerDifSourceOut();
+ static SrcFinfo2< double, double >* outerDifSourceOut();
+ static const Cinfo * initCinfo();
+
+ private:
+
+
+
+};
+
+#endif //_DIFSHELL_BASE_H
diff --git a/moose-core/biophysics/HHGate.cpp b/moose-core/biophysics/HHGate.cpp
index 1f975b660e125eedc38fb29c0e1207165a39fc3d..db9042d251ad1a53fea3100d0b7efc22f7d617fb 100644
--- a/moose-core/biophysics/HHGate.cpp
+++ b/moose-core/biophysics/HHGate.cpp
@@ -1,11 +1,11 @@
/**********************************************************************
-** This program is part of 'MOOSE', the
-** Messaging Object Oriented Simulation Environment.
-** Copyright (C) 2003-2007 Upinder S. Bhalla. and NCBS
-** It is made available under the terms of the
-** GNU Lesser General Public License version 2.1
-** See the file COPYING.LIB for the full notice.
-**********************************************************************/
+ ** This program is part of 'MOOSE', the
+ ** Messaging Object Oriented Simulation Environment.
+ ** Copyright (C) 2003-2007 Upinder S. Bhalla. and NCBS
+ ** It is made available under the terms of the
+ ** GNU Lesser General Public License version 2.1
+ ** See the file COPYING.LIB for the full notice.
+ **********************************************************************/
#include "header.h"
#include "ElementValueFinfo.h"
@@ -15,209 +15,209 @@ static const double SINGULARITY = 1.0e-6;
const Cinfo* HHGate::initCinfo()
{
- ///////////////////////////////////////////////////////
- // Field definitions.
- ///////////////////////////////////////////////////////
- static ReadOnlyLookupValueFinfo< HHGate, double, double > A( "A",
- "lookupA: Look up the A gate value from a double. Usually does"
- "so by direct scaling and offset to an integer lookup, using"
- "a fine enough table granularity that there is little error."
- "Alternatively uses linear interpolation."
- "The range of the double is predefined based on knowledge of"
- "voltage or conc ranges, and the granularity is specified by"
- "the xmin, xmax, and dV fields.",
- &HHGate::lookupA );
- static ReadOnlyLookupValueFinfo< HHGate, double, double > B( "B",
- "lookupB: Look up the B gate value from a double."
- "Note that this looks up the raw tables, which are transformed"
- "from the reference parameters.",
- &HHGate::lookupB );
-
- static ElementValueFinfo< HHGate, vector< double > > alpha( "alpha",
- "Parameters for voltage-dependent rates, alpha:"
- "Set up alpha term using 5 parameters, as follows:"
- "y(x) = (A + B * x) / (C + exp((x + D) / F))"
- "The original HH equations can readily be cast into this form",
- &HHGate::setAlpha,
- &HHGate::getAlpha
- );
-
- static ElementValueFinfo< HHGate, vector< double > > beta( "beta",
- "Parameters for voltage-dependent rates, beta:"
- "Set up beta term using 5 parameters, as follows:"
- "y(x) = (A + B * x) / (C + exp((x + D) / F))"
- "The original HH equations can readily be cast into this form",
- &HHGate::setBeta,
- &HHGate::getBeta
- );
-
- static ElementValueFinfo< HHGate, vector< double > > tau( "tau",
- "Parameters for voltage-dependent rates, tau:"
- "Set up tau curve using 5 parameters, as follows:"
- "y(x) = (A + B * x) / (C + exp((x + D) / F))",
- &HHGate::setTau,
- &HHGate::getTau
- );
-
- static ElementValueFinfo< HHGate, vector< double > > mInfinity(
- "mInfinity",
- "Parameters for voltage-dependent rates, mInfinity:"
- "Set up mInfinity curve using 5 parameters, as follows:"
- "y(x) = (A + B * x) / (C + exp((x + D) / F))"
- "The original HH equations can readily be cast into this form",
- &HHGate::setMinfinity,
- &HHGate::getMinfinity
- );
-
- static ElementValueFinfo< HHGate, double > min( "min",
- "Minimum range for lookup",
- &HHGate::setMin,
- &HHGate::getMin
- );
-
- static ElementValueFinfo< HHGate, double > max( "max",
- "Minimum range for lookup",
- &HHGate::setMax,
- &HHGate::getMax
- );
-
- static ElementValueFinfo< HHGate, unsigned int > divs( "divs",
- "Divisions for lookup. Zero means to use linear interpolation",
- &HHGate::setDivs,
- &HHGate::getDivs
- );
-
- static ElementValueFinfo< HHGate, vector< double > > tableA(
- "tableA",
- "Table of A entries",
- &HHGate::setTableA,
- &HHGate::getTableA
- );
-
- static ElementValueFinfo< HHGate, vector< double > > tableB(
- "tableB",
- "Table of alpha + beta entries",
- &HHGate::setTableB,
- &HHGate::getTableB
- );
-
- static ElementValueFinfo< HHGate, bool > useInterpolation(
- "useInterpolation",
- "Flag: use linear interpolation if true, else direct lookup",
- &HHGate::setUseInterpolation,
- &HHGate::getUseInterpolation
- );
-
- static ElementValueFinfo< HHGate, vector< double > > alphaParms(
- "alphaParms",
- "Set up both gates using 13 parameters, as follows:"
- "setupAlpha AA AB AC AD AF BA BB BC BD BF xdivs xmin xmax"
- "Here AA-AF are Coefficients A to F of the alpha (forward) term"
- "Here BA-BF are Coefficients A to F of the beta (reverse) term"
- "Here xdivs is the number of entries in the table,"
- "xmin and xmax define the range for lookup."
- "Outside this range the returned value will be the low [high]"
- "entry of the table."
- "The equation describing each table is:"
- "y(x) = (A + B * x) / (C + exp((x + D) / F))"
- "The original HH equations can readily be cast into this form",
- &HHGate::setupAlpha,
- &HHGate::getAlphaParms
- );
-
- ///////////////////////////////////////////////////////
- // DestFinfos
- ///////////////////////////////////////////////////////
- static DestFinfo setupAlpha( "setupAlpha",
- "Set up both gates using 13 parameters, as follows:"
- "setupAlpha AA AB AC AD AF BA BB BC BD BF xdivs xmin xmax"
- "Here AA-AF are Coefficients A to F of the alpha (forward) term"
- "Here BA-BF are Coefficients A to F of the beta (reverse) term"
- "Here xdivs is the number of entries in the table,"
- "xmin and xmax define the range for lookup."
- "Outside this range the returned value will be the low [high]"
- "entry of the table."
- "The equation describing each table is:"
- "y(x) = (A + B * x) / (C + exp((x + D) / F))"
- "The original HH equations can readily be cast into this form",
- new EpFunc1< HHGate, vector< double > >( &HHGate::setupAlpha )
- );
- static DestFinfo setupTau( "setupTau",
- "Identical to setupAlpha, except that the forms specified by"
- "the 13 parameters are for the tau and m-infinity curves rather"
- "than the alpha and beta terms. So the parameters are:"
- "setupTau TA TB TC TD TF MA MB MC MD MF xdivs xmin xmax"
- "As before, the equation describing each curve is:"
- "y(x) = (A + B * x) / (C + exp((x + D) / F))",
- new EpFunc1< HHGate, vector< double > >( &HHGate::setupTau )
- );
- static DestFinfo tweakAlpha( "tweakAlpha",
- "Dummy function for backward compatibility. It used to convert"
- "the tables from alpha, beta values to alpha, alpha+beta"
- "because the internal calculations used these forms. Not"
- "needed now, deprecated.",
- new OpFunc0< HHGate >( &HHGate::tweakAlpha )
- );
- static DestFinfo tweakTau( "tweakTau",
- "Dummy function for backward compatibility. It used to convert"
- "the tables from tau, minf values to alpha, alpha+beta"
- "because the internal calculations used these forms. Not"
- "needed now, deprecated.",
- new OpFunc0< HHGate >( &HHGate::tweakTau )
- );
- static DestFinfo setupGate( "setupGate",
- "Sets up one gate at a time using the alpha/beta form."
- "Has 9 parameters, as follows:"
- "setupGate A B C D F xdivs xmin xmax is_beta"
- "This sets up the gate using the equation:"
- "y(x) = (A + B * x) / (C + exp((x + D) / F))"
- "Deprecated.",
- new EpFunc1< HHGate, vector< double > >( &HHGate::setupGate )
- );
- static Finfo* HHGateFinfos[] =
- {
- &A, // ReadOnlyLookupValue
- &B, // ReadOnlyLookupValue
- &alpha, // ElementValue
- &beta, // ElementValue
- &tau, // ElementValue
- &mInfinity, // ElementValue
- &min, // ElementValue
- &max, // ElementValue
- &divs, // ElementValue
- &tableA, // ElementValue
- &tableB, // ElementValue
- &useInterpolation, // ElementValue
- &alphaParms, // ElementValue
- &setupAlpha, // Dest
- &setupTau, // Dest
- &tweakAlpha, // Dest
- &tweakTau, // Dest
- &setupGate, // Dest
- };
-
- static string doc[] =
- {
- "Name", "HHGate",
- "Author", "Upinder S. Bhalla, 2011, NCBS",
- "Description", "HHGate: Gate for Hodkgin-Huxley type channels, equivalent to the "
- "m and h terms on the Na squid channel and the n term on K. "
- "This takes the voltage and state variable from the channel, "
- "computes the new value of the state variable and a scaling, "
- "depending on gate power, for the conductance.",
- };
-
- static Dinfo< HHGate > dinfo;
- static Cinfo HHGateCinfo(
- "HHGate",
- Neutral::initCinfo(),
- HHGateFinfos, sizeof(HHGateFinfos)/sizeof(Finfo *),
- &dinfo,
- doc,
- sizeof(doc)/sizeof(string)
- );
-
- return &HHGateCinfo;
+ ///////////////////////////////////////////////////////
+ // Field definitions.
+ ///////////////////////////////////////////////////////
+ static ReadOnlyLookupValueFinfo< HHGate, double, double > A( "A",
+ "lookupA: Look up the A gate value from a double. Usually does"
+ "so by direct scaling and offset to an integer lookup, using"
+ "a fine enough table granularity that there is little error."
+ "Alternatively uses linear interpolation."
+ "The range of the double is predefined based on knowledge of"
+ "voltage or conc ranges, and the granularity is specified by"
+ "the xmin, xmax, and dV fields.",
+ &HHGate::lookupA );
+ static ReadOnlyLookupValueFinfo< HHGate, double, double > B( "B",
+ "lookupB: Look up the B gate value from a double."
+ "Note that this looks up the raw tables, which are transformed"
+ "from the reference parameters.",
+ &HHGate::lookupB );
+
+ static ElementValueFinfo< HHGate, vector< double > > alpha( "alpha",
+ "Parameters for voltage-dependent rates, alpha:"
+ "Set up alpha term using 5 parameters, as follows:"
+ "y(x) = (A + B * x) / (C + exp((x + D) / F))"
+ "The original HH equations can readily be cast into this form",
+ &HHGate::setAlpha,
+ &HHGate::getAlpha
+ );
+
+ static ElementValueFinfo< HHGate, vector< double > > beta( "beta",
+ "Parameters for voltage-dependent rates, beta:"
+ "Set up beta term using 5 parameters, as follows:"
+ "y(x) = (A + B * x) / (C + exp((x + D) / F))"
+ "The original HH equations can readily be cast into this form",
+ &HHGate::setBeta,
+ &HHGate::getBeta
+ );
+
+ static ElementValueFinfo< HHGate, vector< double > > tau( "tau",
+ "Parameters for voltage-dependent rates, tau:"
+ "Set up tau curve using 5 parameters, as follows:"
+ "y(x) = (A + B * x) / (C + exp((x + D) / F))",
+ &HHGate::setTau,
+ &HHGate::getTau
+ );
+
+ static ElementValueFinfo< HHGate, vector< double > > mInfinity(
+ "mInfinity",
+ "Parameters for voltage-dependent rates, mInfinity:"
+ "Set up mInfinity curve using 5 parameters, as follows:"
+ "y(x) = (A + B * x) / (C + exp((x + D) / F))"
+ "The original HH equations can readily be cast into this form",
+ &HHGate::setMinfinity,
+ &HHGate::getMinfinity
+ );
+
+ static ElementValueFinfo< HHGate, double > min( "min",
+ "Minimum range for lookup",
+ &HHGate::setMin,
+ &HHGate::getMin
+ );
+
+ static ElementValueFinfo< HHGate, double > max( "max",
+ "Minimum range for lookup",
+ &HHGate::setMax,
+ &HHGate::getMax
+ );
+
+ static ElementValueFinfo< HHGate, unsigned int > divs( "divs",
+ "Divisions for lookup. Zero means to use linear interpolation",
+ &HHGate::setDivs,
+ &HHGate::getDivs
+ );
+
+ static ElementValueFinfo< HHGate, vector< double > > tableA(
+ "tableA",
+ "Table of A entries",
+ &HHGate::setTableA,
+ &HHGate::getTableA
+ );
+
+ static ElementValueFinfo< HHGate, vector< double > > tableB(
+ "tableB",
+ "Table of alpha + beta entries",
+ &HHGate::setTableB,
+ &HHGate::getTableB
+ );
+
+ static ElementValueFinfo< HHGate, bool > useInterpolation(
+ "useInterpolation",
+ "Flag: use linear interpolation if true, else direct lookup",
+ &HHGate::setUseInterpolation,
+ &HHGate::getUseInterpolation
+ );
+
+ static ElementValueFinfo< HHGate, vector< double > > alphaParms(
+ "alphaParms",
+ "Set up both gates using 13 parameters, as follows:"
+ "setupAlpha AA AB AC AD AF BA BB BC BD BF xdivs xmin xmax"
+ "Here AA-AF are Coefficients A to F of the alpha (forward) term"
+ "Here BA-BF are Coefficients A to F of the beta (reverse) term"
+ "Here xdivs is the number of entries in the table,"
+ "xmin and xmax define the range for lookup."
+ "Outside this range the returned value will be the low [high]"
+ "entry of the table."
+ "The equation describing each table is:"
+ "y(x) = (A + B * x) / (C + exp((x + D) / F))"
+ "The original HH equations can readily be cast into this form",
+ &HHGate::setupAlpha,
+ &HHGate::getAlphaParms
+ );
+
+ ///////////////////////////////////////////////////////
+ // DestFinfos
+ ///////////////////////////////////////////////////////
+ static DestFinfo setupAlpha( "setupAlpha",
+ "Set up both gates using 13 parameters, as follows:"
+ "setupAlpha AA AB AC AD AF BA BB BC BD BF xdivs xmin xmax"
+ "Here AA-AF are Coefficients A to F of the alpha (forward) term"
+ "Here BA-BF are Coefficients A to F of the beta (reverse) term"
+ "Here xdivs is the number of entries in the table,"
+ "xmin and xmax define the range for lookup."
+ "Outside this range the returned value will be the low [high]"
+ "entry of the table."
+ "The equation describing each table is:"
+ "y(x) = (A + B * x) / (C + exp((x + D) / F))"
+ "The original HH equations can readily be cast into this form",
+ new EpFunc1< HHGate, vector< double > >( &HHGate::setupAlpha )
+ );
+ static DestFinfo setupTau( "setupTau",
+ "Identical to setupAlpha, except that the forms specified by"
+ "the 13 parameters are for the tau and m-infinity curves rather"
+ "than the alpha and beta terms. So the parameters are:"
+ "setupTau TA TB TC TD TF MA MB MC MD MF xdivs xmin xmax"
+ "As before, the equation describing each curve is:"
+ "y(x) = (A + B * x) / (C + exp((x + D) / F))",
+ new EpFunc1< HHGate, vector< double > >( &HHGate::setupTau )
+ );
+ static DestFinfo tweakAlpha( "tweakAlpha",
+ "Dummy function for backward compatibility. It used to convert"
+ "the tables from alpha, beta values to alpha, alpha+beta"
+ "because the internal calculations used these forms. Not"
+ "needed now, deprecated.",
+ new OpFunc0< HHGate >( &HHGate::tweakAlpha )
+ );
+ static DestFinfo tweakTau( "tweakTau",
+ "Dummy function for backward compatibility. It used to convert"
+ "the tables from tau, minf values to alpha, alpha+beta"
+ "because the internal calculations used these forms. Not"
+ "needed now, deprecated.",
+ new OpFunc0< HHGate >( &HHGate::tweakTau )
+ );
+ static DestFinfo setupGate( "setupGate",
+ "Sets up one gate at a time using the alpha/beta form."
+ "Has 9 parameters, as follows:"
+ "setupGate A B C D F xdivs xmin xmax is_beta"
+ "This sets up the gate using the equation:"
+ "y(x) = (A + B * x) / (C + exp((x + D) / F))"
+ "Deprecated.",
+ new EpFunc1< HHGate, vector< double > >( &HHGate::setupGate )
+ );
+ static Finfo* HHGateFinfos[] =
+ {
+ &A, // ReadOnlyLookupValue
+ &B, // ReadOnlyLookupValue
+ &alpha, // ElementValue
+ &beta, // ElementValue
+ &tau, // ElementValue
+ &mInfinity, // ElementValue
+ &min, // ElementValue
+ &max, // ElementValue
+ &divs, // ElementValue
+ &tableA, // ElementValue
+ &tableB, // ElementValue
+ &useInterpolation, // ElementValue
+ &alphaParms, // ElementValue
+ &setupAlpha, // Dest
+ &setupTau, // Dest
+ &tweakAlpha, // Dest
+ &tweakTau, // Dest
+ &setupGate, // Dest
+ };
+
+ static string doc[] =
+ {
+ "Name", "HHGate",
+ "Author", "Upinder S. Bhalla, 2011, NCBS",
+ "Description", "HHGate: Gate for Hodkgin-Huxley type channels, equivalent to the "
+ "m and h terms on the Na squid channel and the n term on K. "
+ "This takes the voltage and state variable from the channel, "
+ "computes the new value of the state variable and a scaling, "
+ "depending on gate power, for the conductance.",
+ };
+
+ static Dinfo< HHGate > dinfo;
+ static Cinfo HHGateCinfo(
+ "HHGate",
+ Neutral::initCinfo(),
+ HHGateFinfos, sizeof(HHGateFinfos)/sizeof(Finfo *),
+ &dinfo,
+ doc,
+ sizeof(doc)/sizeof(string)
+ );
+
+ return &HHGateCinfo;
}
static const Cinfo* hhGateCinfo = HHGate::initCinfo();
@@ -225,22 +225,22 @@ static const Cinfo* hhGateCinfo = HHGate::initCinfo();
// Core class functions
///////////////////////////////////////////////////
HHGate::HHGate()
- : xmin_(0), xmax_(1), invDx_(1),
- originalChanId_(0),
- originalGateId_(0),
- lookupByInterpolation_(0),
- isDirectTable_(0)
+ : xmin_(0), xmax_(1), invDx_(1),
+ originalChanId_(0),
+ originalGateId_(0),
+ lookupByInterpolation_(0),
+ isDirectTable_(0)
{;}
HHGate::HHGate( Id originalChanId, Id originalGateId )
- :
- A_( 1, 0.0 ),
- B_( 1, 0.0 ),
- xmin_(0), xmax_(1), invDx_(1),
- originalChanId_( originalChanId ),
- originalGateId_( originalGateId ),
- lookupByInterpolation_(0),
- isDirectTable_(0)
+ :
+ A_( 1, 0.0 ),
+ B_( 1, 0.0 ),
+ xmin_(0), xmax_(1), invDx_(1),
+ originalChanId_( originalChanId ),
+ originalGateId_( originalGateId ),
+ lookupByInterpolation_(0),
+ isDirectTable_(0)
{;}
///////////////////////////////////////////////////
@@ -249,272 +249,272 @@ HHGate::HHGate( Id originalChanId, Id originalGateId )
double HHGate::lookupTable( const vector< double >& tab, double v ) const
{
- if ( v <= xmin_ ) return tab[0];
- if ( v >= xmax_ ) return tab.back();
- if ( lookupByInterpolation_ ) {
- unsigned int index =
- static_cast< unsigned int >( ( v - xmin_ ) * invDx_ );
- assert( tab.size() > index );
- double frac = ( v - xmin_ - index / invDx_ ) * invDx_;
- return tab[ index ] * ( 1 - frac ) + tab[ index + 1 ] * frac;
- } else {
- return tab[ static_cast< unsigned int >( (v - xmin_) * invDx_ ) ];
- }
+ if ( v <= xmin_ ) return tab[0];
+ if ( v >= xmax_ ) return tab.back();
+ if ( lookupByInterpolation_ ) {
+ unsigned int index =
+ static_cast< unsigned int >( ( v - xmin_ ) * invDx_ );
+ assert( tab.size() > index );
+ double frac = ( v - xmin_ - index / invDx_ ) * invDx_;
+ return tab[ index ] * ( 1 - frac ) + tab[ index + 1 ] * frac;
+ } else {
+ return tab[ static_cast< unsigned int >( (v - xmin_) * invDx_ ) ];
+ }
}
double HHGate::lookupA( double v ) const
{
- return lookupTable( A_, v );
+ return lookupTable( A_, v );
}
double HHGate::lookupB( double v ) const
{
- return lookupTable( B_, v );
+ return lookupTable( B_, v );
}
void HHGate::lookupBoth( double v, double* A, double* B ) const
{
- if ( v <= xmin_ ) {
- *A = A_[0];
- *B = B_[0];
- } else if ( v >= xmax_ ) {
- *A = A_.back();
- *B = B_.back();
- } else {
- unsigned int index =
- static_cast< unsigned int >( ( v - xmin_ ) * invDx_ );
- assert( A_.size() > index && B_.size() > index );
- if ( lookupByInterpolation_ ) {
- double frac = ( v - xmin_ - index / invDx_ ) * invDx_;
- *A = A_[ index ] * ( 1 - frac ) + A_[ index + 1 ] * frac;
- *B = B_[ index ] * ( 1 - frac ) + B_[ index + 1 ] * frac;
- } else {
- *A = A_[ index ];
- *B = B_[ index ];
- }
- }
+ if ( v <= xmin_ ) {
+ *A = A_[0];
+ *B = B_[0];
+ } else if ( v >= xmax_ ) {
+ *A = A_.back();
+ *B = B_.back();
+ } else {
+ unsigned int index =
+ static_cast< unsigned int >( ( v - xmin_ ) * invDx_ );
+ assert( A_.size() > index && B_.size() > index );
+ if ( lookupByInterpolation_ ) {
+ double frac = ( v - xmin_ - index / invDx_ ) * invDx_;
+ *A = A_[ index ] * ( 1 - frac ) + A_[ index + 1 ] * frac;
+ *B = B_[ index ] * ( 1 - frac ) + B_[ index + 1 ] * frac;
+ } else {
+ *A = A_[ index ];
+ *B = B_[ index ];
+ }
+ }
}
vector< double > HHGate::getAlpha( const Eref& e) const
{
- return alpha_;
+ return alpha_;
}
void HHGate::setAlpha( const Eref& e, vector< double > val )
{
- if ( val.size() != 5 ) {
- cout << "Error: HHGate::setAlpha on " << e.id().path() <<
- ": Number of entries on argument vector should be 5, was " <<
- val.size() << endl;
- return;
- }
- if ( checkOriginal( e.id(), "alpha" ) ) {
- alpha_ = val;
- updateTauMinf();
- updateTables();
- }
+ if ( val.size() != 5 ) {
+ cout << "Error: HHGate::setAlpha on " << e.id().path() <<
+ ": Number of entries on argument vector should be 5, was " <<
+ val.size() << endl;
+ return;
+ }
+ if ( checkOriginal( e.id(), "alpha" ) ) {
+ alpha_ = val;
+ updateTauMinf();
+ updateTables();
+ }
}
vector< double > HHGate::getBeta( const Eref& e) const
{
- return beta_;
+ return beta_;
}
void HHGate::setBeta( const Eref& e, vector< double > val )
{
- if ( val.size() != 5 ) {
- cout << "Error: HHGate::setBeta on " << e.id().path() <<
- ": Number of entries on argument vector should be 5, was " <<
- val.size() << endl;
- return;
- }
- if ( checkOriginal( e.id(), "beta" ) ) {
- beta_ = val;
- updateTauMinf();
- updateTables();
- }
+ if ( val.size() != 5 ) {
+ cout << "Error: HHGate::setBeta on " << e.id().path() <<
+ ": Number of entries on argument vector should be 5, was " <<
+ val.size() << endl;
+ return;
+ }
+ if ( checkOriginal( e.id(), "beta" ) ) {
+ beta_ = val;
+ updateTauMinf();
+ updateTables();
+ }
}
vector< double > HHGate::getTau( const Eref& e) const
{
- return tau_;
+ return tau_;
}
void HHGate::setTau( const Eref& e, vector< double > val )
{
- if ( val.size() != 5 ) {
- cout << "Error: HHGate::setTau on " << e.id().path() <<
- ": Number of entries on argument vector should be 5, was " <<
- val.size() << endl;
- return;
- }
- if ( checkOriginal( e.id(), "tau" ) ) {
- tau_ = val;
- updateAlphaBeta();
- updateTables();
- }
+ if ( val.size() != 5 ) {
+ cout << "Error: HHGate::setTau on " << e.id().path() <<
+ ": Number of entries on argument vector should be 5, was " <<
+ val.size() << endl;
+ return;
+ }
+ if ( checkOriginal( e.id(), "tau" ) ) {
+ tau_ = val;
+ updateAlphaBeta();
+ updateTables();
+ }
}
vector< double > HHGate::getMinfinity( const Eref& e) const
{
- return mInfinity_;
+ return mInfinity_;
}
void HHGate::setMinfinity( const Eref& e, vector< double > val )
{
- if ( val.size() != 5 ) {
- cout << "Error: HHGate::setMinfinity on " << e.id().path() <<
- ": Number of entries on argument vector should be 5, was " <<
- val.size() << endl;
- return;
- }
- if ( checkOriginal( e.id(), "mInfinity" ) ) {
- mInfinity_ = val;
- updateAlphaBeta();
- updateTables();
- }
+ if ( val.size() != 5 ) {
+ cout << "Error: HHGate::setMinfinity on " << e.id().path() <<
+ ": Number of entries on argument vector should be 5, was " <<
+ val.size() << endl;
+ return;
+ }
+ if ( checkOriginal( e.id(), "mInfinity" ) ) {
+ mInfinity_ = val;
+ updateAlphaBeta();
+ updateTables();
+ }
}
double HHGate::getMin( const Eref& e) const
{
- return xmin_;
+ return xmin_;
}
void HHGate::setMin( const Eref& e, double val )
{
- if ( checkOriginal( e.id(), "min" ) ) {
- xmin_ = val;
- unsigned int xdivs = A_.size() - 1;
- if ( isDirectTable_ && xdivs > 0 ) {
- // Stuff here to stretch out table using interpolation.
- invDx_ = static_cast< double >( xdivs ) / ( xmax_ - val );
- tabFill( A_, xdivs, val, xmax_ );
- tabFill( B_, xdivs, val, xmax_ );
- } else {
- updateTables();
- }
- }
+ if ( checkOriginal( e.id(), "min" ) ) {
+ xmin_ = val;
+ unsigned int xdivs = A_.size() - 1;
+ if ( isDirectTable_ && xdivs > 0 ) {
+ // Stuff here to stretch out table using interpolation.
+ invDx_ = static_cast< double >( xdivs ) / ( xmax_ - val );
+ tabFill( A_, xdivs, val, xmax_ );
+ tabFill( B_, xdivs, val, xmax_ );
+ } else {
+ updateTables();
+ }
+ }
}
double HHGate::getMax( const Eref& e) const
{
- return xmax_;
+ return xmax_;
}
void HHGate::setMax( const Eref& e, double val )
{
- if ( checkOriginal( e.id(), "max" ) ) {
- xmax_ = val;
- unsigned int xdivs = A_.size() - 1;
- if ( isDirectTable_ && xdivs > 0 ) {
- // Set up using direct assignment of table values.
- invDx_ = static_cast< double >( xdivs ) / ( val - xmin_ );
- tabFill( A_, xdivs, xmin_, val );
- tabFill( B_, xdivs, xmin_, val );
- } else {
- // Set up using functional form. here we just recalculate.
- updateTables();
- }
- }
+ if ( checkOriginal( e.id(), "max" ) ) {
+ xmax_ = val;
+ unsigned int xdivs = A_.size() - 1;
+ if ( isDirectTable_ && xdivs > 0 ) {
+ // Set up using direct assignment of table values.
+ invDx_ = static_cast< double >( xdivs ) / ( val - xmin_ );
+ tabFill( A_, xdivs, xmin_, val );
+ tabFill( B_, xdivs, xmin_, val );
+ } else {
+ // Set up using functional form. here we just recalculate.
+ updateTables();
+ }
+ }
}
unsigned int HHGate::getDivs( const Eref& e) const
{
- return A_.size() - 1;
+ return A_.size() - 1;
}
void HHGate::setDivs( const Eref& e, unsigned int val )
{
- if ( checkOriginal( e.id(), "divs" ) ) {
- if ( isDirectTable_ ) {
- invDx_ = static_cast< double >( val ) / ( xmax_ - xmin_ );
- tabFill( A_, val, xmin_, xmax_ );
- tabFill( B_, val, xmin_, xmax_ );
- } else {
- /// Stuff here to redo sizes.
- A_.resize( val + 1 );
- B_.resize( val + 1 );
- invDx_ = static_cast< double >( val ) / ( xmax_ - xmin_ );
- updateTables();
- }
- }
+ if ( checkOriginal( e.id(), "divs" ) ) {
+ if ( isDirectTable_ ) {
+ invDx_ = static_cast< double >( val ) / ( xmax_ - xmin_ );
+ tabFill( A_, val, xmin_, xmax_ );
+ tabFill( B_, val, xmin_, xmax_ );
+ } else {
+ /// Stuff here to redo sizes.
+ A_.resize( val + 1 );
+ B_.resize( val + 1 );
+ invDx_ = static_cast< double >( val ) / ( xmax_ - xmin_ );
+ updateTables();
+ }
+ }
}
vector< double > HHGate::getTableA( const Eref& e) const
{
- return A_;
+ return A_;
}
void HHGate::setTableA( const Eref& e, vector< double > v )
{
- if ( v.size() < 2 ) {
- cout << "Warning: HHGate::setTableA: size must be >= 2 entries on "
- << e.id().path() << endl;
- return;
- }
- if ( checkOriginal( e.id(), "tableA" ) ) {
- isDirectTable_ = 1;
- A_ = v;
- unsigned int xdivs = A_.size() - 1;
- invDx_ = static_cast< double >( xdivs ) / ( xmax_ - xmin_ );
- }
+ if ( v.size() < 2 ) {
+ cout << "Warning: HHGate::setTableA: size must be >= 2 entries on "
+ << e.id().path() << endl;
+ return;
+ }
+ if ( checkOriginal( e.id(), "tableA" ) ) {
+ isDirectTable_ = 1;
+ A_ = v;
+ unsigned int xdivs = A_.size() - 1;
+ invDx_ = static_cast< double >( xdivs ) / ( xmax_ - xmin_ );
+ }
}
vector< double > HHGate::getTableB( const Eref& e) const
{
- return B_;
+ return B_;
}
void HHGate::setTableB( const Eref& e, vector< double > v )
{
- if ( checkOriginal( e.id(), "tableB" ) ) {
- isDirectTable_ = 1;
- if ( A_.size() != v.size() ) {
- cout << "Warning: HHGate::setTableB: size should be same as table A: " << v.size() << " != " << A_.size() << ". Ignoring.\n";
- return;
- }
- B_ = v;
- }
+ if ( checkOriginal( e.id(), "tableB" ) ) {
+ isDirectTable_ = 1;
+ if ( A_.size() != v.size() ) {
+ cout << "Warning: HHGate::setTableB: size should be same as table A: " << v.size() << " != " << A_.size() << ". Ignoring.\n";
+ return;
+ }
+ B_ = v;
+ }
}
bool HHGate::getUseInterpolation( const Eref& e) const
{
- return lookupByInterpolation_;
+ return lookupByInterpolation_;
}
void HHGate::setUseInterpolation( const Eref& e, bool val )
{
- if ( checkOriginal( e.id(), "useInterpolation" ) )
- lookupByInterpolation_ = val;
+ if ( checkOriginal( e.id(), "useInterpolation" ) )
+ lookupByInterpolation_ = val;
}
void HHGate::setupAlpha( const Eref& e,
- vector< double > parms )
-{
- if ( checkOriginal( e.id(), "setupAlpha" ) ) {
- if ( parms.size() != 13 ) {
- cout << "HHGate::setupAlpha: Error: parms.size() != 13\n";
- return;
- }
- setupTables( parms, false );
- alpha_.resize( 5, 0 );
- beta_.resize( 5, 0 );
- for ( unsigned int i = 0; i < 5; ++i )
- alpha_[i] = parms[i];
- for ( unsigned int i = 5; i < 10; ++i )
- beta_[i - 5] = parms[i];
- }
+ vector< double > parms )
+{
+ if ( checkOriginal( e.id(), "setupAlpha" ) ) {
+ if ( parms.size() != 13 ) {
+ cout << "HHGate::setupAlpha: Error: parms.size() != 13\n";
+ return;
+ }
+ setupTables( parms, false );
+ alpha_.resize( 5, 0 );
+ beta_.resize( 5, 0 );
+ for ( unsigned int i = 0; i < 5; ++i )
+ alpha_[i] = parms[i];
+ for ( unsigned int i = 5; i < 10; ++i )
+ beta_[i - 5] = parms[i];
+ }
}
vector< double > HHGate::getAlphaParms( const Eref& e ) const
{
- vector< double > ret = alpha_;
- ret.insert( ret.end(), beta_.begin(), beta_.end() );
- ret.push_back( A_.size() );
- ret.push_back( xmin_ );
- ret.push_back( xmax_ );
+ vector< double > ret = alpha_;
+ ret.insert( ret.end(), beta_.begin(), beta_.end() );
+ ret.push_back( A_.size() );
+ ret.push_back( xmin_ );
+ ret.push_back( xmax_ );
- return ret;
+ return ret;
}
///////////////////////////////////////////////////
@@ -522,25 +522,25 @@ vector< double > HHGate::getAlphaParms( const Eref& e ) const
///////////////////////////////////////////////////
void HHGate::setupTau( const Eref& e,
- vector< double > parms )
+ vector< double > parms )
{
- if ( checkOriginal( e.id(), "setupTau" ) ) {
- if ( parms.size() != 13 ) {
- cout << "HHGate::setupTau: Error: parms.size() != 13\n";
- return;
- }
- setupTables( parms, true );
- }
+ if ( checkOriginal( e.id(), "setupTau" ) ) {
+ if ( parms.size() != 13 ) {
+ cout << "HHGate::setupTau: Error: parms.size() != 13\n";
+ return;
+ }
+ setupTables( parms, true );
+ }
}
void HHGate::tweakAlpha()
{
- ; // Dummy
+ ; // Dummy
}
void HHGate::tweakTau()
{
- ; // Dummy
+ ; // Dummy
}
/**
@@ -550,96 +550,96 @@ void HHGate::tweakTau()
*/
void HHGate::setupTables( const vector< double >& parms, bool doTau )
{
- assert( parms.size() == 13 );
- static const int XDIVS = 10;
- static const int XMIN = 11;
- static const int XMAX = 12;
- if ( parms[XDIVS] < 1 ) return;
- unsigned int xdivs = static_cast< unsigned int >( parms[XDIVS] );
-
- A_.resize( xdivs + 1 );
- B_.resize( xdivs + 1 );
- xmin_ = parms[XMIN];
- xmax_ = parms[XMAX];
- assert( xmax_ > xmin_ );
- invDx_ = xdivs / (xmax_ - xmin_ );
- double dx = ( xmax_ - xmin_ ) / xdivs;
-
- double x = xmin_;
- double prevAentry = 0.0;
- double prevBentry = 0.0;
- double temp;
- double temp2 = 0.0;
- unsigned int i;
-
- for( i = 0; i <= xdivs; i++ ) {
- if ( fabs( parms[4] ) < SINGULARITY ) {
- temp = 0.0;
- A_[i] = temp;
- } else {
- temp2 = parms[2] + exp( ( x + parms[3] ) / parms[4] );
- if ( fabs( temp2 ) < SINGULARITY ) {
- temp2 = parms[2] + exp( ( x + dx/10.0 + parms[3] ) / parms[4] );
- temp = ( parms[0] + parms[1] * (x + dx/10 ) ) / temp2;
-
- temp2 = parms[2] + exp( ( x - dx/10.0 + parms[3] ) / parms[4] );
- temp += ( parms[0] + parms[1] * (x - dx/10 ) ) / temp2;
- temp /= 2.0;
- // cout << "interpolated temp = " << temp << ", prev = " << prevAentry << endl;
-
- // temp = prevAentry;
- A_[i] = temp;
- } else {
- temp = ( parms[0] + parms[1] * x) / temp2;
- A_[i] = temp;
- }
- }
- if ( fabs( parms[9] ) < SINGULARITY ) {
- B_[i] = 0.0;
- } else {
- temp2 = parms[7] + exp( ( x + parms[8] ) / parms[9] );
- if ( fabs( temp2 ) < SINGULARITY ) {
- temp2 = parms[7] + exp( ( x + dx/10.0 + parms[8] ) / parms[9] );
- temp = (parms[5] + parms[6] * (x + dx/10) ) / temp2;
- temp2 = parms[7] + exp( ( x - dx/10.0 + parms[8] ) / parms[9] );
- temp += (parms[5] + parms[6] * (x - dx/10) ) / temp2;
- temp /= 2.0;
- B_[i] = temp;
- // B_[i] = prevBentry;
- } else {
- B_[i] = (parms[5] + parms[6] * x ) / temp2;
- // B_.table_[i] = ( parms[5] + parms[6] * x ) / temp2;
- }
- }
- // There are cleaner ways to do this, but this keeps
- // the relation to the GENESIS version clearer.
- // Note the additional SINGULARITY check, to fix a bug
- // in the earlier code.
- if ( doTau == 0 && fabs( temp2 ) > SINGULARITY )
- B_[i] += temp;
-
- prevAentry = A_[i];
- prevBentry = B_[i];
- x += dx;
- }
-
- prevAentry = 0.0;
- prevBentry = 0.0;
- if ( doTau ) {
- for( i = 0; i <= xdivs; i++ ) {
- temp = A_[i];
- temp2 = B_[i];
- if ( fabs( temp ) < SINGULARITY ) {
- A_[i] = prevAentry;
- B_[i] = prevBentry;
- } else {
- A_[i] = temp2 / temp;
- B_[i] = 1.0 / temp;
- }
- prevAentry = A_[i];
- prevBentry = B_[i];
- }
- }
+ assert( parms.size() == 13 );
+ static const int XDIVS = 10;
+ static const int XMIN = 11;
+ static const int XMAX = 12;
+ if ( parms[XDIVS] < 1 ) return;
+ unsigned int xdivs = static_cast< unsigned int >( parms[XDIVS] );
+
+ A_.resize( xdivs + 1 );
+ B_.resize( xdivs + 1 );
+ xmin_ = parms[XMIN];
+ xmax_ = parms[XMAX];
+ assert( xmax_ > xmin_ );
+ invDx_ = xdivs / (xmax_ - xmin_ );
+ double dx = ( xmax_ - xmin_ ) / xdivs;
+
+ double x = xmin_;
+ double prevAentry = 0.0;
+ double prevBentry = 0.0;
+ double temp;
+ double temp2 = 0.0;
+ unsigned int i;
+
+ for( i = 0; i <= xdivs; i++ ) {
+ if ( fabs( parms[4] ) < SINGULARITY ) {
+ temp = 0.0;
+ A_[i] = temp;
+ } else {
+ temp2 = parms[2] + exp( ( x + parms[3] ) / parms[4] );
+ if ( fabs( temp2 ) < SINGULARITY ) {
+ temp2 = parms[2] + exp( ( x + dx/10.0 + parms[3] ) / parms[4] );
+ temp = ( parms[0] + parms[1] * (x + dx/10 ) ) / temp2;
+
+ temp2 = parms[2] + exp( ( x - dx/10.0 + parms[3] ) / parms[4] );
+ temp += ( parms[0] + parms[1] * (x - dx/10 ) ) / temp2;
+ temp /= 2.0;
+ // cout << "interpolated temp = " << temp << ", prev = " << prevAentry << endl;
+
+ // temp = prevAentry;
+ A_[i] = temp;
+ } else {
+ temp = ( parms[0] + parms[1] * x) / temp2;
+ A_[i] = temp;
+ }
+ }
+ if ( fabs( parms[9] ) < SINGULARITY ) {
+ B_[i] = 0.0;
+ } else {
+ temp2 = parms[7] + exp( ( x + parms[8] ) / parms[9] );
+ if ( fabs( temp2 ) < SINGULARITY ) {
+ temp2 = parms[7] + exp( ( x + dx/10.0 + parms[8] ) / parms[9] );
+ temp = (parms[5] + parms[6] * (x + dx/10) ) / temp2;
+ temp2 = parms[7] + exp( ( x - dx/10.0 + parms[8] ) / parms[9] );
+ temp += (parms[5] + parms[6] * (x - dx/10) ) / temp2;
+ temp /= 2.0;
+ B_[i] = temp;
+ // B_[i] = prevBentry;
+ } else {
+ B_[i] = (parms[5] + parms[6] * x ) / temp2;
+ // B_.table_[i] = ( parms[5] + parms[6] * x ) / temp2;
+ }
+ }
+ // There are cleaner ways to do this, but this keeps
+ // the relation to the GENESIS version clearer.
+ // Note the additional SINGULARITY check, to fix a bug
+ // in the earlier code.
+ if ( doTau == 0 && fabs( temp2 ) > SINGULARITY )
+ B_[i] += temp;
+
+ prevAentry = A_[i];
+ prevBentry = B_[i];
+ x += dx;
+ }
+
+ prevAentry = 0.0;
+ prevBentry = 0.0;
+ if ( doTau ) {
+ for( i = 0; i <= xdivs; i++ ) {
+ temp = A_[i];
+ temp2 = B_[i];
+ if ( fabs( temp ) < SINGULARITY ) {
+ A_[i] = prevAentry;
+ B_[i] = prevBentry;
+ } else {
+ A_[i] = temp2 / temp;
+ B_[i] = 1.0 / temp;
+ }
+ prevAentry = A_[i];
+ prevBentry = B_[i];
+ }
+ }
}
/**
@@ -649,97 +649,97 @@ void HHGate::setupTables( const vector< double >& parms, bool doTau )
*/
void HHGate::tweakTables( bool doTau )
{
- unsigned int i;
- unsigned int size = A_.size();
- assert( size == B_.size() );
- if ( doTau ) {
- for ( i = 0; i < size; i++ ) {
- double temp = A_[ i ];
- double temp2 = B_[ i ];
- if ( fabs( temp ) < SINGULARITY ) {
- if ( temp < 0.0 )
- temp = -SINGULARITY;
- else
- temp = SINGULARITY;
- }
- A_[ i ] = temp2 / temp;
- B_[ i ] = 1.0 / temp;
- }
- } else {
- for ( i = 0; i < size; i++ )
- B_[i] = A_[i] + B_[i];
- }
+ unsigned int i;
+ unsigned int size = A_.size();
+ assert( size == B_.size() );
+ if ( doTau ) {
+ for ( i = 0; i < size; i++ ) {
+ double temp = A_[ i ];
+ double temp2 = B_[ i ];
+ if ( fabs( temp ) < SINGULARITY ) {
+ if ( temp < 0.0 )
+ temp = -SINGULARITY;
+ else
+ temp = SINGULARITY;
+ }
+ A_[ i ] = temp2 / temp;
+ B_[ i ] = 1.0 / temp;
+ }
+ } else {
+ for ( i = 0; i < size; i++ )
+ B_[i] = A_[i] + B_[i];
+ }
}
void HHGate::setupGate( const Eref& e,
- vector< double > parms )
-{
- // The nine arguments are :
- // A B C D F size min max isbeta
- // If size == 0 then we check that the gate has already been allocated.
- // If isbeta is true then we also have to do the conversion to
- // HHGate form of alpha, alpha+beta, assuming that the alpha gate
- // has already been setup. This uses tweakTables.
- // We may need to resize the tables if they don't match here.
- if ( !checkOriginal( e.id(), "setupGate" ) )
- return;
-
- if ( parms.size() != 9 ) {
- cout << "HHGate::setupGate: Error: parms.size() != 9\n";
- return;
- }
-
- double A = parms[0];
- double B = parms[1];
- double C = parms[2];
- double D = parms[3];
- double F = parms[4];
- int size = static_cast< int > (parms[5] );
- double min = parms[6];
- double max = parms[7];
- bool isBeta = static_cast< bool >( parms[8] );
-
- vector< double >& ip = ( isBeta ) ? B_ : A_;
-
- if ( size <= 0 ) { // Look up size, min, max from the interpol
- size = ip.size() - 1;
- if ( size <= 0 ) {
- cout << "Error: setupGate has zero size\n";
- return;
- }
- } else {
- ip.resize( size + 1 );
- }
-
- double dx = ( max - min ) / static_cast< double >( size );
- double x = min + dx / 2.0;
- for ( int i = 0; i <= size; i++ ) {
- if ( fabs ( F ) < SINGULARITY ) {
- ip[i] = 0.0;
- } else {
- double temp2 = C + exp( ( x + D ) / F );
- if ( fabs( temp2 ) < SINGULARITY )
- ip[i] = ip[i-1];
- else
- ip[i] = ( A + B * x ) / temp2;
- }
- }
-
- if ( isBeta ) {
- assert( A_.size() > 0 );
- // Here we ensure that the tables are the same size
- if ( A_.size() != B_.size() ) {
- if ( A_.size() > B_.size() ) {
- // Note that the tabFill expects to allocate the
- // terminating entry, so we put in size - 1.
- tabFill( B_, A_.size() - 1, xmin_, xmax_ );
- } else {
- tabFill( A_, B_.size() - 1, xmin_, xmax_ );
- }
- }
- // Then we do the tweaking to convert to HHChannel form.
- tweakTables( 0 );
- }
+ vector< double > parms )
+{
+ // The nine arguments are :
+ // A B C D F size min max isbeta
+ // If size == 0 then we check that the gate has already been allocated.
+ // If isbeta is true then we also have to do the conversion to
+ // HHGate form of alpha, alpha+beta, assuming that the alpha gate
+ // has already been setup. This uses tweakTables.
+ // We may need to resize the tables if they don't match here.
+ if ( !checkOriginal( e.id(), "setupGate" ) )
+ return;
+
+ if ( parms.size() != 9 ) {
+ cout << "HHGate::setupGate: Error: parms.size() != 9\n";
+ return;
+ }
+
+ double A = parms[0];
+ double B = parms[1];
+ double C = parms[2];
+ double D = parms[3];
+ double F = parms[4];
+ int size = static_cast< int > (parms[5] );
+ double min = parms[6];
+ double max = parms[7];
+ bool isBeta = static_cast< bool >( parms[8] );
+
+ vector< double >& ip = ( isBeta ) ? B_ : A_;
+
+ if ( size <= 0 ) { // Look up size, min, max from the interpol
+ size = ip.size() - 1;
+ if ( size <= 0 ) {
+ cout << "Error: setupGate has zero size\n";
+ return;
+ }
+ } else {
+ ip.resize( size + 1 );
+ }
+
+ double dx = ( max - min ) / static_cast< double >( size );
+ double x = min + dx / 2.0;
+ for ( int i = 0; i <= size; i++ ) {
+ if ( fabs ( F ) < SINGULARITY ) {
+ ip[i] = 0.0;
+ } else {
+ double temp2 = C + exp( ( x + D ) / F );
+ if ( fabs( temp2 ) < SINGULARITY )
+ ip[i] = ip[i-1];
+ else
+ ip[i] = ( A + B * x ) / temp2;
+ }
+ }
+
+ if ( isBeta ) {
+ assert( A_.size() > 0 );
+ // Here we ensure that the tables are the same size
+ if ( A_.size() != B_.size() ) {
+ if ( A_.size() > B_.size() ) {
+ // Note that the tabFill expects to allocate the
+ // terminating entry, so we put in size - 1.
+ tabFill( B_, A_.size() - 1, xmin_, xmax_ );
+ } else {
+ tabFill( A_, B_.size() - 1, xmin_, xmax_ );
+ }
+ }
+ // Then we do the tweaking to convert to HHChannel form.
+ tweakTables( 0 );
+ }
}
///////////////////////////////////////////////////////////////////////
@@ -753,54 +753,54 @@ void HHGate::setupGate( const Eref& e,
* subdivisions that the table represents.
*/
void HHGate::tabFill( vector< double >& table,
- unsigned int newXdivs, double newXmin, double newXmax )
-{
- if ( newXdivs < 3 ) {
- cout << "Error: tabFill: # divs must be >= 3. Not filling table.\n";
- return;
- }
-
- vector< double > old = table;
- double newDx = ( newXmax - newXmin ) / newXdivs;
- table.resize( newXdivs + 1 );
- bool origLookupMode = lookupByInterpolation_;
- lookupByInterpolation_ = 1;
-
- for( unsigned int i = 0; i <= newXdivs; ++i ) {
- table[i] = lookupTable( table, newXmin + i * newDx );
- }
+ unsigned int newXdivs, double newXmin, double newXmax )
+{
+ if ( newXdivs < 3 ) {
+ cout << "Error: tabFill: # divs must be >= 3. Not filling table.\n";
+ return;
+ }
+
+ vector< double > old = table;
+ double newDx = ( newXmax - newXmin ) / newXdivs;
+ table.resize( newXdivs + 1 );
+ bool origLookupMode = lookupByInterpolation_;
+ lookupByInterpolation_ = 1;
+
+ for( unsigned int i = 0; i <= newXdivs; ++i ) {
+ table[i] = lookupTable( table, newXmin + i * newDx );
+ }
- lookupByInterpolation_ = origLookupMode;
+ lookupByInterpolation_ = origLookupMode;
}
bool HHGate::checkOriginal( Id id, const string& field ) const
{
- if ( id == originalGateId_ )
- return 1;
+ if ( id == originalGateId_ )
+ return 1;
- cout << "Warning: HHGate: attempt to set field '" << field << "' on " <<
- id.path() << "\nwhich is not the original Gate element. Ignored.\n";
- return 0;
+ cout << "Warning: HHGate: attempt to set field '" << field << "' on " <<
+ id.path() << "\nwhich is not the original Gate element. Ignored.\n";
+ return 0;
}
bool HHGate::isOriginalChannel( Id id ) const
{
- return ( id == originalChanId_ );
+ return ( id == originalChanId_ );
}
bool HHGate::isOriginalGate( Id id ) const
{
- return ( id == originalGateId_ );
+ return ( id == originalGateId_ );
}
Id HHGate::originalChannelId() const
{
- return originalChanId_;
+ return originalChanId_;
}
Id HHGate::originalGateId() const
{
- return originalGateId_;
+ return originalGateId_;
}
void HHGate::updateAlphaBeta()
@@ -813,13 +813,13 @@ void HHGate::updateTauMinf()
void HHGate::updateTables()
{
- if ( alpha_.size() == 0 || beta_.size() == 0 )
- return;
- vector< double > parms = alpha_;
- parms.insert( parms.end(), beta_.begin(), beta_.end() );
- parms.push_back( A_.size() );
- parms.push_back( xmin_ );
- parms.push_back( xmax_ );
+ if ( alpha_.size() == 0 || beta_.size() == 0 )
+ return;
+ vector< double > parms = alpha_;
+ parms.insert( parms.end(), beta_.begin(), beta_.end() );
+ parms.push_back( A_.size() );
+ parms.push_back( xmin_ );
+ parms.push_back( xmax_ );
- setupTables( parms, 0 );
+ setupTables( parms, 0 );
}
diff --git a/moose-core/biophysics/MMPump.cpp b/moose-core/biophysics/MMPump.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f9ed90ef61899052687647517323cfe0a1e51115
--- /dev/null
+++ b/moose-core/biophysics/MMPump.cpp
@@ -0,0 +1,141 @@
+#include "header.h"
+#include "MMPump.h"
+#include "ElementValueFinfo.h"
+#include "../utility/numutil.h"
+#include <cmath>
+
+SrcFinfo2< double,double >* MMPump::PumpOut()
+{
+ static SrcFinfo2< double,double > pumpOut( "PumpOut",
+ "Sends out MMPump parameters.");
+ return &pumpOut;
+}
+
+const Cinfo * MMPump::initCinfo()
+{
+ static DestFinfo process( "process",
+ "Handles process call",
+ new ProcOpFunc< MMPump>(&MMPump::process ) );
+ static DestFinfo reinit( "reinit",
+ "Reinit happens only in stage 0",
+ new ProcOpFunc< MMPump>( &MMPump::reinit ));
+
+ static Finfo* processShared[] = {
+ &process, &reinit
+ };
+
+ static SharedFinfo proc(
+ "proc",
+ "Shared message to receive Process message from scheduler",
+ processShared, sizeof( processShared ) / sizeof( Finfo* ));
+
+
+ ////////////////////////////
+ // Field defs
+ ////////////////////////////
+
+
+
+ static ElementValueFinfo<MMPump, double> Vmax("Vmax",
+ "maximum pump velocity, scaled by mebrane"
+ "surface area. i.e., max ion flux in moles/sec",
+ &MMPump::setVmax,
+ &MMPump::getVmax);
+
+ static ElementValueFinfo<MMPump, double> Kd("Kd",
+ "half-maximal activating concentration in mM",
+ &MMPump::setKd,
+ &MMPump::getKd);
+
+ ////
+ // DestFinfo
+ ////
+ static Finfo * difMMPumpFinfos[] = {
+ //////////////////////////////////////////////////////////////////
+ // Field definitions
+ //////////////////////////////////////////////////////////////////
+
+
+ &Vmax,
+ &Kd,
+ //////////////////////////////////////////////////////////////////
+ // SharedFinfo definitions
+ /////////////////////////////////////////////////////////////////
+
+ &proc,
+ PumpOut(),
+ //////////////////////////////////////////////////////////////////
+ // DestFinfo definitions
+ //////////////////////////////////////////////////////////////////
+
+ };
+
+ static string doc[] = {
+ "Name", "MMPump",
+ "Author", "Subhasis Ray (ported from GENESIS2)",
+ "Description", "Models Michaelis-Menten pump. It is coupled with a DifShell.",
+ };
+ static ZeroSizeDinfo<int> dinfo;
+ static Cinfo MMPumpCinfo(
+ "MMPump",
+ Neutral::initCinfo(),
+ difMMPumpFinfos,
+ sizeof(difMMPumpFinfos)/sizeof(Finfo*),
+ &dinfo,
+ doc,
+ sizeof(doc)/sizeof(string));
+
+ return &MMPumpCinfo;
+}
+
+static const Cinfo * MMpumpCinfo = MMPump::initCinfo();
+
+
+////////////////////////////////////////////////////////////////////////////////
+// Class functions
+////////////////////////////////////////////////////////////////////////////////
+
+MMPump::MMPump()
+{ ; }
+
+
+
+
+double MMPump::getVmax(const Eref& e) const
+{
+
+ return Vmax_;
+}
+void MMPump::setVmax(const Eref& e,double value)
+{
+if ( value < 0.0 ) {
+ cerr << "Error: MMPump: Vmax cannot be negative!\n";
+ return;
+ }
+ Vmax_ = value;
+}
+
+
+double MMPump::getKd(const Eref& e) const
+{
+ return Kd_;
+}
+void MMPump::setKd(const Eref& e,double value)
+{
+if ( value < 0.0 ) {
+ cerr << "Error: MMPump: Kd cannot be negative!\n";
+ return;
+ }
+ Kd_ = value;
+}
+
+
+void MMPump::process(const Eref& e, ProcPtr p)
+{
+ PumpOut()->send(e,Vmax_,Kd_);
+}
+
+void MMPump::reinit(const Eref& e, ProcPtr p)
+{
+ return;
+}
diff --git a/moose-core/biophysics/MMPump.h b/moose-core/biophysics/MMPump.h
new file mode 100644
index 0000000000000000000000000000000000000000..911fb1d4b118cb49e4b1b1a3ab5d3fe06ada55c7
--- /dev/null
+++ b/moose-core/biophysics/MMPump.h
@@ -0,0 +1,45 @@
+/**********************************************************************
+ ** This program is part of 'MOOSE', the
+ ** Multiscale Object Oriented Simulation Environment.
+ ** copyright (C) 2003-2008
+ ** Upinder S. Bhalla, Niraj Dudani and NCBS
+ ** It is made available under the terms of the
+ ** GNU Lesser General Public License version 2.1
+ ** See the file COPYING.LIB for the full notice.
+ **********************************************************************/
+
+#ifndef _MMPUMP_H
+#define _MMPUMP_H
+/*This is base class for MMPump*/
+class MMPump
+{
+public:
+ MMPump();
+ void process( const Eref & e, ProcPtr p);
+ void reinit(const Eref &e, ProcPtr p);
+
+ double getVmax(const Eref& e) const; //Vmax of the pump
+ void setVmax(const Eref& e,double value);
+
+ int getVal(const Eref& e) const; //Valence
+ void setVal(const Eref& e,int value);
+
+ double getKd(const Eref& e) const; // Pump's Kd
+ void setKd(const Eref& e,double value);
+
+
+
+ static SrcFinfo2< double,double >* PumpOut(); //Pump parameters;
+
+ static const Cinfo * initCinfo();
+
+private:
+
+ double Vmax_;
+ double Kd_;
+};
+
+
+
+
+#endif //_MMPUMP_BASE_H
diff --git a/moose-core/biophysics/Makefile b/moose-core/biophysics/Makefile
index ae8391efc9bb97a14a709d0f892885c0ddbda182..bdd2a109a84b5f517f7c07a01c159ce971d3dce4 100644
--- a/moose-core/biophysics/Makefile
+++ b/moose-core/biophysics/Makefile
@@ -38,7 +38,11 @@ OBJ = \
NMDAChan.o \
testBiophysics.o \
IzhikevichNrn.o \
+ DifShellBase.o \
DifShell.o \
+ DifBufferBase.o \
+ DifBuffer.o \
+ MMpump.o \
Leakage.o \
VectorTable.o \
MarkovRateTable.o \
@@ -94,7 +98,11 @@ ReadCell.o: CompartmentBase.h Compartment.h SymCompartment.h ReadCell.h ../shell
SwcSegment.o: SwcSegment.h ../utility/Vec.h
ReadSwc.o: CompartmentBase.h Compartment.h SymCompartment.h SwcSegment.h ReadSwc.h ../shell/Shell.h ../utility/Vec.h
IzhikevichNrn.o: IzhikevichNrn.h
-DifShell.o: DifShell.h
+DifShellBase.o: DifShellBase.h
+DifShell.o: DifShellBase.h DifShell.h
+DifBufferBase.o: DifBufferBase.h
+DifBuffer.o: DifBufferBase.h DifBuffer.h
+MMPump.o: MMPump.h
testBiophysics.o: IntFire.h CompartmentBase.h Compartment.h HHChannel.h HHGate.h
VectorTable.o : VectorTable.h
MarkovGslSolver.o : MarkovGslSolver.h
diff --git a/moose-core/scheduling/Clock.cpp b/moose-core/scheduling/Clock.cpp
index e48b6a343506df5110c2dfc9320e7b5c5847dd99..434d2f62e94e54674310b3d363d1ef9a3733a6d3 100644
--- a/moose-core/scheduling/Clock.cpp
+++ b/moose-core/scheduling/Clock.cpp
@@ -358,6 +358,10 @@ const Cinfo* Clock::initCinfo()
" CaConc 1 50e-6\n"
" CaConcBase 1 50e-6\n"
" DifShell 1 50e-6\n"
+ " DifShellBase 1 50e-6\n"
+ " MMPump 1 50e-6\n"
+ " DifBuffer 1 50e-6\n"
+ " DifBufferBase 1 50e-6\n"
" MgBlock 1 50e-6\n"
" Nernst 1 50e-6\n"
" RandSpike 1 50e-6\n"
@@ -837,6 +841,10 @@ void Clock::buildDefaultTick()
defaultTick_["CaConc"] = 1;
defaultTick_["CaConcBase"] = 1;
defaultTick_["DifShell"] = 1;
+ defaultTick_["DifShellBase"] = 1;
+ defaultTick_["MMPump"] = 1;
+ defaultTick_["DifBuffer"] = 1;
+ defaultTick_["DifBufferBase"] = 1;
defaultTick_["MgBlock"] = 1;
defaultTick_["Nernst"] = 1;
defaultTick_["RandSpike"] = 1;