C+
C	Title:	JULIA
C	Author:	T. R. Wyant
C	Date:	27-Oct-1987
C	Modified:
C
C	Remarks:
C		Computes the Julia set around any point. The
C		Julia set is defined as the set of points for
C		which the complex equation
C			Z <- Z*Z + C
C		remains finite no matter how often it is iterated.
C		The user enters a value for the complex constant C,
C		and a range of starting values for Z. Points which
C		still have an absolute value less than 4 after 100
C		iterations are considered part of the set.
C-
	INTEGER*2	LUNTO		! Terminal output.
	PARAMETER	(LUNTO = 5)	! Terminal output to LUN 5.
	INTEGER*2	LUNTI		! Terminal input.
	PARAMETER	(LUNTI = 5)	! Terminal input from LUN 5.
	INTEGER*2	LUNOUT		! Graphics output LUN.
	PARAMETER	(LUNOUT = 4)	! Graphics output to LUN 4.
	INTEGER*4	OUTMAX		! Maximum bits in output.
	PARAMETER	(OUTMAX = 8*144)! Will handle 8 inches @144 bits
	INTEGER*2	HUEMAX		! Maximum colors.
	PARAMETER	(HUEMAX = 15)	! Max. is 16 (counting 0).
	REAL*4		OUTWID		! Output width, inches.
	PARAMETER	(OUTWID = 8.)	! Output is 8 inches wide.
	INTEGER*2	BITMAX		! Max. bit to use in mask.
	PARAMETER	(BITMAX = 6)	! 6 bits in a sixel.
	LOGICAL*1	MASKMX		! Highest set bit.
	PARAMETER	(MASKMX = 2**(BITMAX-1))

	CHARACTER*1	ANSWER		! Answer to questions.
	REAL*4		ASPECT		! Pixel aspect ratio (vert/hor)
	INTEGER*2	BANDWD		! Coutour band width.
	INTEGER*2	BITHI		! Highest bit.
	INTEGER*2	BITLO		! Lowest bit.
	COMPLEX*8	CCONS		! Complex constant.
	COMPLEX*8	CITERA		! Complex iterate
	INTEGER*2	COLCU		! Current column.
	INTEGER*2	COLHI		! High column number for this point.
	INTEGER*2	COLIX		! Step index for columns.
	INTEGER*2	COLLO		! Low column number for this point.
	INTEGER*4	COLMAX		! Maximum columns displayed.
	INTEGER*2	COLOCC		! Highest occupied column.
	CHARACTER*40	FILNAM		! Output file name.
	INTEGER*2	HORHI		! Horizontal end.
	INTEGER*2	HORLO		! Horizontal start.
	REAL*4		HORSIZ		! Horizontal size.
	INTEGER*2	HORST		! Horizontal steps.
	INTEGER*2	HUECUR		! Current color.
	REAL*4		IMAGCN		! Imaginary center part.
	REAL*4		IMAGCU		! Current value of imaginary part.
	INTEGER*2	ITER8		! Iteration counter.
	INTEGER*2	ITERMX		! Max. value of iterate.
	INTEGER*2	LANDSC		! .TRUE. if landscape mode.
	INTEGER*2	LDGFF		! .TRUE. for leading formfeed.
	INTEGER*2	MANDEL		! .TRUE. if finding Mandelbrot set.
	LOGICAL*1	MASKRG (9)	! Masks for setting range of bits.
	LOGICAL*1	MASKSX		! Sixel mask for current row.
	INTEGER*2	MXCMAP		! Max. color map entries (0=B&W)
	INTEGER*2	NEGFLG		! .TRUE. for negative image.
	REAL*4		PIXHOR		! Pixels per inch, horizontal.
	REAL*4		RCONS (2)	! Real constant.
	REAL*4		REALCN		! Real center part.
	REAL*4		REALCU		! Current value of real part.
	REAL*4		RITERA (2)	! Real iterate.
	INTEGER*2	ROWHI		! High row number for this point.
	INTEGER*2	ROWIX		! Step index for rows.
	INTEGER*2	ROWLO		! Low row number for this point.
	INTEGER*4	ROWMAX		! Maximum row number.
	REAL*4		SCLFAC		! Scale factor.
	LOGICAL*1	SIXELS (OUTMAX, HUEMAX)	! Sixel array.
	INTEGER*2	VERHI		! Vertical end.
	INTEGER*2	VERLO		! Vertical start.
	REAL*4		VERSIZ		! Vertical size.
	INTEGER*2	VERST		! Vertical steps.


	EQUIVALENCE	(CCONS, RCONS)
	EQUIVALENCE	(CITERA, RITERA)

	DATA	MASKRG	/'377'O, '376'O, '374'O, '370'O,
     1			'360'O, '340'O, '300'O, '200'O, '0'O/

	ITERMX = 100

	WRITE (LUNTO, 1010) ' '
