function fosc(action, subcase) % FOSC Gui to demonstrate a forced oscillator. % The equation of motion is % % 2 % d x dx % m --- + r -- = F % 2 dt % dt % % where % % m = mass of oscillator (block) % % r = coefficient of friction % % F = -k e - alpha e^2 tan(e/2) = force exerted by spring on m % % k = spring constant % % e = x + R sin( omega t ) = extension of spring % % R = amplitude of driving oscillation % % omega = angular frequency of driving oscillation % % alpha = coefficient of nonlinear term % % Bill Mclean, UNSW, 5/2000. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Persistent variables % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% persistent control animation curves ... % figure handles go stop examples exmenu ... timestep parameter initcond ... curve_type cumenu axmenu ... % ui handles m r k omega R alpha ... % parameter values paramhelp ... % tooltip string exno exdata ranges ... % example data play ... % stop/go flag rod dot spr mass ... % animation handles s w h xmax xdotmax Rmax a b ... % animation parameters xvt_axes xvt_line ... phase_axes phase_line ... force_axes force_line ... % plotting handles steps_per_plot ... % scale for t axis x0 xdot0 ... % initial conditions x xdot t ... % current state prevx prevxdot prevt ... % previous state j dt delay ... % time stepping %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if nargin == 0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Initialisations % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Figures % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% if isempty(control) control = figure('Resize', 'off', 'Name', ... 'Forced Oscillator - control panel', ... 'NumberTitle', 'off', 'Pointer', 'watch', ... 'Color', [1 .85 .5], ... 'Position', [90 60 300 440], ... 'MenuBar', 'none'); else fprintf(1, '\a') errordlg('fosc GUI already exists') return end animation = figure('Name', 'Forced Oscillator - animation', ... 'NumberTitle', 'off', 'MenuBar', 'none', ... 'Pointer', 'watch', ... 'Position', [400 560 560 200]); s = .1; w = .2; h = .3; a = 3; b = 4; xmax = 4; xdotmax = 4; Rmax = 1; curves = figure('Name', 'Forced Oscillator - curves', ... 'NumberTitle', 'off', 'MenuBar', 'none', ... 'Pointer', 'watch', ... 'Position', [400 60 560 440]); figure(control) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Examples % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% examples = uimenu(control, 'Label', 'Examples'); exmenu(1) = uimenu(examples, 'Label', 'Simple harmonic motion', ... 'Interruptible', 'off', ... 'Checked', 'on', 'CallBack', 'fosc examples'); exmenu(2) = uimenu(examples, 'Label', 'Resonance', ... 'Interruptible', 'off', ... 'CallBack', 'fosc examples'); exmenu(3) = uimenu(examples, 'Label', 'Near resonance (beats)', ... 'Interruptible', 'off', ... 'CallBack', 'fosc examples'); exmenu(4) = uimenu(examples, 'Label', 'Underdamped motion', ... 'Interruptible', 'off', ... 'CallBack', 'fosc examples'); exmenu(5) = uimenu(examples, 'Label', 'Overdamped motion', ... 'Interruptible', 'off', ... 'CallBack', 'fosc examples'); exmenu(6) = uimenu(examples, 'Label', 'Forced damped motion', ... 'Interruptible', 'off', ... 'CallBack', 'fosc examples'); exmenu(7) = uimenu(examples, 'Label', 'Nonlinear oscillation', ... 'Interruptible', 'off', ... 'CallBack', 'fosc examples'); exno = 1; exdata = [1.5 0 0.5 0 0 0 2 0 1 0 1 1 0.1 0 0 0 1 0 1 0.85 0.5 0 0 0 1 0.15 1 0 0 0 2 -2 1 0.21 0.01 0 0 0 3 -1.5 1.3 0.1 0.8 0.6 0.8 0 0 0 1 0.2 -0.36 1.04 0.5 0.09 0 0]; ranges = [0 0 -1 0 0 0 -xmax -xdotmax 2 2 1 2 Rmax 1 xmax xdotmax]; m = exdata(exno,1); r = exdata(exno,2); k = exdata(exno,3); omega = exdata(exno,4); R = exdata(exno,5); alpha = exdata(exno,6); x0 = exdata(exno,7); xdot0 = exdata(exno,8); steps_per_plot = 1000; dt = 0.1; delay = 0.02; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Help % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% uimenu(control, 'Label', 'Help', 'Interruptible', 'off', ... 'CallBack', 'helpwin fosc') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Go, Stop, Close % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% go = uicontrol('Style', 'pushbutton', ... 'Position', [20 360 60 60], ... 'String', 'Go', ... 'Interruptible', 'on', ... 'CallBack', 'fosc go'); stop = uicontrol('Style', 'pushbutton', ... 'Position', [120 360 60 60], ... 'String', 'Stop', 'Enable', 'off', ... 'Interruptible', 'off', ... 'CallBack', 'fosc stop'); uicontrol('Style', 'pushbutton', ... 'Position', [220 360 60 60], ... 'String', 'Close', ... 'Interruptible', 'off', ... 'CallBack', 'fosc close'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Time step % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% uicontrol('Style', 'frame', 'Position', [20 280 260 60]) uicontrol('Style', 'text', 'String', 'Time step', ... 'Position', [40 290 40 30]) timestep(1) = uicontrol('Style', 'edit', 'String', num2str(dt), ... 'Position', [80 300 50 20], ... 'Interruptible', 'off', ... 'CallBack', 'fosc setdt'); uicontrol('Style', 'text', 'String', 'Pause', ... 'Position', [160 290 50 30], ... 'TooltipString', 'Time delay at each step') timestep(2) = uicontrol('Style', 'edit', 'String', num2str(delay), ... 'Position', [210 300 50 20], ... 'Interruptible', 'off', ... 'CallBack', 'fosc setpause'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Parameters % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% uicontrol('Style', 'frame', 'Position', [20 150 260 110]) uicontrol('Style', 'text', 'String', 'Parameters', ... 'Position', [40 210 60 30]) menustr = str2mat('m', 'r', 'k', 'omega', 'R', 'alpha'); paramhelp = {'mass of oscillator (block)', ... 'coefficient of friction', ... 'spring constant (coefficient of linear term)', ... 'angular frequency of driving oscillation', ... 'amplitude of driving oscillation', ... 'strength of nonlinearity'}; parameter(1) = uicontrol('Style', 'popupmenu', ... 'String', menustr, 'Value', 1, ... 'TooltipString', paramhelp{1}, ... 'Position', [110 220 150 30], ... 'Interruptible', 'off', ... 'CallBack', 'fosc(''parameter'', ''select'')'); parameter(2) = uicontrol('Style', 'slider', ... 'Position', [40 190 220 15], ... 'Value', m, 'Min', ranges(1,1), 'Max', ranges(2,1), ... 'Interruptible', 'off', ... 'CallBack', 'fosc(''parameter'', ''slider'')'); parameter(3) = uicontrol('Style', 'text', ... 'Position', [40 160 40 20], ... 'HorizontalAlignment', 'left', ... 'String', num2str(get(parameter(2), 'Min'))); parameter(4) = uicontrol('Style', 'edit', ... 'Position', [130 160 60 20], ... 'String', num2str(get(parameter(2), 'Value')), ... 'Interruptible', 'off', ... 'CallBack', 'fosc(''parameter'', ''edit'')'); parameter(5) = uicontrol('Style', 'text', ... 'Position', [220 160 40 20], ... 'HorizontalAlignment', 'right', ... 'String', num2str(get(parameter(2), 'Max'))); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Initial conditions % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% uicontrol('Style', 'frame', 'Position', [20 20 260 110]) uicontrol('Style', 'text', 'String', 'Initial Conditions', ... 'Position', [40 80 60 30]) menustr = str2mat('x0', 'xdot0'); initcond(1) = uicontrol('Style', 'popupmenu', ... 'String', menustr, 'Value', 1, ... 'Position', [110 90 150 30], ... 'Interruptible', 'off', ... 'CallBack', 'fosc(''initcond'', ''select'')'); initcond(2) = uicontrol('Style', 'slider', ... 'Position', [40 60 220 15], ... 'Value', m, 'Min', ranges(1,7), 'Max', ranges(2,7), ... 'Interruptible', 'off', ... 'CallBack', 'fosc(''initcond'', ''slider'')'); initcond(3) = uicontrol('Style', 'text', ... 'Position', [40 30 40 20], ... 'HorizontalAlignment', 'left', ... 'String', num2str(get(initcond(2), 'Min'))); initcond(4) = uicontrol('Style', 'edit', ... 'Position', [130 30 60 20], ... 'String', num2str(get(initcond(2), 'Value')), ... 'Interruptible', 'off', ... 'CallBack', 'fosc(''initcond'', ''edit'')'); initcond(5) = uicontrol('Style', 'text', ... 'Position', [220 30 40 20], ... 'HorizontalAlignment', 'right', ... 'String', num2str(get(initcond(2), 'Max'))); j = 0; t = 0; prevt = t; x = x0; prevx = x; xdot = xdot0; prevxdot = xdot; fosc draw setup %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Curve menus % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% curve_type = uimenu(curves, 'Label', 'Select'); cumenu(1) = uimenu(curve_type, 'Label', 'x versus t', ... 'Interruptible', 'off', ... 'Checked', 'on', 'CallBack', 'fosc curve xvt'); cumenu(2) = uimenu(curve_type, 'Label', 'phase plane', ... 'Interruptible', 'off', ... 'CallBack', 'fosc curve phase'); cumenu(3) = uimenu(curve_type, 'Label', 'Spring response', ... 'Interruptible', 'off', ... 'CallBack', 'fosc curve response'); axmenu = uimenu(curves, 'Label', 'Axes', ... 'Interruptible', 'off', ... 'CallBack', 'fosc setaxes'); uimenu(curves, 'Label', 'Refresh', 'Interruptible', 'off', ... 'CallBack', 'refresh') fosc curve xvt watchoff(control) watchoff(animation) watchoff(curves) else %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Process callback % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% switch(action) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % case 'examples' % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% set(exmenu(exno), 'Checked', 'off') exno = get(gcbo, 'Position'); m = exdata(exno,1); r = exdata(exno,2); k = exdata(exno,3); omega = exdata(exno,4); R = exdata(exno,5); alpha = exdata(exno,6); x0 = exdata(exno,7); xdot0 = exdata(exno,8); x = x0; prevx = x; xdot = xdot0; prevxdot = xdot; j = 0; t = 0; prevt = t; fosc('parameter', 'select') fosc('initcond', 'select') set(exmenu(exno), 'Checked', 'on') fosc('draw', 'update') refresh(curves) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % case 'go' % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% set(go, 'Enable', 'off') set(stop, 'Enable', 'on') set(examples, 'Enable', 'off') set(timestep, 'Enable', 'off') set(parameter([1 2 4]), 'Enable', 'off') set(initcond([1 2 4]), 'Enable', 'off') set(axmenu, 'Enable', 'off') set(curve_type, 'Enable', 'off') fosc play %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % case 'play' % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% play = 1; while play prevx = x; prevxdot = xdot; prevt = t; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % Classical 4-stage Runge-Kutta step % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% z = [x; xdot]; Phi1 = rhs(t, z, m, r, k, R, omega, alpha); Phi2 = rhs(t+dt/2, z+dt*Phi1/2, m, r, k, R, omega, alpha); Phi3 = rhs(t+dt/2, z+dt*Phi2/2, m, r, k, R, omega, alpha); Phi4 = rhs(t+dt, z+dt*Phi3, m, r, k, R, omega, alpha); z = z + (dt/6) * ( Phi1 + 2*Phi2 + 2*Phi3 + Phi4 ); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% j = j + 1; t = j * dt; x = z(1); xdot = z(2); fosc('draw', 'update') if abs(x) > xmax set(go, 'Enable', 'on') errordlg('Model breakdown: x exceeds physical limits') fosc stop end pause(delay) end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % case 'stop' % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% play = 0; set(go, 'String', 'Proceed', 'CallBack', 'fosc proceed', ... 'Enable', 'on') set(stop, 'String', 'Reset', 'CallBack', 'fosc reset') set(curve_type, 'Enable', 'on') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % case 'proceed' % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% set(go, 'Enable', 'off') set(stop, 'String', 'Stop', 'CallBack', 'fosc stop') set(curve_type, 'Enable', 'off') fosc play %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % case 'reset' % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% set(stop, 'String', 'Stop', 'CallBack', 'fosc stop', ... 'Enable', 'off') set(go, 'String', 'Go', 'CallBack', 'fosc go') j = 0; t = j * dt; prevt = t; x = x0; prevx = x; xdot = xdot0; prevxdot = xdot; fosc('draw', 'update') refresh(curves) drawnow set(examples, 'Enable', 'on') set(timestep, 'Enable', 'on') set(parameter([1 2 4]), 'Enable', 'on') set(initcond([1 2 4]), 'Enable', 'on') set(axmenu, 'Enable', 'on') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % case 'close' % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% close(control) close(animation) close(curves) clear fosc %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % case 'setdt' % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% dtstr = get(timestep(1), 'String'); try dtval = eval( dtstr ); catch errordlg(lasterr) return end dt = dtval; set(xvt_axes, 'XLim', [0 steps_per_plot*dt]) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % case 'setpause' % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% delaystr = get(timestep(2), 'String'); try delayval = eval( delaystr ); catch errordlg(lasterr) return end delay = delayval; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % case 'parameter' % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% selection = get(parameter(1), 'Value'); set(parameter(1), 'TooltipString', paramhelp{selection}) if strcmp(subcase, 'select') switch(selection) case 1 v = m; case 2 v = r; case 3 v = k; case 4 v = omega; case 5 v = R; case 6 v = alpha; end else switch(subcase) case 'edit' vstr = get(parameter(4), 'String'); try v = eval( vstr ); catch errordlg(lasterr) return end if ( v < ranges(1, selection) | ... v > ranges(2, selection) ) fprintf(1, '\a') errordlg('Value outside allowed range') return end case 'slider' v = get(parameter(2), 'Value'); end switch(selection) case 1 m = v; case 2 r = v; case 3 k = v; case 4 omega = v; case 5 R = v; case 6 alpha = v; end end update_sldr(parameter, selection, ranges, v) set(exmenu, 'Checked', 'off') fosc draw update %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % case 'initcond' % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% selection = get(initcond(1), 'Value'); if strcmp(subcase, 'select') switch(selection) case 1 v = x0; case 2 v = xdot0; end else switch(subcase) case 'edit' vstr = get(initcond(4), 'String'); try v = eval( vstr ); catch errordlg(lasterr) return end if ( v < ranges(1, 6+selection) | ... v > ranges(2, 6+selection) ) fprintf(1, '\a') errordlg('Value outside allowed range') return end case 'slider' v = get(initcond(2), 'Value'); end switch(selection) case 1 x0 = v; case 2 xdot0 = v; end end update_sldr(initcond, 6+selection, ranges, v) x = x0; xdot = xdot0; fosc('draw', 'update') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % case 'curve' % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% set(xvt_axes, 'Visible', 'off') set(get(xvt_axes, 'Children'), 'Visible', 'off') set(phase_axes, 'Visible', 'off') set(get(phase_axes, 'Children'), 'Visible', 'off') set(force_axes, 'Visible', 'off') set(get(force_axes, 'Children'), 'Visible', 'off') set(cumenu, 'Checked', 'off') switch(subcase) case 'xvt' set(xvt_axes, 'Visible', 'on') set(get(xvt_axes, 'Children'), 'Visible', 'on') set(cumenu(1), 'Checked', 'on') set(curves, 'CurrentAxes', xvt_axes) case 'phase' set(phase_axes, 'Visible', 'on') set(get(phase_axes, 'Children'), 'Visible', 'on') set(cumenu(2), 'Checked', 'on') set(curves, 'CurrentAxes', phase_axes) case 'response' set(force_axes, 'Visible', 'on') set(get(force_axes, 'Children'), 'Visible', 'on') set(cumenu(3), 'Checked', 'on') set(curves, 'CurrentAxes', force_axes) otherwise error('Unknown subcase') end fosc draw update %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % case 'setaxes' % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure(curves) xlim = [ '[' num2str( get(gca, 'XLim') ) ']' ]; ylim = [ '[' num2str( get(gca, 'YLim') ) ']' ]; switch(gca) case xvt_axes prompt = {'Time steps per plot:', 'Vertical axis range:'}; answer = inputdlg(prompt, 'Axes dialog box', 1, ... {num2str(steps_per_plot) ylim}); if ~isempty(answer) steps_per_plot = eval(answer{1}); ylim = eval(answer{2}); set(gca, 'XLim', [0 steps_per_plot*dt], 'YLim', ylim); end case {phase_axes, force_axes} prompt = {'Horizontal axis range:', 'Vertical axis range:'}; answer = inputdlg(prompt, 'Axes dialog box', 1, {xlim ylim}); if ~isempty(answer) xlim = eval(answer{1}); ylim = eval(answer{2}); set(gca, 'XLim', xlim, 'YLim', ylim); end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % case 'draw' % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% theta = linspace(0,2*pi); slot = -a-b-R*sin(omega*t); rod1 = [slot-2*s slot-2*s slot+2*s slot+2*s slot+a NaN slot+2*s ... slot+2*s slot-2*s]; rod2 = [-Rmax-2*s Rmax+2*s Rmax+2*s 0 0 NaN 0 ... -Rmax-2*s -Rmax-2*s]; dot1 = -a-b-R*sin(omega*t)+s*[-1 -1 1 1 -1]; dot2 = R*cos(omega*t)+s*[-1 1 1 -1 -1]; tau = linspace(0, 16*pi,80); spr1 = slot + a + tau * ( x - slot - a - h ) / tau(end); spr2 = w * sin(tau); mass1 = [x-h x-h x+h x+h x-h]; mass2 = [-h h h -h -h]; switch(subcase) case 'setup' figure(animation) clf axes('XLim', [-a-b-1.5*Rmax xmax+Rmax/2], 'YLim', [-1.5*Rmax 1.5*Rmax], ... 'DataAspectRatio', [1 1 1], ... 'PlotBoxAspectRatioMode', 'auto', 'DrawMode', 'fast') xlabel('x') unitcircle = [cos(theta); sin(theta)]; line(-a-b+(Rmax+2*s)*unitcircle(1,:), (Rmax+2*s)*unitcircle(2,:), ... 'Color', 'g'); rod = line(rod1, rod2, 'Color', 'k'); dot = patch(dot1, dot2, 'g', 'EdgeColor', 'none'); spr = line(spr1, spr2, 'Color', 'k'); mass = patch(mass1, mass2, 'k', 'EdgeColor', 'none'); line([-b-Rmax xmax], [-h -h], 'Color', 'k') set(get(gca, 'Children'), 'EraseMode', 'xor'); figure(curves) xvt_axes = axes('XLim', [0 steps_per_plot*dt], ... 'YLim', [-b-Rmax xmax], 'DrawMode', 'fast'); xlabel('t'), ylabel('x') xvt_line = line(0, x0, 'EraseMode', 'none'); phase_axes = axes('XLim', [-xmax xmax], ... 'YLim', [-xdotmax xdotmax], ... 'DrawMode', 'fast'); title('Phase Plane'), xlabel('x'), ylabel('xdot') phase_line = line(x0, xdot0, 'Color', 'r', 'EraseMode', 'none'); force_axes = axes('XLim', response('erange'), ... 'YLim', response('Frange'), ... 'XGrid', 'on', 'YGrid', 'on'); title('Spring response'), xlabel('extension'), ylabel('force') e = response('erange'); e = linspace(e(1), e(2), 50); F = response(e, k, alpha); force_line = line(e, F); set(curves, 'CurrentAxes', xvt_axes) case 'update' set(rod, 'XData', rod1, 'YData', rod2); set(dot, 'XData', dot1, 'YData', dot2); set(spr, 'XData', spr1, 'YData', spr2); set(mass, 'XData', mass1, 'YData', mass2); switch(get(curves, 'CurrentAxes')) case xvt_axes trange = get(xvt_axes, 'XLim'); if t < trange(1) | t > trange(2) t_axis_range = steps_per_plot * dt; l = floor( t / t_axis_range ); set(xvt_axes, 'XLim', [l*t_axis_range (l+1)*t_axis_range]) end set(xvt_line, 'XData', [prevt t], 'YData', [prevx x]) case phase_axes set(phase_line, 'XData', [prevx x], 'YData', [prevxdot xdot]) case force_axes e = response('erange'); e = linspace(e(1), e(2), 50); F = response(e, k, alpha); set(force_line, 'XData', e, 'YData', F); otherwise errordlg('Unknown axes handle') end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % end % switch(action) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % end % if else % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function update_sldr(hndl, selection, ranges, v) set(hndl(2), 'Value', v); set(hndl(2), 'Min', ranges(1, selection)) set(hndl(2), 'Max', ranges(2, selection)) set(hndl(3), 'String', num2str(get(hndl(2), 'Min'))) set(hndl(4), 'String', num2str(v)) set(hndl(5), 'String', num2str(get(hndl(2), 'Max'))) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function zdot = rhs(t, z, m, r, k, R, omega, alpha); x = z(1); xdot = z(2); e = x + R * sin(omega*t); zdot = [ xdot; m \ ( response(e, k, alpha) - r*xdot ) ]; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function F = response(varargin) switch(nargin) case 1 flag = varargin{1}; switch(flag) case 'erange' F = [-3 3]; case 'Frange' F = [-2 2]; otherwise error('unknown flag') end case 3 [e, k, alpha] = deal(varargin{:}); F = -k*e-alpha*e.^2.*tan(e/2); otherwise error('illegal number of arguments') end