|
|
マンデルブロートを計算し、ボリュームレンダリングによって
立体的に表示させます。3次元なので計算に時間がかかります。
また、ファイルサイズも大きくなります。
マシンパワーによって格子のサイズ及び繰り返し数(MAXITER)を調整してください。
LET XSIZE=100 !'各格子の大きさ
LET YSIZE=100
LET ZSIZE=100
OPTION BASE 0
DIM CC(3),ZZ(3)
READ XS,XE,YS,YE,ZS,ZE,MAXITER
DATA -2.2,0.5,-1.35,1.35,-1.35,1.35,50
FILE GETSAVENAME F$,"vtkファイル|*.vtk"
IF F$="" THEN STOP
IF POS(UCASE$(F$),".VTK")=0 THEN LET F$=F$&".vtk"
OPEN #1:NAME F$
ERASE #1
PRINT #1:"# vtk DataFile Version 2.0"
PRINT #1:"mandel"
PRINT #1:"ASCII"
PRINT #1:"DATASET STRUCTURED_POINTS"
PRINT #1:"DIMENSIONS";XSIZE;YSIZE;ZSIZE
PRINT #1:"ORIGIN 0 0 0"
PRINT #1:"SPACING 1 1 1"
PRINT #1:"POINT_DATA";XSIZE*YSIZE*ZSIZE
PRINT #1:"SCALARS scalar float"
PRINT #1:"LOOKUP_TABLE default"
FOR Z=0 TO ZSIZE-1
FOR Y=0 TO YSIZE-1
FOR X=0 TO XSIZE-1
LET CC(0)=XS+X*(XE-XS)/XSIZE
LET CC(1)=YS+Y*(YE-YS)/YSIZE
LET CC(2)=ZS+Z*(ZE-ZS)/ZSIZE
LET CC(3)=0
MAT ZZ=ZER
FOR I=1 TO MAXITER
CALL QMUL(ZZ,ZZ,ZZ) !' CALL HMUL(ZZ,ZZ,ZZ)
CALL ADD(ZZ,ZZ,CC)
IF ZZ(0)^2+ZZ(1)^2+ZZ(2)^2>4 THEN EXIT FOR
NEXT I
IF I>=MAXITER THEN LET DA=255 ELSE LET DA=0
PRINT #1:DA
NEXT X
NEXT Y
NEXT Z
CLOSE #1
END
EXTERNAL SUB QMUL(S(),A(),B()) !'QUATERNION 掛け算 S=A*B
OPTION BASE 0
DIM SS(3)
LET SS(0)=A(0)*B(0)-A(1)*B(1)-A(2)*B(2)-A(3)*B(3)
LET SS(1)=A(0)*B(1)+A(1)*B(0)+A(2)*B(3)-A(3)*B(2)
LET SS(2)=A(0)*B(2)-A(1)*B(3)+A(2)*B(0)+A(3)*B(1)
LET SS(3)=A(0)*B(3)+A(1)*B(2)-A(2)*B(1)+A(3)*B(0)
MAT S=SS
END SUB
EXTERNAL SUB HMUL(S(),A(),B()) !'HYPER COMPLEX 掛け算 S=A*B
OPTION BASE 0
DIM SS(3)
LET P1=SQR(A(0)^2+A(1)^2)
LET P2=SQR(B(0)^2+B(1)^2)
IF P1<>0 AND P2<>0 THEN
LET AA=1-A(2)*B(2)/(P1*P2)
LET SS(0)=AA*(A(0)*B(0)-A(1)*B(1))
LET SS(1)=AA*(B(0)*A(1)+A(0)*B(1))
END IF
LET SS(2)=P1*B(2)+P2*A(2)
MAT S=SS
END SUB
EXTERNAL SUB ADD(S(),A(),B()) !'足し算 S=A+B
OPTION BASE 0
DIM SS(3)
MAT SS=A+B
MAT S=SS
END SUB
|
|