  % Simulate Oven Temperature Control
  %
  % This Octave/MATLAB script was written to mimic a corresponding
  % Labview module (Module 4.1) that accompanies the free text 
  % "Signals and Systems: Theory and Applications" by Ulaby and Yagle
  % (ss2.eecs.umich.edu).
  %
  % In almost all cases the compuation code used in these scripts
  % is the same as that used in the Labview module.  Credit for the
  % computation code should go to the textbook authors.
  
  % These scripts have been tested under Octave 4.2.2, Octave 4.4.0
  % and MATLAB R2015a.  The Octave GUI implementation is still a
  % little buggy.  Some of strange plotting code used here attempts
  % to work around a few Octave bugs.

  % Copyright (C) 2018 Anthony Richardson <richardson.tony@gmail.com>
  %
  % 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/>.
  %
  
  % R is thermal resistance of the walls
  % C is the heat capacity of the interior space
  % T is the desired temp change
  % K is the feedback loop constant
  RC = 50; T = 3; K = 0.05; 

  % Period of the sinusoidal input
  T0 = 200;
  
  % Slope of ramp input
  r = 0.0025;
  
  % Save initial simulation state
  simState('RC', RC); simState('T', T); simState('K', K);
  simState('T0', T0); simState('r', r);
  
  clf;
  createGUI();
  createPlot();
end

function createGUI
  % Set up plot windows
  h = subplot(2, 4, 1:3);
  guiState('hfigo', h);
  h = subplot(2, 4, 5:7);
  guiState('hfigi', h);
  
  % The slider handle and text handle for each slider are saved in GUI state
  pos = [.75 .0 .25 1];
  panel = uipanel ('position', pos, 'visible', 'off');
  guiState('cpanel', panel);
  
  %% BUG workaround - may be able to remove this in a
  %% future version of octave
  if (exist('OCTAVE_VERSION', 'builtin'))
    % The panel position seems to get messed up when the
    % window is resized in Octave. This was an attempt to fix it.
    %set(gcf, 'sizechangedfcn', @resetControlPanel);
    % That didn't work but this seems too.
    dellistener(gcf, 'position');
    addlistener (gcf, 'position', {@resetControlPanel, panel, pos});
  end
  
  h = sliderControl(panel, 'RC', 10, ...
                      100, simState('RC'), [2 5], ...
                      [.1 .9 1 0.05], @guiCallback );
  guiState('rch', h); 
  h = sliderControl(panel, 'T', 0.1, ...
                      5, simState('T'), [0.1 0.2] , ...
                      [.1 .8 1 0.05], @guiCallback );
  guiState('th', h);
  h = sliderControl(panel, 'K', 0.01, 0.2, ...
                      simState('K'), [0.01 0.05], ...
                      [.1 .7 1 0.05], @guiCallback );
  guiState('kh', h);
  h = sliderControl(panel, 'Period', 100, 300, ...
                      simState('T0'), [10 50], ...
                      [.1 .6 1 0.05], @guiCallback );
  guiState('t0h', h);
  h = sliderControl(panel, 'Ramp', 0, 0.005, ...
                      simState('r'), [0.00025 0.001], ...
                      [.1 .5 1 0.05], @guiCallback );
  guiState('rh', h);

  % Create a subpanel for the input type selection buttons
  cbpanel1 = uipanel (panel,'position', [0 .3 1 .1]);
  h = uicontrol(cbpanel1, 'style', 'checkbox', ...
               'units', 'normalized',...
                'string', 'step', ...
                'value', 1, ...
                'callback',@guiCallback,...
                'position', [0.1 0.5 0.2 0.2]);
  guiState('st_cb_h', h);
   
  h = uicontrol(cbpanel1, 'style', 'checkbox', ...
               'units', 'normalized',...
                'string', 'ramp', ...
                'value', 1, ...
                'callback',@guiCallback,...
                'position', [0.4 0.5 0.2 0.2]);
  guiState('ra_cb_h', h);

  h = uicontrol(cbpanel1, 'style', 'checkbox', ...
               'units', 'normalized',...
                'string', 'sine', ...
                'value', 1, ...
                'callback',@guiCallback,...
                'position', [0.7 0.5 0.2 0.2]);
  guiState('si_cb_h', h);
 
  % Create a subpanel for the feedback selection buttons
  cbpanel2 = uipanel (panel,'position', [0 .1 1 .1]);

  h = uicontrol(cbpanel2, 'style', 'checkbox', ...
               'units', 'normalized',...
                'string', 'W/O Feedback', ...
                'value', 1, ...
                'callback',@guiCallback,...
                'position', [0.1 0.5 0.4 0.2]);
  guiState('wof_cb_h', h);
  h = uicontrol(cbpanel2, 'style', 'checkbox', ...
               'units', 'normalized',...
                'string', 'With Feedback', ...
                'value', 1, ...
                'callback',@guiCallback,...
                'position', [0.6 0.5 0.4 0.2]);
  guiState('wif_cb_h', h);
