diff --git a/matlab/brew_eval.m b/matlab/brew_eval.m
index 321bae1..d2e35da 100755
--- a/matlab/brew_eval.m
+++ b/matlab/brew_eval.m
@@ -20,10 +20,11 @@ function brew_eval(varargin)
 %
 % brew_eval('M', 20.0, 'cont_ctrl', 0, 'kp', 0.1, 'ki', 0.1, 'kd', 100, 'T', 1800, 'mass_kleak', 0.2, 'P_q', 100, 'mass_delay', 100, [50 60 70 80]);
 
-WITH_LOWPASS = 1;
+WITH_LOWPASS = 0;
 WITH_RUNNING_MEAN = 0;
+WITH_RATE_CONTROLLER = 1;
 
-params = struct ('cont_ctrl',0, 'C',4.19e3, 'M',10, 'mass_kleak',0.1, 'mass_delay',50, 'P_min',0, 'P_max',3000, 'P_q',100, 'T',100, 'dt',1.0, 'Td',10.0, 'theta_lock',0.5, 'theta_amb',20.0, 'kp',10, 'ki',0.04, 'kd',500, 'pid_kleak',0.1, 'k_noise', 0.01, 'k_noise2', 0.0, 'f_noise2', 0.01);
+params = struct ('cont_ctrl',0, 'C',4.19e3, 'M',10, 'mass_kleak',0.5, 'mass_delay',50, 'P_min',0, 'P_max',3000, 'P_q',100, 'T',100, 'dt',1.0, 'Td',10.0, 'theta_lock',0.5, 'theta_amb',20.0, 'kp',10, 'ki',0.04, 'kd',500, 'pid_kleak',0.1, 'k_noise', 0.0, 'k_noise2', 0.0, 'f_noise2', 0.01);
 units = struct ('cont_ctrl', '', 'C', 'J/(kg*K)', 'M', 'kg', 'mass_kleak', '', 'mass_delay', 's', 'P_min', 'W', 'P_max', 'W', 'P_q','W', 'T', 's', 'dt', 's', 'Td', 's', 'theta_lock', '°C', 'theta_amb', '°C', 'kp', '', 'ki', '', 'kd', '', 'pid_kleak', '', 'k_noise', '', 'k_noise2', '', 'f_noise2', 'Hz');
 
 % Parse parameters
@@ -46,17 +47,22 @@ theta_ist = params.theta_amb;
 theta_last = params.theta_amb;
 thetas = varargin{nargin}
 heatRate = 0;
-heatRateMeasureInterval = 30.0; % seconds
+heatRate_f = 0;
+heatRateMeasureInterval = 10.0; % seconds
 heatRateMeasureCount = 0; % seconds
 n = 0;
 nSteps = params.P_max / params.P_q;
 
 N = 2*fix(length(thetas)*params.T/params.dt);
 
+maz = 17;
 lp_N = 17;
 [lp_b, lp_a] = butter (lp_N-1, 0.125);
 lp_z = zeros(1, length(lp_b)-1);
-maz = 17;
+
+lp_rate_N = 17;
+[lp_rate_b, lp_rate_a] = butter (lp_rate_N-1, 0.125);
+lp_rate_z = zeros(1, length(lp_rate_b)-1);
 
 for k=1:length(thetas)
 	theta_soll = thetas(k);
@@ -64,17 +70,23 @@ for k=1:length(thetas)
 	temp_ok = 0;
 	while n < N,
 		err = theta_ist - theta_soll;
-		[y, yc, s_pid] = pid(s_pid, params.kp, params.ki*params.dt, params.kd*params.dt, 1-params.pid_kleak, -err);
-
+		[y, yc, s_pid] = jpid(s_pid, params.kp, params.ki*params.dt, params.kd*params.dt, 1-params.pid_kleak, -err);
+		s_pid.y_max = 50;
+		
+		P_rate = params.P_max;
 		% Rate controller
 		err_rate = heatRate - 1.0;
