diff --git a/Main.cpp b/Main.cpp index d38e8c3..5346ed3 100644 --- a/Main.cpp +++ b/Main.cpp @@ -22,19 +22,18 @@ /****************************************************************************************************/ /* Main file for compilation tests */ -/* The Simulation requires the following boost libraries: Preprocessor */ -/* Random */ +/* The Simulation requires the following boost libraries: Random */ /****************************************************************************************************/ #include -#include +#include #include "Thalamic_Column.h" /****************************************************************************************************/ /* Fixed simulation settings */ /****************************************************************************************************/ +typedef std::chrono::high_resolution_clock::time_point timer; extern const int T = 30; /* Simulation length s */ extern const int res = 1E4; /* number of iteration steps per s */ -extern const int red = res/100; /* number of iterations that is saved */ extern const double dt = 1E3/res; /* duration of a timestep in ms */ extern const double h = sqrt(dt); /* squareroot of dt for SRK iteration */ /****************************************************************************************************/ @@ -43,43 +42,27 @@ extern const double h = sqrt(dt); /* squareroot of dt for SRK iteration */ /****************************************************************************************************/ -/* Constants for SRK4 iteration */ -/****************************************************************************************************/ -extern const vector B1 = {0, - 0.626708569400000081728308032325, - 1.7296310295000001389098542858846, - 1.2703689705000000831347506391467}; -extern const vector B2 = {0, - 0.78000033203198970710445792065002, - 1.28727807507536762265942797967, - 0.44477273249350995909523476257164}; -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - - -/****************************************************************************************************/ /* Main simulation routine */ /****************************************************************************************************/ int main(void) { /* Initialize the populations */ - Thalamic_Column Thalamus; + Thalamic_Column Thalamus = Thalamic_Column(); - /* Takes the time of the simulation */ - time_t start,end; - time (&start); + /* Take the time of the simulation */ + timer start,end; - /* Simulation */ - for (int t=0; t< T*res; ++t) { - Thalamus.iterate_ODE(); - } - time (&end); + /* Simulation */ + start = std::chrono::high_resolution_clock::now(); + for (int t=0; t< T*res; ++t) { + Thalamus.iterate_ODE(); + } + end = std::chrono::high_resolution_clock::now(); - /* Time consumed by the simulation */ - double dif = difftime(end,start); - std::cout << "simulation done!\n"; - std::cout << "took " << dif << " seconds" << "\n"; - std::cout << "end\n"; + /* Time consumed by the simulation */ + double dif = 1E-3*std::chrono::duration_cast( end - start ).count(); + std::cout << "simulation done!\n"; + std::cout << "took " << dif << " seconds" << "\n"; + std::cout << "end\n"; } /****************************************************************************************************/ /* end */ diff --git a/Plots.m b/Plots.m index 394a3e7..b8da657 100644 --- a/Plots.m +++ b/Plots.m @@ -4,10 +4,10 @@ function Plots(T) if nargin == 0 - Con = [ 0.062; % g_h + Con = [ 0.063; % g_h 0.02; % g_LK_t 3; % N_tr - 4; % N_rt + 5; % N_rt 30]; % N_rr @@ -25,20 +25,20 @@ function Plots(T) 1]; % time until stimuli after min in ms T = 30; % duration of the simulation end -[Vt, Vr] = Thalamus(T, Con, var_stim); +[Vt, Vr, ah] = Thalamus(T, Con, var_stim); L = max(size(Vt)); timeaxis = linspace(0,T,L); - -fs = L/T; -[Pxx,f] = pwelch(Vt-mean(Vt), [], [], [], fs,'onesided'); -n = find(f<=60, 1, 'last' ); +% +% fs = L/T; +% [Pxx,f] = pwelch(Vt-mean(Vt), 20*FS, 4*FS, [], fs,'onesided'); +% n = find(f<=60, 1, 'last' ); figure(1) subplot(311), plot(timeaxis,Vt) -title('Thalamic relay membrane voltage'), xlabel('time in s'), ylabel('Vt in mV') +title('Thalamic relay membrane voltage'), xlabel('time in s'), ylabel('Vt [mV]') subplot(312), plot(timeaxis,Vr) -title('Thalamic reticular membrane voltage'), xlabel('time in s'), ylabel('Vr in mV') -subplot(313), plot(f(1:n),Pxx(1:n)) -title('Powerspectrum of Steyn-Ross model with pwelch'), xlabel('frequency in Hz'), ylabel('Power') -save('Thalamus.mat','Vt','Vr') +title('Thalamic reticular membrane voltage'), xlabel('time in s'), ylabel('Vr [mV]') +subplot(313), plot(timeaxis,ah) +title('Thalamic relay I_h activation'), xlabel('time in s'), ylabel('m_h') +%save('Thalamus.mat','Vt','Vr') \ No newline at end of file diff --git a/T_model.pro b/T_model.pro new file mode 100644 index 0000000..2c79a62 --- /dev/null +++ b/T_model.pro @@ -0,0 +1,17 @@ +TEMPLATE = app +CONFIG += console +CONFIG -= app_bundle +CONFIG -= qt + +SOURCES += Main.cpp \ + Thalamic_Column.cpp \ + Thalamus.cpp + +HEADERS += ODE.h \ + saves.h \ + Thalamic_Column.h \ + Stimulation.h + +QMAKE_CXXFLAGS += -std=c++11 -O3 + +SOURCES -= Thalamus.cpp diff --git a/Thalamic_Column.cpp b/Thalamic_Column.cpp index e1a8f19..8b33f82 100644 --- a/Thalamic_Column.cpp +++ b/Thalamic_Column.cpp @@ -26,20 +26,33 @@ #include "Thalamic_Column.h" /****************************************************************************************************/ -/* Initialization of RNG */ +/* Parameters for SRK4 iteration */ +/****************************************************************************************************/ +extern const vector A = {0.5, 0.5, 1.0, 1.0}; +extern const vector B = {0.75, 0.75, 0.0, 0.0}; +/****************************************************************************************************/ +/* end */ /****************************************************************************************************/ -void Thalamic_Column::set_RNG(void) { - /* Number of independent streams */ - int N = 2; - /* Create RNG for each stream */ - for (int i=0; i B1, B2; - double n = 1 / h * (B1[N-1] * Rand_vars[0] + B2[N-1] * Rand_vars[1]); - return n; +double Thalamic_Column::noise_xRK(int N, int M) const{ + return gamma_e * gamma_e * (Rand_vars[2*M] + Rand_vars[2*M+1]/std::sqrt(3))*B[N]; +} + +double Thalamic_Column::noise_aRK(int M) const{ + return gamma_e * gamma_e * (Rand_vars[2*M] - Rand_vars[2*M+1]*std::sqrt(3))/4; } /****************************************************************************************************/ /* end */ @@ -253,25 +251,21 @@ double Thalamic_Column::noise_xRK(int N) const{ /****************************************************************************************************/ void Thalamic_Column::set_RK (int N) { extern const double dt; - _SWITCH((Ca) - (Phi_tt)(Phi_tr)(Phi_rt)(Phi_rr) - (x_tt) (x_tr) (x_rt) (x_rr) - (h_T_t) (h_T_r) (m_h) (m_h2)) - Vt [N] = dt*(-(I_L_t(N) + I_et(N) + I_it(N))/tau_t - (I_LK_t(N) + I_T_t(N) + I_h(N))); - Vr [N] = dt*(-(I_L_r(N) + I_er(N) + I_ir(N))/tau_r - (I_LK_r(N) + I_T_r(N))); - Ca [N] = dt*(alpha_Ca * I_T_t(N) - (var_Ca - Ca_0)/tau_Ca); - h_T_t [N] = dt*(h_inf_T_t(N) - var_h_T_t)/tau_h_T_t(N); - h_T_r [N] = dt*(h_inf_T_r(N) - var_h_T_r)/tau_h_T_r(N); - m_h [N] = dt*((m_inf_h(N) * (1 - var_m_h2) - var_m_h)/tau_m_h(N) - k3 * P_h(N) * var_m_h + k4 * var_m_h2); - m_h2 [N] = dt*(k3 * P_h(N) * var_m_h - k4 * var_m_h2); - Phi_tt [N] = dt*(var_x_tt); - Phi_tr [N] = dt*(var_x_tr); - Phi_rt [N] = dt*(var_x_rt); - Phi_rr [N] = dt*(var_x_rr); - x_tt [N] = dt*(pow(gamma_e, 2) * (noise_xRK(N) - var_Phi_tt) - 2 * gamma_e * var_x_tt); - x_tr [N] = dt*(pow(gamma_e, 2) * (N_tr * get_Qt(N) - var_Phi_tr) - 2 * gamma_e * var_x_tr); - x_rt [N] = dt*(pow(gamma_i, 2) * (N_rt * get_Qr(N) - var_Phi_rt) - 2 * gamma_i * var_x_rt); - x_rr [N] = dt*(pow(gamma_i, 2) * (N_rr * get_Qr(N) - var_Phi_rr) - 2 * gamma_i * var_x_rr); + Vt [N+1] = Vt [0] + A[N]*dt*(-(I_L_t(N) + I_et(N) + I_it(N))/tau_t - (I_LK_t(N) + I_T_t(N) + I_h(N))); + Vr [N+1] = Vr [0] + A[N]*dt*(-(I_L_r(N) + I_er(N) + I_ir(N))/tau_r - (I_LK_r(N) + I_T_r(N))); + Ca [N+1] = Ca [0] + A[N]*dt*(alpha_Ca * I_T_t(N) - (Ca[N] - Ca_0)/tau_Ca); + h_T_t [N+1] = h_T_t[0] + A[N]*dt*(h_inf_T_t(N) - h_T_t[N])/tau_h_T_t(N); + h_T_r [N+1] = h_T_r[0] + A[N]*dt*(h_inf_T_r(N) - h_T_r[N])/tau_h_T_r(N); + m_h [N+1] = m_h [0] + A[N]*dt*((m_inf_h(N) * (1 - m_h2[N]) - m_h[N])/tau_m_h(N) - k3 * P_h(N) * m_h[N] + k4 * m_h2[N]); + m_h2 [N+1] = m_h2 [0] + A[N]*dt*(k3 * P_h(N) * m_h[N] - k4 * m_h2[N]); + y_tt [N+1] = y_tt [0] + A[N]*dt*(x_tt[N]); + y_tr [N+1] = y_tr [0] + A[N]*dt*(x_tr[N]); + y_rt [N+1] = y_rt [0] + A[N]*dt*(x_rt[N]); + y_rr [N+1] = y_rr [0] + A[N]*dt*(x_rr[N]); + x_tt [N+1] = x_tt [0] + A[N]*dt*(pow(gamma_e, 2) * ( - y_tt[N]) - 2 * gamma_e * x_tt[N]) + noise_xRK(N,0); + x_tr [N+1] = x_tr [0] + A[N]*dt*(pow(gamma_e, 2) * (N_tr * get_Qt(N) - y_tr[N]) - 2 * gamma_e * x_tr[N]); + x_rt [N+1] = x_rt [0] + A[N]*dt*(pow(gamma_i, 2) * (N_rt * get_Qr(N) - y_rt[N]) - 2 * gamma_i * x_rt[N]); + x_rr [N+1] = x_rr [0] + A[N]*dt*(pow(gamma_i, 2) * (N_rr * get_Qr(N) - y_rr[N]) - 2 * gamma_i * x_rr[N]); } /****************************************************************************************************/ /* end */ @@ -282,22 +276,21 @@ void Thalamic_Column::set_RK (int N) { /* Function that adds all SRK terms */ /****************************************************************************************************/ void Thalamic_Column::add_RK(void) { - extern const double h; - Vt [0] += (Vt [1] + Vt [2] * 2 + Vt [3] * 2 + Vt [4])/6; - Vr [0] += (Vr [1] + Vr [2] * 2 + Vr [3] * 2 + Vr [4])/6; - Ca [0] += (Ca [1] + Ca [2] * 2 + Ca [3] * 2 + Ca [4])/6; - Phi_tt [0] += (Phi_tt [1] + Phi_tt [2] * 2 + Phi_tt [3] * 2 + Phi_tt [4])/6; - Phi_tr [0] += (Phi_tr [1] + Phi_tr [2] * 2 + Phi_tr [3] * 2 + Phi_tr [4])/6; - Phi_rt [0] += (Phi_rt [1] + Phi_rt [2] * 2 + Phi_rt [3] * 2 + Phi_rt [4])/6; - Phi_rr [0] += (Phi_rr [1] + Phi_rr [2] * 2 + Phi_rr [3] * 2 + Phi_rr [4])/6; - x_tt [0] += (x_tt [1] + x_tt [2] * 2 + x_tt [3] * 2 + x_tt [4])/6 + pow(gamma_e, 2) * h * Rand_vars[0]; - x_tr [0] += (x_tr [1] + x_tr [2] * 2 + x_tr [3] * 2 + x_tr [4])/6; - x_rt [0] += (x_rt [1] + x_rt [2] * 2 + x_rt [3] * 2 + x_rt [4])/6; - x_rr [0] += (x_rr [1] + x_rr [2] * 2 + x_rr [3] * 2 + x_rr [4])/6; - h_T_t [0] += (h_T_t [1] + h_T_t [2] * 2 + h_T_t [3] * 2 + h_T_t [4])/6; - h_T_r [0] += (h_T_r [1] + h_T_r [2] * 2 + h_T_r [3] * 2 + h_T_r [4])/6; - m_h [0] += (m_h [1] + m_h [2] * 2 + m_h [3] * 2 + m_h [4])/6; - m_h2 [0] += (m_h2 [1] + m_h2 [2] * 2 + m_h2 [3] * 2 + m_h2 [4])/6; + Vt [0] =(-3*Vt [0] + 2*Vt [1] + 4*Vt [2] + 2* Vt [3] + Vt [4])/6; + Vr [0] =(-3*Vr [0] + 2*Vr [1] + 4*Vr [2] + 2* Vr [3] + Vr [4])/6; + Ca [0] =(-3*Ca [0] + 2*Ca [1] + 4*Ca [2] + 2* Ca [3] + Ca [4])/6; + y_tt [0] =(-3*y_tt [0] + 2*y_tt [1] + 4*y_tt [2] + 2* y_tt [3] + y_tt [4])/6; + y_tr [0] =(-3*y_tr [0] + 2*y_tr [1] + 4*y_tr [2] + 2* y_tr [3] + y_tr [4])/6; + y_rt [0] =(-3*y_rt [0] + 2*y_rt [1] + 4*y_rt [2] + 2* y_rt [3] + y_rt [4])/6; + y_rr [0] =(-3*y_rr [0] + 2*y_rr [1] + 4*y_rr [2] + 2* y_rr [3] + y_rr [4])/6; + x_tt [0] =(-3*x_tt [0] + 2*x_tt [1] + 4*x_tt [2] + 2* x_tt [3] + x_tt [4])/6 + noise_aRK(0); + x_tr [0] =(-3*x_tr [0] + 2*x_tr [1] + 4*x_tr [2] + 2* x_tr [3] + x_tr [4])/6; + x_rt [0] =(-3*x_rt [0] + 2*x_rt [1] + 4*x_rt [2] + 2* x_rt [3] + x_rt [4])/6; + x_rr [0] =(-3*x_rr [0] + 2*x_rr [1] + 4*x_rr [2] + 2* x_rr [3] + x_rr [4])/6; + h_T_t [0] =(-3*h_T_t[0] + 2*h_T_t[1] + 4*h_T_t[2] + 2* h_T_t [3] + h_T_t [4])/6; + h_T_r [0] =(-3*h_T_r[0] + 2*h_T_r[1] + 4*h_T_r[2] + 2* h_T_r [3] + h_T_r [4])/6; + m_h [0] =(-3*m_h [0] + 2*m_h [1] + 4*m_h [2] + 2* m_h [3] + m_h [4])/6; + m_h2 [0] =(-3*m_h2 [0] + 2*m_h2 [1] + 4*m_h2 [2] + 2* m_h2 [3] + m_h2 [4])/6; /* Generate noise for the next iteration */ for (unsigned i=0; i #include #include -#include "macros.h" using std::vector; /****************************************************************************************************/ @@ -44,17 +43,28 @@ typedef boost::variate_generator GEN; /* Variate generator */ /****************************************************************************************************/ +/* Macro for vector initialization */ +/****************************************************************************************************/ +#ifndef _INIT +#define _INIT(x) {x, 0.0, 0.0, 0.0, 0.0} +#endif +/****************************************************************************************************/ +/* end */ +/****************************************************************************************************/ + + +/****************************************************************************************************/ /* Implementation of the thalamic module */ /****************************************************************************************************/ class Thalamic_Column { public: /* Constructors for compile check */ - Thalamic_Column(void) + Thalamic_Column(void) {set_RNG();} /* Constructor for simulation */ Thalamic_Column(double* Par) - : g_LK_t (Par[1]), g_LK_r (Par[1]), g_h (Par[0]), + : g_LK_t (Par[1]), g_LK_r (Par[1]), g_h (Par[0]), N_tr (Par[2]), N_rt (Par[3]), N_rr (Par[4]) {set_RNG();} @@ -80,6 +90,7 @@ class Thalamic_Column { double m_inf_h (int) const; double tau_m_h (int) const; double P_h (int) const; + double act_h (void)const; /* Deactivation functions */ double h_inf_T_t (int) const; @@ -97,7 +108,8 @@ class Thalamic_Column { double I_h (int) const; /* Noise functions */ - double noise_xRK (int) const; + double noise_xRK (int,int) const; + double noise_aRK (int) const; /* ODE functions */ void set_RK (int); @@ -105,21 +117,21 @@ class Thalamic_Column { void iterate_ODE (void); /* Data storage access */ - friend void get_data (int, Thalamic_Column&, _REPEAT(double*, 2)); + friend void get_data (int, Thalamic_Column&, double*, double*, double*); private: /* Population variables */ - vector Vt = _INIT(E_L_t), /* TC membrane voltage */ + vector Vt = _INIT(E_L_t), /* TC membrane voltage */ Vr = _INIT(E_L_r), /* RE membrane voltage */ Ca = _INIT(Ca_0), /* Calcium concentration of TC population */ - Phi_tt = _INIT(0.0), /* PostSP from TC population to TC population */ - Phi_tr = _INIT(0.0), /* PostSP from TC population to RE population */ - Phi_rt = _INIT(0.0), /* PostSP from RE population to TC population */ - Phi_rr = _INIT(0.0), /* PostSP from RE population to RE population */ - x_tt = _INIT(0.0), /* derivative of Phi_tt */ - x_tr = _INIT(0.0), /* derivative of Phi_tr */ - x_rt = _INIT(0.0), /* derivative of Phi_rt */ - x_rr = _INIT(0.0), /* derivative of Phi_rr */ + y_tt = _INIT(0.0), /* PostSP from TC population to TC population */ + y_tr = _INIT(0.0), /* PostSP from TC population to RE population */ + y_rt = _INIT(0.0), /* PostSP from RE population to TC population */ + y_rr = _INIT(0.0), /* PostSP from RE population to RE population */ + x_tt = _INIT(0.0), /* derivative of y_tt */ + x_tr = _INIT(0.0), /* derivative of y_tr */ + x_rt = _INIT(0.0), /* derivative of y_rt */ + x_rr = _INIT(0.0), /* derivative of y_rr */ h_T_t = _INIT(0.0), /* inactivation of T channel */ h_T_r = _INIT(0.0), /* inactivation of T channel */ m_h = _INIT(0.0), /* activation of h channel */ @@ -129,7 +141,7 @@ class Thalamic_Column { vector MTRands; /* Container for noise */ - vector Rand_vars; + vector Rand_vars; /* Declaration and Initialization of parameters */ /* Membrane time in ms */ @@ -169,7 +181,7 @@ class Thalamic_Column { const double g_T_r = 2.3; /* h current */ - const double g_h = 0.07; + const double g_h = 0.062; /* Reversal potentials in mV */ /* Synaptic */ @@ -190,9 +202,9 @@ class Thalamic_Column { const double E_h = -40; /* Calcium parameters */ - const double alpha_Ca = -50E-6; /* influx per spike in nmol */ + const double alpha_Ca = -51.8E-6; /* influx per spike in nmol */ const double tau_Ca = 10; /* calcium time constant in ms */ - const double Ca_0 = 2.4E-4; /* resting concentration */ + const double Ca_0 = 2.4E-4; /* resting concentration */ /* I_h activation parameters */ const double k1 = 2.5E7; @@ -204,12 +216,12 @@ class Thalamic_Column { /* Noise parameters in ms^-1 */ const double mphi = 0E-3; - const double dphi = 10E-3; + const double dphi = 120E-3; double input = 0.0; /* Connectivities (dimensionless) */ const double N_tr = 3; - const double N_rt = 3; + const double N_rt = 5; const double N_rr = 30; friend class Stim; diff --git a/Thalamus.cpp b/Thalamus.cpp index 797e17b..b096aaa 100644 --- a/Thalamus.cpp +++ b/Thalamus.cpp @@ -45,23 +45,6 @@ extern const double h = sqrt(dt); /* squareroot of dt for SRK iteration */ /* end */ /****************************************************************************************************/ - -/****************************************************************************************************/ -/* Constants for SRK4 iteration */ -/****************************************************************************************************/ -extern const vector B1 = {0, - 0.626708569400000081728308032325, - 1.7296310295000001389098542858846, - 1.2703689705000000831347506391467}; -extern const vector B2 = {0, - 0.78000033203198970710445792065002, - 1.28727807507536762265942797967, - 0.44477273249350995909523476257164}; -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - - /****************************************************************************************************/ /* Simulation routine */ /* lhs defines outputs */ @@ -86,10 +69,12 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { /* Create data containers */ mxArray* Vt = SetMexArray(1, T*res/red); mxArray* Vr = SetMexArray(1, T*res/red); + mxArray* ah = SetMexArray(1, T*res/red); /* Pointer to the actual data block */ double* Pr_Vt = mxGetPr(Vt); double* Pr_Vr = mxGetPr(Vr); + double* Pr_ah = mxGetPr(ah); /* Simulation */ int count = 0; @@ -97,7 +82,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { Thalamus.iterate_ODE(); Stimulation.check_stim(t); if(t>=onset*res && t%red==0){ - get_data(count, Thalamus, Pr_Vt, Pr_Vr); + get_data(count, Thalamus, Pr_Vt, Pr_Vr, Pr_ah); ++count; } } @@ -105,6 +90,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { /* Output of the simulation */ plhs[0] = Vt; plhs[1] = Vr; + plhs[2] = ah; return; } /****************************************************************************************************/ diff --git a/macros.h b/macros.h deleted file mode 100644 index 1baa7a1..0000000 --- a/macros.h +++ /dev/null @@ -1,81 +0,0 @@ -/* -* Copyright (c) 2014 Michael Schellenberger Costa -* -* Permission is hereby granted, free of charge, to any person obtaining a copy -* of this software and associated documentation files (the "Software"), to deal -* in the Software without restriction, including without limitation the rights -* to use, copy, modify, merge, publish, distribute, sublicense, and/or sell -* copies of the Software, and to permit persons to whom the Software is -* furnished to do so, subject to the following conditions: -* -* The above copyright notice and this permission notice shall be included in -* all copies or substantial portions of the Software. -* -* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR -* IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, -* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE -* AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER -* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, -* OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN -* THE SOFTWARE. -*/ - -/****************************************************************************************************/ -/* Definition of all macros used */ -/****************************************************************************************************/ -#pragma once -#include -#include -#include -#include -#include -#include - - -/****************************************************************************************************/ -/* Macro for vector initialization */ -/****************************************************************************************************/ -#define _INIT(x) {x, 0.0, 0.0, 0.0, 0.0} -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - - -/****************************************************************************************************/ -/* Macros for calculation of nth RK term */ -/****************************************************************************************************/ -#define _GET_VARNAMES(r, data, i, elem) BOOST_PP_COMMA_IF( i ) BOOST_PP_CAT(data, elem) -#define _SET_RK1(r, data, elem) BOOST_PP_CAT(data, elem) = (elem[0]); -#define _SET_RK2(r, data, elem) BOOST_PP_CAT(data, elem) = (elem[0] + elem[1] * 0.5); -#define _SET_RK3(r, data, elem) BOOST_PP_CAT(data, elem) = (elem[0] + elem[2] * 0.5); -#define _SET_RK4(r, data, elem) BOOST_PP_CAT(data, elem) = (elem[0] + elem[3]); - -#define _SWITCH(...) \ - double BOOST_PP_SEQ_FOR_EACH_I(_GET_VARNAMES, var_, __VA_ARGS__); \ - switch(N) { \ - default: \ - BOOST_PP_SEQ_FOR_EACH(_SET_RK1, var_, __VA_ARGS__) \ - break; \ - case 2: \ - BOOST_PP_SEQ_FOR_EACH(_SET_RK2, var_, __VA_ARGS__) \ - break; \ - case 3: \ - BOOST_PP_SEQ_FOR_EACH(_SET_RK3, var_, __VA_ARGS__) \ - break; \ - case 4: \ - BOOST_PP_SEQ_FOR_EACH(_SET_RK4, var_, __VA_ARGS__) \ - break; \ - } -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ - - -/****************************************************************************************************/ -/* Macros for repeated strings */ -/****************************************************************************************************/ -#define _REPEAT1(z, num, elem) BOOST_PP_COMMA_IF(num) elem -#define _REPEAT(elem, num) BOOST_PP_REPEAT(num , _REPEAT1 , elem) -/****************************************************************************************************/ -/* end */ -/****************************************************************************************************/ diff --git a/saves.h b/saves.h index 5d846ff..ff7c3a3 100644 --- a/saves.h +++ b/saves.h @@ -29,9 +29,10 @@ /****************************************************************************************************/ /* Save data */ /****************************************************************************************************/ -inline void get_data(int counter, Thalamic_Column& Col, double* Vt, double* Vr) { +inline void get_data(int counter, Thalamic_Column& Col, double* Vt, double* Vr, double* ah) { Vt [counter] = Col.Vt [0]; Vr [counter] = Col.Vr [0]; + ah [counter] = Col.act_h (); } /****************************************************************************************************/ /* end */