funcprot(0);
clear;
//exec('bode2.sci');
global axes1_h axes2_h axes3_h axes4_h  g;

x = 0; y = 0 ;
handles.axes1  = 1;
handles.axes2  = 1;
handles.axes3  = 1;
handles.axes4  = 1;
h=figure('position',[200 50 1000 550],'units','pixels');

//*************
//Initialize Axes
//*************
subplot(4,1,1);
handles.axes1=gca();
handles.axes1.axes_bounds = [0.15 -0.1 0.42,1.2];
plot(x(1),y(1));
handles.axes1.axes_visible = ['off','off','off']; 
    
subplot(4,1,2);
handles.axes2=gca();
handles.axes2.axes_bounds = [0.52 0 0.52,0.3];
plot(x(1),y(1));
handles.axes2.axes_visible = ['off','off','off']; 

subplot(4,1,3);
handles.axes3=gca(); 
handles.axes3.axes_bounds = [0.52 0.25 0.52,0.3];
plot(x(1),y(1));
handles.axes3.axes_visible = ['off','off','off']; 

subplot(4,1,4);
handles.axes4=gca(); 
handles.axes4.axes_bounds = [0.52 0.5 0.52,0.5];
plot(x(1),y(1));
handles.axes4.axes_visible = ['off','off','off']; 


//*************
// Plant Panel
//*************
uicontrol(h, "units","pixels", "style", "frame","string","Plant","position",[5 465 150 80],'fontsize',12);
  uicontrol(h, "units","pixels", "style", "text","string","Plant","position",[60 520 80 20],'fontsize',12);
  uicontrol(h, "units","pixels", "style", "text","string","Num","position",[10 500 100 20],'fontsize',12);
  uicontrol(h, "units","pixels", "style", "text","string","Den","position",[10 470 100 20],'fontsize',12);
  handles.num= uicontrol(h, "units","pixels", "style", "edit","string","3","position",[40 500 100 20],'fontsize',12);
  handles.den= uicontrol(h, "units","pixels", "style", "edit","string","3*s^2 + s + 2","position",[40 470 100 20],'fontsize',12);
  
//*************    
// Control Panel
//*************
handles.update= uicontrol(h, "units","pixels", "style", "pushbutton","string","Update","position",[40 430 100 30],"Callback","pbUpdate(handles,h)");
handles.gain= uicontrol(h, "units","pixels", "style", "slider","string","Update","position",[40 400 100 20],"Min",1,"Max",100,"SliderStep",[1 10],"Callback","sldGain(handles,h)");
handles.sldval = uicontrol(h, "units","pixels", "style", "text","string","1","position",[80 380 80 20],'fontsize',12);

//*************    
// Controller Panel
//*************    
uicontrol(h, "units","pixels", "style", "frame","string","Controller","position",[5 300 150 82],'fontsize',12);
  uicontrol(h, "units","pixels", "style", "text","string","Controller","position",[60 360 80 20],'fontsize',12);
  uicontrol(h, "units","pixels", "style", "text","string","Num","position",[10 340 100 20],'fontsize',12);
  uicontrol(h, "units","pixels", "style", "text","string","Den","position",[10 310 100 20],'fontsize',12);
  handles.numC= uicontrol(h, "units","pixels", "style", "edit","string","1","position",[40 340 100 20],'fontsize',12);
  handles.denC= uicontrol(h, "units","pixels", "style", "edit","string","1","position",[40 310 100 20],'fontsize',12);

handles.closef= uicontrol(h, "units","pixels", "style", "pushbutton","string","Close","position",[40 100 100 30],"Callback","pbClose(handles,h)");

// Initialization
axes1_h = 0;
axes2_h = 0;
axes3_h = 0;
axes4_h = 0;
g = 0;

//*************    
// Update Button
//*************    
function pbUpdate(handles,h)
  
// Preparing handles for plot updates 
    clear axes1_h; global axes1_h; axes1_h = 0;  
    clear axes2_h; global axes2_h; axes2_h = 0;   
    clear axes3_h; global axes3_h; axes3_h = 0;   
    clear axes4_h; global axes4_h; axes4_h = 0;
    val = get(handles.gain,'value');    
    