-		[y_rate, yc_rate, s_pid_rate] = pid(s_pid_rate, 1.0*params.dt, 0.001*params.dt, 0*params.dt, 1-params.pid_kleak, -err_rate);
-		P_rate = max(params.P_min+1, min(params.P_max, y_rate));
+		[y_rate, yc_rate, s_pid_rate] = jpid(s_pid_rate, 0.4*params.dt, 0.001*params.dt, 8.0*params.dt, 0.999, -err_rate);
+		s_pid_rate.y_max = 4000;
+		if WITH_RATE_CONTROLLER
+			if (hold_theta == 0)
+				P_rate = min(params.P_max, params.P_max*y_rate);
+			end
+		end
 		% Rate controller
 
-		P_norm = max(params.P_min/params.P_max, min(1.0, yc));
-		P_cont = P_norm * params.P_max;
-		P_cont_q = fix(P_norm * nSteps)*params.P_q;
+		P_cont = max(params.P_min, min(P_rate, P_rate*y));
+		P_cont_q = P_cont; %fix(P_norm * nSteps)*params.P_q;
 		if params.cont_ctrl
 			P = P_cont_q;
 		else
@@ -92,22 +104,24 @@ for k=1:length(thetas)
 			[theta_ist, lp_z] = filter(lp_b, lp_a, theta_ist, lp_z);
 		end
 		if heatRateMeasureCount <= 0
-			heatRate = 60*(theta_ist-theta_last)/params.dt/heatRateMeasureInterval;
+			heatRate = 60*(theta_ist-theta_last)/heatRateMeasureInterval;
 			heatRateMeasureCount = heatRateMeasureInterval;
 			theta_last = theta_ist ;
 		else
 			heatRateMeasureCount = heatRateMeasureCount - params.dt;
 		end	
+		[heatRate_f, lp_rate_z] = filter(lp_rate_b, lp_rate_a, heatRate, lp_rate_z);
 		n = n + 1;
 		theta_ist = theta_ist;
-		heatRate_(n) = heatRate;
+		heatRate_(n) = heatRate_f;
 		theta_(n) = theta_ist;
 		err_(n) = err;
 		p_(n) = P;
 		pc_(n) = P_cont;
 		pcq_(n) = P_cont_q;
 		y_(n) = y;
-		P_rate_(n) = y_rate;
+		P_rate_(n) = P_rate;
+		err_rate_(n) = err_rate;
 		if hold_theta > 0
 			hold_theta = hold_theta - 1;
 			if hold_theta == 0
@@ -133,11 +147,11 @@ plot(t, err_); grid; xlabel('t/s'); ylabel('err/K'); title('Error'); legend('Err
 
 figure;
 subplot(3,1,1)
-plot(t, heatRate_); grid; xlabel('t/s'); ylabel('Theta/°'); title('Aufheizrate'); axis ([0 length(t)-1 0 5])
-subplot(3,1,2)
-plot(t, P_rate_); grid; xlabel('t/s'); ylabel('Theta/°'); title('Aufheizrate');
-subplot(3,1,3)
 plot(t, heatRate_); grid; xlabel('t/s'); ylabel('Theta/°'); title('Aufheizrate');
+subplot(3,1,2)
+plot(t, P_rate_); grid; xlabel('t/s'); ylabel('Theta/°'); title('P_rate');
+subplot(3,1,3)
+plot(t, err_rate_); grid; xlabel('t/s'); ylabel('Theta/°'); title('Error');
 
 function [theta, s] = mass(si, dt, C, M, L, P, Td, dTheta)
 s = si;
@@ -146,11 +160,11 @@ if ~isstruct(si)
 	if Td > 0
 		alpha = dt/Td;
 	end
-	s = struct('e', 0, 'a', alpha, 'x', 0);
+	s = struct('e', 0, 'a', alpha, 'x', 0, 'gain', 0.7);
 end
 
 s.e = s.e * (1-((L*dTheta)*dt)/(M*C));
-s.x = (1-s.a)*s.x + s.a*P*dt;
+s.x = (1-s.a)*s.x + s.gain*s.a*P*dt;
 s.e = s.e + s.x;
 theta = s.e/(M*C);
 endfunction