end

function resetControlPanel(s, e, p, pos)
  %p = guiState('cpanel');
  set(p, 'visible', 'off');
  set(p, 'position', [0 0 1 1]);
  set(p, 'position', pos);
  set(p, 'visible', 'on');
end

function guiCallback(s, e)
  replot = false;
  g = guiState();
  if (s == g.rch)
    % RC slider
    v = round(get(s, 'value'));
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(v));
    simState('RC', v);
    replot = true;
  elseif (s == g.th)
    % temp change slider
    v = round(get(s, 'value')/0.1)*0.1;
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(v));
    simState('T', v);
    replot = true;
  elseif (s == g.kh)
    % K coefficient
    v = round(get(s, 'value')/0.01)*0.01;
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(v));
    simState('K', v);
    replot = true;
  elseif (s == g.t0h)
    % Period
    v = round(get(s, 'value')/10)*10;
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(v));
    simState('T0', v);
    replot = true;
  elseif (s == g.rh)
    % Ramp Slope
    v = round(get(s, 'value')/0.00025)*0.00025;
    set(s, 'value', v);
    set(get(s, 'userdata'), 'string', num2str(v));
    simState('r', v);
    replot = true;
  elseif (s == g.st_cb_h)
    % Step input check box
    if(get(s, 'value') == 0)            % If step input is 'off'
      set(g.stih, 'visible', 'off');    % hide step input
      set(g.stowih, 'visible', 'off');  % hide step output with fdback
      set(g.stowoh, 'visible', 'off');  % hide step output without fdback
    else                                % Else step input is 'on'
      set(g.stih, 'visible', 'on');     % Show step input
      if(get(g.wif_cb_h, 'value') == 1) % Is with fdback selected?
        set(g.stowih, 'visible', 'on'); % Show step fdback output
      end
      if(get(g.wof_cb_h, 'value') == 1) % Is without fdback selected?
        set(g.stowoh, 'visible', 'on'); % Show step wo fdback output
      end
    end
  elseif (s == g.ra_cb_h)
    % Ramp input check box
    if(get(s, 'value') == 0)            % Is ramp input 'off'
      set(g.raih, 'visible', 'off');
      set(g.raowih, 'visible', 'off');
      set(g.raowoh, 'visible', 'off');
    else
      set(g.raih, 'visible', 'on');
      if(get(g.wif_cb_h, 'value') == 1)
        set(g.raowih, 'visible', 'on');
      end
      if(get(g.wof_cb_h, 'value') == 1)
        set(g.raowoh, 'visible', 'on');
      end
    end
  elseif (s == g.si_cb_h)
    % Sine input check box
    if(get(s, 'value') == 0)            % Is sine input 'off'
      set(g.siih, 'visible', 'off');
      set(g.siowih, 'visible', 'off');
      set(g.siowoh, 'visible', 'off');
    else
      set(g.siih, 'visible', 'on');
      if(get(g.wif_cb_h, 'value') == 1)
        set(g.siowih, 'visible', 'on');
      end
      if(get(g.wof_cb_h, 'value') == 1)
        set(g.siowoh, 'visible', 'on');
      end
    end
  elseif (s == g.wof_cb_h)
    % Without feedbox check box
    if(get(s, 'value') == 0)            % Is wo fdback selection off?
      set(g.stowoh, 'visible', 'off');  % Hide all wo fdback outputs
      set(g.raowoh, 'visible', 'off');
      set(g.siowoh, 'visible', 'off');
    else                                % Show wo fdback output if input is
      if (get(g.st_cb_h, 'value') == 1) % also on.
        set(g.stowoh, 'visible', 'on');
      end
      if (get(g.ra_cb_h, 'value') == 1)
        set(g.raowoh, 'visible', 'on');
      end
      if (get(g.si_cb_h, 'value') == 1)
        set(g.siowoh, 'visible', 'on');
      end
    end
  elseif (s == g.wif_cb_h)              % Is with fdback selection off
    % With feedback check box
    if(get(s, 'value') == 0)
      set(g.stowih, 'visible', 'off');
      set(g.raowih, 'visible', 'off');
      set(g.siowih, 'visible', 'off');
    else                                % Show with fdback output plots
      if (get(g.st_cb_h, 'value') == 1)
        set(g.stowih, 'visible', 'on');
      end
      if (get(g.ra_cb_h, 'value') == 1)
        set(g.raowih, 'visible', 'on');
      end
      if (get(g.si_cb_h, 'value') == 1)
        set(g.siowih, 'visible', 'on');
      end
    end
  end
  
  % Don't replot on button selection/deselection
  % Plots are hidden/unhidden according to button selection
  % Replot if there is a change in a simulation value though.
  if (replot)
    createPlot();
  end
