Rcjp's Weblog

June 19, 1992

Fractals

Filed under: physics — rcjp @ 11:27 am

This is a historical blog entry of some old Fortran files I found on my HDD…

I think this was my first attempt to draw fractals, unfortunately now (in 2007) I no longer have the compiler to try reproducing the graphs.


          PROGRAM FRACTAL
          INTEGER NTRANS,XSCALE,YSCALE,XOFFSET,YOFFSET,IY,IX
          REAL A,B,C,D,E,F,P,X,Y,NEWX,NEWY
          DOUBLE PRECISION PK,N
          CHARACTER*9 NAME
          DIMENSION A(100),B(100),C(100),D(100),E(100),F(100),
         * P(100)
          
          PRINT *,'ENTER FILENAME FOR DATA 5CHAR.DAT'     
          READ 100,NAME
          OPEN (5,FILE=NAME,STATUS='OLD')

          CALL ZEROISE(A,B,C,D,E,F,P)
          CALL GETDATA(A,B,C,D,E,F,P,NTRANS,XSCALE,YSCALE)
          CALL SETSCREEN(X,Y,XOFFSET,YOFFSET)
          CALL DATE_TIME_SEED@
         
        1 DO 10 N=1,200000
             PK=RANDOM()
             DO 5 M=NTRANS,1,-1
                IF (PK.LE.P(M)) K=M
        5    CONTINUE
             NEWX=A(K)*X+B(K)*Y+E(K)
             NEWY=C(K)*X+D(K)*Y+F(K)
             X=NEWX
             Y=NEWY
             IX=INT(X*XSCALE+XOFFSET)
             IY=INT(Y*YSCALE+YOFFSET)
             CALL GET_KEY1@(J)
             IF (J.NE.0) GOTO 12
             IF (N.LT.10) GOTO 10
             CALL SET_PIXEL@(IX,IY,14)
       10 CONTINUE
       12 IF (J.EQ.73) THEN
            XSCALE=XSCALE+10
            YSCALE=YSCALE+10
          ENDIF
          IF (J.EQ.79) THEN
            XSCALE=XSCALE-10
            YSCALE=YSCALE-10
          ENDIF
          IF (J.EQ.328) YOFFSET=YOFFSET-10
          IF (J.EQ.336) YOFFSET=YOFFSET+10
          IF (J.EQ.331) XOFFSET=XOFFSET-10
          IF (J.EQ.333) XOFFSET=XOFFSET+10
          CALL CLEAR_SCREEN@
          IF (J.NE.81) GOTO1
          CALL TEXT_MODE@

      100 FORMAT (A) 
          END

          SUBROUTINE ZEROISE(A,B,C,D,E,F,P)
          REAL A,B,C,D,E,F,P,PK
          DIMENSION A(100),B(100),C(100),D(100)
         * ,E(100),F(100),P(100)
          DO 1000 I=1,100
              A(I)=0.0
              B(I)=0.0
              C(I)=0.0
              D(I)=0.0
              E(I)=0.0
              F(I)=0.0
              P(I)=0.0
     1000 CONTINUE
          END

          SUBROUTINE GETDATA(A,B,C,D,E,F,P,NTRANS,XSCALE,YSCALE)
          INTEGER NTRANS,XSCALE,YSCALE
          REAL A,B,C,D,E,F,P,PK,PT
          DIMENSION A(100),B(100),C(100),D(100)
         * ,E(100),F(100),P(100)
          
          READ (5,*) NTRANS
          READ (5,*) XSCALE,YSCALE
          PT=0
          DO 2000 I=1,NTRANS
             READ (5,*) A(I),B(I),C(I),D(I),E(I),F(I),PK
             PT=PT+PK
             P(I)=PT
     2000 CONTINUE
          END

          SUBROUTINE SETSCREEN(X,Y,XOFFSET,YOFFSET)
          INTEGER XOFFSET,YOFFSET
          REAL X,Y
          X=0.0
          Y=0.0
          XOFFSET=220
          YOFFSET=0
          CALL VGA@
          CALL CLEAR_SCREEN@
          END

and there were several datafiles of parameters which I think drew fern patterns etc.


    4
    50 40
    0,0,0,0.16,0,0,0.01
    0.2,-0.26,0.23,0.22,0,1.6,0.07
    -0.15,0.28,0.26,0.24,0,0.44,0.07
    0.85,0.04,-0.04,0.85,0,1.6,0.85

    4
    800 800
    0    0  0  0.5  0  0   0.05
    0.1  0  0  0.1  0  0.2 0.15
    0.42 -0.42 0.42 0.42 0 0.2 0.4
    0.42 0.42 -0.42 0.42 0 0.2 0.4

    3
    500 500
    0.5  0   0   0.5   0   0   0.33
    0.5  0   0   0.5   1   0   0.33
    0.5  0   0   0.5   0.5 0.5 0.34

There was also a very primitive mandelbrot program which I managed to make work under wine and take a screenshot!


          PROGRAM MANDEL
          DIMENSION ZRL(60),ZPX(60)

          SCALEX=266.666
          SCALEY=200.0
          LOFFX=320
          LOFFY=240

          CALL VGA@
          CALL CLEAR_SCREEN@

          DO 30 CPX=1.2,-1.2,-5.0E-3
             DO 20 CRL=-1.2,1.2,3.75E-3
               ZRL(1)=0.0
               ZPX(1)=0.0
               ZRL(2)=ZRL(1)*ZRL(1)-ZPX(1)*ZPX(1)+CRL
               ZPX(2)=2.0*ZRL(1)*ZPX(1)+CPX
               IF (ABS(ZRL(2)).GT.1E10) GOTO 11
               IF (ABS(ZPX(2)).GT.1E10) GOTO 11
               ZRL(3)=ZRL(2)*ZRL(2)-ZPX(2)*ZPX(2)+CRL
               ZPX(3)=2.0*ZRL(2)*ZPX(2)+CPX
               IF (ABS(ZRL(3)).GT.1E10) GOTO 11
               IF (ABS(ZPX(3)).GT.1E10) GOTO 11
                  DO 10 I=3,50,1
                     ZRL(I+1)=ZRL(I)*ZRL(I)-ZPX(I)*ZPX(I)+CRL
                     ZPX(I+1)=2.0*ZRL(I)*ZPX(I)+CPX
                     DET=((ZRL(I)-ZRL(I-1))-(ZRL(I-1)-ZRL(I-2)))/2.0
                     DET1=((ZPX(I)-ZPX(I-1))-(ZPX(I-1)-ZPX(I-2)))/2.0
    C               PRINT *,I,DET,DET1
                     IF (ABS(DET).GT.1.0 .OR. ABS(DET1).GT.1.0) THEN
                       CALL SET_PIXEL@(INT(LOFFX+CRL*SCALEX)
         *             ,INT(LOFFY-CPX*SCALEY),15)
                     GOTO 20
                     ENDIF
       10         CONTINUE
       11       CALL SET_PIXEL@(INT(LOFFX+CRL*SCALEX)
         *     ,INT(LOFFY-CPX*SCALEY),0)
       20     CONTINUE
       30 CONTINUE
          END

Basic Mandelbrot

Advertisements

Leave a Comment »

No comments yet.

RSS feed for comments on this post. TrackBack URI

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

Blog at WordPress.com.

%d bloggers like this: