1010 OPTION ARITHMETIC NATIVE
1020 OPTION ANGLE DEGREES
1030 DECLARE EXTERNAL SUB Plot3d.PLOT, Plot3d.PenUp, PlotLabel
1040 DECLARE PICTURE Lorenz
1050 DECLARE NUMERIC theta0, tx, ty
1060 LET theta0=45      ! Rotation angle about the z-axis
1070 LET tx=-60         ! Rotation angle about the x-axis
1080 LET ty=0           ! Rotation angle about the y-axis
1150 DIM m(4,4),mm(4,4)
1160 MAT m=ROTATE(theta0)
1170 SET POINT STYLE 1
1180 SET WINDOW -120, 120, -60, 180
1200 DIM rotx(4,4)      ! rotation about x-axis
1210 MAT rotx=IDN
1220 LET rotx(2,2)=COS(tx)
1230 LET rotx(2,3)=SIN(tx)
1240 LET rotx(3,2)=-SIN(tx)
1250 LET rotx(3,3)=COS(tx)
1260 DIM roty(4,4)      ! rotation about y-axis
1270 MAT roty=IDN
1280 LET roty(1,1)=COS(ty)
1290 LET roty(1,3)=-SIN(ty)
1300 LET roty(3,1)=SIN(ty)
1310 LET roty(3,3)=COS(ty)
1320 MAT m=m *  rotx * roty 
1330 DRAW Lorenz WITH m
1500 DECLARE NUMERIC s,b,r,t,dt,x,y,z,xx,yy,zz
1510 PICTURE Lorenz
1530    LET s=11
1540    LET b=8/4
1550    LET r=88
1560    LET x=15
1570    LET y=25
1580    LET z=75 
1590    LET dt=1/256
1600    FOR t=0 TO 75 STEP dt
1610       LET xx=x+(-s*x+s*y)*dt
1620       LET yy=y+(r*x-y-x*z)*dt
1630       LET zz=z+(-b*z+x*y)*dt
1640       LET x=xx
1650       LET y=yy
1660       LET z=zz
1670       CALL PLOT(x,y,z)
1680    NEXT t
1690    REM x-axis
1700    CALL PenUp
1710    CALL PLOT(0,0,0)
1720    CALL PLOT(50,0,0)
1730    CALL PLOTLABEL(50,0,0,"x") 
1740    REM y-axis
1750    CALL PenUp
1760    CALL PLOT(0,0,0)
1770    CALL PLOT(0,50,0)
1780    CALL PLOTLABEL(0,50,0,"y") 
1790    REM z-axis
1800    CALL PenUp
1810    CALL PLOT(0,0,0)
1820    CALL PLOT(0,0,180)
1830    CALL PLOTLABEL(0,0,180,"z") 
1840    CALL PenUp
1850 END PICTURE  
1860 END
2000 MODULE PLot3D
2010 MODULE OPTION ARITHMETIC NATIVE
2020 PUBLIC SUB PLOT
2030 SHARE NUMERIC x0,y0
2040 SHARE NUMERIC PenState
2050 LET PenState=0
2060 EXTERNAL SUB PLOT(x,y,z)
2070    DECLARE NUMERIC POINT(4) 
2080    DECLARE NUMERIC m(4,4)
2090    DECLARE NUMERIC xx,yy
2100    MAT m=TRANSFORM
2110    LET POINT(1)=x
2120    LET POINT(2)=y
2130    LET POINT(3)=z
2140    LET POINT(4)=1
2150    MAT POINT=POINT*m
2160    LET xx=POINT(1)/POINT(4)
2170    LET yy=POINT(2)/POINT(4)
2180    IF PenState<>0 THEN GRAPH LINES:x0,y0; xx,yy
2190    LET x0=xx
2200    LET y0=yy
2210    LET PenState=1
2220 END SUB
2230 EXTERNAL SUB PenUp
2240    LET PenState=0
2250 END SUB
2260 END MODULE 
3000 EXTERNAL SUB PLOTLABEL(x,y,z,s$)
3010 OPTION ARITHMETIC NATIVE
3020 DIM m(4,4)
3030 MAT m=TRANSFORM
3040 DIM POINT(4) 
3050 LET POINT(1)=x
3060 LET POINT(2)=y
3070 LET POINT(3)=z
3080 LET POINT(4)=1
3090 MAT POINT=POINT*m
3100 GRAPH TEXT ,AT POINT(1)/POINT(4),POINT(2)/POINT(4):s$
3110 END SUB
   