end

function createPlot()
  % The t vector is only computed once for efficiency.
  persistent t;
  t = linspace(0, 500, 100);
  persistent xstep;
  xstep = ones(size(t));
  
  % Get the current state
  s = simState();
  a = 1/s.RC;
  b = a + s.K;
  T = s.T;
  W0 = 2*pi/s.T0;
  g = guiState();
  
  xramp = s.r *t;
  xcos = cos(W0*t);

  % Setting subplot by use of saved handle does not seem to work in Octave
  % when window is resized
  %subplot(g.hfigi);
  subplot(2, 4, 5:7);
  h = plot(t, xstep, 'r', 'linewidth', 2);
  guiState('stih', h);
  val = get(g.st_cb_h, 'value');
  if(val == 0)
    set(h, 'visible', 'off');
  end
  hold on;
  h = plot(t, xramp, 'b', 'linewidth', 2);
  guiState('raih', h);
  val = get(g.ra_cb_h, 'value');
  if(val == 0)
    set(h, 'visible', 'off');
  end
  
  h = plot(t, xcos, 'k', 'linewidth', 2);
  guiState('siih', h);
  val = get(g.si_cb_h, 'value');
  if(val == 0)
    set(h, 'visible', 'off');
  end
  hold off;
  grid on;
  xlabel('Time (s)');
  ylabel('Q (normalized by T/R)');
  title('Power Input');
  % MATLAB ignores trailing space
  legend('step                   .', 'ramp', 'sine', ...
          'location', 'northeastoutside');
  legend('boxoff');
  
  % Setting subplot by use of saved handle does not seem to work in Octave
  %subplot(g.hfigo);
  subplot(2, 4, 1:3);
  Xstep = T*(1 - exp(-a.*t));
  Ystep = T*(1 - exp(-b.*t));
  
  Xramp = s.r*T*(t - 1/a + 1/a*exp(-a*t));
  Yramp = s.r*T*(t - 1/b + 1/b*exp(-b.*t));

  C = -a^2/(a^2 + W0^2);
  A = -C; B = a*(1 + C);
  Xcos = T*(C*exp(-a*t) + A*cos(W0*t) + (B/W0)*sin(W0*t));
  
  C = -b^2/(b^2 + W0^2);
  A = -C; B = b*(1 + C);
  Ycos = T*(C*exp(-b*t) + A*cos(W0*t) + (B/W0)*sin(W0*t));
  
  h = plot(t, Xstep, 'r', 'linewidth', 2);
  guiState('stowoh', h);
  if(or(get(g.st_cb_h, 'value') == 0, get(g.wof_cb_h, 'value') == 0))
    set(h, 'visible', 'off');
  end
  hold on;
  h = plot(t, Ystep, '--r', 'linewidth', 2);
  guiState('stowih', h);
  if(or(get(g.st_cb_h, 'value') == 0, get(g.wif_cb_h, 'value') == 0))
    set(h, 'visible', 'off');
  end
  
  h = plot(t, Xramp, 'b', 'linewidth', 2);
  guiState('raowoh', h);
  if(or(get(g.ra_cb_h, 'value') == 0, get(g.wof_cb_h, 'value') == 0))
    set(h, 'visible', 'off');
  end
  h = plot(t, Yramp, '--b', 'linewidth', 2);
  guiState('raowih', h);
  if(or(get(g.ra_cb_h, 'value') == 0, get(g.wif_cb_h, 'value') == 0))
    set(h, 'visible', 'off');
  end

  h = plot(t, Xcos, 'k', 'linewidth', 2);
  guiState('siowoh', h);
  if(or(get(g.si_cb_h, 'value') == 0, get(g.wof_cb_h, 'value') == 0))
    set(h, 'visible', 'off');
  end
  h = plot(t, Ycos, '--k', 'linewidth', 2);
  guiState('siowih', h);
  if(or(get(g.si_cb_h, 'value') == 0, get(g.wif_cb_h, 'value') == 0))
    set(h, 'visible', 'off');
  end

  hold off;
  grid on;
  xlabel('Time (s)');
  ylabel('Temperature Change');
  title('Oven Response');
  legend('step wo fdback', 'step with fdback', ...
         'ramp wo fdback', 'step with fdback', ...
         'sine wo fdback', 'sine with fdback', ...
         'location', 'northeastoutside');
   legend('boxoff');
   
   set(guiState('cpanel'), 'visible', 'on');
