%% Initialisation
initial;
uuref=[-0.0103   -0.0084   -0.0048   -0.0172    0.0353];
uF3=0;
uunew=uuref;
uu=uuref;

sim('al');
Tref=Ttot(end);
Tnew=Tref;
Told=0;

j=0; %j is used to identify the first pass of the loop
delta=0.01; %delta is used to enhance the precision from passes to passes

% Gradient computation
stepu=0.00000000001;
for i=1:5
    uu(i)=uunew(i)+stepu;
    sim('al');
    Tdelta(i)=Ttot(end);
    Tgrad(i)=(Tdelta(i)-Tnew)/stepu;
    uu(i)=uunew(i);
end

%% Loop
while max(abs(Tgrad))>0.1 %when <= 0.1, the optimum is found with temperature precision of less than 0.001°C

%One pass is used to find the optimum in the gradient direction

linu(1,:)=uunew;
linT(1)=Tnew;


% If it is the first pass (j=0), the search scale is very large so we are sure to go
% behind the maximum temperature in the gradient direction
if j==0
linstep=0.05/max(abs(Tgrad));
uu=uunew+Tgrad*linstep;
else
linstep=0.005/max(abs(Tgrad));
uu=uunew+Tgrad*linstep;
end

%Three points are always considered by the algorithm, the maximum is always
%between point 1 and point 3
sim('al');
linu(3,:)=uu;
linT(3)=Ttot(end);
linu(2,:)=(linu(1,:)+linu(3,:))/2;
uu=linu(2,:);
sim('al');
linT(2)=Ttot(end);
linstep=linstep/2; %the point 2 is initialised at equal distance from 1 to 3

%safety condition, if the search scale was not big enough
if linT(3)<linT(1) 
    j=1;
else
    j=0;
end

%the distance is reduced between the points until a higher value is
%obtained on point 2
while linT(2)<linT(1)
    linT(3)=linT(2);
    linu(3,:)=linu(2,:);
    linstep=linstep/2; 
    uu=linu(3,:)-Tgrad*linstep;
    sim('al');
    linT(2)=Ttot(end);
    linu(2,:)=uu;
end




% then, the point 1 or 3 are replaced by new inside points
% the highest temperature point always becomes the point 2 

while linT(2)>max(linT(1),linT(3))+delta
[Tmin,index]=min(linT);
uu=(linu(index,:)+linu(2,:))/2;
sim('al');
if index==3
    if    Ttot(end)>linT(2)
          linT(1)=linT(2);
          linu(1,:)=linu(2,:);
          linT(2)=Ttot(end);
          linu(2,:)=uu;
    else  linT(3)=Ttot(end);
          linu(3,:)=uu;
    end
else
    if    Ttot(end)>linT(2)
          linT(3)=linT(2);
          linu(3,:)=linu(2,:);
          linT(2)=Ttot(end);
          linu(2,:)=uu;
    else  linT(1)=Ttot(end);
          linu(1,:)=uu;
    end    
end
end

%update
Told=Tnew;
uunew=linu(2,:);
Tnew=linT(2);

uu=uunew;
sim('al');
% the eventual new optimal F* split is then computed 
optimF2F3;

sim('al');
Tnew=Ttot(end);

% computation of the new local gradient
stepu=0.00000000001;
for i=1:5
    uu(i)=uunew(i)+stepu;
    sim('al');
    Tdelta(i)=Ttot(end);
    Tgrad(i)=(Tdelta(i)-Tnew)/stepu;
    uu(i)=uunew(i);
end

%Display the new point and gradient on workspace for supervision of the
%work being done by the algorithm
'New'
Tnew
uunew
Tgrad
delta=delta/2;
end

%% Display the result
'Finale'
uA=uAin+uunew(2)
pause
uB=uBin+uunew(3)
pause
uC=uCin+uunew(4)
pause
uD=uDin+uunew(5)
pause
uE=uEin+uunew(1)
pause
Toptim=Tnew
'End'
