function ss3d(action)
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%                   ss3d.m
% 
%    Plot the trajectory of the 3D Swingin Spring.
%
%    Written by Peter Lynch (Peter.Lynch@met.ie  )
%
%    Adapted from MATLAB demonstration program
%    lorenz, written for Expo by Ned Gulley, 6-21-93.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%   Set the initial conditions.
    x0 = 0.01;      xdot0 = 0.00;
    y0 = 0.00;      ydot0 = 0.02;
    zprime0 = 0.1;  zdot0 = 0.00;

%   Change to alternative initial conditions as required.
%   These ones give a nice six-pointed star for the projection.
    x0 = 0.04;      xdot0 = 0.00;
    y0 = 0.00;      ydot0 = 0.034272;
    zprime0 = 0.08; zdot0 = 0.00;

%   Change scale to alter character of the motion.
    skale = 1.00; 


  yinit1=x0*skale;          yinit2=xdot0*skale;
  yinit3=y0*skale;          yinit4=ydot0*skale;
  yinit5=-1+zprime0*skale;  yinit6=zdot0*skale;    
  yinit = [ yinit1 yinit2 yinit3 yinit4 yinit5 yinit6 ];

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Information regarding the play status will be held in
% the axis user data according to the following table:
play= 1;
stop=-1;

if nargin<1,
    action='initialize';
end;

if strcmp(action,'initialize'),
    oldFigNumber=watchon;

    figNumber=figure( ...
        'Name','Swingin Spring', ...
        'NumberTitle','off', ...
        'Color',[0.15 0.15 0.30], ...
	'Visible','off');

    axes( ...
        'Units','normalized', ...
        'Position',[0.10 0.10 0.7 0.85], ...
        'Visible','off');

    text(0,0,'Press the "Start" button to see the Swingin Spring', ...
	'HorizontalAlignment','center');
    axis([-1 1 -1 1]);

    %===================================
    % Information for all buttons
    labelColor=[0.8 0.8 1.0];
    yInitPos=0.90;
    xPos=0.85;
    btnLen=0.10;
    btnWid=0.10;
    % Spacing between the button and the next command's label
    spacing=0.05;

    %====================================
    % The CONSOLE frame
    frmBorder=0.02;
    yPos=0.05-frmBorder;
    frmPos=[xPos-frmBorder yPos btnLen+2*frmBorder 0.9+2*frmBorder];
    h=uicontrol( ...
        'Style','frame', ...
        'Units','normalized', ...
        'Position',frmPos, ...
	'BackgroundColor',[0.70 0.70 0.80]);        

    %====================================
    % The START button
    btnNumber=1;
    yPos=0.90-(btnNumber-1)*(btnWid+spacing);
    labelStr='Start';
    cmdStr='start';
    callbackStr='ss3d(''start'');';
                      
    % Generic button information
    btnPos=[xPos yPos-spacing btnLen btnWid];
    startHndl=uicontrol( ...
        'Style','pushbutton', ...
        'BackGroundColor','b', ...
        'Units','normalized', ...
        'Position',btnPos, ...
        'String',labelStr, ...
        'Interruptible','yes', ...
        'Callback',callbackStr);

    %====================================
    % The STOP button
    btnNumber=2;
    yPos=0.90-(btnNumber-1)*(btnWid+spacing);
    labelStr='Stop';
    % Setting userdata to -1 (=stop) will stop the program.
    callbackStr='set(gca,''Userdata'',-1)';
        
    % Generic  button information
    btnPos=[xPos yPos-spacing btnLen btnWid];
    stopHndl=uicontrol( ...
        'Style','pushbutton', ...
        'Units','normalized', ...
        'Position',btnPos, ...
        'Enable','off', ...
        'String',labelStr, ...
        'Callback',callbackStr);

    %====================================
    % The INFO button
    labelStr='Info';
    callbackStr='ss3d(''info'')';
    infoHndl=uicontrol( ...
        'Style','push', ...
        'Units','normalized', ...
        'position',[xPos 0.20 btnLen 0.10], ...
        'string',labelStr, ...
        'call',callbackStr);

    %====================================
    % The CLOSE button
    labelStr='Close';
    callbackStr='close(gcf)';
    closeHndl=uicontrol( ...
        'Style','push', ...
        'Units','normalized', ...
        'position',[xPos 0.05 btnLen 0.10], ...
        'string',labelStr, ...
        'call',callbackStr);
        
    % Uncover the figure
    hndlList=[startHndl stopHndl infoHndl closeHndl];
    set(figNumber,'Visible','on', ...
	'UserData',hndlList);
    
    watchoff(oldFigNumber);
    figure(figNumber);

