diff --git a/Plots.m b/Plots.m
index af7c3b9..7226504 100644
--- a/Plots.m
+++ b/Plots.m
@@ -1,13 +1,13 @@
 % mex command is given by: 
-% mex CXXFLAGS="\$CXXFLAGS -std=gnu++0x -fpermissive -O3" Thalamus.cpp Thalamic_Column.cpp
+% mex CXXFLAGS="\$CXXFLAGS -std=c++11" Thalamus.cpp Thalamic_Column.cpp
 
 function Plots(T)
 
 if nargin == 0
-    Con     	= [ 4;      % sigma
-                    6;		% N_tr
-                    5;		% N_rt
-                    100];		% N_rr 
+    Con     	= [ 6;      % sigma
+                    4;		% N_tr
+                    4;		% N_rt
+                    20];		% N_rr 
                
 
     % stimulation parameters
@@ -18,19 +18,19 @@ function Plots(T)
     % 3 == phase dependend down state
     
     var_stim    = [ 0;           % mode of stimulation
-                    200;          % strength of the stimulus      in Hz (spikes per second)
-                    100;       	% duration of the stimulus      in ms
-                    6;          % time between stimuli          in s    
+                    50;          % strength of the stimulus      in Hz (spikes per second)
+                    120;       	% duration of the stimulus      in ms
+                    7;          % time between stimuli          in s    
                     1];         % time until stimuli after min 	in ms
     T       	= 30;  		% duration of the simulation
 end
-[Vt, Vr] = Thalamus(T, Con, var_stim);
+[Vt, Vr, I_t, I_h] = 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');
+Fs      = L/T;
+[Pxx,f] = pwelch(Vt-mean(Vt), hamming(4*Fs), 2*Fs, [], Fs,'ConfidenceLevel', 0.9, 'power');
 n       = find(f<=60, 1, 'last' );
 
 figure(1)
@@ -40,4 +40,12 @@ function Plots(T)
 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')