1010	FORMAT (16A)
1020	FORMAT (F20.0)
	WRITE (LUNTO, 1010) '$Do you want to compute the ',
     1		'Mandelbrot set (no=Julia set)? [Y/N d:N]: '
	READ (LUNTI, 1010,END=9990) ANSWER
	MANDEL = INDEX('TtYy', ANSWER) .GT. 0
	WRITE (LUNTO, 1010) '$Do you want a negative image? [Y/N d:N]: '
	READ (LUNTI, 1010,END=9990) ANSWER
	NEGFLG = (INDEX('TtYy', ANSWER) .GT. 0)
	WRITE (LUNTO, 1010) '$Enter width of contour bands (0=none): '
	READ (LUNTI, 1030, END=9990) BANDWD
	IF (BANDWD .LE. 0 .OR. BANDWD .GT. ITERMX) THEN
	    MXCMAP = 0
	    BANDWD = ITERMX
	  ELSE
	    WRITE (LUNTO, 1010) '$Enter number of color map entries',
     1			' (0=B&W): '
	    READ (LUNTI, 1030, END=9990) MXCMAP
	  END IF
	MXCMAP = MIN (MAX (MXCMAP, 2), HUEMAX + 1)
	IF (.NOT. MANDEL) THEN
	    WRITE (LUNTO, 1010) '$Enter real part of constant: '
	    READ (LUNTI, 1020 ,END=9990) RCONS(1)
	    WRITE (LUNTO, 1010) '$Enter imaginary part of constant: '
	    READ (LUNTI, 1020 ,END=9990) RCONS(2)
	  END IF
	WRITE (LUNTO, 1010) '$Do you want a landscape image? [Y/N d:N]: '
	READ (LUNTI, 1010,END=9990) ANSWER
	LANDSC = (INDEX('TtYy', ANSWER) .GT. 0)
	WRITE (LUNTO, 1010) '$Enter center real part: '
	READ (LUNTI, 1020, END=9990) REALCN
	WRITE (LUNTO, 1010) '$Enter center imaginary part: '
	READ (LUNTI, 1020, END=9990) IMAGCN
	WRITE (LUNTO, 1010) '$Enter scale factor: '
	READ (LUNTI, 1020, END=9990) SCLFAC
	WRITE (LUNTO, 1010) '$Enter horizontal size (d: 8.5): '
	READ (LUNTI, 1020, END=9990) HORSIZ
	IF (HORSIZ .LE. 0.) HORSIZ = 8.5
	WRITE (LUNTO, 1010) '$Enter vertical size (d: 11.0): '
	READ (LUNTI, 1020, END=9990) VERSIZ
	IF (VERSIZ .LE. 0.) VERSIZ = 11.
	WRITE (LUNTO, 1010) '$Enter number of horizontal steps: '
	READ (LUNTI, 1030, END=9990) HORST
1030	FORMAT (I)
	IF (HORST .LE. 0) GO TO 9990

	WRITE (LUNTO, 1010) '$Enter output file name: '
	READ (LUNTI, 1010, END=9990) FILNAM
	WRITE (LUNTO, 1010) '$Do you want leading and trailing ',
     1		'formfeeds? [Y/N d:N]: '
	READ (LUNTI, 1010,END=9990) ANSWER
	LDGFF = INDEX('TtYy', ANSWER) .GT. 0
	WRITE (LUNTO, 1010) '$Enter horizontal pixels/inch ',
     1		'(default: 144.): '
	READ (LUNTI, 1020, END=9990) PIXHOR
	IF (PIXHOR .LE. 0.) PIXHOR = 144.
	WRITE (LUNTO, 1010) '$Enter pixel aspect ratio (vert/horiz)',
     1		'(default: 2.): '
	READ (LUNTI, 1020, END=9990) ASPECT
	IF (ASPECT .LE. 0.) ASPECT = 2.

	VERST = VERSIZ/HORSIZ * HORST	! Find vertical steps.
	HORLO = -HORST/2		! Find horizontal start.
	HORHI = HORST + HORLO		! Find horizontal end.
	VERLO = -VERST/2		! Find vertical start.
	VERHI = VERST + VERLO		! Find vertical end.
	COLMAX = HORSIZ*PIXHOR		! Find maximum column.
	ROWMAX = VERSIZ*PIXHOR/ASPECT	! Find maximum row.

	OPEN (UNIT=LUNOUT, NAME=FILNAM, TYPE='NEW',
     1		CARRIAGECONTROL='LIST')

	IF (LDGFF) WRITE (LUNOUT, 4420) 12
	IF (MANDEL) THEN
	    WRITE (LUNOUT, 4070) REALCN, IMAGCN, SCLFAC, HORST, VERST
4070	    FORMAT ('Mandelbrot set - Iterates of Z <- Z*Z + C', /,
     1		    '"C" Centered on (', F12.6, ',', F12.6, ')', /,
     2		    'Scale factor is', F12.2, /,
     3		    'Horizontal steps:', I5, 5X, 'Vertical steps:',
     4				I5, /, /, /)
	  ELSE
	    WRITE (LUNOUT, 4080) RCONS, REALCN, IMAGCN, SCLFAC,
     1			HORST, VERST
