home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / basic / library / trubasic / chaos / xbasinpo.tru < prev    next >
Encoding:
Text File  |  1980-01-01  |  4.8 KB  |  145 lines
  1. !PROGRAM XBASINPO(INCARE)
  2. CLEAR
  3. PRINT"                      ***PENDULUM - BASINS OF ATTRACTION***"
  4. PRINT"This program calculates the average angular velocity for"
  5. PRINT"an array of initial conditions:"
  6. PRINT"velocity (-3,3) and angle (-pi ,pi).  If the average angular velocity"
  7. PRINT"from a given initial point is positive then circle is placed at the"
  8. PRINT"point, otherwise the average is negative.  This program can also "
  9. PRINT"superpose the appropriate Poincare section on the graph."
  10. PRINT"The program can also save the data sets to two different files."
  11. PRINT
  12. LIBRARY "sglib.trc"
  13. DIM a(1),b(1)
  14.  
  15. DECLARE DEF accel
  16.  
  17. INPUT prompt"Input driving force strength: ":g
  18. INPUT prompt"input damping :":q
  19. INPUT Prompt"Input min. and max. time of averaging:":tmin,tmax
  20. INPUT Prompt"Poincare attractor yes (1), no(2) ":Poincare
  21. INPUT prompt"Save data to a file, yes(1), no(2):":sv
  22. IF sv=1 then
  23.    PRINT"Basin of attraction File name includes:"
  24.    PRINT"  1)First 2 digits for q value"
  25.    PRINT"  2)Next needed digits for g value"
  26.    PRINT"  3)Last digits as 0's"
  27.    INPUT prompt"file name  - format ex. 20150000:":filename
  28.    INPUT prompt"data file drive a,b,c,etc.?":b$
  29.    IF poincare=1 then
  30.       PRINT"Superposed Poincare section file name is similar to above except"
  31.       PRINT"that the numeral '0' is added, for example, 20015000."
  32.       INPUT prompt"Input corresponding Poincare file name:":poinfile
  33.    END IF
  34.    LET FILENAME$=STR$(FILENAME)
  35.    LET POINFILE$=STR$(POINFILE)
  36. END IF
  37. !
  38. CALL PARAMS(W,EPS,TSTEP)
  39. CALL SETXSCALE(-3,3)
  40. CALL SETYSCALE(-3,3)
  41. CALL SETTEXT("PENDULUM - BASINS OF ATTRACTION","INIT. ANGLE","INIT. ANG. VEL.")
  42. CALL RESERVELEGEND
  43. IF sv=1 then
  44.    OPEN #1:name b$&":"&filename$,organization record,create newold
  45.    ASK #1: FILESIZE length
  46.    IF length=0 then SET#1:rECSIZE 10
  47.    SET #1: POINTER end
  48.    IF poincare=1 then
  49.       OPEN #2:name b$&":"&poinfile$,organization record,create newold
  50.       ASK #2: FILESIZE length
  51.       IF length=0 then SET#2:rECSIZE 10
  52.       SET #2:pOINTER end
  53.    END IF
  54. END IF
  55.  
  56. DATA 0,0
  57. CALL DATAGRAPH(A,B,1,0,"WHITE")
  58. CALL gotocanvas
  59. LET phi=0
  60. FOR xint=-pi to pi step .1
  61.     FOR vint=-3 to 3 step .15
  62.         LET x= xint
  63.         LET v=vint
  64.         LET t=0
  65.         LET s=0
  66.         FOR i=1 to 1000000
  67.             CALL rk4(x,v,tstep,xnew,vnew,t,w,g,q)     ! Take a step of size tstep
  68.             LET tshalf=tstep/2
  69.             CALL rk4(x,v,tshalf,xnh,vnh,t,w,g,q)      !Take two half steps
  70.             CALL rk4(xnh,vnh,tshalf,xn,vn,t+tshalf,w,g,q)
  71.             LET d1=abs(xn-xnew)
  72.             LET d2=abs(vn-vnew)
  73.             LET delta=max(d1,d2)
  74.             IF delta<eps then
  75.                IF t>tmin then
  76.                   LET tnew=t+tstep
  77.                   LET s=s+vnew*tstep
  78.                   IF Poincare=1 then
  79.                      LET w1=mod(phi-w*t,2*pi)
  80.                      LET w2=mod(w*tnew-phi,2*pi)
  81.                      IF w1<w*tstep then
  82.                         IF w2<w*tstep then
  83.                            LET ts=w1/w
  84.                            CALL rk4(x,v,ts,xp,vp,t,w,g,q)
  85.                            IF abs(xp)>pi then LET xp=xp-2*pi*abs(xp)/xp
  86.                            CALL GRAPHPOINT(xp,vp,1)
  87.                            IF sv=1 then WRITE #2:xp,vp
  88.                         END IF
  89.                      END IF
  90.                   END IF
  91.                END IF
  92.                LET x=xnew
  93.                LET v=vnew
  94.                LET t=t+tstep      !Expand step size
  95.                LET tstep=tstep*.95*(eps/delta)^.25
  96.                IF abs(x)>pi then  !bring theta back in range
  97.                   LET x=x-2*pi*abs(x)/x
  98.                END IF
  99.             ELSE                  !else reduce step size
  100.                LET tstep=tstep*.95*(eps/delta)^.2
  101.             END IF
  102.             IF t>tmax then
  103.                LET average=s/(tmax-tmin)
  104.                IF average>0 then
  105.                   CALL GRAPHPOINT(XINT,VINT,4)
  106.                   IF sv=1 then WRITE #1:xint,vint
  107.                END IF
  108.                LET i=1000001
  109.             END IF
  110.         NEXT i
  111.     NEXT vint
  112. NEXT xint
  113. CALL addlegend("g="&str$(g)&"  q="&str$(q)&"  circle=positive",0,1,"white")
  114. CALL drawlegend
  115. GET KEY VARIABLE
  116. CLEAR
  117. print"press <esc> key to finish"
  118. END
  119. !
  120. !
  121. SUB rk4(x,v,tstep,xnew,vnew,t,w,g,q)
  122.     DECLARE DEF accel
  123.     LET xk1=tstep*v
  124.     LET vk1=tstep*accel(x,v,t,w,g,q)
  125.     LET xk2=tstep*(v+vk1/2)
  126.     LET vk2=tstep*accel(x+xk1/2,v+vk1/2,t+tstep/2,w,g,q)
  127.     LET xk3=tstep*(v+vk2/2)
  128.     LET vk3=tstep*accel(x+xk2/2,v+vk2/2,t+tstep/2,w,g,q)
  129.     LET xk4=tstep*(v+vk3)
  130.     LET vk4=tstep*accel(x+xk3,v+vk3,t+tstep,w,g,q)
  131.     LET vnew=v+(vk1+2*vk2+2*vk3+vk4)/6
  132.     LET xnew=x+(xk1+2*xk2+2*xk3+xk4)/6
  133. END SUB
  134. !
  135. !
  136. DEF accel(x,v,t,w,g,q)
  137.     LET accel=-sin(x)-(1/q)*v+g*cos(w*t)
  138. END DEF
  139. !
  140. SUB PARAMS(W,EPS,TSTEP)
  141.     LET W=.6666666666
  142.     LET EPS=1E-6
  143.     LET TSTEP=.5
  144. END SUB
  145.