シェルピンスキー曲線

 投稿者:しばっち  投稿日:2017年 5月27日(土)09時50分35秒
  OPTION BASE 0
PUBLIC NUMERIC KK,XX(8500),YY(8500)
DIM FR(8500),FI(8500),TX(8500),TY(8500)
CALL GINIT(800,700)
LET N=6
LET B$="LF"
LET L=INT(750/2^N)
LET X0=50
LET Y0=650
LET TH=0
SET LINE COLOR 7
PLOT LINES:X0,Y0;
LET KK=KK+1
LET XX(KK)=X0
LET YY(KK)=Y0
FOR LEV=1 TO N
   LET A$=B$
   LET B$=""
   FOR K=1 TO LEN(A$)
      SELECT CASE MID$(A$,K,1)
      CASE "F"
         LET B$=B$&""
      CASE "L"
         LET B$=B$&"-RF+LF+RF-"
      CASE "R"
         LET B$=B$&"+LF-RF-LF+"
      CASE "+"
         LET B$=B$&"+"
      CASE "-"
         LET B$=B$&"-"
      CASE ELSE
      END SELECT
   NEXT K
NEXT LEV
FOR I=1 TO LEN(B$)
   SELECT CASE MID$(B$,I,1)
   CASE "+"
      LET TH=TH-60
   CASE "-"
      LET TH=TH+60
   CASE "F"
      LET X1=X0+L*COS(TH*PI/180)
      LET Y1=Y0-L*SIN(TH*PI/180)
      PLOT LINES:X1,Y1;
      LET X0=X1
      LET Y0=Y1
      LET KK=KK+1
      LET XX(KK)=X1
      LET YY(KK)=Y1
   CASE "L","R"
   END SELECT
NEXT I
LET LV=INT(LOG2(KK))+1
LET N=2^LV
MAT TX=XX
MAT TY=YY
DO
   MAT XX=TX
   MAT YY=TY
   MAT FR=ZER
   MAT FI=ZER
   FOR I=1 TO KK-1
      LET FR(I+1)=XX(I+1)-XX(I)
      LET FI(I+1)=YY(I+1)-YY(I)
   NEXT I
   CALL CDFT(2*N,COS(PI/N),SIN(PI/N),FR,FI)
   LOCATE VALUE ,RANGE 2 TO 2^(LV-1)-1:NCUT
   LET NCUT=INT(NCUT)
   FOR I=NCUT+1 TO 2^LV-NCUT-1
      LET FR(I)=0
      LET FI(I)=0
   NEXT I
   CALL CDFT(2*N,COS(PI/N),-SIN(PI/N),FR,FI)
   CLEAR
   PLOT LINES
   FOR I=1 TO KK
      LET XX(I+1)=XX(I)+FR(I+1)/N
      LET YY(I+1)=YY(I)+FI(I+1)/N
      IF I=1 THEN PLOT LINES:XX(I),YY(I);
      PLOT LINES:XX(I+1),YY(I+1);
   NEXT I
LOOP
END

EXTERNAL SUB GINIT(XSIZE,YSIZE)
SET BITMAP SIZE XSIZE,YSIZE
SET WINDOW 0,XSIZE-1,YSIZE-1,0
SET POINT STYLE 1
SET COLOR MODE "REGULAR"
SET COLOR MIX(0) 0,0,0
SET COLOR MIX(1) 0,0,1
SET COLOR MIX(2) 1,0,0
SET COLOR MIX(3) 1,0,1
SET COLOR MIX(4) 0,1,0
SET COLOR MIX(5) 0,1,1
SET COLOR MIX(6) 1,1,0
SET COLOR MIX(7) 1,1,1
CLEAR
END SUB

EXTERNAL SUB CDFT(N,WR,WI,AR(),AI())
OPTION BASE 0
DIM A(N)
FOR I=0 TO N/2-1
   LET A(2*I)=AR(I)
   LET A(2*I+1)=AI(I)
NEXT I
LET WMR=WR
LET WMI=WI
LET M=N
DO WHILE M>4
   LET L=M/2
   LET WKR=1
   LET WKI=0
   LET WDR=1-2*WMI*WMI
   LET WDI=2*WMI*WMR
   LET SS=2*WDI
   LET WMR=WDR
   LET WMI=WDI
   FOR J=0 TO N-M STEP M
      LET I=J+L
      LET XR=A(J)-A(I)
      LET XI=A(J+1)-A(I+1)
      LET A(J)=A(J)+A(I)
      LET A(J+1)=A(J+1)+A(I+1)
      LET A(I)=XR
      LET A(I+1)=XI
      LET XR=A(J+2)-A(I+2)
      LET XI=A(J+3)-A(I+3)
      LET A(J+2)=A(J+2)+A(I+2)
      LET A(J+3)=A(J+3)+A(I+3)
      LET A(I+2)=WDR*XR-WDI*XI
      LET A(I+3)=WDR*XI+WDI*XR
   NEXT J
   FOR K=4 TO L-4 STEP 4
      LET WKR=WKR-SS*WDI
      LET WKI=WKI+SS*WDR
      LET WDR=WDR-SS*WKI
      LET WDI=WDI+SS*WKR
      FOR J=K TO N-M+K STEP M
         LET I=J+L
         LET XR=A(J)-A(I)
         LET XI=A(J+1)-A(I+1)
         LET A(J)=A(J)+A(I)
         LET A(J+1)=A(J+1)+A(I+1)
         LET A(I)=WKR*XR-WKI*XI
         LET A(I+1)=WKR*XI+WKI*XR
         LET XR=A(J+2)-A(I+2)
         LET XI=A(J+3)-A(I+3)
         LET A(J+2)=A(J+2)+A(I+2)
         LET A(J+3)=A(J+3)+A(I+3)
         LET A(I+2)=WDR*XR-WDI*XI
         LET A(I+3)=WDR*XI+WDI*XR
      NEXT J
   NEXT K
   LET M=L
LOOP
IF M>2 THEN
   FOR J=0 TO N-4 STEP 4
      LET XR=A(J)-A(J+2)
      LET XI=A(J+1)-A(J+3)
      LET A(J)=A(J)+A(J+2)
      LET A(J+1)=A(J+1)+A(J+3)
      LET A(J+2)=XR
      LET A(J+3)=XI
   NEXT J
END IF
IF N>4  THEN CALL BITRV2(N,A)
FOR I=0 TO N/2-1
   LET AR(I)=A(2*I)
   LET AI(I)=A(2*I+1)
NEXT I
END SUB

EXTERNAL SUB BITRV2(N,A())
LET M=N/4
LET M2=2*M
LET N2=N-2
LET K=0
FOR J=0 TO M2-4 STEP 4
   IF J<K THEN
      LET XR=A(J)
      LET XI=A(J+1)
      LET A(J)=A(K)
      LET A(J+1)=A(K+1)
      LET A(K)=XR
      LET A(K+1)=XI
   ELSEIF J>K THEN
      LET J1=N2-J
      LET K1=N2-K
      LET XR=A(J1)
      LET XI=A(J1+1)
      LET A(J1)=A(K1)
      LET A(J1+1)=A(K1+1)
      LET A(K1)=XR
      LET A(K1+1)=XI
   END IF
   LET K1=M2+K
   LET XR=A(J+2)
   LET XI=A(J+3)
   LET A(J+2)=A(K1)
   LET A(J+3)=A(K1+1)
   LET A(K1)=XR
   LET A(K1+1)=XI
   LET L=M
   DO WHILE K>=L
      LET K=K-L
      LET L=L/2
   LOOP
   LET K=K+L
NEXT J
END SUB

 

戻る