%
%  manif.m -- 21 jan 07 -- djm
%
%  - ODE solver for inertial manifold
%  - NOTE:  cntl-c to exit IV loop
%
%- - - - - - - - - - - - - - - - - - - - - - - - 

%  parameters (damping)
global rr epsi
%epsi = (1/1);  rr = (0.25)^2;  BX = 1.0;
epsi = (1/7);  rr = (0.25)^2;  BX = 1.0;

%  ODE, timing & accuracy
Xode = 'Xmanif';  T = 4;  cycles = 15;  tol = 1e-3;

figure(1);  clf;  
x1 = -BX:BX/20:BX;  plot(x1,-x1.^2,'r')
title('\bf slow manifold phase plane','fontsize',12)
xlabel('\bf u_1(t)','fontsize',12)
ylabel('\bf u_2(t)','fontsize',12)

figure(2);  clf

%  multi-IV loop
while(1==1)
    %  too fancy IV input (discretized for repeatability)
    disp(' ');  disp(' click to choose initial point');  disp(' ');

    in = [];  IV = [2 1]*BX;
    figure(1);  hold on;  box 
    axis equal;  axis([-1 1 -1 1]*BX)
    while (isempty(in))
      in = floor(ginput(1)/(BX/25))*(BX/25); 

      if(isempty(in))
        break
      else
        plot(IV(1),IV(2),'w.')
        disp(['  IV:  ' num2str(in) '  (rtn to confirm)']); 
        IV = in;  in = [];
        plot(IV(1),IV(2),'r.')
      end
    end
    plot(IV(1),IV(2),'w.');  plot(IV(1),IV(2),'rx');

    %  ODE solver options
    options = odeset('reltol',tol*1e-3,'abstol',tol*1e-6, ...
            'maxstep',T/5,'stats','off','events','off');

    %  plot IV
    figure(2);  clf;  hold on
    plot(0,IV(1),'b.');  plot(0,IV(2),'r.')

    %  ODE solve & plot
    for C = 1:cycles
      [t y] = ode15s(Xode,[(C-1)*T C*T],IV,options);

      %  update IV by last point
      IV = y(end,:);  

      figure(1);
      plot(y(:,1),y(:,2),'b');

      figure(2);
      plot(t,y(:,1),'b')
      plot(t,y(:,2),'r')
      plot(t(end),y(end,1),'b.');
      plot(t(end),y(end,2),'r.');
      title('\bf slow manifold solutions','fontsize',12)
      xlabel('\bf time','fontsize',12)
      ylabel('\bf u_1(t) = blue ; u_2(t) = red','fontsize',12)

      drawnow
    end

    figure(1)
    plot(y(end,1),y(end,2),'k.');
end