function edlorenz(action,kind),
% This demo animates the integration of the three coupled nonlinear
% differential equations that define the "Lorenz Attractor", a chaotic
% system first described by Edward Lorenz of the MIT. 
% As the integration proceeds you will see:
% Picture on the left: a point moving around in a curious orbit in 3-D
%    space known as a strange attractor. The orbit is bounded, but not
%    periodic and not convergent. C represents horizontally averaged
%    temperature mode.                   
% Picture on the right: point is in time-parameter axes associated with the
%    stream function mode A. It can be switched to display the tent map.
% Picture in the center: the temperature difference between actual 
%    temperature and temperature in absence of convection. 
%
% References:
%    Lorenz, E.N., Deterministic non-periodic flow,
%	J. Atmos. Sci., 20, 130-141, 1963
%    Ott, Edward, Chaos in dynamical systems,
%	Cambridge University Press, 1993.  

% File edlorenz.m modified from lorenz.m Matlab-file 
%    by Svetlana Panasyuk, for educational purposes of MIT class 
%    on GeoSystems
%    last modified: August 2000

% Related questions can be addressed to:
% Svetlana Panasyuk
%    e-mail:  panasyuk@alum.mit.edu

if nargin<1,action = 'initialize';kind = 1;end;
%======================================================================
if strcmp(action,'initialize'),% set up the input graphical window
%======================================================================
  oldFigNumber = watchon;
  pos = [50 50 560 420;50 100 700 700];
  figNumber = figure('Position',pos(kind,:),'Name','Lorenz Attractor', ...
    'NumberTitle','off','Sharecolors','off','Renderer','painters', ...
    'Visible','off');
  MaaxHndl = axes('Units','normalized','Position',[.05 .10 .75 .95], ...
    'Visible','off');
    
  %------------------------------------------------------------------------  
  % Information for general descriptive and immediate action buttons: 
  %------------------------------------------------------------------------
  %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  % The CONSOLE frames, their titles and the desciptive texts:
  %
  if kind == 1,
    text(0,0,'Set initial conditions and Press the "Start" button ', ...
      'HorizontalAlignment','center');
    text(0,-0.1,'to begin the integration', ...
      'HorizontalAlignment','center');
    axis([-1 1 -1 1]);
    frmHndl = uicontrol('Style','frame','Units','normalized', ...
      'Position',[.78 .08 .21 .4],'BackgroundColor',[1 1 1]*.5);
    t = uicontrol('Style','text','Units','normalized', ...
      'Position',[.8 .425 .13 .04],'String','Input','BackgroundColor',[1 1 1]*.5);

    %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    % The INFO/Start/Repeat/Close buttons
    %
    strings = str2mat('Info','Start','Repeat','Close');
    btnpos = [.05 .02 .18 .05;.795 .02 .18 .05;.55 .02 .18 .05;.3 .02 .18 .05];
    clbacks = str2mat(['edlorenz(''info'',1)'],['edlorenz(''start'',1)'],...
      ['edlorenz(''repeat'',1)'],'close(gcf)');
    for n = 1:length(strings(:,1)),
      var =  uicontrol('Style','push','Units','normalized', ...
        'Position',btnpos(n,:),'String',strings(n,:),'CallBack',clbacks(n,:));
      eval([strings(n,1:4) 'Hndl = var;']);
    end;

    %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    % The InputInfo (II) bottons  and  the Input bottons
    %
    strings = str2mat('A','B','C','t','r');
    A = 5;B = -5;C = 20;t = 10;r = 28;
    for n = 1:5,
      InInHndl(n) = uicontrol('Style','push','Units','normalized', ...
        'Position',[.79 .45-.07*n .09 .05],'String',strings(n,:), ...
        'CallBack',['edlorenz(''InputInfo' num2str(n) ''')']);
      InpuHndl(n) = uicontrol('Style','Edit','Units','normalized', ...
        'Position',[.88 .45-.07*n .1 .05],'String',num2str(eval(strings(n,:))));
    end;

  else InfoHndl=-1;InInHndl=-ones(1,5);InpuHndl=-ones(1,5);
    StarHndl=-1;RepeHndl=-1;ClosHndl=-1;axis([-1 1 -1 1]);
  end;  
  
  %------------------------------------------------------------------------
  % The axes for the plots: BaaxHndl displays parameter space trajectory
  % TmaxHndl displays the temperature variations vs time,
  % AdaxHndl displays the advection temperature variations
  %------------------------------------------------------------------------
  ABC = uicontrol('Style','Text','Units','normalized',...
    'Position',[.4 .93 .4 .05],'HorizontalAlignment','left');

  strings = str2mat('Baax','Tmax','Adax','Teax');
  btnpos = [.42 .6 .3 .3;.07 .6 .3 .3;.02 .07 0.8 0.4;.78 .6 .21 .3];
  labels = str2mat('C','A','time','temperature variations (B)',' ',' ',...
    'max(i)','max(i+1)');
  for n = 1:length(strings(:,1)),
    var =  axes('Position',btnpos(n,:),'XTick',[],'YTick',[], ...
      'DrawMode','normal','NextPlot','add','Visible','off');
    xlabel(labels(2*n-1,:));ylabel(labels(2*n,:));
    eval([strings(n,1:4) 'Hndl = var;']);
  end;
  set(BaaxHndl,'XTick',0,'YTick',0,'ZTick',0,'View',[-37.5,30]);
  set(AdaxHndl,'Box','off');
  
  %----------------------------------------------------------------------
  % Uncover the figure: 1-14 BtnHndl, 15-19 AxesHndl , 20 ABC output   
  %----------------------------------------------------------------------
  hndlList = [InInHndl InpuHndl InfoHndl StarHndl RepeHndl ClosHndl ...
    MaaxHndl BaaxHndl TmaxHndl AdaxHndl TeaxHndl ABC];
  set(figNumber,'Visible','on','UserData',hndlList);
  watchoff(oldFigNumber);
  figure(figNumber);

  if kind == 2,edlorenz('start',2);end;

%===========================================================================
elseif strcmp(action,'start'), % To start the demo
%===========================================================================
  global SIGMA BETA RHO;if kind == 2,global tfinal y0 ANIM;end;
  %----------------------------------------------------------------------
  % gather all the input information and prepare for the run    
  %----------------------------------------------------------------------
  strings = str2mat('InInHnd','InpuHndl','InfoHndl','StarHndl','RepeHndl', ...
    'ClosHndl','MaaxHndl','BaaxHndl','TmaxHndl','AdaxHndl','TeaxHndl','ABC');
  hndlList = get(gcf,'UserData');
  for n = 3:length(strings(:,1)),
    eval([strings(n,:) ' = hndlList(8+n);']);
  end;
  InInHnd = hndlList(1:5);InpuHndl = hndlList(6:10);

  set(ABC,'String',' ');
  axes(MaaxHndl);cla;
  text(-.5,-.25,'Advection Temperature Field');
  set([BaaxHndl TmaxHndl TeaxHndl],'Visible','off','Visible','on');

  %---------------------------------------------------------------------
  % Start of Demo
  % The orbit ranges chaotically back and forth around two different points,
  % or attractors.  It is bounded, but not periodic and not convergent.
  %----------------------------------------------------------------------

  %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  % The initial conditions & run parameters
  %
  %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  SIGMA = 10;BETA = 8/3;
  if kind == 1,% nice results are for y0 = [20 5 -5]
    y0(1,1) = str2num(get(InpuHndl(3),'String'));
    y0(2,1) = str2num(get(InpuHndl(1),'String'));
    y0(3,1) = str2num(get(InpuHndl(2),'String'));
    RHO = str2num(get(InpuHndl(5),'String'));
    tfinal = str2num(get(InpuHndl(4),'String'));
  elseif kind == 2,% given: RHO tfinal y0 ANIM
  end;
  pow = 1/3;
  t0 = 0;t = t0;
  hmax = (tfinal - t)/5;hmin = (tfinal - t)/200000;h = (tfinal - t)/100;
  y = y0;
  tol = 0.001;tau = tol*max(norm(y,'inf'),1);
  CABT = [y0;t0];

  %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  % The main loop, do while until hmin
  % 
  while h >= hmin,
    if t + h > tfinal, h = tfinal - t; end
    % Compute the slopes
    s1 = [-BETA 0  y(2);0 -SIGMA SIGMA;-y(2)  RHO -1]*y;ys = y+h*s1;
    s2 = [-BETA 0 ys(2);0 -SIGMA SIGMA;-ys(2) RHO -1]*ys;ys= y+h*(s1+s2)/4;
    s3 = [-BETA 0 ys(2);0 -SIGMA SIGMA;-ys(2) RHO -1]*ys;

    % Estimate the error and the acceptable error
    delta = norm(h*(s1 - 2*s3 + s2)/3,'inf');
    tau = tol*max(norm(y,'inf'),1.0);
 
    ts = t;ys = y;
    if delta <=  tau,% Update the solution only if the error is acceptable
      t = t + h;y = y + h*(s1 + 4*s3 + s2)/6;
      CABT = [CABT [y;t]];
    end
    if delta ~=  0.0,% Update the step size
      h = min(hmax, 0.9*h*(tau/delta)^pow);
    end
  end;

  %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  % End of Demo
  %
  set(MaaxHndl,'UserData',CABT);
  ys = round(CABT(1:3,length(CABT))*10)/10;
  set(ABC,'String',['Final:  A=' num2str(ys(2)) '  B=' num2str(ys(3)) ...
    '  C=' num2str(ys(1)) ]);

  edlorenz('repeat',kind);
%===========================================================================
elseif strcmp(action,'repeat');
%===========================================================================
  global SIGMA BETA RHO;if kind == 2,global tfinal y0 ANIM;end;
  strings = str2mat('InInHnd','InpuHndl','InfoHndl','StarHndl','RepeHndl', ...
    'ClosHndl','MaaxHndl','BaaxHndl','TmaxHndl','AdaxHndl','TeaxHndl','ABC');
  hndlList = get(gcf,'UserData');
  for n = 3:length(strings(:,1)),
    eval([strings(n,:) ' = hndlList(8+n);']);
  end;
  InInHnd = hndlList(1:5);InpuHndl = hndlList(6:10);

  CABT = get(MaaxHndl,'UserData');len = length(CABT);
  L = (len<11)*round(len/2) + (len>=11)*10;% Save L steps and plot like a comet tail.
  mima = [min(CABT.');max(CABT.')];nums = find(diff(mima)==0);
  for ii = 1:length(nums),mima(:,nums(ii)) = [-1;1]*1e-14 + mima(1,nums(ii));end;

  %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  % Define the comet on both upper axes and calculate&plot the tent map
  %
  axes(BaaxHndl);cla;
  set(BaaxHndl,'Drawmode','fast','Visible','off','Visible','on', ...
    'XLim',mima([1 2],1),'YLim',mima([1 2],2),'ZLim',mima([1 2],3), ...
    'XGrid','on','YGrid','on','ZGrid','on','nextplot','add');
  xlabel('C');ylabel('A');zlabel('B');
  head = line('color','r','Marker','.','MarkerSize',25,'EraseMode','xor',...
    'xdata',CABT(1,1),'ydata',CABT(2,1),'zdata',CABT(3,1));
  body = line('color','r','linestyle','-','EraseMode','none',...
    'LineWidth',1,'xdata',[],'ydata',[],'zdata',[]);
  tail = line('color','b','linestyle','-','EraseMode','none',...
    'LineWidth',1,'xdata',[],'ydata',[],'zdata',[]);

axes(TmaxHndl);cla;
set(TmaxHndl,'Drawmode','fast','YLim',mima([1 2],3),'XLim',[0 mima(2,4)],...
   'YTick',[0],'YGrid','on','nextplot','add');
  xlabel('time');ylabel('temperature variations (B)');
  thead = line('color','r','Marker','.','MarkerSize',20,'EraseMode',...
    'xor','xdata',CABT(4,1),'ydata',CABT(3,1));
  tbody = line('color','r','linestyle','-','EraseMode','none',... 
    'xdata',[],'ydata',[]);
  ttail = line('color','b','linestyle','-','EraseMode','none',...
    'xdata',[],'ydata',[]);

axes(TeaxHndl);cla;
set(TeaxHndl,'Drawmode','fast','YTick',[]);
  xlabel('max(i)');ylabel('max(i+1)');title('1D tent map');
  y = abs(CABT(3,len:-1:1).');y_ = y.*(y>mean(y));
  yy = y_(1:length(y)-1)-y_(2:length(y));
  num1 = (yy>0);num2 = (yy<0);
  ln = length(num1)-1;NN = find(num2(1:ln).*num1(2:ln+1));ln = length(NN);
  if ln>2,
    plot(y(NN(2:ln)),y(NN(1:ln-1)),'ro');
    mimay=[min(y(NN(2:ln))) min(y(NN(1:ln-1)));max(y(NN(2:ln))) max(y(NN(1:ln-1)))];
    set(TeaxHndl,'AspectRatio',[NaN 1],'XLim',[mimay(1,1) mimay(2,1)], ...
      'YLim',[mimay(1,2) mimay(2,2)]);
    if RHO > 40,
      text(mimay(1,1)+1,round(mimay(2,1)/1.2),'tent-map algorithm might');
      text(mimay(1,1)+1,round(mimay(2,1)/3),'be unstable for this Ra#');
    end;
    axis([min(y(NN)) max(y(NN)) min(y(NN)) max(y(NN))]);
  else axis([0 CABT(4,len) min(CABT(3,:)) max(CABT(3,:))]);
    text(1,0,'not enough data');
  end;
  
  %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  % The parameters to display advection field and define colormap
  %
  l = 2/sqrt(4/BETA-1);coef1 = (4+l*l)/sqrt(2)/l;coef2 = (pi*coef1).^3/sqrt(2)/l;
  xx = [0:7]/7/8;yy = ([0:15]/15).';
  lenx = length(xx);leny = length(yy);lenxy = lenx*leny;
  s = coef2*sin(2*pi*yy)*ones(size(xx));
  cs = sqrt(2)*coef2*sin(pi*yy)*cos(2*pi*l*xx);

  axes(AdaxHndl);cla;
  junk = zeros(leny,lenx);Theta = junk;Ksi = junk;junk(:) = [1:lenxy].';
  pc=pcolor([junk junk(:,8:-1:1) junk junk(:,8:-1:1)]);
  set(pc,'EraseMode','none','EdgeColor','none','cdata',[]);
  Theta = CABT(1,1)*s - CABT(3,1)*cs;
  mimaT = [min(min(Theta)) max(max(Theta))];
  if diff(mimaT) == 0,Theta = ones(size(Theta));
  else Theta = Theta-mimaT(1);Theta = floor(Theta/max(max(Theta))*(lenxy-1))+1;
  end;

  colormap('default');colormap(hot);
  MyColorMap = (interp1([1:64]/64,colormap,[3:130]/130));
  set(gcf,'MinColormap',length(MyColorMap),'colormap',MyColorMap);

  %^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  % The drawing loop
  % 
  junk = [CABT(:,1)*ones(1,L) CABT];
  if kind == 1,ANIM = 1;elseif kind == 2, ANIM = (len-1)*(ANIM~=1)+1;end
  
  for ii = L + [ANIM:len],
    if abs(junk(1,ii))>0.23 & abs(junk(3,ii))>0.23,% to ensure stable colormap
      Theta = junk(1,ii)*s - junk(3,ii)*cs;
      Theta = Theta-min(min(Theta));
      Theta = floor(Theta/max(max(Theta))*(lenxy-1))+1;
      set(pc,'cdata',[Theta Theta(:,8:-1:1) Theta Theta(:,8:-1:1)]);
    end;
    vec1 = ii-[1 0];vec2 = ii-[L L-1];
    set(thead,'xdata',junk(4,ii),'ydata',junk(3,ii));
    set(tbody,'xdata',junk(4,vec1),'ydata',junk(3,vec1))
    set(ttail,'xdata',junk(4,vec2),'ydata',junk(3,vec2))
    set(head,'xdata',junk(1,ii),'ydata',junk(2,ii),'zdata',junk(3,ii))
    set(body,'xdata',junk(1,vec1),'ydata',junk(2,vec1),'zdata',junk(3,vec1))
    set(tail,'xdata',junk(1,vec2),'ydata',junk(2,vec2),'zdata',junk(3,vec2))
    drawnow;
 end;
 
  vec1 = len + [-L:-1];vec2 = [1:len-L];% plot all over again, for print option
  set(tail,'xdata',CABT(1,vec2),'ydata',CABT(2,vec2),'zdata',CABT(3,vec2));
  set(body,'xdata',CABT(1,vec1),'ydata',CABT(2,vec1),'zdata',CABT(3,vec1))
  set(head,'xdata',CABT(1,len),'ydata',CABT(2,len),'zdata',CABT(3,len));
  set(ttail,'xdata',CABT(4,vec2),'ydata',CABT(3,vec2));
  set(tbody,'xdata',CABT(4,vec1),'ydata',CABT(3,vec1));
  set(thead,'xdata',CABT(4,len),'ydata',CABT(3,len));

  axes(MaaxHndl);
  texts = str2mat('t_0=','t=','A=','B=','C=','A=','B=','C=','Ra=');
  strings = num2str([CABT(4,[1 len]).';CABT([2 3 1],1);CABT([2 3 1],len);RHO],3);
  pos = [.95 -.4;.95 -.7;1.2 -.4;1.2 -.5;1.2 -.6;1.2 -.7;1.2 -.8;1.2 -.9;1.2 -1];
  for n=1:length(pos(:,1)),text(pos(n,1),pos(n,2),[texts(n,:) strings(n,:)]);end;

%===========================================================================
elseif strcmp(action,'info');
%===========================================================================
  ttlStr = 'Lorenz Attractor Info';
  hlpStr =  ...                                             
        [' This demo animates the integration of the three  '  
         ' coupled nonlinear differential equations that    '  
         ' define the "Lorenz Attractor", a chaotic system  '  
         ' first described by Edward Lorenz of the MIT.     '  
         '    As the integration proceeds you will see:     '
         ' Picture on the left: a point moving around in a  '  
         ' curious orbit in 3-D space known as a strange    '  
         ' attractor.The orbit is bounded, but not periodic '  
         ' and not convergent. C represents horizontally    '  
         ' averaged temperature mode.                       '
         ' Picture on the right: point is in time-parameter '
         ' associated with the stream function mode A.      '
         ' Picture in the center: the temperature difference' 
         ' between actual temperature and temperature in    '
         ' absence of convection.                           '
         '                                                  '
         'File edlorenz.m modified from lorenz.m Matlab-file'
         '            by SvetLana Panasyuk while at MIT, for'
         'educational purposes of MIT class on GeoSystems   '
         '                          panasyuk@alum.mit.edu   '];
  helpfun(ttlStr,hlpStr);                                 

%===========================================================================
elseif strcmp(action,'InputInfo1');
%===========================================================================
  RttlStr = 'A-parameter Info';
  RhlpStr =  ...           
        [' A = 5 is a good example                         '
         '                                                 '];
  helpfun(RttlStr,RhlpStr);                                 
%===========================================================================
elseif strcmp(action,'InputInfo2');
%===========================================================================
  RttlStr = 'B-parameter Info';
  RhlpStr =  ...           
        [' B = -5 is a good example                        '
         '                                                 '];
  helpfun(RttlStr,RhlpStr);                                 
%===========================================================================
elseif strcmp(action,'InputInfo3');
%===========================================================================
  RttlStr = 'C-parameter Info';
  RhlpStr =  ...           
        [' C = 20 is a good example                        '
         '                                                 '];
  helpfun(RttlStr,RhlpStr);                                 
%===========================================================================
elseif strcmp(action,'InputInfo4');
%===========================================================================
  RttlStr = 'Maximal time Info';
  RhlpStr =  ...           
        [' t = 10 is a good example                        '
         '                                                 '];
  helpfun(RttlStr,RhlpStr);                                 
%===========================================================================
elseif strcmp(action,'InputInfo5');
%===========================================================================
  RttlStr = 'Rayleigh Number Info';
  RhlpStr =  ...           
        [' r = Ra/Ra_c is Rayleigh number relative to its  '
         ' critical value (which equals to 657.5 for the   '
         ' onset of thermal convection in a fluid layer    '  
         ' heated from below.                              '
         '                                                 '
         '             important r values are:             '  
         ' 0<r<1         stable 0 attractor                '  
         ' 1<r<13.96     unstable 0, stable C1 & C2        ' 
         ' ----------------------------------------------- ' 
         ' r ~ 13.96     homoclinic orbit appears          '
         ' 13.96<r<24.06 transient chaos, stable C1 & C2   '  
         ' r ~ 24.06     transient chaos->chaotic attractor'  
         ' 24.06<r<24.74 three possible stable attractors: '
         '                        chaotic attractor, C1, C2' 
         ' ----------------------------------------------- ' 
         ' r ~ 24.74     C1 & C2 become unstable,          '
         '                                 Hopf bifurcation'  
         ' r > 24.74     Lorenz Chaotic Attractor stable   '];
  helpfun(RttlStr,RhlpStr);                                 
end;
