function [g_msg,g_status]=process_eme_121_lab1(data) % called by the Matlab GUI to solve equation of motion Lab 1 MAE 121 % This function is called by the Matlab GUI to solve % the equation of motion for Lab 1 problem, EME 121, % spring 2011, UC Davis. % % INPUT: data, a record which contains the problem parameters % % The function will call ode45 solver and plot the solution % % by Nasser M. Abbasi % 4/12/2011 % Reset axes for plotting cla(data.handles.axes,'reset'); cla(data.handles.time_angle_axes,'reset'); cla(data.handles.time_r_axes,'reset'); cla(data.handles.angle_r_axes,'reset'); %Check if using auto step size of an explicit step size. see %GUI for input choice. if data.time_step==0 time_span=[0 data.max_t]; else time_span=0:data.time_step/1000:data.max_t; end %set the y-axis for plotting max_y = data.L; %set the initial conditions for ode45 IC = [data.r_zero;data.r_speed_zero;data.angle_zero;data.angle_speed_zero]; options = odeset('OutputFcn',@output); g_status = true; g_msg=''; %call the solver [t,x] = feval(data.solver,@rhs,time_span,IC,options,data); %plot the result N=length(t); picks=1:N; set(data.handles.figure1, 'CurrentAxes',data.handles.time_angle_axes); angle=x(:,3); plot(t(picks),angle(picks)*180/pi,'-'); title('pendulum angle vs time'); xlabel('Time (sec)'); ylabel('angle (degree)'); grid; drawnow; set(data.handles.figure1, 'CurrentAxes',data.handles.time_r_axes); r=x(:,1); stairs(t(picks),r(picks),'-'); title('pendulum length vs time'); xlabel('Time (sec)'); ylabel('length (m)'); grid; drawnow; set(data.handles.figure1, 'CurrentAxes',data.handles.angle_r_axes); plot(angle(picks)*180/pi,r(picks),'-'); title('pendulum angle vs length'); xlabel('Angle (degree)'); ylabel('length (m)'); grid; drawnow; %----------------------------------- % %----------------------------------- function dx=rhs(~,x,data) %the ode45 function temp = x(1)-data.L; if abs(temp<2*eps) term = 0; else term = data.k/data.mass*temp; end dx=[x(2); x(1)*x(4)^2+data.g*cos(x(3))-term; x(4); -2*x(2)*x(4)/x(1)- data.g/x(1)*sin(x(3)) ]; if x(1)>max_y max_y = x(1)+0.2*x(1); end end %----------------------------------- % %----------------------------------- function status= output(t,x,~,data) %called by ode45 after each step. Plot the current %pendulum position for simulation userData = get(data.handles.figure1, 'UserData'); if userData.stop status=true; g_status =true; else status = false; if not(isempty(t)) plot_solution(x(1,end),x(3,end),t(1),data); if x(1,end)>10*data.L %spring stretched too far status = true; g_status = false; g_msg ='spring over stretched'; end end end end %----------------------------------- % %----------------------------------- function plot_solution(r,theta,t,data) %plot the final solution set(0,'CurrentFigure',data.handles.figure1); set(data.handles.figure1, 'CurrentAxes',data.handles.axes); [x_local,y] = pol2cart(theta-pi/2,r); plot(x_local,y,'o','MarkerSize',20,'MarkerEdgeColor','r',... 'MarkerFaceColor','r','LineWidth',2); [x_local,y] = nma_spring.make(r,theta-(pi/2),30); hold on; line(x_local,y); plot(0,0,'o','MarkerSize',5,'MarkerEdgeColor','k',... 'MarkerFaceColor','k','LineWidth',1); title(sprintf(... 'time = %3.4f, angle(degree) = %3.3f, r = %3.3f, k=%3.4f',... t,theta*180/pi,r,data.k)); xlim([-max_y max_y]); ylim([-max_y max_y]); %refreshdata grid; drawnow; hold off; end end