end

function v = guiState(key, val)
  persistent state;
  if (nargin == 0)  
    v = state;
  elseif (nargin == 1)
    v = state.(key);
  elseif (nargin == 2)
    state.(key) = val;
    v = val;
  end
end

function v = simState(key, val)
  persistent state;
  if (nargin == 0)  
    % Return the entire struct
    v = state;
  elseif (nargin == 1)
    % Return a particular field from the struct
    v = state.(key);
  elseif (nargin == 2)
    % Set a key value (key must be a string)
    state.(key) = val;
    v = val;
  end
end

function h = sliderControl(parent, label, min, max, val, step, pos, cback)
  % sliderstep = [minstep maxstep]
  % where minstep is fraction of range on arrow click
  % and maxstep is fraction of range on trough click
  h = uicontrol(parent, 'style' ,'slider',...
           'units', 'normalized',...
           'position', [pos(1)+0.1*pos(3) pos(2) 0.7*pos(3) pos(4)],...
           'min', min, 'max', max,...
           'sliderstep', [step(1)/(max-min) step(2)/(max-min)],...
           'value', val, 'callback', cback);
  uicontrol(parent, 'style', 'text',...
           'units', 'normalized',...
           'position', [pos(1) pos(2) 0.1*pos(3) pos(4)],...
           'String', num2str(min));
  %addproperty('my_slider_min', h, 'handle', hc);
  uicontrol(parent, 'style', 'text',...
           'units', 'normalized',...
           'position', [pos(1)+0.8*pos(3) pos(2) 0.1*pos(3) pos(4)],...
           'String', num2str(max));         
  %addproperty('my_slider_max', h, 'handle', hc);
  uicontrol(parent, 'style', 'text',...
           'units', 'normalized',...
           'position', [pos(1)+0.15*pos(3) pos(2)-pos(4) 0.4*pos(3) pos(4)],...
           'horizontalalignment','left',...
           'String', label);        
  hc = uicontrol(parent, 'style', 'text',...
           'units', 'normalized',...
           'position', [pos(1)+.5*pos(3) pos(2)-pos(4) 0.3*pos(3) pos(4)],...
           'String', num2str(val));
  %addproperty('my_slider_val', h, 'handle', hc);
  set(h, 'userdata', hc);
end