// Calculation    
 
  s=poly(0,'s');
  
  // Get the Plant Tranfer Functions
  num = eval(get(handles.num,'string'));
  den = eval(get(handles.den,'string'));
  Gol=syslin('c',(num),den);
  
  // Get the Controller Tranfer Functions
  numC = eval(get(handles.numC,'string'));
  denC = eval(get(handles.denC,'string'));
  Gc=syslin('c',(numC),denC);
  
  // Open Loop TF with Controller
  g = Gol * Gc;
  
// 1. Root Locus Plot
  sca(handles.axes1);
  delete(handles.axes1.children);
  evans(g);
  
  // Calculate close loop gain
  gl=1+val*g;
  numl=numer(gl);
  r=roots(numl);
  x=real(r);
  y=imag(r);
  plot(x,y,'o');
  axes1_h = gce();
  
// 2. Bode Mag Plot  
  sca(handles.axes2);
  //wrect=axes.axes_bounds;
  //axes.axes_bounds=[wrect(1)+0,wrect(2)+0,wrect(3)*1.0,wrect(4)*hx*0.95]
  [fdata,ddata,pdata,gdata]=bode2(val*g,0.01,100);
  handles.axes2.data_bounds = [min(fdata),min(ddata);max(fdata),max(ddata)];
  handles.axes2.log_flags = "lnn" ;
  handles.axes2.grid=color('lightgrey')*ones(1,3);
  handles.axes2.axes_visible="on";
  handles.axes2.clip_state = "clipgrf";
  mtlb_grid on;
  delete(handles.axes2.children);
  xpolys(fdata,ddata,gdata);
  axes2_h =  gce();
  
// 3. Bode Phase Plot  
  sca(handles.axes3);
  //wrect=axes.axes_bounds;
  //axes.axes_bounds=[wrect(1)+0,wrect(2)+0,wrect(3)*1.0,wrect(4)*hx*0.95]
  handles.axes3.data_bounds = [min(fdata),min(pdata);max(fdata),max(pdata)];
  handles.axes3.log_flags = "lnn" ;
  handles.axes3.grid=color('lightgrey')*ones(1,3);
  handles.axes3.axes_visible="on";
  handles.axes3.clip_state = "clipgrf";
  mtlb_grid on;
  delete(handles.axes3.children);
  xpolys(fdata,pdata,gdata);
  axes3_h =  gce();
  
// 4. Step Response
  sca(handles.axes4);
  mtlb_grid on;
  delete(handles.axes4.children);
  g2 = syslin('c',1,1);
  gcl = (val*g)/.g2;
  t=0:0.001:10;
  y=csim('step',t,gcl);  
  plot(t,y);
  axes4_h = gce();
  
// Resuming Sstates: ??? Is it correct way???
  [g,axes1_h,axes2_h,axes3_h,axes4_h] = resume(g,axes1_h,axes2_h,axes3_h,axes4_h);
  
endfunction


function sldGain(handles,h)
try
   if axes4_h ~= 0 
    val1 = get(handles.gain,'value');
    sleep(50);
    val2 = get(handles.gain,'value');
    sleep(50);
    while val1 ~= val2 do
      abort;//disp('Moving');
   end
    
    
    
  end
catch
  disp(lasterror);
end
UpdatePlot(handles,h)
  //sldGain(handles,h)
endfunction


function pbClose(handles,h)
  clear;
  close;
endfunction