elseif strcmp(action,'start'),
    axHndl=gca;    
    figNumber=gcf;
    hndlList=get(figNumber,'UserData');
    startHndl=hndlList(1);
    stopHndl=hndlList(2);
    infoHndl=hndlList(3);
    closeHndl=hndlList(4);
    set([startHndl closeHndl infoHndl],'Enable','off');
    set(stopHndl,'Enable','on');


    % ====== Start of Integration
    set(figNumber,'Backingstore','off');
    % The graphics axis limits are set to values known
    % to contain the solution.
    
      XLimn=-0.2*skale; XLimx=+0.2*skale; 
      YLimn=-0.2*skale; YLimx=+0.2*skale; 
      ZLimn=-1+(1-1.3)*skale; ZLimx=-1+(1-0.9)*skale; 
      tickx=linspace(XLimn,XLimx,5);
      ticky=linspace(YLimn,YLimx,5);
      tickz=linspace(ZLimn,ZLimx,5);
      set(axHndl, ...
        'XLim',[XLimn XLimx],'YLim',[YLimn YLimx],'ZLim',[ZLimn ZLimx], ...
        'XColor','b','YColor','b','ZColor','b',  ...
	'Userdata',play, ...
        'XTick',tickx,'YTick',ticky,'ZTick',tickz, ...
	'Drawmode','fast', ...
        'Visible','on', ...
	'NextPlot','add', ...
        'Userdata',play, ...
        'View',[-30,30]     );
         %%% 'View',[-37.5,30]);
         xlabel('X');
         ylabel('Y');
         zlabel('Z');   
         grid on

    % Values of the  physical spring parameters.
      m = 1;        %%% mass on spring
      g = pi^2;     %%% gravity
      k = 4*pi^2;   %%% stiffness 
      l = 1;        %%% length at equilibrium.

    epssq=(m*g)/(k*l);
    epsilon=sqrt(epssq);
    l_0=(1-epssq)*l;

    omega_R = sqrt(g/l);  omega_Z = sqrt(k/m);
    tau_R = 2*pi/omega_R; tau_Z = 2*pi/omega_Z;

    % The orbit ranges (regularly or chaotically) back and forth around
    % the stable equilibrium point.  It is bounded, but not periodic and
    % not convergent. The numerical integration, and the display of the 
    % evolving solution, are handled by the method of ODE23P (inline).


    t0=0;
    tfinal=1000;
    pow = 1/3;
    tol = 0.0001;  %%% reduced by factor of 10. 
    
    t = t0;
    hmax = (tfinal - t)/5;
    hmin = (tfinal - t)/200000;
    h = (tfinal - t)/100;
    y = yinit(:);
    tau = tol * max(norm(y,'inf'),1);
    
    % Save L steps and plot like a comet tail.
    L = 50;
    Y = y*ones(1,L);

    cla;
    head = line( ...
	'color','r', ...
	'linestyle','.', ...
	'Markersize',50, ...
	'erase','xor', ...
      	'xdata',y(1),'ydata',y(3),'zdata',y(5));
    cord = line( ...
	'color','r', ...
	'linestyle','-', ...
	'linewidth',2.0, ...
	'erase','xor', ...
      	'xdata',[0;y(1)],'ydata',[0;y(3)],'zdata',[0;y(5)]);
    body = line( ...
	'color','g', ...
	'linestyle','-', ...
	'erase','none', ...
	'xdata',[],'ydata',[],'zdata',[]);
    tail=line( ...
	'color','m', ...
	'linestyle','-', ...
	'erase','none', ...
	'xdata',[],'ydata',[],'zdata',[]);
    base = line( ...
	'color','y', ...
	'linestyle','-', ...
	'erase','none', ...
	'xdata',[],'ydata',[],'zdata',[]);
    
    % The main loop
    while (get(axHndl,'Userdata')==play) & (h >= hmin)
   	if t + h > tfinal, h = tfinal - t; end
   	% Compute the slopes
        %%%%%%%%%   FunFcn='sprroc';
