home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / basic / library / trubasic / chaos / xphase2d.tru < prev    next >
Encoding:
Text File  |  1980-01-01  |  2.7 KB  |  90 lines
  1. !PROGRAM TITLE - "XPHASE2D"
  2. !THIS PROGRAM DISPLAYS THE 2-DIMENSIONAL PHASE DIAGRAM 
  3. !FOR THE DRIVEN AND UNDRIVEN PENDULUM.
  4. !
  5. LIBRARY "SGLIB.TRC"
  6. DECLARE DEF ACCEL
  7. DIM A(1), B(1)
  8. !INPUT STATEMENTS
  9.  INPUT PROMPT"INPUT DRIVING FORCE STRENGTH:":g
  10.  INPUT PROMPT"INPUT DAMPING (IF NO DAMPING THEN INPUT 9999999):":q
  11.  INPUT PROMPT"INPUT: INITIAL ANGLE , ANGULAR VELOCITY:":XINT,VINT
  12.  INPUT PROMPT"INPUT: MINUMUM TIME , MAXIMUM TIME:":TMIN,TMAX
  13.  INPUT PROMPT"AVERAGE VELOCTY CALCULATION; YES(1) , NO(2):":AVC
  14. CALL PARAMS(W,EPS,TSTEP,XMIN,XMAX,YMIN,YMAX)   !SETS MISC AND GRAPH PARAMETERS 
  15. CALL SETXSCALE(XMIN,XMAX)  !FROM SGLIB
  16. CALL SETYSCALE(YMIN,YMAX)  !FROM SGLIB
  17. CALL SETTEXT("PENDULUM - 2-D PHASE DIAGRAM","ANGLE","ANGULAR VELOCITY")
  18. CALL RESERVELEGEND  !FROM SGLIB , SAVES SPACE FOR LEGENDS
  19.  
  20. DATA 0,0
  21. CALL DATAGRAPH(A,B,1,0,"WHITE")  !FROM SGLIB - PLOTS INITIAL POINT
  22. LET T=0
  23. LET X=XINT
  24. LET V=VINT
  25. CALL GOTOCANVAS   !SETS SCREEN FOR GRAPH
  26. !
  27. !CALCULATION AND GRAPHNG BLOCK
  28. FOR I=1 TO 10000000
  29.   CALL RK4(X,V,TSTEP,XNEW,VNEW,T,W,G,Q)   !CALL RUNGE-KUTTA, STEP = TSTEP
  30.    LET TSHALF=TSTEP/2    ! SPLIT INTERVAL
  31.   CALL RK4(X,V,TSHALF,XNH,VNH,T,W,G,Q)  ! DO TWO HALF STEPS
  32.   CALL RK4(XNH,VNH,TSHALF,XN,VN,T+TSHALF,W,G,Q)
  33.   LET D1=ABS(XN-XNEW)
  34.   LET D2=ABS(VN-VNEW)
  35.   LET DELTA=MAX(D1,D2)
  36.   IF DELTA<EPS THEN
  37.     IF T>TMIN THEN
  38.      IF ABS(X)>PI THEN LET X=X-2*PI*ABS(X)/X
  39.      CALL GRAPHPOINT(X,V,1)
  40.      LET SUMVEL=SUMVEL+V*TSTEP   !UPDATE AVERAGE
  41.     END IF
  42.      LET X=XNEW
  43.      LET V=VNEW
  44.      LET T=T+TSTEP
  45.      LET TSTEP=TSTEP*.95*(EPS/DELTA)^.25
  46.      IF ABS(X)>PI THEN LET X=X-2*PI*ABS(X)/X
  47.   ELSE
  48.      LET TSTEP=TSTEP*.95*(EPS/DELTA)^.2   !REDUCE STEP SIZE
  49.   END IF
  50.   IF T>TMAX THEN LET I=10000001
  51.   NEXT I
  52.   LET MEANVEL=SUMVEL/(TMAX-TMIN)
  53.   CALL ADDLEGEND("G="&STR$(G)&"   Q="&STR$(Q),0,1,"WHITE")  
  54.   IF AVC=1 THEN CALL ADDLEGEND("AV. VEL. = "&STR$(MEANVEL),0,1,"WHITE")
  55.   CALL DRAWLEGEND    !ADDS G AND Q VALUES TO LEGEND
  56. get key variable
  57. clear
  58. print"press <esc> key to finish"
  59. END 
  60. !
  61. SUB RK4(X,V,TSTEP,XNEW,VNEW,T,W,G,Q)   !RUNGE-KUTTA INTEGRATOR
  62.    DECLARE DEF ACCEL
  63.    LET XK1=TSTEP*V
  64.    LET VK1=TSTEP*ACCEL(X,V,T,W,G,Q)
  65.    LET XK2=TSTEP*(V+VK1/2)
  66.    LET VK2=TSTEP*ACCEL(X+XK1/2,V+VK1/2,T+TSTEP/2,W,G,Q)
  67.    LET XK3=TSTEP*(V+VK2/2)
  68.    LET VK3=TSTEP*ACCEL(X+XK2/2,V+VK2/2,T+TSTEP/2,W,G,Q)
  69.    LET XK4=TSTEP*(V+VK3)
  70.    LET VK4=TSTEP*ACCEL(X+XK3,V+VK3,T+TSTEP,W,G,Q)
  71.    LET VNEW=V+(VK1+2*VK2+2*VK3+VK4)/6
  72.    LET XNEW=X+(XK1+2*XK2+2*XK3+XK4)/6
  73. END SUB
  74. !
  75. DEF ACCEL(X,V,T,W,G,Q)
  76.    LET DAMP=1/Q
  77.    LET ACCEL=-SIN(X)-DAMP*V+G*COS(W*T)
  78. END DEF
  79. !
  80. SUB PARAMS(W,EPS,TSTEP,XMIN,XMAX,YMIN,YMAX)
  81.  LET W=0.66666666
  82.  LET EPS=1.0E-6
  83.  LET TSTEP=0.5
  84.  LET XMIN=-3
  85.  LET XMAX=3
  86.  LET YMIN=-3
  87.  LET YMAX=3
  88. END SUB
  89.  
  90.