4080	    FORMAT ('Julia set - Iterates of Z <- Z*Z + C', /,
     1		    '"C" = (', F12.6, ',', F12.6, ')', /,
     2		    '"Z" Centered on (', F12.6, ',', F12.6, ')', /,
     3		    'Scale factor is', F12.2, /,
     4		    'Horizontal steps:', I5, 5X, 'Vertical steps:',
     5				I5, /, /, /)
	  END IF
	SCLFAC = SCLFAC * HORST / COLMAX * PIXHOR
	WRITE (LUNOUT, 4420) 27, 'P', 'q'

	ROWHI = 0
	DO 4690 ROWIX = VERLO, VERHI
	    ROWLO = ROWHI + 1
	    ROWHI = (ROWIX - VERLO + 1) * ROWMAX / VERST
	    BITLO = MOD (ROWLO-1, BITMAX) + 1
	    BITHI = ROWHI - ROWLO + BITLO
	    MASKSX = MASKRG (BITLO) .AND. .NOT. MASKRG (BITMAX+1)
	    IF (BITHI .LT. BITMAX)
     1		MASKSX = MASKSX .AND. .NOT. MASKRG (BITHI+1)
	    IF (LANDSC) THEN
		IMAGCU = - ROWIX / SCLFAC + IMAGCN
	      ELSE
		REALCU = ROWIX / SCLFAC + REALCN
	      END IF
	    COLHI = 0
	    DO 4290 COLIX = HORLO, HORHI
		COLLO = COLHI + 1
		COLHI = (COLIX - HORLO + 1) * COLMAX / HORST
		IF (LANDSC) THEN
		    REALCU = COLIX / SCLFAC + REALCN
		  ELSE
		    IMAGCU = COLIX / SCLFAC + IMAGCN
		  END IF
		IF (MANDEL) THEN
		    RCONS (1) = REALCU
		    RCONS (2) = IMAGCU
		    RITERA (1) = 0.
		    RITERA (2) = 0.
		  ELSE
		    RITERA (1) = REALCU
		    RITERA (2) = IMAGCU
		  END IF
		ITER8 = ITERMX
4220		CITERA = CITERA*CITERA + CCONS
		ITER8 = ITER8 - 1
		IF (RITERA(1)*RITERA(1) + RITERA(2)*RITERA(2) .LT. 4.
     1			.AND. ITER8 .GT. 0) GO TO 4220

		HUECUR = MOD ((ITER8 + BANDWD - 1)/BANDWD, MXCMAP)
		IF (.NOT. NEGFLG) HUECUR = MXCMAP - HUECUR - 1
		IF (HUECUR .GT. 0) THEN
		    DO 4240 COLCU = COLLO, COLHI
			SIXELS (COLCU, HUECUR) =
     1				SIXELS (COLCU, HUECUR) .OR. MASKSX
4240			CONTINUE
		  END IF
4290		CONTINUE

4400		CONTINUE
		IF (BITHI .GE. BITMAX) THEN
		    DO 4450 HUECUR = 1, MXCMAP - 1
			COLOCC = 0
			DO 4410 COLCU = 2, COLMAX, 2
			    SIXELS (COLCU, HUECUR) = 0
			    IF (SIXELS(COLCU-1, HUECUR) .NE. 0)
     1				COLOCC = COLCU - 1
4410			    CONTINUE
			IF (COLOCC .GT. 0) THEN
			    IF (MXCMAP .GT. 2)
     1				WRITE (LUNOUT, 4415) HUECUR
4415			    FORMAT ('#', I2.2)
			    DO 4440 COLLO = 1, COLOCC, 80
				COLHI = MIN (COLLO+79, COLOCC)
				WRITE (LUNOUT, 4420)
     1				        (SIXELS(COLCU, HUECUR)+'77'O,
     2					COLCU = COLLO, COLHI)
4420				FORMAT (80A1)
4440				CONTINUE
			  END IF
			IF (MXCMAP .GT. 2 .AND. HUECUR .LT. MXCMAP-1)
     1				THEN
			    IF (COLOCC .GT. 0)
     1				WRITE (LUNOUT, 4420) '$'
			  ELSE
			    WRITE (LUNOUT, 4420) '-'
			  END IF
4450			CONTINUE
		    BITLO = 1
		    BITHI = BITHI - BITMAX
		    MASKSX = MASKRG (BITLO) .AND.
     1			.NOT. MASKRG (BITMAX+1)
		    IF (BITHI .LT. BITMAX)
     1			MASKSX = MASKSX .AND. .NOT. MASKRG (BITHI+1)
		    DO 4470 HUECUR = 1, MXCMAP - 1
			DO 4460 COLCU = 1, COLMAX
			    IF ((SIXELS (COLCU, HUECUR) .AND. MASKMX)
     1					.NE. 0) THEN
				SIXELS (COLCU, HUECUR) = MASKSX
			      ELSE
				SIXELS (COLCU, HUECUR) = 0
			      END IF
4460			    CONTINUE
4470			CONTINUE
		    GO TO 4400
		  END IF
4690	    CONTINUE

	WRITE (LUNOUT, 4420) 27, '\'
	IF (LDGFF) WRITE (LUNOUT, 4420) 12
	CLOSE (UNIT=LUNOUT)

9990	CALL EXIT
	END