home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / basic / library / trubasic / chaos / xinitflo.tru < prev    next >
Encoding:
Text File  |  1980-01-01  |  3.8 KB  |  112 lines
  1.  
  2. !PROGRAM TITLE "XINITFLO"
  3. CLEAR
  4. PRINT"             ***PENDULUM - INITIAL FLOW OF BLOCK OF COORDS***"
  5. PRINT"THIS PROGRAM DISPLAYS THE EVOLUTION OF A BLOCK OF PHASE POINTS PENDULUM"
  6. PRINT"AT VARIOUS TIMES DURING THE DRIVE CYCLE.  FOR EXAMPLE, IF PHI=0.5 THEN"
  7. PRINT"THE BLOCK WILL BE PRINTED EVERY 1/2 DRIVE PERIOD UNTIL THE MAXIMUM TIME."
  8. PRINT"AS A FIRST TRIAL LET THE MINIMUM TIME BE ZERO AND THE MAXIMUM BE 5,"
  9. PRINT" WHICH IS ABOUT 1/2 A DRIVE PERIOD."
  10. LIBRARY "SGLIB.TRC"
  11. !
  12. DECLARE DEF accel
  13. DIM A(1),B(1)
  14. INPUT prompt"Input driving force strength: ":g
  15. INPUT prompt"Input damping (If no damping then input 9999999):":q
  16. !INPUT prompt"Input initial angle, angular velocity: ": xint,vint
  17.  
  18. INPUT prompt"input lower and upper xint:":lowerxint, upperxint
  19. INPUT prompt"input lower and upper vint:":lowervint, uppervint
  20. INPUT Prompt"Input min. and max. time:":tmin,tmax
  21. INPUT prompt"Input phase angle/(2*pi): ":phi
  22. CALL PARAMS(W,EPS,TSTEP,XMIN,XMAX,YMIN,YMAX)
  23. CALL SETXSCALE(XMIN,XMAX)
  24. CALL SETYSCALE(YMIN,YMAX)
  25. CALL SETTEXT("PENDULUM INITIAL FLOW","ANGLE","ANGULAR VELOCITY")
  26. CALL RESERVELEGEND
  27. !
  28. DATA O,O
  29. CALL DATAGRAPH(A,B,1,O,"WHITE")
  30. CALL gotocanvas
  31. LET phi=phi*2*pi
  32. FOR xint=lowerxint to upperxint step 0.03
  33.     FOR vint=lowervint to uppervint step 0.03
  34.         LET x=xint
  35.         LET v=vint
  36.         LET t=0
  37.         CALL graphpoint(xint,vint,1)
  38.         !
  39.         !CALCULATION AND GRAPHING BLOCK
  40.         FOR i=1 to 1000000
  41.             CALL rk4(x,v,tstep,xnew,vnew,t,w,g,q)     ! Take a step of size tstep
  42.             LET tshalf=tstep/2
  43.             CALL rk4(x,v,tshalf,xnh,vnh,t,w,g,q)      !Take two half steps
  44.  
  45.             CALL rk4(xnh,vnh,tshalf,xn,vn,t+tshalf,w,g,q)
  46.             LET d1=abs(xn-xnew)
  47.             LET d2=abs(vn-vnew)
  48.             LET delta=max(d1,d2)
  49.             IF delta<eps then
  50.                IF t>tmin then
  51.                   LET tnew=t+tstep
  52.                   LET w1=mod(phi-w*t,2*pi)  !Check for Poincare section
  53.                   LET w2=mod(w*tnew-phi,2*pi)
  54.                   IF w1<w*tstep then
  55.                      IF w2<w*tstep then
  56.                         LET ts=w1/w
  57.                         CALL rk4(x,v,ts,xp,vp,t,w,g,q)     !CALCULATES POINT AT SECTION
  58.                         !                        IF abs(xp)>pi then LET xp=xp-2*pi*abs(xp)/xp
  59.                         CALL GRAPHPOINT(XP,VP,1)
  60.                      END IF
  61.                   END IF
  62.                END IF
  63.                LET x=xnew
  64.                LET v=vnew
  65.                LET t=t+tstep      !Expand step size
  66.                LET tstep=tstep*.95*(eps/delta)^.25
  67.                !              IF abs(x)>pi then  !bring theta back into range
  68.                !                LET x=x-2*pi*abs(x)/x
  69.                !            END IF
  70.             ELSE                  !else reduce step size
  71.                LET tstep=tstep*.95*(eps/delta)^.2
  72.             END IF
  73.             IF t>tmax then LET i=1000001
  74.         NEXT i
  75.     NEXT vint
  76. NEXT xint
  77. LET G$=STR$(G)
  78. LET Q$=STR$(Q)
  79. CALL ADDLEGEND("G="&STR$(G)&"   Q="&STR$(Q)&"   TIME="&STR$(PHI*3/2),0,1,"WHITE")
  80. CALL DRAWLEGEND
  81. get key variable
  82. clear
  83. END 
  84. !
  85. !
  86. SUB rk4(x,v,tstep,xnew,vnew,t,w,g,q)
  87.     DECLARE DEF accel
  88.     LET xk1=tstep*v
  89.     LET vk1=tstep*accel(x,v,t,w,g,q)
  90.     LET xk2=tstep*(v+vk1/2)
  91.     LET vk2=tstep*accel(x+xk1/2,v+vk1/2,t+tstep/2,w,g,q)
  92.     LET xk3=tstep*(v+vk2/2)
  93.     LET vk3=tstep*accel(x+xk2/2,v+vk2/2,t+tstep/2,w,g,q)
  94.     LET xk4=tstep*(v+vk3)
  95.     LET vk4=tstep*accel(x+xk3,v+vk3,t+tstep,w,g,q)
  96.     LET vnew=v+(vk1+2*vk2+2*vk3+vk4)/6
  97.     LET xnew=x+(xk1+2*xk2+2*xk3+xk4)/6
  98. END SUB
  99. DEF accel(x,v,t,w,g,q)
  100.     LET accel=-sin(x)-(1/q)*v+g*cos(w*t)
  101. END DEF
  102. !
  103. SUB PARAMS(W,EPS,TSTEP,XMIN,XMAX,YMIN,YMAX)
  104.     LET W=0.66666666
  105.     LET EPS=1.0E-6
  106.     LET TSTEP=0.5
  107.     LET XMIN=-6
  108.     LET XMAX=6
  109.     LET YMIN=-4
  110.     LET YMAX=4
  111. END SUB
  112.