expeirmental campaing 1: data analysis complete

.gitignore

@@ -2,4 +2,4 @@
 **/*.pyc
 **/*.bak
 **/*.dat
-
+wiki

experimental_setup/data_analysis/air_density.m

+## Copyright (C) 2018 Juan Pablo Carbajal
+##
+## 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 this program. If not, see <http://www.gnu.org/licenses/>.
+
+## Author: Juan Pablo Carbajal <ajuanpi+dev@gmail.com>
+## Created: 2018-01-12
+
+## -*- texinfo -*-
+## @defun {@var{} =} air_density (@var{T}, @var{P}, @var{H})
+## 
+## T: Temperature in Celcius
+## P: Presure in Pascals
+## H: Relative humidity
+##
+## Reference: @url{https://en.wikipedia.org/wiki/Density_of_air}
+## @seealso{}
+## @end defun
+
+function [rho Psat] = air_density (Tc, P, H)
+
+  ## Constants
+  R  = 8.314;     #[J/K/mol] Universal gas constant
+  Rd = 287.058;   #[J/kg/K] Specific gas constant for dry air
+  Rv = 461.495;   #[J/kg/K] Specific gas constant for water vapor
+
+  # Conversion Pvsat --> Pv
+  kPv = 6.1078e2;
+  aPv = 7.5;
+  bPv = 237.3;
+
+  T = Tc + 273.15; # Kelvin
+
+  Psat = kPv * 10 .^ (aPv * Tc ./ (Tc + bPv)); # Pa
+  Pv   = H .* Psat;                            # Pa
+  Pd   = P - Pv;                               # Pa
+  
+  rho = Pd ./ (Rd * T) + Pv ./ (Rv * T);
+
+endfunction
+
+%!test
+%! T = 35:-5:-25; # C
+%! P = 101.325e3; # Pa
+%! H = 0;
+%! rho = [1.1455, 1.1644, 1.1839, 1.2041, 1.2250, 1.2466, 1.2690, 1.2922, ...
+%!        1.3163, 1.3413, 1.3673, 1.3943, 1.4224];
+%! assert (air_density(T,P,H), rho, 1e-4);
+

experimental_setup/sandbox/maxonRE50_model.m → experimental_setup/data_analysis/maxonRE50_model.m

@@ -25,14 +25,18 @@
 ## @var{V} should be in volts and @var{I} in amperes.
 ## @var{T} is given in Nm and @var{W} in rad/s.
 ##
-## If a gearbox is attached to the motor, the reduction ratio  and its efficinecy
-## should be given in the structure @var{gear}.
-## The reduced speed and torque are returned in @var{Tg} and @var{Wg}. 
+## If a gearbox is attached to the motor, the reduction ratio  and its efficiency
+## should be given in the structure @var{gear}, i.e. 
+## @code{@var{gear} = struct ('ratio', @var{r}, 'efficiency', @var{f}})}; where
+## @var{r} and @var{f} ar ethe ratio and efficiancy of the gearbox, resp.
+## The reduced speed and torque are returned in @var{Tg} and @var{Wg}.
+##
+## Run @code{demo maxonRE50_model} for examples of use.
 ##
 ## Reference: @url{https://www.maxonmotor.com/maxon/view/product/370354}
 ## @end defun
 
-function [T W Tg Wg] = maxonRE50_model (V, I, gear = [])
+function [T W varargout] = maxonRE50_model (V, I, gear = [])
 
   ## Datasheet values
   kM   = 242;   # Torque constant [mNm/A].
@@ -53,10 +57,12 @@ function [T W Tg Wg] = maxonRE50_model (V, I, gear = [])
   T *= fT;
   W *= fW;
 
-  Tg = T; Wg = W;
   if !isempty (gear)
+    Tg = T; Wg = W;
     Wg /= gear.ratio;
     Tg *= gear.ratio * gear.efficiency;
+    if nargout > 2; varargout{1} = Tg; endif
+    if nargout > 3; varargout{2} = Wg; endif
   endif
 
 endfunction

experimental_setup/data_analysis/pickplot2d.m

+## Copyright (C) 2008-2018 Juan Pablo Carbajal
+##
+## 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 this program. If not, see <http://www.gnu.org/licenses/>.
+
+## Author: Juan Pablo Carbajal <ajuanpi+dev@gmail.com>
+
+## -*- texinfo -*-
+## @defun {@var{idx} =} pickplot2d (@var{x}, @var{y})
+## @defunx {@var{idx} =} pickplot2d (@var{h}, @dots{})
+## 
+## @seealso{ginput}
+## @end defun
+
+function idx = pickplot2d (fig, x, y, fmt)
+
+  if !ishandle (fig)
+    y = x;
+    x = fig;
+    fig = figure ();
+  else
+    figure (fig);
+  endif
+  plot (x, y, fmt);
+  axis tight;
+
+  basetxt = ['Click the %s button to define\nthe %s of the selection box.\n' ...
+             'Click %s button to cancel'];
+  %corner = {'upper-left','botton-right'};
+  corner = {'1st corner','2nd corner'};
+  button = {'left','middle','right'};
+
+  % box
+  xrange = [min(x), max(x)];
+  yrange = [min(y), max(y)];
+
+  hold on
+  hp = plot (x, y,'sr');
+  set(hp, 'visible', 'off')
+  axis tight
+  hold off
+
+  htxt = text (0.05,0.9,'','Units','normalized');
+
+  idx = [];
+  done = 0;
+  while done ~=1
+    %% Upper-left
+    txt = sprintf (basetxt, button{1}, corner{1}, button{3});
+    set (htxt, 'String', txt);
+    btn = 0;
+
+    while btn ~= 1
+     [pickx picky btn] = ginput (1);
+     if btn == 3
+        disp ('Aborting...')
+        return
+     endif
+    endwhile
+    xrange(1) = pickx;
+    yrange(1) = picky;
+
+    %% Bottom-right
+    txt = sprintf (basetxt, corner{2}, button{1}, button{3});
+    set (htxt, 'String', txt);
+    btn = 0;
+
+    while btn ~= 1
+      [pickx picky btn] = ginput(1);
+      if btn == 3
+         disp ('Aborting...')
+         return
+      endif
+    endwhile
+    xrange(2) = pickx;
+    yrange(2) = picky;
+
+    xrange = sort (xrange, 'ascend');
+    yrange = sort (yrange, 'ascend');
+    
+    indx = x >= xrange(1) & x <= xrange(2);
+    indy = y >= yrange(1) & y <= yrange(2);
+    idx  = sort (find (indx & indy));
+
+    set (hp, 'xdata', x(idx), 'ydata', y(idx))
+    set (hp, 'visible', 'on')
+    txt = sprintf (['Happy with the selection?\nPress %s button to restart or'...
+                    ' %s button to finish.'], button{2}, button{1});
+    set (htxt, 'String', txt)
+
+    [pickx picky btn] = ginput (1);
+    if btn == 2
+       done = 0;
+       set (hp, 'visible', 'off')
+    else
+       done = 1;
+       set (htxt, 'visible', 'off')
+       set(hp, 'visible', 'off')
+    endif
+
+  endwhile
+
+endfunction

experimental_setup/data_logger/logger_single.py

@@ -173,7 +173,7 @@ def zmq_logger (outputfile, sensor, port = PORT, freq=100):
         fout.write(u'# Wind turbine HSR looger\n')
         fout.write(u'# Automatically created on {}\n'.format(datetime.now()))
         fout.write(u'# Desired sampling rate {} Hz, data is not uniformly sampled\n'.format(freq))
-        fout.write(u'# Time {}\n'.format(sensor.label))
+        fout.write(u'# Time [s], {}\n'.format(sensor.label))
 
         vals = [''] * 2
         live = ''
@@ -229,11 +229,11 @@ if __name__ == "__main__":
             options['freq'] = float(arg)
           elif opt in ("-s", "--sensor"):
             if arg == "T":
-              options['sensor'] = Torque (label='Torque')
+              options['sensor'] = Torque (label='Torque [mNm]')
             elif arg == "T2":
-              options['sensor'] = TorqueBi (label='Torque(Bi)')
+              options['sensor'] = TorqueBi (label='Torque(Bi) [mNm]')
             elif arg == "W":
-              options['sensor'] = RotSpeed (label='Ang. speed')
+              options['sensor'] = RotSpeed (label='Ang. speed [Hz]')
               if 'freq' not in options.keys():
                 options['freq'] = 50
 

experimental_setup/experiments/steadyrotation_fixedwind/README.md

@@ -109,7 +109,7 @@ fedw = [100  45 13.99 11.39 8.809 6.209 3.609 1.109].';
 
 ### Relaxation of rotational speed from high initial value and no wind, 1st modality
 
-Data stored in `friction_Urelax_*.dat`
+Data stored in `friction_Urelax_0_*.dat`
 
 Set maximum voltage and let the turbine go down to rest.
 The power supply controller is set to off but the motor is still
@@ -128,7 +128,7 @@ fe = [19 120 100 90].';
 
 ### Relaxation of rotational speed from high initial value and no wind, 2nd modality
 
-Data stored in `friction_Urelax_01_*.dat`
+Data stored in `friction_Urelax_1_*.dat`
 
 Set maximum voltage and let the turbine go down to rest.
 The power supply controller is set to off but the motor is still
@@ -145,7 +145,7 @@ I = [220 232; NA NA] * 1e-3;
 fe = [20; NA].';
 
 ### Relaxation of rotational speed from high initial value and no wind, 3rd modality
-Data stored in `friction_Urelax_02_*.dat`
+Data stored in `friction_Urelax_2_*.dat`
 
 Set maximum voltage and let the turbine go down to rest.
 The power supply controller is disconnected.

experimental_setup/experiments/steadyrotation_fixedwind/s_Cp_TSR.m

+## Copyright (C) 2018 Juan Pablo Carbajal
+##
+## 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 this program. If not, see <http://www.gnu.org/licenses/>.
+
+## Author: Juan Pablo Carbajal <ajuanpi+dev@gmail.com>
+## Created: 2018-01-10
+
+## Script to preprocess measured data
+#
+
+## Friction
+fname_F = 'tauF.dat';
+load (fname_F)
+
+fname_vp  = @(w) sprintf ("vp%s.dat", strrep (num2str (w), '.', '_'));
+fname_raw = @(w) sprintf ("rawdata_vp%s.dat", strrep (num2str (w), '.', '_'));
+
+vp = [12.5, 15];
+A = 2 * radiusR * 0.7;
+for i=1:2
+  load (fname_vp(vp(i)));
+  v_w{i} = v;
+  tsr{i} = vtR_data / v_w{i};
+  #tsr_std{i} = vtR_data_std / v;
+  tau{i} = - tauR - tauF(omegaR);
+  #tau_std{i} = tauR_std;
+
+  load (fname_raw(vp(i)),'Temp', 'Hum', 'Pres');
+  rho{i} = air_density (Temp, Pres, Hum);
+
+  Pw{i}  = 0.5 * A  * v^3 * rho{i};
+  C_p{i} = tau{i} .* omegaR ./ Pw{i};
+
+endfor
+
+
+## Plots
+#
+txt_label = @(v) {sprintf('%.1f m/s down',v), sprintf('%.1f m/s up',v)};
+figure (1)
+clf
+#h1  = errorbar (tsr{1}, tau{1}, tsr_std{1}, tau_std{1}, '~>o');
+h1  = plot (tsr{1}, tau{1}, 'o');
+hold on
+#h2  = errorbar (tsr{2}, tau{2}, tsr_std{2}, tau_std{2}, '~>^');
+h2  = plot (tsr{2}, tau{2}, '^');
+hold off
+axis tight
+xlabel ('TSR [a.u.]')
+ylabel ('\tau_W [Nm]')
+lbl = horzcat (txt_label(v_w{1}), txt_label(v_w{2}));
+legend ([h1; h2], lbl, ...
+'Location', 'NorthOutside','Orientation', 'horizontal');
+set (h1(1), 'color', get (h2(1), 'color'))
+set (h1(2), 'color', get (h2(2), 'color'))
+
+figure (2)
+clf
+h1  = plot (tsr{1}, C_p{1}, 'o');
+hold on
+h2  = plot (tsr{2}, C_p{2}, '^');
+hold off
+axis tight
+xlabel ('TSR [a.u.]')
+ylabel ('C_p [a.u.]')
+lbl = horzcat (txt_label(v_w{1}), txt_label(v_w{2}));
+legend ([h1; h2], lbl, ...
+'Location', 'NorthOutside','Orientation', 'horizontal');
+set (h1(1), 'color', get (h2(1), 'color'))
+set (h1(2), 'color', get (h2(2), 'color'))
+

experimental_setup/experiments/steadyrotation_fixedwind/s_friction.m

+## Copyright (C) 2018 Juan Pablo Carbajal
+##
+## 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 this program. If not, see <http://www.gnu.org/licenses/>.
+
+## Author: Juan Pablo Carbajal <ajuanpi+dev@gmail.com>
+## Created: 2018-01-11
+
+fname_data = "freq_friction.dat";
+load (fname_data);
+
+fname = 'tauF.dat';
+vars  = {'tauF', 'J', 'J_du', 'pwdot', 'pt'};
+
+##
+# synchronize
+
+np = numel (T);
+fmin = min(arrayfun(@(i)freqR{i}(1), 1:np));
+t = w = [];
+for i=1:np
+  j = find (freqR{i} <= fmin, 1, 'first');
+  t = [t; T{i}-T{i}(j)]; 
+  w = [w; freqR{i}];
+endfor
+w *= 2 * pi;
+[t,o] = sort (t);
+w = w(o);
+w = w;
+pw     = splinefit (t, w, 6, 'robust', true);
+pwdot  = ppder (pw);
+pt     = splinefit (w, t, 6, 'robust', true);
+tauF_J = @(w) ppval (pwdot, ppval(pt, w));
+
+##
+# properties of gearbox installed in motor
+gearGP52 = struct ('ratio', 4.3, 'efficiency', 0.91);
+radiusR  = 0.35; # in meters
+
+U(:,1) = [2 6 10 14 18 22 26 30].';
+U(:,2) = [30 26 22 18 14 10 6 2];
+
+I_lu(:,:,1) = [29 42; 45 65; 78 86; 104 113; 133 145; 155 187; 196 211; 234 244] * 1e-3;
+I_lu(:,:,2) = [221 232; 195 205; 158 173; 129 140; 104 110; 77 85; 51 64; 32 42] * 1e-3;
+I = squeeze (mean (I_lu, 2));
+
+fe(:,1) = [1.1 3.7 6.3   8.8   11.4  14.0  16.5  23].';
+fe(:,2) = [100  45 13.99 11.39 8.809 6.209 3.609 1.109].';
+
+# Error of the mean from interval
+# Measurements are consider 0.99 confidence intervals
+coeff = 1 ./ erfinv (0.99);
+lu2std = @(LU) squeeze (std (LU, 0, 2)) * coeff;
+
+# Motor/rotor torque and rotation
+[TM, WM, TR, WR]         = maxonRE50_model (U, I, gearGP52);
+[TMlu, WMlu, TRlu, WRlu] = maxonRE50_model (U, lu2std(I_lu), gearGP52);
+
+tauM = TM; tauM_std = TMlu;
+tauR = TR; tauR_std = TRlu;
+omegaM = WM; omegaM_std = WMlu;
+omegaR = WR; omegaR_std = WRlu;
+
+# moment of inertia
+
+J_du(1) = -tauF_J(omegaR(:,1)) \ tauR(:,1);
+J_du(2) = -tauF_J(omegaR(:,2)) \ tauR(:,2);
+J       = mean (J_du);
+
+tauF  = @(w) J * ppval (pwdot, ppval(pt, w));
+
+## Plots
+wf    = linspace (0, max(w)*1.1, 100);
+figure (1);
+plot (t, w, '.r', t, ppval(pw, t), 'k-', ppval(pt, wf), wf, 'g-');
+axis tight
+xlabel ('Time [s]')
+ylabel ('\omega [rad/s]')
+
+figure(2)
+plot (wf, tauF(wf), 'k-', omegaR, -tauR, 'o');
+axis tight
+ylabel ('\tau_F [Nm]')
+xlabel ('\omega [rad/s]')
+
+save (fname, vars{:});
+fprintf ("Saved to %s\n", fname);
+

experimental_setup/experiments/steadyrotation_fixedwind/s_preprocess_rawdata.m

@@ -16,64 +16,143 @@
 ## Author: Juan Pablo Carbajal <ajuanpi+dev@gmail.com>
 ## Created: 2018-01-10
 
-## Voltage ramp up and down
-# Last index is down (=1) and up (=2)
-U      = zeros (8,2);
-U(:,1) = [30 26 22 18 14 10 6 2].';
-U(:,2) = flipud (U(:,1));
+## Script to preprocess measured data
+#
+# add folder with data analysis functions
 
-## Wind power 12.5 %
-# Wind speed from anemometer @ 12.5 %
-v12p5 = 5.0;
+fname_raw = @(w) sprintf ("rawdata_vp%s.dat", strrep (num2str (w), '.', '_'));
+fname     = @(w) sprintf ("vp%s.dat", strrep (num2str (w), '.', '_'));
+vars      = {'gearGP52', 'radiusR', ...
+             'tauM', 'tauR', 'tauM_std', 'tauR_std', ...
+             'omegaM', 'omegaR', 'omegaM_std', 'omegaR_std', ...
+             'omegaR_data', 'omegaR_data_std', ...
+             'vtR', 'vtR_std', 'vtR_data', 'vtR_data_std', 'v'};
 
-## 
-# Current @ 12.5 %
-I12p5_lu     = zeros (8,2,2);
-I12p5_lu(:,:,1) = [7 11; 11 16; 32 55; 60 70; 61 68; 49 55; 28 38; 8 16] * 1e-3;
-I12p5_lu(:,:,2) = [8 14; 34 41; 48 59; 59 70; 59 73; 32 59; 16 22; 9 15] * 1e-3;
+txt_label = @(v) {sprintf('%.1f m/s down',v), sprintf('%.1f m/s up',v)};
 
-I12p5 = squeeze (mean (I12p5_lu, 2));
+#vp = 15;
+load (fname_raw(vp))
+
+if ~exist ('vtR', 'var')
+  # properties of gearbox installed in motor
+  gearGP52 = struct ('ratio', 4.3, 'efficiency', 0.91);
+  radiusR  = 0.35; # in meters
+
+  # Error of the mean from interval
+  # Measurements are consider 0.99 confidence intervals
+  coeff = 1 ./ erfinv (0.99);
+  lu2std = @(LU) squeeze (std (LU, 0, 2)) * coeff;
+
+  ##
+  # Motor/rotor torque and rotation
+  [TM, WM, TR, WR]         = maxonRE50_model (U, I, gearGP52);
+  [TMlu, WMlu, TRlu, WRlu] = maxonRE50_model (U, lu2std(I_lu), gearGP52);
+
+  tauM = TM; tauM_std = TMlu;
+  tauR = TR; tauR_std = TRlu;
+  omegaM = WM; omegaM_std = WMlu;
+  omegaR = WR; omegaR_std = WRlu;
+  
+  ##
+  # Segmentation of measured data, see s_segment_omega.m
+  fn = sprintf ("freq_vp%s.dat", strrep (num2str (vp), '.', '_'));
+  load (fn,'freqR_data', 'freqR_data_std');
+  omegaR_data      = 2 * pi * freqR_data;
+  omegaR_data_std  = 2 * pi * freqR_data_std;
+  vtR_data         = omegaR_data * radiusR;
+  vtR_data_std     = omegaR_data_std * radiusR;
+
+endif
+
+## Plots
+#
 
 ##
-# Encoder rotation frequency
-fe12p5      = NA (8,2);
-fe12p5(:,2) = [NA NA NA 8.97 11.5 14.2 17.0 19.7].';
+# Motor voltage and current at stead state
+# The voltage was scanned with increasing (up) and decreasing values (down)
+# for two values of the wind speed.
+figure (1)
+clf
+h = plot (U,I,'o-');
+axis tight
+xlabel ('Voltage [V]')
+ylabel ('Current [A]')
+legend (h, txt_label(v));
 
 ##
-# Angular speed from sensor @ 12.5 %
-# Use s_extract_idx_angspeed.m to extract the indexes manually
-idx12p5 = zeros (8,2,2);
+# Motor rotational frequency as measured from the encoder and
+# as caculated using the manufacturer's model.
+# The plots show that the error propagated through the model 
+# is an over estimation of the uncertainty in the angular speed.
+# Therefore we assign the erro of th emeasured speed to the model output.
+
+figure (2)
+clf
+fm  = omegaM / 2 / pi; 
+if ~exist ('vtR', 'var')
+  subplot (1,2,1)
+  errf = omegaM_std / 2 / pi;
+  h   = errorbar (fm, fe, errf, fe_err*ones(size(errf)),'~>o');
+  ff  = [fm(:) fe(:)]; ff(!isfinite(ff(:,2)),:) = [];
+  vax = [min(ff(:)) max(ff(:))];
+  axis (vax([1 2 1 2]));
+  line (vax, vax);
+  xlabel ('Freq. motor (model) [Hz]')
+  ylabel ('Freq. encoder (measured) [Hz]')
+  legend (h, txt_label(v), ...
+  'Location', 'NorthOutside','Orientation', 'horizontal');
 
-## Wind power 15 %
-v15 = 6.0;
+  omegaM_std  = fe_err * 2 * pi;
+  omegaR_std = omegaM_std / gearGP52.ratio;
 
-## 
-# Current @ 15 %
-I15_lu        = zeros (8,2,2);
-I15_lu(:,:,1) = [0 0;  1  8; 45 60; 54 66; 53 62; 34 44; 16 25; 0 0] * 1e-3;
-I15_lu(:,:,2) = [0 1; 19 26; 36 42; 55 64; 56 69; 37 49;  8 14; 0 0] * 1e-3;
 
-I15 = squeeze (mean (I15_lu, 2));
+  subplot (1,2,2)
+endif
 
-##
-# Encoder rotation frequency
-fe15      = NA (8,2);
-fe15(:,1) = [20.0 17.1 14.2 11.5 8.86 6.36 3.76 1.36].';
-fe15(:,2) = [1.3 3.8 6.4 9.08 11.6 14.3 17.1 19.9].';
+ferr = fe_err * ones(size(fm));
+h = errorbar (fm, fe, ferr, ferr,'~>o');
+axis tight
+vax = axis ();
+line (vax(1:2), vax(3:4));
+xlabel ('Freq. motor (model) [Hz]')
+ylabel ('Freq. encoder (measured) [Hz]')
+legend (h, txt_label(v), ...
+'Location', 'NorthOutside','Orientation', 'horizontal');
 
 ##
-# Angular speed from sensor @ 15 %
-# Use s_extract_idx_angspeed.m to extract the indexes manually
-idx15 = zeros (8,2,2);
+# Torque of the motor applied to the turbine rotor vs. tangential speed of
+# turbine rotor
+# Tangential speed
+vtR     = omegaR * radiusR;
+vtR_std = omegaR_std * radiusR;
 
+figure (3)
+clf
+subplot(2,1,1)
+errV = vtR_std * ones (size (vtR));
+h    = errorbar (vtR, tauR, errV, tauR_std, '~>o');
+axis tight
+xlabel ('Tangetial speed rotor (model) [rad/s]')
+ylabel ('\tau_M [Nm]')
+legend (h, txt_label(v), ...
+'Location', 'NorthOutside','Orientation', 'horizontal');
 
-## Plots
-figure (1)
+subplot(2,1,2)
+h    = errorbar (vtR_data, tauR, vtR_data_std, tauR_std, '~>o');
+axis tight
+xlabel ('Tangetial speed rotor (sensor)[rad/s]')
+ylabel ('\tau_M [Nm]')
+legend (h, txt_label(v), ...
+'Location', 'NorthOutside','Orientation', 'horizontal');
+
+figure (4)
 clf
-h = plot(U,I15,'^-',U,I12p5,'o-');
-set (h(3), 'color', get (h(1),'color'))
-set (h(4), 'color', get (h(2),'color'))
+h  = errorbar (vtR_data / v, tauR, vtR_data_std / v, tauR_std, '~>o');
 axis tight
-xlabel ('Voltage [V]')
-ylabel ('Current [A]')
-legend (h, {'15% down', '15% up', '12.5% down', '12.5% up'})
+xlabel ('TSR [a.u.]')
+ylabel ('\tau_M [Nm]')
+legend (h, txt_label(v), ...
+'Location', 'NorthOutside','Orientation', 'horizontal');
+
+save (fname(vp), vars{:});
+

experimental_setup/experiments/steadyrotation_fixedwind/s_rawdata.m

+## Copyright (C) 2018 Juan Pablo Carbajal
+##
+## 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 this program. If not, see <http://www.gnu.org/licenses/>.
+
+## Author: Juan Pablo Carbajal <ajuanpi+dev@gmail.com>
+## Created: 2018-01-10
+
+## Script to collect and save measured data
+#
+fname = @(w) sprintf ("rawdata_vp%s.dat", strrep (num2str (w), '.', '_'));
+vars  = {'vp', 'v', 'U', 'I', 'I_lu', 'fe', 'fe_err', 'Temp', 'Pres', 'Hum'};
+
+## Voltage ramp up and down
+# Last index is down (=1) and up (=2)
+U      = zeros (8,2);
+U(:,1) = [30 26 22 18 14 10 6 2].';
+U(:,2) = flipud (U(:,1));
+
+## Wind power 12.5 %
+# Wind speed from anemometer @ 12.5 %
+vp = 12.5;
+v = 5.0;  # m/s
+
+## 
+# Current @ 12.5 %
+I_lu        = zeros (8,2,2);
+I_lu(:,:,1) = [7 11; 11 16; 32 55; 60 70; 61 68; 49 55; 28 38; 8 16] ...
+                  * 1e-3;
+I_lu(:,:,2) = [8 14; 34 41; 48 59; 59 70; 59 73; 32 59; 16 22; 9 15] ...
+                  * 1e-3;
+I           = squeeze (mean (I_lu, 2));
+
+##
+# Encoder rotation frequency
+fe      = NA (8,2);
+fe(:,2) = [NA NA NA 8.97 11.5 14.2 17.0 19.7].';
+# Resolution in frequency measurement [Hz]
+fe_err = 0.1; 
+
+##
+# Temprature, Pressure and Relative humidity
+Pres = 963.0e2 * [1 1];  #Pa
+Hum  = [0.45 0.47];
+Temp = [14.9 13.5];     #C
+
+save (fname (vp), vars{:});
+
+## Wind power 15 %
+# Wind speed from anemometer @ 15 %
+vp  = 15;
+v = 6.0;  # m/s
+
+## 
+# Current @ 15 %
+I_lu        = zeros (8,2,2);
+I_lu(:,:,1) = [0 0;  1  8; 45 60; 54 66; 53 62; 34 44; 16 25; 0 0] ...
+                 * 1e-3;
+I_lu(:,:,2) = [0 1; 19 26; 36 42; 55 64; 56 69; 37 49;  8 14; 0 0] ...
+                 * 1e-3;
+I           = squeeze (mean (I_lu, 2));
+
+##
+# Encoder rotation frequency
+fe      = NA (8,2);
+fe(:,1) = [20.0 17.1 14.2 11.5 8.86 6.36 3.76 1.36].';
+fe(:,2) = [1.3 3.8 6.4 9.08 11.6 14.3 17.1 19.9].';
+# Resolution in frequency measurement [Hz]
+fe_err = 0.1;
+
+##
+# Temprature, Pressure and Relative humidity
+Pres = 963.0e2 * [1 1];  #Pa
+Hum  = [0.49 0.51];
+Temp = [12.2 11.4];     #C
+
+save (fname (vp), vars{:});

experimental_setup/experiments/steadyrotation_fixedwind/s_segment_freq.m

+## Copyright (C) 2018 Juan Pablo Carbajal
+##
+## 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 this program. If not, see <http://www.gnu.org/licenses/>.
+
+## Author: Juan Pablo Carbajal <ajuanpi+dev@gmail.com>
+## Created: 2018-01-11
+
+fname_raw = @(vp,d) sprintf ("steadyW_%s_wind%s_speed.dat", d, ...
+                          strrep (num2str(vp),'.','p'));
+fname = @(vp) sprintf ("freq_vp%s.dat", strrep (num2str (vp), '.', '_'));
+vars  = {'idx', 'freqR_data', 'freqR_data_std'};
+
+vp     = 12.5;             # wind generator power percentage
+direc  = {"down","up"};
+
+np  = 8;  # numer of measurements
+idx = zeros (np, 2, 2);
+freqR_data = freqR_data_std = zeros (np, 2);
+
+for j = 1:2
+  d = direc{j};  # direction
+  fn = fname_raw(vp,d);
+  printf ("Load from %s\n", fn); fflush (stdout);
+  data = load (fn);
+  t = data(:,1);
+  w = data(:,2);
+  idx_o = 1:length(t);
+
+  ##
+  # Filter NA values or too large
+  tf    = !isfinite (w) | w > 10;
+  t(tf) = [];
+  w(tf) = [];
+
+  f = figure (1,'Name', num2str(vp));
+  clf
+  for i=1:np
+    tmp = pickplot2d (f, t, w, '.');
+    ww                  = w(tmp);
+    freqR_data(i,j)     = mean (ww);
+    freqR_data_std(i,j) = std (ww) / sqrt (length (ww));
+    idx(i,:,j) = idx_o(tmp([1 end]));
+  endfor
+
+endfor
+close (f);
+fn = fname(vp);
+save (fn, vars{:});
+printf ("Saved to %s\n", fn); fflush (stdout);

experimental_setup/experiments/steadyrotation_fixedwind/s_segment_freq_friction.m

+## Copyright (C) 2018 Juan Pablo Carbajal
+##
+## 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 this program. If not, see <http://www.gnu.org/licenses/>.
+
+## Author: Juan Pablo Carbajal <ajuanpi+dev@gmail.com>
+## Created: 2018-01-11
+
+fname = "freq_friction.dat";
+vars  = {'idx', 'freqR', 'T'};
+
+if !exist(fname, 'file')
+  files = ls("friction_Urelax_*_speed.dat");
+  nf  = size(files, 1);  # numer of files
+  freqR = T = {};
+  idx = [];
+  for i = 1:nf
+    fn = files(i,:);
+    printf ("Load from %s\n", fn); fflush (stdout);
+    data = load (fn);
+    t = data(:,1);
+    w = data(:,2);
+    idx_o = 1:length(t);
+
+    ##
+    # Filter NA values or too large
+    tf    = !isfinite (w) | w > 5;
+    t(tf) = [];
+    w(tf) = [];
+
+    f = figure (1);
+    clf
+    axis tight
+    tmp = 0;
+    while true
+      tmp          = pickplot2d (f, t, w, '.');
+      if isempty (tmp); break; endif
+
+      T{end+1}     = t(tmp) - t(tmp(1));
+      freqR{end+1} = w(tmp);
+      idx(end+1,:) = idx_o(tmp([1 end]));
+    endwhile
+  endfor
+  close (f);
+else
+  load (fname);
+endif
+np = numel (T);
+fmin = min(arrayfun(@(i)freqR{i}(1), 1:np));
+## Plots
+# Detect outliers
+outlier = 6;
+
+figure(1)
+clf 
+hold on
+for i=1:np 
+  j = find(freqR{i}<=fmin,1,'first');
+  plot (T{i}-T{i}(j),freqR{i}, sprintf('o;%d;',i));
+endfor
+axis tight
+xlabel ('Time [s]')
+ylabel ('Freq [Hz]')
+
+to_keep = setdiff (1:np, outlier);
+T = T(to_keep);
+freqR = freqR(to_keep);
+idx = idx(to_keep,:);
+
+save (fname, vars{:});
+printf ("Saved to %s\n", fname); fflush (stdout);

experimental_setup/experiments/steadyrotation_fixedwind/startup.m

+## Copyright (C) 2018 Juan Pablo Carbajal
+##
+## 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 this program. If not, see <http://www.gnu.org/licenses/>.
+
+## Author: Juan Pablo Carbajal <ajuanpi+dev@gmail.com>
+## Created: 2018-01-11
+
+addpath ('../../data_analysis'); 
+