%  *******  LEVEL CONTROL AND EXAMPLES 18.2-4  ****
%  *******        from  MARLIN                 ****
%  COPYRIGHT McMASTER UNIVERSITY 1994
%  FILENAME : LVLSIM.M
%  LAST UPDATED: DECEMBER 1994
%  VERSION 2.0
%  DESCRIPTION
%  dynamic simulation of single level with feedback control
%  called from main menu
% All files in this program:
%    LVLINIT .M  -Initial data
%    LVLPLOT .M  -Simulation is executed
%
%  ********************************************************
%   INITIALIZATION
%  ********************************************************
% monitors
frate=0;
trate=0;

lvlinit ; % call program with initial data ;
% initialization
t=0.001     ;
omega = 3.14159*2/t ; % frequency of sine input ;
tdist = 0.05 * tend ;
l = Linit     ;  lm1 = l ;
fout = fouti; foutm1=fout; fin = fouti;
sp = l; spm1 = sp ;
cnt = 0 ;

%      initialize arrays used for plots
stepss = tend/deltat;
     Level = zeros(1,stepss);
     finp = Level; foutp = Level; time = Level ; 
     Levelsp = Level;
%  *******************************************************
   disp(' ******   SIMULATION BEING CALCULATED  ****')
%  *******************************************************
%  SIMULATION
%  *******************************************************

while t<tend

     % disturbance

     fin1 = amp*sin (omega*t)  ;
     fin2=0;
     if t>tdist
          fin2 = delf ;
     end 
     fin=fin1+fin2 +fouti;

     foutm1=fout;  % store for rate check
     %set gain
     kc=kcs;
     if abs(sp-l) >= brk
         kc=kcl;
     end
                    % p-only if ti=0
     if ti > 0
        fout = fout + kc*( (sp - l)-(spm1-lm1) + deltat*(sp - l)/ti);
     else
        fout = fout + kc*( (sp - l)-(spm1-lm1) );
     end
        
     lm1=l  ; spm1 = sp;

     l = l + deltat * (fin - fout)/A ;
     t=t+deltat;
     cnt=cnt+1;

     if abs(fout-foutm1)/deltat>frate   % find max rate of change
          frate=abs(fout-foutm1)/deltat;% time when max occurs
          trate=t;
     end

     Level(cnt)=l;
     finp(cnt) = fin;
     foutp(cnt)= fout;
     time(cnt) = t  ;
     Levelsp(cnt) = Linit;
end

     disp('  ***    SIMULATION COMPLETE   ***')
     pause (3) 
disp('Calculating...    Please wait several moments.') ;
disp('') ;
pause (3) ;
disp('When finished viewing graph, press ENTER.') ;
pause (5) ;
lvlplot % call M-file to execute the simulation ;
end