'Program FIG10.BAS calculates quasiperiodic regions in a 3-D stereogram 'by J. C. Sprott 'Compile with the PowerBASIC Console Compiler DECLARE FUNCTION gauss AS EXT DEFEXT a-z 'do all calculations in extended (80-bit) precision FUNCTION PBMAIN() n&=3 'number of nonlinear odes h=0.05 'stepsize of integration tmax=1e3 'time to spend at each initial condition pi=4#*ATN(1##) s=pi 'scale factor on graph CONSOLE NAME "Fig 10" CONSOLE SET LOC 808,0 GRAPHIC WINDOW "Fig 10", 0, 0, 800, 400 TO hWin1& GRAPHIC ATTACH hWin1&, 0 GRAPHIC WIDTH 1 GRAPHIC COLOR RGB(0,0,0), RGB(255,255,255) GRAPHIC CLEAR GRAPHIC ATTACH hWin1&, 0 GRAPHIC BOX(0,0)-(800,400) GRAPHIC LINE(400,0)-(400,800) GRAPHIC BITMAP LOAD "fig10.bmp", 800, 400 TO hBmp??? GRAPHIC COPY hBmp???, hWin1& DIM xo(n&), x(n&), xb(n&), dxdt(n&),zb&(800,400) RANDOMIZE TIMER xb(2)=pi/2 eps=.01 WHILE INKEY\$<>CHR\$(27) 'try different random initial conditions FOR i&=1 TO n& 'explore the neighborhood of the last solution xo(i&)=xb(i&)+eps*gauss WHILE xo(i&)>pi: xo(i&)=xo(i&)-2##*pi: WEND WHILE xo(i&)<-pi: xo(i&)=xo(i&)+2##*pi: WEND x(i&)=xo(i&) NEXT i& eps=1.01##*eps IF eps>2## THEN eps=2## xp&=200*xo(1) zp&=200*xo(3) xp1&=200+(xp&+.05*zp&)/s xp2&=600+(xp&-.05*zp&)/s yp&=200-200*xo(2)/s IF xp1&<=1 THEN ITERATE IF xp2&>=799 THEN ITERATE IF yp&<=1 THEN ITERATE IF yp&>=399 THEN ITERATE zb1&=zb&(xp1&,yp&) 'this is the z-buffer zb2&=zb&(xp2&,yp&) IF zp&0 THEN ITERATE zb&(xp1&,yp&)=zp& zb&(xp2&,yp&)=zp& FOR t=0 TO tmax/h CALL rk4(x(),dxdt(),n&,h,x()) NEXT t d2=0: FOR i&=1 TO n&: d2=d2+x(i&)^2: NEXT i& x2=0: FOR i&=1 TO n&: x2=MAX(x2,x(i&)^2): NEXT i& IF d2>0.16*tmax^2 THEN 'orbit escaped IF x(1)^2>x(2)^2 AND x(1)^2>x(3)^2 THEN c&=RGB(255,0,0): IF x(1)<0 THEN c&=RGB(127,127,0) IF x(2)^2>x(3)^2 AND x(2)^2>x(1)^2 THEN c&=RGB(0,255,0): IF x(2)<0 THEN c&=RGB(0,255,255) IF x(3)^2>x(1)^2 AND x(3)^2>x(2)^2 THEN c&=RGB(0,0,255): IF x(3)<0 THEN c&=RGB(255,0,255) PRINT CLNG(d2),CLNG(d2-x2),ROUND(eps,3) IF zp&>zb1& THEN GRAPHIC SET PIXEL(xp1&,yp&),c& IF zp&>zb2& THEN GRAPHIC SET PIXEL(xp2&,yp&),c& FOR i&=1 TO n&: xb(i&)=xo(i&): NEXT i& eps=.01## END IF WEND GRAPHIC SAVE"fig10.bmp" BEEP WAITKEY\$ END FUNCTION SUB DERIVS (x(), dxdt(), n&) 'Returns the time derivatives dxdt(i&) of x(i&) FOR i&=1 TO n& i1&=1+i& MOD n& dxdt(i&)=SIN(x(i1&)) NEXT i& END SUB SUB RK4 (X(), DXDT(), N&, H, XNEW()) 'Fourth-order Runge-Kutta integrator DIM XT(N&), DXT(N&), DXM(N&) HH = H * .5## H6 = H / 6## CALL DERIVS(X(), DXDT(), N&) FOR I& = 1 TO N& XT(I&) = X(I&) + HH * DXDT(I&) NEXT I& CALL DERIVS(XT(), DXT(), N&) FOR I& = 1 TO N& XT(I&) = X(I&) + HH * DXT(I&) NEXT I& CALL DERIVS(XT(), DXM(), N&) FOR I& = 1 TO N& XT(I&) = X(I&) + H * DXM(I&) DXM(I&) = DXT(I&) + DXM(I&) NEXT I& CALL DERIVS(XT(), DXT(), N&) FOR I& = 1 TO N& XNEW(I&) = X(I&) + H6 * (DXDT(I&) + DXT(I&) + 2## * DXM(I&)) NEXT I& ERASE DXM, DXT, XT END SUB FUNCTION gauss 'Returns a normally (Gaussian) distributed random deviate 'with mean of zero and standard deviation of 1.0 DO V1 = 2 * RND - 1 V2 = 2 * RND - 1 R = V1 * V1 + V2 * V2 LOOP WHILE R >= 1 OR R = 0 GAUSS = V2 * SQR(-2 * LOG(R) / R) END FUNCTION