home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / basic / library / trubasic / chaos / xpoincar.tru < prev    next >
Encoding:
Text File  |  1980-01-01  |  3.3 KB  |  113 lines
  1. !PROGRAM TITLE "XPOINCAR(E)"
  2. CLEAR
  3. PRINT"               ***PENDULUM - POINCARE SECTION***"
  4. PRINT
  5. PRINT"THIS PROGRAM DISPLAYS THE POINCARE SECTION OF THE PENDULUM"
  6. PRINT"AND CAN SAVE THE DATA TO A FILE."
  7. LIBRARY "SGLIB.TRC"
  8. !
  9. DECLARE DEF accel
  10. DIM A(1),B(1)
  11. INPUT prompt"Input driving force strength: ":g
  12. INPUT prompt"Input damping (If no damping then input 9999999):":q
  13. INPUT prompt"Input initial angle, angular velocity: ": xint,vint
  14. INPUT Prompt"Input min. and max. time:":tmin,tmax
  15. INPUT prompt"Input phase angle/(2*pi): ":phi
  16. INPUT PROMPT" SAVE DATA TO A FILE? YES(1), NO(2):":SAVEFILE
  17. IF SAVEFILE=1 THEN
  18.  INPUT PROMPT"FILE NAME FORMAT EX. 14954020 :":FILENAME
  19.  INPUT PROMPT"DRIVE FOR FILE DISK A,B,C,ETC.:":DISK$
  20.  LET NAME$=DISK$&":"&STR$(FILENAME)
  21. END IF
  22. !
  23. CALL PARAMS(W,EPS,TSTEP,XMIN,XMAX,YMIN,YMAX)
  24. CALL SETXSCALE(XMIN,XMAX)
  25. CALL SETYSCALE(YMIN,YMAX)
  26. CALL SETTEXT("PENDULUM POINCARE SECTION","THETA","OMEGA")
  27. CALL RESERVELEGEND
  28. !
  29. DATA O,O
  30. CALL DATAGRAPH(A,B,1,O,"RED")
  31. LET t=0
  32. LET x=xint
  33. LET v=vint
  34. CALL GOTOCANVAS
  35. !
  36. !CALCULATION AND GRAPHING BLOCK
  37. LET phi=phi*2*pi
  38. IF SAVEFILE=1 THEN
  39. OPEN #1:NAME NAME$, ORGANIZATION RECORD, CREATE NEWOLD
  40. ASK #1:FILESIZE LENGTH
  41. IF LENGTH=0 THEN SET#1:RECSIZE 10
  42. SET #1: POINTER END
  43. END IF
  44. FOR i=1 to 1000000
  45.     CALL rk4(x,v,tstep,xnew,vnew,t,w,g,q)     ! Take a step of size tstep
  46.     LET tshalf=tstep/2
  47.     CALL rk4(x,v,tshalf,xnh,vnh,t,w,g,q)      !Take two half steps
  48.  
  49.     CALL rk4(xnh,vnh,tshalf,xn,vn,t+tshalf,w,g,q)
  50.     LET d1=abs(xn-xnew)
  51.     LET d2=abs(vn-vnew)
  52.     LET delta=max(d1,d2)
  53.     IF delta<eps then
  54.        IF t>tmin then
  55.           LET tnew=t+tstep
  56.           LET w1=mod(phi-w*t,2*pi)     !Check for Poincare section
  57.           LET w2=mod(w*tnew-phi,2*pi)
  58.           IF w1<w*tstep then
  59.              IF w2<w*tstep then
  60.                 LET ts=w1/w
  61.                 CALL rk4(x,v,ts,xp,vp,t,w,g,q)  !CALCULATES POINT AT SECTION
  62.                 IF abs(xp)>pi then LET xp=xp-2*pi*abs(xp)/xp
  63.                 CALL GRAPHPOINT(XP,VP,1)
  64.                 IF SAVEFILE=1 THEN WRITE #1:XP,VP
  65.              END IF
  66.           END IF
  67.        END IF
  68.        LET x=xnew
  69.        LET v=vnew
  70.        LET t=t+tstep              !Expand step size
  71.        LET tstep=tstep*.95*(eps/delta)^.25
  72.        IF abs(x)>pi then          !bring theta back into range
  73.           LET x=x-2*pi*abs(x)/x
  74.        END IF
  75.     ELSE                          !else reduce step size
  76.        LET tstep=tstep*.95*(eps/delta)^.2
  77.     END IF
  78.     IF t>tmax then LET i=1000001
  79. NEXT i
  80. LET G$=STR$(G)
  81. LET Q$=STR$(Q)
  82. CALL ADDLEGEND("G="&STR$(G)&"   Q="&STR$(Q)&"   PHI="&STR$(PHI),0,1,"WHITE")
  83. CALL DRAWLEGEND
  84. END
  85. !
  86. !
  87. SUB rk4(x,v,tstep,xnew,vnew,t,w,g,q)
  88.     DECLARE DEF accel
  89.     LET xk1=tstep*v
  90.     LET vk1=tstep*accel(x,v,t,w,g,q)
  91.     LET xk2=tstep*(v+vk1/2)
  92.     LET vk2=tstep*accel(x+xk1/2,v+vk1/2,t+tstep/2,w,g,q)
  93.     LET xk3=tstep*(v+vk2/2)
  94.     LET vk3=tstep*accel(x+xk2/2,v+vk2/2,t+tstep/2,w,g,q)
  95.     LET xk4=tstep*(v+vk3)
  96.     LET vk4=tstep*accel(x+xk3,v+vk3,t+tstep,w,g,q)
  97.     LET vnew=v+(vk1+2*vk2+2*vk3+vk4)/6
  98.     LET xnew=x+(xk1+2*xk2+2*xk3+xk4)/6
  99. END SUB
  100. DEF accel(x,v,t,w,g,q)
  101.     LET accel=-sin(x)-(1/q)*v+g*cos(w*t)
  102. END def
  103. !
  104. SUB PARAMS(W,EPS,TSTEP,XMIN,XMAX,YMIN,YMAX)
  105.     LET W=0.66666666
  106.     LET EPS=1.0E-6
  107.     LET TSTEP=0.5
  108.     LET XMIN=-3
  109.     LET XMAX=3
  110.     LET YMIN=-3
  111.     LET YMAX=3
  112. END SUB
  113.