home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / basic / library / trubasic / chaos / xinitial.tru < prev    next >
Encoding:
Text File  |  1980-01-01  |  3.3 KB  |  110 lines
  1. !PROGRAM TITLE "XINITIAL( EVOLUTION)"
  2. !THIS PROGRAM DISPLAYS THE EVOLUTION OF A BLOCK OF PHASE POINTS FOR THE PENDULUM
  3. LIBRARY "SGLIB.TRC"
  4. !
  5. DECLARE DEF accel
  6. DIM A(1),B(1)
  7. INPUT prompt"Input driving force strength: ":g
  8. INPUT prompt"Input damping (If no damping then input 9999999):":q
  9. PRINT"THE NEXT SET OF PROMPTS ASK FOR THE BOUNDARIES OF THE AREA  OF PHASE"
  10. PRINT"POINTS WHOSE INITIAL EVOLUTION WILL BE GRAPHED."
  11. INPUT PROMPT"LOWEST X VALUE:":LOWX
  12. INPUT PROMPT"HIGHEST X VALUE:":HIGHX
  13. INPUT PROMPT"LOWEST Y VALUE:":LOWY
  14. INPUT PROMPT"HIGHEST Y VALUE:":HIGHY  
  15. INPUT Prompt"Input min. and max. time:":tmin,tmax
  16. INPUT prompt"Input TIME STEP AS A FRACTION OF THE FORCING PERIOD: ":phi
  17. LET PHISTEP=PHI*2*PI
  18. CALL PARAMS(W,EPS,TSTEP,XMIN,XMAX,YMIN,YMAX)
  19. CALL SETXSCALE(XMIN,XMAX)
  20. CALL SETYSCALE(YMIN,YMAX)
  21. CALL SETTEXT("EVOLUTION OF INITIAL PENDULUM STATES","THETA","OMEGA")
  22. CALL RESERVELEGEND
  23. !
  24. DATA O,O
  25. CALL DATAGRAPH(A,B,1,O,"RED")
  26. CALL GOTOCANVAS
  27. !
  28. !LOOP FOR BLOCK OF POINT
  29. FOR XINT = LOWX TO HIGHX STEP .1
  30.  FOR VINT = LOWY TO HIGHY STEP .1
  31.  CALL GRAPHPOINT(XINT,VINT,1)
  32.   LET T = 0
  33. LET X=XINT
  34. LET V=VINT
  35.  
  36.  
  37.  
  38. !CALCULATION AND GRAPHING BLOCK
  39. LET phi=phi*2*pi
  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(-W*T,PHISTEP)     !Check for TIME STEP TO PLOT POINT
  53.           LET w2=mod(W*TNEW,PHISTEP)
  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 TIME STEP
  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 STEP="&STR$(PHISTEP*1.5),0,1,"WHITE")
  80. CALL DRAWLEGEND
  81. END
  82. !
  83. !
  84. SUB rk4(x,v,tstep,xnew,vnew,t,w,g,q)
  85.     DECLARE DEF accel
  86.     LET xk1=tstep*v
  87.     LET vk1=tstep*accel(x,v,t,w,g,q)
  88.     LET xk2=tstep*(v+vk1/2)
  89.     LET vk2=tstep*accel(x+xk1/2,v+vk1/2,t+tstep/2,w,g,q)
  90.     LET xk3=tstep*(v+vk2/2)
  91.     LET vk3=tstep*accel(x+xk2/2,v+vk2/2,t+tstep/2,w,g,q)
  92.     LET xk4=tstep*(v+vk3)
  93.     LET vk4=tstep*accel(x+xk3,v+vk3,t+tstep,w,g,q)
  94.     LET vnew=v+(vk1+2*vk2+2*vk3+vk4)/6
  95.     LET xnew=x+(xk1+2*xk2+2*xk3+xk4)/6
  96. END SUB
  97. DEF accel(x,v,t,w,g,q)
  98.     LET accel=-sin(x)-(1/q)*v+g*cos(w*t)
  99. END def
  100. !
  101. SUB PARAMS(W,EPS,TSTEP,XMIN,XMAX,YMIN,YMAX)
  102.     LET W=0.66666666
  103.     LET EPS=1.0E-6
  104.     LET TSTEP=0.5
  105.     LET XMIN=-PI
  106.     LET XMAX=PI
  107.     LET YMIN=-4
  108.     LET YMAX=4
  109. END SUB
  110.