function nma_steady_state() %simulation of steady state single degree of freedom system with harmonic input %by Nasser M. Abbasi, %version March 1, 2013 %free to copy and use as long as this note remains in place. %initial values. Can be changed using sliders maxTime = 14; k = 1000; zeta = 0.01; m = 915; omega = 1.11; %hz forcingPeriod = 1; timeStep = forcingPeriod/20; data = zeros(length(0:timeStep:maxTime),3); omegaRadians = omega*2*pi; omegaN = sqrt(k/m); r = 0; phase = 0; currentStep = 0; stateOfSimulation = 'PAUSED'; magnificationFactor = 0; phase = 0; forceVectorLastChangedIndex=1; set(0,'Units','Normalized'); f=figure('Units','normalized','Position',[.3 .3 .42 .42]); set(f,'name','steady state simulation, version 030113',... 'numbertitle','off'); axesCommonOptions={'Units','normalized','Parent',f,'FontName',... 'FixedWidth','FontUnits','points','FontSize',9}; h1=axes(axesCommonOptions{:},'Position',[.54 .52 .40 .38]); set(h1,'Title',text('String','response','Color','k')); h2=axes(axesCommonOptions{:},'Position',[.54 .07 .40 .32],... 'XTick', [],'YTick', []); set(h2,'Title',text('String','phase','Color','k')); hMaxForceDispAxes=axes(axesCommonOptions{:},'Position',... [.23 .04 .28 .28]); set(hMaxForceDispAxes,'Title',text('String','f(t) vs. normalize response',... 'Color','k','FontSize',8),... 'XTick', [],'YTick', []); h3 = uipanel('Units','normalized','Parent',f,'Title',... 'system parameters','FontSize',9,... 'BackgroundColor','white',... 'Position',[.01 .40 .44 .58]); %zeta controls slidersWidth = 0.52; slidersHeight = 0.11; textRightOfSliderHeight = 0.10; sliderSep = 0.02; valueFieldWidth = 0.2; currentSliderH = .85; leftSliderTextWidth = 0.11; sliderLeftStart = 0.13; textRightSliderFontSize=9; h4 = uicontrol('Units','normalized','Parent',h3,'style',... 'slider','position',... [sliderLeftStart currentSliderH slidersWidth slidersHeight],... 'callback',@zetaCallback,'Min',0,'Max',2,'Value',... zeta,'SliderStep',[0.01/2 0.1/2]); h41=uicontrol('Units','normalized','Parent',h3,'position',... [.01 currentSliderH leftSliderTextWidth slidersHeight],... 'style','text','String',char(hex2dec('7A')),'fontname','symbol',... 'FontSize',12); h42=uicontrol('Units','normalized','Parent',h3,'position',... [slidersWidth+leftSliderTextWidth+.03 currentSliderH valueFieldWidth... textRightOfSliderHeight],... 'style','text','String','0.01','HorizontalAlignment',... 'left','FontSize',textRightSliderFontSize); %k controls currentSliderH=currentSliderH-(slidersHeight+sliderSep); h5 = uicontrol('Units','normalized','Parent',h3,'style',... 'slider','position',... [sliderLeftStart currentSliderH slidersWidth slidersHeight ],... 'callback',@sliderStiffness,'Min',10,... 'Max',10010,'Value',k,... 'SliderStep',[10/10000 100/10000]); h51=uicontrol('Units','normalized','Parent',h3,... 'position',[.01 currentSliderH leftSliderTextWidth slidersHeight],... 'style','text','String',char(hex2dec('6B')),'fontname','symbol',... 'FontSize',12); h52=uicontrol('Units','normalized','Parent',h3,... 'position',[slidersWidth+leftSliderTextWidth+.03 ... currentSliderH valueFieldWidth textRightOfSliderHeight],... 'style','text','String','1000','HorizontalAlignment','left',... 'FontSize',textRightSliderFontSize); h53=uicontrol('Units','normalized','Parent',h3,'position',... [slidersWidth+.35 currentSliderH .1 slidersHeight],... 'style','text','String','N/m','FontSize',8); %mass controls currentSliderH=currentSliderH-(slidersHeight+sliderSep); h6 = uicontrol('Units','normalized','Parent',h3,... 'style','slider','position',... [sliderLeftStart currentSliderH slidersWidth slidersHeight],... 'callback',@sliderMass,'Min',1,'Max',... 1001,'Value',m,... 'SliderStep',[1/1000 10/1000]); h61=uicontrol('Units','normalized','Parent',h3,'position',... [.01 currentSliderH leftSliderTextWidth slidersHeight],... 'style','text','String',char(hex2dec('4D')),... 'fontname','symbol','FontSize',12); h62=uicontrol('Units','normalized','Parent',h3,... 'position',[slidersWidth+leftSliderTextWidth+.03 currentSliderH... valueFieldWidth textRightOfSliderHeight],... 'style','text','String','915','HorizontalAlignment','left',... 'FontSize',textRightSliderFontSize); h63=uicontrol('Units','normalized','Parent',h3,'position',... [slidersWidth+.35 currentSliderH .1 slidersHeight],... 'style','text','String','kg','FontSize',8); %forcing frequency controls currentSliderH=currentSliderH-(slidersHeight+sliderSep); h70 = uicontrol('Units','normalized','Parent',h3,'style','slider',... 'position',... [sliderLeftStart currentSliderH slidersWidth slidersHeight],... 'callback',@sliderFreq,'Min',0.0001,'Max',3,'Value',omega,... 'SliderStep',[0.0001/(3-0.0001) 0.01/(3-0.0001)]); h71=uicontrol('Units','normalized','Parent',h3,'position',... [.01 currentSliderH leftSliderTextWidth slidersHeight],... 'style','text','String',char(hex2dec('76')),'fontname',... 'symbol','FontSize',12); h72=uicontrol('Units','normalized','Parent',h3,'position',... [slidersWidth+leftSliderTextWidth+.03 currentSliderH ... valueFieldWidth textRightOfSliderHeight],... 'style','text','String','1.11','HorizontalAlignment','left',... 'FontSize',textRightSliderFontSize); h73=uicontrol('Units','normalized','Parent',h3,'position',... [slidersWidth+.35 currentSliderH .1 slidersHeight],... 'style','text','String','Hz','FontSize',8); currentSliderH=currentSliderH-(slidersHeight+sliderSep); hForceFrequencyBtn=uicontrol('Units','Normalized','Parent',h3,... 'Style', 'pushbutton',... 'String',[char(hex2dec('76')) char(hex2dec('3D')) ... char(hex2dec('77'))],'fontname','symbol',... 'Position', [.3 currentSliderH .2 .12],... 'Callback', @forceFrequencyCallback,'FontSize',12); %max time currentSliderH=currentSliderH-(slidersHeight+sliderSep); h80 = uicontrol('Units','normalized','Parent',h3,'style',... 'slider','position',... [sliderLeftStart currentSliderH slidersWidth slidersHeight],... 'callback',@maxTimeCallback,'Min',1,... 'Max',1000,'Value',maxTime,... 'SliderStep',[1/(1000-1) 5/(1000-1)]); h81=uicontrol('Units','normalized','Parent',h3,'position',... [.01 currentSliderH leftSliderTextWidth slidersHeight],... 'style','text','String','time',... 'FontSize',9,'HorizontalAlignment','left'); h82=uicontrol('Units','normalized','Parent',h3,'position',... [slidersWidth+leftSliderTextWidth+.03 currentSliderH ... valueFieldWidth textRightOfSliderHeight],... 'style','text','String','1000','HorizontalAlignment',... 'left','FontSize',textRightSliderFontSize); h83=uicontrol('Units','normalized','Parent',h3,'position',... [slidersWidth+.35 currentSliderH .1 slidersHeight],... 'style','text','String','sec','FontSize',8); %run button currentSliderH=currentSliderH-(1.5*slidersHeight+sliderSep); h9=uicontrol('Units','Normalized','Parent',h3,'Style', 'pushbutton',... 'String', 'RUN',... 'Position', [.01 currentSliderH .2 .16],... 'Callback', @runCallback); %pause button h10=uicontrol('Units','Normalized','Parent',h3,'Style', 'pushbutton',... 'String', 'PAUSE',... 'Position', [.22 currentSliderH .2 .16],... 'Callback', @pauseCallback); %reset button h11=uicontrol('Units','Normalized','Parent',h3,'Style', 'pushbutton',... 'String', 'RESET',... 'Position', [.43 currentSliderH .2 .16],... 'Callback', @resetCallback); %step button h12=uicontrol('Units','Normalized','Parent',h3,'Style', 'pushbutton',... 'String', 'STEP',... 'Position', [.64 currentSliderH .2 .16],... 'Callback', @stepCallback); %information panel h13 = uipanel('Units','normalized','Parent',f,'Title',... 'calculated','FontSize',9,... 'BackgroundColor','white',... 'Position',[.01 .02 .20 .36]); %natural fequency leftSliderTextWidth=0.16; currentSliderH=0.8; slidersHeight=0.18; textRightSliderFontSize=10; hNaturalFrequency=uicontrol('Units','normalized','Parent',h13,... 'position',[0.01 currentSliderH leftSliderTextWidth slidersHeight ],... 'style','text','String',char(hex2dec('77')),'fontname',... 'symbol','FontSize',12,... 'HorizontalAlignment','left','FontSize',textRightSliderFontSize); hNaturalFrequency1=uicontrol('Units','normalized','Parent',h13,'position',... [0.18 currentSliderH .55 slidersHeight],... 'style','text','String','0','FontSize',textRightSliderFontSize,... 'HorizontalAlignment','left'); hNaturalFrequency2=uicontrol('Units','normalized','Parent',h13,... 'position',[0.76 currentSliderH .15 slidersHeight],... 'style','text','String','Hz','FontSize',8); %r currentSliderH=currentSliderH-slidersHeight-sliderSep; hr=uicontrol('Units','normalized','Parent',h13,... 'position',[0.01 currentSliderH leftSliderTextWidth slidersHeight ],... 'style','text','String','r','FontSize',12,... 'HorizontalAlignment','left','FontSize',textRightSliderFontSize); hr1=uicontrol('Units','normalized','Parent',h13,'position',... [0.18 currentSliderH .55 slidersHeight],... 'style','text','String','0','FontSize',textRightSliderFontSize,... 'HorizontalAlignment',... 'left'); %magnification factor currentSliderH=currentSliderH-slidersHeight-sliderSep; hMagnificationFactor=uicontrol('Units','normalized','Parent',h13,... 'position',[0.01 currentSliderH leftSliderTextWidth slidersHeight ],... 'style','text','String',char(hex2dec('62')),'fontname','symbol',... 'HorizontalAlignment','left','FontSize',textRightSliderFontSize); hMagnificationFactor1=uicontrol('Units','normalized','Parent',h13,... 'position',... [0.18 currentSliderH .55 slidersHeight],... 'style','text','String','0','FontSize',textRightSliderFontSize,... 'HorizontalAlignment','left'); %phase angle currentSliderH=currentSliderH-slidersHeight-sliderSep; hPhaseAngle=uicontrol('Units','normalized','Parent',h13,... 'position',[0.01 currentSliderH leftSliderTextWidth slidersHeight ],... 'style','text','String',char(hex2dec('66')),'fontname','symbol',... 'HorizontalAlignment','left','FontSize',textRightSliderFontSize); hPhaseAngle1=uicontrol('Units','normalized','Parent',h13,'position',... [0.18 currentSliderH .55 slidersHeight],... 'style','text','String','0','FontSize',textRightSliderFontSize,... 'HorizontalAlignment','left'); hPhaseAngle2=uicontrol('Units','normalized','Parent',h13,... 'position',[0.76 currentSliderH .18 slidersHeight],... 'style','text','String','deg','FontSize',8); resetAll(); %--------------------------------------------- function resetAll() omegaN = sqrt(k/m); set(hNaturalFrequency1,'String',num2str(omegaN/(2*pi))); r = omegaRadians/omegaN; set(hr1,'String',num2str(r)); updateMagnificationFactor(); set(hMagnificationFactor1,'String',num2str(magnificationFactor)); updatePhaseAngle(); set(hPhaseAngle1,'String',num2str(-phase*180/pi)); forceVectorLastChangedIndex=currentStep+1; makePlot(); cla(hMaxForceDispAxes); end %--------------------------------------------- function zetaCallback(hObject,eventdata) zeta = round(get(hObject,'Value')*10^(2))/10^(2); resetAll(); set(h42,'String', num2str(zeta),'HorizontalAlignment','left') end %--------------------------------------------- function sliderStiffness(hObject,eventdata) k = round(get(hObject,'Value')); resetAll(); set(h52,'String',num2str(k),... 'HorizontalAlignment','left') end %--------------------------------------------- function sliderMass(hObject,eventdata) m = round(get(hObject,'Value')); resetAll(); set(h62,'String', num2str(m),... 'HorizontalAlignment','left') end %--------------------------------------------- function sliderFreq(hObject,eventdata) omega = round(get(hObject,'Value')*10^4)/10^4; omegaRadians = omega*2*pi; forcingPeriod = 1/omega; timeStep = forcingPeriod/20; resetAll(); if currentStep>0 && currentStep < size(data,1) additionalNumberOfStepsNeeded = ... round((maxTime - data(currentStep,1))/timeStep); tmp = zeros(additionalNumberOfStepsNeeded,3); data = [data(1:currentStep,:) ; tmp]; end set(h72,'String', num2str(omega),... 'HorizontalAlignment','left'); omegaRadians = omega*2*pi; end %--------------------------------------------- function maxTimeCallback(hObject,eventdata) maxTime = round(get(hObject,'Value')); if currentStep==0 resetData(); else tmp = zeros(currentStep+length(data(currentStep,1)+... timeStep:timeStep:maxTime),3); tmp(1:currentStep,:) = data(1:currentStep,:); data = tmp; end; makePlot(); set(h82,'String', num2str(maxTime ),... 'HorizontalAlignment','left') end %--------------------------------------------- function forceFrequencyCallback(hObject,eventdata) omega=omegaN/(2*pi); omegaRadians=omegaN; set(h72,'String', num2str(omega),... 'HorizontalAlignment','left'); set(h70,'Value',omega); resetAll() makePlot(); end %--------------------------------------------- function runCallback(hObject,eventdata) stateOfSimulation='RUNNING'; if currentStep == size(data,1) resetData(); end while (currentStep < size(data,1)&& ... strcmp(stateOfSimulation,'RUNNING')) makeOneStep(); makePlot(); end makePlot(); end %--------------------------------------------- function makePlot() if currentStep>0 set(f,'CurrentAxes',h1); hold(h1,'on'); grid(h1,'on'); xlim(h1,[0 maxTime]); plot(h1,data(1:currentStep,1),data(1:currentStep,2),'r-'); title(h1,{'y(t) at time',sprintf('%3.2f sec',... data(currentStep,1))},'FontSize',9); set(f,'CurrentAxes',h2); hold(h2,'off'); quiver(0,0,cos(omegaRadians*data(currentStep,1)),... sin(omegaRadians*data(currentStep,1)),... 'color',[1,0,0],'linewidth',2); hold(h2,'on'); quiver(0,0,cos(omegaRadians*data(currentStep,1)-phase),... sin(omegaRadians*data(currentStep,1)-phase),... 'color',[0,0,0],'linewidth',2); set(gca, 'XLim', [-1.1 1.1], 'YLim', [-1.1 1.1],'FontSize',6); set(gca,'PlotBoxAspectRatio',[1 1 1]) ; title(h2,{'load/displacement phase diagram, Force is red arrow'},... 'FontSize',8); line([0 0],[-1 1]); line([-1 1],[0 0]); grid(h2,'on'); set(f,'CurrentAxes',hMaxForceDispAxes); plot(hMaxForceDispAxes,... data(forceVectorLastChangedIndex:currentStep,2)/... (1/k*magnificationFactor),... data(forceVectorLastChangedIndex:currentStep,3),'r-'); title(hMaxForceDispAxes,'f(t) vs. normalize response','FontSize',8); set(hMaxForceDispAxes,'PlotBoxAspectRatio',[1 1 1],... 'XLim', [-1.1 1.1], 'YLim', [-1.1 1.1],'FontSize',6) ; line([0 0],[-1 1]); line([-1 1],[0 0]); grid(hMaxForceDispAxes,'on'); end drawnow; end %--------------------------------------------- function makeOneStep() currentStep = currentStep+1; if currentStep>1 data(currentStep,1) = data(currentStep-1,1)+timeStep; end data(currentStep,2) = 1/k * magnificationFactor * ... cos(omegaRadians*data(currentStep,1)-phase); data(currentStep,3) = cos(omegaRadians*data(currentStep,1)); end %--------------------------------------------- function resetData() data = zeros(length(0:timeStep:maxTime),3); currentStep = 0; cla(h1,'reset') end %--------------------------------------------- function resetCallback(hObject,eventdata) resetData(); stateOfSimulation = 'RESET'; cla(h1,'reset') end %--------------------------------------------- function pauseCallback(hObject,eventdata) stateOfSimulation = 'PAUSED'; end %--------------------------------------------- function stepCallback(hObject,eventdata) stateOfSimulation = 'STEP'; if currentStep1) phase = atan2(2*zeta*r,1-r^2); else phase = pi/2; end end end end