+%save('Thalamus.mat','Vt','Vr')
+% 
+% figure(2)
+% subplot(311), plot(timeaxis,Vt)
+% title('Thalamic relay membrane voltage'), xlabel('time in s'), ylabel('Vt in mV')
+% subplot(312), plot(timeaxis,I_t)
+% title('Thalamic reticular membrane voltage'), xlabel('time in s'), ylabel('Vr in mV')
+% subplot(313), plot(timeaxis,I_h)
+% title('Thalamic reticular membrane voltage'), xlabel('time in s'), ylabel('Vr in mV')
\ No newline at end of file
diff --git a/Thalamic_Column.cpp b/Thalamic_Column.cpp
index 1ca2007..f25e461 100644
--- a/Thalamic_Column.cpp
+++ b/Thalamic_Column.cpp
@@ -114,7 +114,7 @@ double Thalamic_Column::m_inf_T_t	(int N) const{
 /* Activation in RE population after Destexhe 1996 */
 double Thalamic_Column::m_inf_T_r	(int N) const{
 	_SWITCH((Vr))
-	double m = 1/(1+exp(-(var_Vr+52)/7.2));
+	double m = 1/(1+exp(-(var_Vr+52)/7.4));
 	return m;
 }
 
@@ -132,7 +132,7 @@ double Thalamic_Column::h_inf_T_r	(int N) const{
 	return h;
 }
 
-/* deactivation time in RE population after Destexhe 1996 */
+/* Deactivation time in RE population after Destexhe 1996 */
 double Thalamic_Column::tau_h_T_t	(int N) const{
 	_SWITCH((Vt))
 	double tau =  (30.8 + (211.4 + exp((var_Vt+115.2)/5))/(1 + exp((var_Vt+86)/3.2)))/pow(3, 1.2);
@@ -163,7 +163,8 @@ double Thalamic_Column::m_inf_h	(int N) const{
 /* Activation time for slow components in TC population after Destexhe 1993 */
 double Thalamic_Column::tau_m_h	(int N) const{
 	_SWITCH((Vt))
-	double tau = 1. / (exp(-14.59 - 0.086 * var_Vt) + exp(-1.87 + 0.07 * var_Vt));
+	//double tau = 1. / (exp(-14.59 - 0.086 * var_Vt) + exp(-1.87 + 0.07 * var_Vt));
+	double tau = (20 + 1000/(exp((var_Vt + 71.5)/14.2) + exp(-(var_Vt + 89)/11.6)));
 	return tau;
 }
 /****************************************************************************************************/
diff --git a/Thalamic_Column.h b/Thalamic_Column.h
index 02b2371..11d1a78 100644
--- a/Thalamic_Column.h
+++ b/Thalamic_Column.h
@@ -104,7 +104,7 @@ 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&, _REPEAT(double*, 4));
 
 private:
 	/* Population variables */
@@ -145,8 +145,8 @@ class Thalamic_Column {
 	const double 	theta_r		= -58.6;
 
 	/* Sigmoid gain in mV */
-	const double 	sigma_t		= 4;
-	const double 	sigma_r		= 4;
+	const double 	sigma_t		= 6;
+	const double 	sigma_r		= 6;
 
 	/* Scaling parameter for sigmoidal mapping (dimensionless) */
 	const double 	C1          = (3.14159265/sqrt(3));
@@ -166,10 +166,10 @@ class Thalamic_Column {
 
 	/* T current */
 	const double	g_T_t		= 3;
-	const double	g_T_r		= 2;
+	const double	g_T_r		= 2.3;
 
 	/* h current */
-	const double	g_h			= 0.08;
+	const double	g_h			= 0.05;
 
 	/* Reversal potentials in mV */
 	/* Synaptic */
@@ -187,16 +187,16 @@ class Thalamic_Column {
 	const double 	E_Ca    	= 120;
 
 	/* I_h current */
-	const double 	E_h    		= -43;
+	const double 	E_h    		= -40;
 
 	/* Calcium parameters */
-	const double	alpha_Ca	= -50E-6;			/* influx per spike in nmol		*/
+	const double	alpha_Ca	= -52E-6;			/* influx per spike in nmol		*/
 	const double	tau_Ca		= 10;				/* calcium time constant in ms	*/
-	const double	Ca_0		= 2E-4;				/* resting concentration 		*/
+	const double	Ca_0		= 2.4E-4;				/* resting concentration 		*/
 
 	/* I_h activation parameters */
 	const double 	k1			= 2.5E7;
-	const double 	k2			= 5E-4;
+	const double 	k2			= 4E-4;
 	const double 	k3			= 1E-1;
 	const double 	k4			= 1E-3;
 	const double 	n_P			= 4;
@@ -208,9 +208,9 @@ class Thalamic_Column {
 	double			input		= 0.0;
 
 	/* Connectivities (dimensionless) */
-	const double 	N_tr		= 6;
-	const double 	N_rt		= 5;
-	const double 	N_rr		= 100;
+	const double 	N_tr		= 4;
+	const double 	N_rt		= 4;
+	const double 	N_rr		= 20;
 
 	friend class 	Stim;
 };
diff --git a/Thalamus.cpp b/Thalamus.cpp
index 797e17b..f355a25 100644
--- a/Thalamus.cpp
+++ b/Thalamus.cpp
@@ -86,10 +86,14 @@ 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* I_T	= SetMexArray(1, T*res/red);
+	mxArray* I_h	= SetMexArray(1, T*res/red);
 
 	/* Pointer to the actual data block */
 	double* Pr_Vt	= mxGetPr(Vt);
 	double* Pr_Vr	= mxGetPr(Vr);
+	double* Pr_I_T	= mxGetPr(I_T);
+	double* Pr_I_h	= mxGetPr(I_h);
 
 	/* Simulation */
 	int count = 0;
@@ -97,7 +101,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_I_T, Pr_I_h);
 			++count;
 		}
 	}
@@ -105,6 +109,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 	/* Output of the simulation */
 	plhs[0] = Vt;
 	plhs[1] = Vr;
+	plhs[2] = I_T;
+	plhs[3] = I_h;
 return;
 }
 /****************************************************************************************************/
diff --git a/saves.h b/saves.h
index 5d846ff..d089033 100644
--- a/saves.h
+++ b/saves.h
@@ -29,9 +29,12 @@
 /****************************************************************************************************/
 /*											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* I_T,
+														double* I_h) {
 	Vt 		[counter] = Col.Vt 		[0];
 	Vr 		[counter] = Col.Vr 		[0];
+	I_T		[counter] = -Col.I_T_t	(0);
+	I_h		[counter] = -Col.I_h	(0);
 }
 /****************************************************************************************************/
 /*										 		end													*/