rem INVERTED PENDULUM Simulation .... Rev 2.0 rem A J Tooth / 22nd Dec 2002 rem Modified 31st December 2003 : Phase Plot Option added on error if (err=17) then quit proc_fullscreen(22) : rem Set up use of Full Screen proc_pi : rem Defines a symbol for "pi" proc_theta : rem Defines a symbol for theta proc_thetadot : rem Defines a symbol for theta-dot proc_sqrd : rem Defines a superfix symbol for "squared" : 100 rem Choose colour scheme dependant on whether printing required cls : print "Press -i- to use a WHITE background" inv$=get$ if inv$="i" or inv$="I" then colour 128,15 : colour 7,0 : colour 3,1 : colour 6,0 else colour 128,0 : colour 7,7 : colour 3,3 : colour 6,6 endif clg rem Format of Main Differential Function phi$="((V*W*W*COS(W*T))-1)*(SIN(P))*D" rem Enter Parametric Options print " Suggested trial values appear in brackets." input " How many steps for the Time parameter T (5000)",stps% if stps%=0 then stps%=5000 input " What value for initial angle (radians) from vertical (1.3)",P0S if P0S=0 then P0S=1.3 input " What initial value for angular velocity of pendulum, if any (0)",Q0S input " What value for number of oscillations per second of the pivot (25)",w if w=0 then w=25 input " What value for max amplitude of pivot osciallations. Range 0 to 1. (0.1)",v if v=0 then v=0.1 input " How many cycles (30) ",lim% if lim%=0 then lim%=30 print " See actual motion or phase-space representation?" print " Press -p- for phase space, or any other key for actual motion." p$=get$ scal=600 : md=pi clg : cls xmax%=1024 : ymax%=768 rem Place Origin and draw Axes setx%=xmax% : sety%=ymax% origin setx%,sety% : gcol 3 move -xmax%,0 : draw xmax%,0 move 0,-ymax% : draw 0,ymax% rem Displays Expression for Interest / Reference xdis$="(vw"+chr$(131)+"cos(wt)-1).sin"+chr$(129)+".dt" print tab(5,4);"The INVERTED PENDULUM" print tab(5,2);" d";chr$(130);"=";xdis$ rem Curve Drawing Section rem Sets initial and incremental values dtim=1/stps% P0=pi-P0S : Q0=Q0S rem Plots Initial Value if p$="p" then colour 2: print tab(5,7);"PHASE-SPACE PLOT" colour 7: print tab(65,3);chr$(130);" Axis" print tab(115,24);chr$(129);" Axis" gcol 5: move 150*md,700 : draw 150*md,-700 move (-150*md),700 : draw (-150*md),-700 print tab(32,24);"-";chr$(128);tab(94,24);chr$(128) else proc_plot(P0,scal,v) print tab(5,7);"ACTUAL MOTION" endif gcol 6 for par%= 0 to lim% for tim%=0 to stps% timlit=par%+(dtim*tim%) : timlit2=timlit + dtim rem Computational Error Trapping on error local goto 1210 rem Main calculations - Improved Euler Method dphi=fn_dfunc(phi$,timlit,Q0,P0,dtim,w,v) : Q1=Q0 + dphi dtheta=Q0*dtim : P1=P0 + dtheta dphi2=fn_dfunc(phi$,timlit2,Q1,P1,dtim,w,v) dtheta2=Q1*dtim Q1=Q0 + ((dphi+dphi2)/2) : P1=P0 + ((dtheta+dtheta2)/2) rem Plot formatting section vact = v*cos(w*timlit) if p$="p" then proc_phas(P1,Q1,scal) else proc_plot(P1,scal,vact) endif rem Resets current values for Phi and Theta at time T Q0=Q1 : P0=P1 960 next tim% next par% print tab(5,10);"DONE" rem Repeat Options rep$=get$ clg print tab(10,10);" Re-run, or QUIT the programme? ( r or q )" rep$=get$ : clg if rep$="r" or rep$="R" then flag%=1 else flag%=2 on flag% goto 100,1080 1080 quit : rem END OF PROGRAMME end rem ===================================================================== rem Various Functions and Procedures Section rem Main Calculation for x and y def fn_dfunc(main$,T,Q,P,D,W,V) =eval(main$) rem Traps Computational Range Errors ONLY 1210 if err=21 or err=23 or err=18 or err=20 or err=22 or err=24 then goto 960 else restore error endif goto 960 rem Plot Formatting Procedure def proc_plot(Pval,Scal,V) local XP,YP,VP,sx%,sy%,x%,y%,delay%,xl%,yl% VP=Scal*V XP=Scal*sin(Pval) : YP=VP-Scal*cos(Pval) sx%=sgn(XP) : sy%=sgn(YP) x%=sx%*int(abs(XP)) y%=sy%*int(abs(YP)) if abs(y%)>1600 then goto 1430 if abs(x%)>2100 then goto 1430 rem Prints pendulum shaft, and locus of end-point xl%= int(0.6*x%) : yl%= int(0.6*y%) move 0,VP : gcol 1 : plot 5,xl%,yl% : gcol 6 : plot 69,x%,y% sys "UpdateWindow",@hwnd% rem Un-prints pendulum shaft, to simulate movement move 0,VP : gcol 0 : plot 5,xl%,yl% 1430 endproc rem Set up use of Full Screen def proc_fullscreen(N%) sys "GetSystemMetrics", 0 to xscreen% sys "GetSystemMetrics", 1 to yscreen% sys "SetWindowLong",@hwnd%,-16,&16000000 sys "SetWindowPos",@hwnd%,-1,0,0,xscreen%,yscreen%,0 mode N% mouse off : off : rem Turns off the Mouse Pointer and the Cursor endproc rem Defines ASCII character 128 as a symbol for Pi def proc_pi local r1%,r2%,r3%,r4%,r5%,r6%,r7%,r8% r1%= %00000000 r2%= %00000000 r3%= %11111111 r4%= %01101100 r5%= %01101100 r6%= %01100110 r7%= %00000000 r8%= %00000000 vdu 23,128,r1%,r2%,r3%,r4%,r5%,r6%,r7%,r8% endproc rem Defines ASCII character 129 as a symbol for Theta def proc_theta local r1%,r2%,r3%,r4%,r5%,r6%,r7%,r8% r1%= %00000000 r2%= %00111100 r3%= %01100110 r4%= %01111110 r5%= %01100110 r6%= %00111100 r7%= %00000000 r8%= %00000000 vdu 23,129,r1%,r2%,r3%,r4%,r5%,r6%,r7%,r8% endproc rem Defines ASCII character 130 as a symbol for ThetaDot def proc_thetadot local r1%,r2%,r3%,r4%,r5%,r6%,r7%,r8% r1%= %00011000 r2%= %00000000 r3%= %00111100 r4%= %01100110 r5%= %01111110 r6%= %01100110 r7%= %00111100 r8%= %00000000 vdu 23,130,r1%,r2%,r3%,r4%,r5%,r6%,r7%,r8% endproc rem Defines ASCII character 131 as a superfix symbol for "squared" def proc_sqrd local r1%,r2%,r3%,r4%,r5%,r6%,r7%,r8% r1%= %00110000 r2%= %01011000 r3%= %00110000 r4%= %01111100 r5%= %00000000 r6%= %00000000 r7%= %00000000 r8%= %00000000 vdu 23,131,r1%,r2%,r3%,r4%,r5%,r6%,r7%,r8% endproc rem Displays Phase-Space representation def proc_phas(P,Q,Scal) local Pref,dv% P=md-P : Pref=P rem Count how many times round dv%=int((1+((abs(Pref))/md))/2) if Pref<-md then P=P+(dv%*2*md) if Pref>md then P=P-(dv%*2*md) plot 69,150*P,100*Q endproc