%Make figure: F0maxVBmax.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Steady-state optimization - 1: F0max fri, VBmax begrensende
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if 1==0

options=optimset;
options.DerivativeCheck='off';
options.Diagnostics='on';
options.Display='iter';
options.GradConstr='on';
options.GradObj='on';
options.LargeScale='off';
%options.TolCon=1e-9;
%options.TolFun=1e-9;
%options.TolX=1e-9;

%Disturbances
p1=[1 10:10:2000]; p1=p1(:);

%Initial values
load results/opt0
F0max=10000; p=p0; p(20)=F0max;
ux=ux0;
np1=size(p1,1);
YY=[];

%Optimization
tcpu=cputime;
for i=np1:-1:1
   p(1)=p1(i);
   [ux,fval,exitflag,output,lambda,grad,hessian]=fmincon('fun',ux,[],[],[],[],[],[],'nonlcon',options,p);
   [ux,fval,exitflag,output,lambda,grad,hessian]=fmincon('fun',ux,[],[],[],[],[],[],'nonlcon',options,p);
   VBmax=p(1); F0max=p(20); Mr=Y(7); xb=Y(12); V=Y(2); F0=Y(6);
   disp(sprintf('VBmax=%4.2f F0max=%4.2f Mr=%6.4f xb=%6.4f V=%6.4f F0=%6.4f exitflag=%d max(g)=%6.4f',VBmax,F0max,Mr,xb,V,F0,exitflag,max(g(3:8))));   
   YY=[YY;VBmax F0max Mr xb V F0 exitflag max(g(3:8))]; 
end
cputime-tcpu
plot(YY(:,1),YY(:,6))


%VBmax=YY(:,1);F0max=YY(:,6);EXITFLAG=YY(:,7);save results/F0maxVBmax F0max VBmax YY EXITFLAG
%VBmax=YY(:,1);F0max=YY(:,6);EXITFLAG=YY(:,7);save results/F0maxVBmax2 F0max VBmax YY EXITFLAG
%VBmax=YY(:,1);F0max=YY(:,6);EXITFLAG=YY(:,7);save results/F0maxVBmax3 F0max VBmax YY EXITFLAG

end %if 1==0

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Steady-state optimization - 1: F0max fri, VBmax begrensende
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if 1==0

options=optimset;
options.DerivativeCheck='off';
options.Diagnostics='on';
options.Display='iter';
options.GradConstr='on';
options.GradObj='on';
options.LargeScale='off';
%options.TolCon=1e-9;
%options.TolFun=1e-9;
%options.TolX=1e-9;

%Disturbances
%p1=[1 10:10:2000]; p1=p1(:);
p20=[1 10:10:500]; p20=p20(:);

%Initial values
load results/opt0
%F0max=10000; p=p0; p(20)=F0max;
Vmax=10000; p=p0; p(1)=Vmax;
ux=ux0;
np20=size(p20,1);
YY=[];

%Optimization
tcpu=cputime;
for i=np20:-1:1
   p(20)=p20(i);
   [ux,fval,exitflag,output,lambda,grad,hessian]=fmincon('fun',ux,[],[],[],[],[],[],'nonlcon',options,p);
   [ux,fval,exitflag,output,lambda,grad,hessian]=fmincon('fun',ux,[],[],[],[],[],[],'nonlcon',options,p);
   VBmax=p(1); F0max=p(20); Mr=Y(7); xb=Y(12); V=Y(2); F0=Y(6);
   disp(sprintf('VBmax=%4.2f F0max=%4.2f Mr=%6.4f xb=%6.4f V=%6.4f F0=%6.4f exitflag=%d max(g)=%6.4f',VBmax,F0max,Mr,xb,V,F0,exitflag,max(g(3:8))));   
   YY=[YY;VBmax F0max Mr xb V F0 exitflag max(g(3:8))]; 
end
cputime-tcpu
plot(YY(:,5),YY(:,6))

VBmax=YY(:,1);F0max=YY(:,6);EXITFLAG=YY(:,7);save results/F0maxVBmax3b F0max VBmax YY EXITFLAG

end %if 1==0

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Figure: F0maxVBmax.eps
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if 1==0

load results/F0maxVBmax3
VBmax=YY(:,1); F0max=YY(:,6);
load results/F0maxVBmax3b
VBmaxb=YY(:,1); F0maxb=YY(:,6);
plot(VBmax,F0max,'b-',VBmaxb,F0maxb,'r-');
hold on
xlabel('V_{max}');
ylabel('F_{0,max}');
h=gca;set(h,'FontSize',20);set(h,'LineWidth',1); 
h2=get(h,'XLabel');set(h2,'FontSize',20);
h3=get(h,'YLabel');set(h3,'FontSize',20);
%text(500,500,'Active: V_{max}');
text(500,500,'Active: V, x_B, M_r');
h4=get(h,'Children');set(h4(1),'FontSize',20);
%text(500,100,'Active: F_{0,max}');
text(500,100,'Active: F_0, x_{B}, M_r');
h4=get(h,'Children');set(h4(1),'FontSize',20);
set(h4(3),'LineWidth',2);


%print -deps /home/msg/LaTeX/larsson2002/figure/F0maxVBmax.eps %J=F0
%print -deps /home/msg/LaTeX/larsson2002/figure/F0maxVBmax2.eps%J=F0-0.1V
%print -deps /home/msg/LaTeX/larsson2002/figure/F0maxVBmax3.eps%J=F0-0.1D-0.01V

end %if 1==0

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Table: F0maxVBmax
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Vmax=YY(:,1); Mr=YY(:,3); xb=YY(:,4); F0=YY(:,6);
disp(sprintf('Vmax Fmax Mr xb V F0 exitflag max(g(3:8))'))
for i=1:size(YY,1)
   disp(sprintf('%9.4f %6.4f %6.4f %6.4f %9.4f %9.4f %6.4f %9.4f',YY(i,:)));
end  

Vmax=YY(:,1); Mr=YY(:,3); xb=YY(:,4); F0=YY(:,6);
disp(sprintf('$V_{max}$&$M_r$&$x_b$&$F_0$\\\\'))
for i=size(YY,1)-1:-1:1
   disp(sprintf('$%9.4f$&$%6.4f$&$%6.4f$&$%6.4f$\\\\',Vmax(i),Mr(i),xb(i),F0(i)))
end  







