Added a project file for QT Creator.

Reworked the iterative solver for Stochastic RK after Rößler2010.
The new scheme ist much cleaner and easier. However, it seems the old one was slightly off, so new parameter tunings might be needed.
Removed unneeded use of variadic macros.

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 <iostream>
-#include <ctime>
+#include <chrono>
 #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<double> B1 = {0,
-								  0.626708569400000081728308032325,
-								  1.7296310295000001389098542858846,
-		 	 	 	 	 	 	  1.2703689705000000831347506391467};
-extern const vector<double> 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<std::chrono::milliseconds>( end - start ).count();
+    std::cout << "simulation done!\n";
+    std::cout << "took " << dif 	<< " seconds" << "\n";
+    std::cout << "end\n";
 }
 /****************************************************************************************************/
 /*										 		end			 										*/

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')

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