% s1 = (feval(FunFcn, t, y))';
    yy = y;
    r=sqrt(yy(1)^2+yy(3)^2+yy(5)^2);  fac=omega_Z^2*(r-l_0)/r;
    ydot(1) = yy(2);  ydot(2) = - fac*yy(1);
    ydot(3) = yy(4);  ydot(4) = - fac*yy(3);
    ydot(5) = yy(6);  ydot(6) = - fac*yy(5) - g;
    s1 = ydot';
% s2 = (feval(FunFcn, t+h, y+h*s1))';
    yy = y+h*s1;
    r=sqrt(yy(1)^2+yy(3)^2+yy(5)^2);  fac=omega_Z^2*(r-l_0)/r;
    ydot(1) = yy(2);  ydot(2) = - fac*yy(1);
    ydot(3) = yy(4);  ydot(4) = - fac*yy(3);
    ydot(5) = yy(6);  ydot(6) = - fac*yy(5) - g;
    s2 = ydot';
% s3 = (feval(FunFcn, t+h/2, y+h*(s1+s2)/4))'; 
    yy = y+h*(s1+s2)/4;
    r=sqrt(yy(1)^2+yy(3)^2+yy(5)^2);  fac=omega_Z^2*(r-l_0)/r;
    ydot(1) = yy(2);  ydot(2) = - fac*yy(1);
    ydot(3) = yy(4);  ydot(4) = - fac*yy(3);
    ydot(5) = yy(6);  ydot(6) = - fac*yy(5) - g;
    s3 = ydot';

        % Estimate the error and the acceptable error
      	delta = norm(h*(s1 - 2*s3 + s2)/3,'inf');
      	tau = tol*max(norm(y,'inf'),1.0);
       	
        % Update the solution only if the error is acceptable
      	ts = t;
      	ys = y;
      	if delta <= tau
            t = t + h;
            y = y + h*(s1 + 4*s3 + s2)/6;
            
            % Update the plot
            Y = [y Y(:,1:L-1)];
            zbase=ZLimn;

            set(head,'xdata',Y(1,1),'ydata',Y(3,1),'zdata',Y(5,1))
            set(cord,'xdata',[0;Y(1,1)],'ydata',[0;Y(3,1)],'zdata',[0;Y(5,1)])
            set(body,'xdata',Y(1,1:2),'ydata',Y(3,1:2),'zdata',Y(5,1:2))
            set(tail,'xdata',Y(1,L-1:L),'ydata',Y(3,L-1:L),'zdata',Y(5,L-1:L))
            set(base,'xdata',Y(1,1:2),'ydata',Y(3,1:2),'zdata',[zbase,zbase])
            drawnow;
        end

       	% Update the step size
      	if delta ~= 0.0
            h = min(hmax, 0.9*h*(tau/delta)^pow);
      	end

    end;  % Main loop ... 
    % ====== End of Program
    set([startHndl closeHndl infoHndl],'Enable','on');
    set(stopHndl,'Enable','off');

elseif strcmp(action,'info');
    ttlStr='Swingin Spring Info';
    hlpStr= ...
        ['                                               '
         ' This program animates the integration of the  '
         ' six coupled nonlinear differential equations  '
         ' that determine the motion of the elastic      '
         ' pendulum, more colloquially called the        ' 
         '             "Swingin Spring"                  ' 
         ' As the integration proceeds you will see a    '
         ' point moving around in a curious orbit in     '
         ' 3-D space. The red blob is the pendulum bob.  '
         ' The green streak shows the recent track.      '
         ' The magenta shows the complete 3D trajectory. '
         ' The yellow is a projection of the motion on   '
         ' a horizontal plane.                           '
         '                                               '
         ' The initial conditions are in the array yinit.'
         ' They can easily be changed by changing the    '
         ' values x0, y0, z0, xdot0, ydot0,zdot0 or the  '
         ' value of the parameter "skale", defined       '
         ' near the top of the program. For small values '
         ' of skale the motion is quasi-regular. For     '
         ' large values it is chaotic.                   '
         ' [ for references, see publications in         '
         '     http://www.maths.tcd.ie/~plynch  ]        '
         '                                               '
         ' File name: ss3d.m                             '];

        helpfun(ttlStr,hlpStr);

end;    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%