function UpdatePlot(handles,h)
  
    val = get(handles.gain,'value');    
    set(handles.sldval,'string',string(val));
    
    gl=1+val*g;
    numl=numer(gl);
    r=roots(numl);
    x=real(r);
    y=imag(r);
    
    g2 = syslin('c',1,1);
    gcl = (val*g)/.g2;
    t=0:0.001:10;
    y2=csim('step',t,gcl);  
    [fdata,ddata,pdata,gdata]=bode2(val*g,0.01,100);
    
    // 1. Root Locus
    sca(handles.axes1);
    set(axes1_h.children,'data',[x,y]);
    
    // 2. Bode Mag Plot  
    sca(handles.axes2);
    mtlb_grid on;
    set(axes2_h.children,'data',[fdata,ddata]);
    handles.axes2.data_bounds = [mini(fdata),mini(ddata);maxi(fdata),maxi(ddata)];
    
    // 3. Phase Mag Plot  
    sca(handles.axes3);
    mtlb_grid on;
    set(axes3_h.children,'data',[fdata,pdata]);
    handles.axes3.data_bounds = [mini(fdata),mini(pdata);maxi(fdata),maxi(pdata)];
        
    // 4. Step Response
    sca(handles.axes4);  
    mtlb_axis([0 10 0 max(y2)]);
    set(axes4_h.children,'data',[t',y2']);

    
   [g,axes1_h,axes2_h,axes3_h,axes4_h] = resume(g,axes1_h,axes2_h,axes3_h,axes4_h);
   
endfunction


// Do note that bode2 is mmodified from the Scilab Function "bode"
// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
// Copyright (C) INRIA Serge Steer
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution.  The terms
// are also available at    
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt

function [fdata,ddata,pdata,gdata]=bode2(varargin)
  rhs=size(varargin)
  if type(varargin($))==10 then
    comments=varargin($),rhs=rhs-1;
  else
    comments=[];
  end
  fname="bode";//for error messages
  fmax=[]
  if or(typeof(varargin(1))==['state-space' 'rational']) then 
    //sys,fmin,fmax [,pas] or sys,frq
    refdim=1 //for error message
    if rhs==1 then //sys
      [frq,repf]=repfreq(varargin(1),1d-3,1d3)
    elseif rhs==2 then //sys,frq
      if size(varargin(2),2)<2 then
	error(msprintf(_("%s: Wrong size for input argument #%d: A row vector with length>%d expected.\n"),..
		       fname,2,1))
      end
      [frq,repf]=repfreq(varargin(1:rhs))
    elseif or(rhs==(3:4)) then //sys,fmin,fmax [,pas]
      [frq,repf]=repfreq(varargin(1:rhs))
    else
      error(msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"),fname,1,5))
    end
    [phi,d]=phasemag(repf)
    if rhs>=3 then fmax=varargin(3),end
  elseif  type(varargin(1))==1 then 
    //frq,db,phi [,comments] or frq, repf [,comments]
    refdim=2
    select rhs
    case 2 then //frq,repf
      frq=varargin(1);
      if size(frq,2)<2 then
	error(msprintf(_("%s: Wrong size for input argument #%d: A row vector with length>%d expected.\n"),..
		       fname,1,1))
      end
      if size(frq,2)<>size(varargin(2),2) then
	error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),..
			 fname,1,2))
      end
      [phi,d]=phasemag(varargin(2))
    case 3 then  //frq,db,phi
      [frq,d,phi]=varargin(1:rhs)
      if size(frq,2)<>size(d,2) then
	error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),..
			 fname,1,2))
      end
      if size(frq,2)<>size(phi,2) then
	error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same column dimensions expected.\n"),..
			 fname,1,3))
      end
    else 
       error(msprintf(_("%s: Wrong number of input arguments: %d to %d expected.\n"),fname,2,4))
    end
  else
    error(msprintf(_("%s: Wrong type for input argument #%d: Linear dynamical system or row vector of floats expected.\n"),fname,1))
  end;
   frq=frq';d=d',phi=phi'
  [n,mn]=size(d)

  if comments==[] then
    hx=0.48
  else
    if size(comments,'*')<>mn then
      error(msprintf(_("%s: Incompatible input arguments #%d and #%d: Same number of elements expected.\n"),...
		     fname,refdim,rhs+1))
    end
    hx=0.43
  end;

  if size(d,2)>1&size(frq,2)==1 then
    fdata = frq(:,ones(1,mn));
    ddata = d;
    gdata = 1:mn;
  else
    fdata = frq;
    ddata = d;
    gdata = 1:mn;
  end
 
  if size(phi,2)>1&size(frq,2)==1 then
    pdata = phi;
  else
    pdata = phi;
  end


endfunction


