function SingCtrlExample63
% Figs. A.5a, A.5b Book Illustration for Singular Control Example (3/2007): 
% Linear State Dynamics, Power Utility.
% 
clc % clear variables, but must come before globals, 
    %    else clears globals too.
clf % clear figures
fprintf('\nExample:  Singular Control OutPut:');
%
global mu0 p0 c0 delta0 x0 tf xsing umax TF 
%
nfig = 0;
scrsize = get(0,'ScreenSize');
ss = [3.0,2.8,2.6,2.4,2.2,2.0];
mu0 = 0.08; delta0 = 1.8*mu0; umax = 2*mu0; umin = 0;
p0 = 5.0; c0dp0 = 4.0; c0 = c0dp0*p0;
tf = 15.0; %t0 = 0; 
xsing = c0dp0/(1 - mu0/delta0);
using = mu0;
TF = tf + log(1-mu0*(delta0+umax-mu0)/(delta0*umax))/(delta0+umax-mu0);
% Assuming TF \geq T0max or TF \geq T0min;
x0v = [xsing+1,xsing-1]; nx0 = 2;
% cannot predeclare dynamic t,x,u vectors : ignore mlint warnings.
for ix0 = 1:nx0
    x0 = x0v(ix0);
    T0max = -log(xsing/x0)/(umax-mu0);
    T0min = +log(xsing/x0)/mu0;
    fprintf('\n x0=%f; c0/p0=%f; xsing=%f; XF(tf)=%f;',x0,c0dp0,xsing,XF(tf));
    fprintf('\n T0min=%f, T0max=%f; TF=%f; tf=%f;',T0min,T0max,TF,tf);
    dt = 0.01; % Comparison time step size;
    if x0 > xsing
        fprintf('\n Max0Case: x0 = %f > xsing = %f;',x0,xsing);
        n0 = ceil(T0max/dt); dt0 = T0max/n0;
        % t , x, u are dynamic vectors
        for i = 1:n0+1
            t(i) = (i-1)*dt0;
            x(i) = X0max(t(i));
            u(i) = umax;
        end
        ns = ceil((TF-T0max)/dt); dts = (TF-T0max)/ns; 
        for i = 2:ns+1
            it = i + n0;
            t(it) = T0max+(i-1)*dts;
            x(it) = xsing;
            u(it) = using;
        end
        nf = ceil((tf-TF)/dt); dtf = (tf-TF)/nf; 
        for i = 2:nf+1
            it = i + n0 + ns;
            t(it) = TF+(i-1)*dtf;
            x(it) = XF(t(it));
            u(it) = umax;
        end
        nt = n0 + ns + nf + 1;
    end
    if x0 < xsing
        fprintf('\n Min0Case: x0 = %f < xsing = %f;',x0,xsing);
        n0 = ceil(T0min/dt); dt0 = T0min/n0; 
        for i = 1:n0+1
            t(i) = (i-1)*dt0;
            x(i) = X0min(t(i));
            u(i) = umin;
        end
        ns = ceil((TF-T0min)/dt); dts = (TF-T0min)/ns; 
        for i = 2:ns+1
            it = i + n0;
            t(it) = T0min+(i-1)*dts;
            x(it) = xsing;
            u(it) = using;
        end
        nf = ceil((tf-TF)/dt); dtf = (tf-TF)/nf; 
        for i = 2:nf+1
            it = i + n0 + ns;
            t(it) = TF+(i-1)*dtf;
            x(it) = XF(t(it));
            u(it) = umax;
        end
        nt = n0 + ns + nf + 1;
    end
    t(nt) = tf;
%%%%%%%%%%%%%%%
    nfig = nfig + 1;
    fprintf('\n\nFigure(%i): Singular Control Example\n',nfig)
    figure(nfig)
    plot(t,x,'k-','LineWidth',4);
    if ix0 == 1, htitle=title('Singular Control Example (a)'); end
    if ix0 == 2, htitle=title('Singular Control Example (b)'); end
    if ix0 == 3, htitle=title('Singular Control Example (c)'); end
    if ix0 == 4, htitle=title('Singular Control Example (d)'); end
    set(htitle,'Fontsize',44,'FontWeight','Bold')
    ylabel('X^*(t), Optimal State'...
    ,'Fontsize',44,'FontWeight','Bold');
    xlabel('t, time'...
    ,'Fontsize',44,'FontWeight','Bold');
    hlegend=legend('X^*(t) Optimal State','Location','Best');
    set(hlegend,'Fontsize',36,'FontWeight','Bold');
    axis([0 tf 0 xsing+2]);
    set(gca,'Fontsize',36,'FontWeight','Bold','linewidth',3);
    set(gcf,'Color','White','Position' ...
        ,[scrsize(3)/ss(nfig) 70 scrsize(3)*0.60 scrsize(4)*0.80]);
end
%
% End Optimal Example Code
%
% function fv = F(x,u,t)
% global mu0 p0 c0 delta0 x0 tf xsing umax TF
% fv = (mu0-u)*x;
%
% function cv = C(x,u,t)
% global mu0 p0 c0 delta0 x0 tf xsing umax TF
% cv = exp(-delta0*t)*max(p0*x-c0,0)*u;
%
% function sv = S(x,t)
% global mu0 p0 c0 delta0 x0 tf xsing umax TF
% sv = 0;
%
function xv = XF(t)
global mu0 xsing umax TF
xv = xsing*exp(-(umax-mu0)*(t-TF));
%
function xv = X0max(t)
global mu0 x0 umax
xv = x0*exp(-(umax-mu0)*t);
%
function xv = X0min(t)
global mu0 x0
xv = x0*exp(mu0*t);
%
% End SingCtrlExample63