	PROGRAM D10R4
C	Driver for routine DBRENT
	PARAMETER(TOL=1.0E-6,EQL=1.E-4)
	EXTERNAL BESSJ0,DERIV
	DIMENSION AMIN(20)
	NMIN=0
	WRITE(*,'(/1X,A)') 'Minima of the function BESSJ0'
	WRITE(*,'(/1X,T6,A,T19,A,T27,A,T39,A,T53,A/)') 'Min. #','X',
     *		'BESSJ0(X)','BESSJ1(X)','DBRENT'
	DO 12 I=1,100
		AX=I
		BX=I+1.0
		CALL MNBRAK(AX,BX,CX,FA,FB,FC,BESSJ0)
		DBR=DBRENT(AX,BX,CX,BESSJ0,DERIV,TOL,XMIN)
		IF (NMIN.EQ.0) THEN
			AMIN(1)=XMIN
			NMIN=1
			WRITE(*,'(1X,5X,I2,3X,4F12.6)') NMIN,XMIN,
     *				BESSJ0(XMIN),DERIV(XMIN),DBR
		ELSE
			IFLAG=0
			DO 11 J=1,NMIN
				IF (ABS(XMIN-AMIN(J)).LE.EQL*XMIN) IFLAG=1
11			CONTINUE
			IF (IFLAG.EQ.0) THEN
				NMIN=NMIN+1
				AMIN(NMIN)=XMIN
				WRITE(*,'(1X,5X,I2,3X,4F12.6)') NMIN,XMIN,
     *					BESSJ0(XMIN),DERIV(XMIN),DBR
			ENDIF
		ENDIF
12	CONTINUE
	END

	FUNCTION DERIV(X)
	DERIV=-BESSJ1(X)
	END
	PROGRAM D5R2
C	Driver for routine DDPOLY
C	Polynomial (X-1)**5
	PARAMETER(NC=6,NCM1=5,NP=20)
	DIMENSION C(NC),PD(NCM1),D(NCM1,NP)
	CHARACTER A(NCM1)*15
	DATA A/'polynomial:','first deriv:','second deriv:',
     *		'third deriv:','fourth deriv:'/
	DATA C/-1.0,5.0,-10.0,10.0,-5.0,1.0/
	DO 12 I=1,NP
		X=0.1*I
		CALL DDPOLY(C,NC,X,PD,NC-1)
		DO 11 J=1,NC-1
			D(J,I)=PD(J)
11		CONTINUE
12	CONTINUE
	DO 14 I=1,NC-1
		WRITE(*,'(1X,T7,A)') A(I)
		WRITE(*,'(1X,T13,A,T25,A,T40,A)') 'X','DDPOLY','actual'
		DO 13 J=1,NP
			X=0.1*J
			WRITE(*,'(1X,3F15.6)') X,D(I,J),
     *  		    FACTRL(NC-1)/FACTRL(NC-I)*((X-1.0)**(NC-I))
13		CONTINUE
		WRITE(*,*) 'press ENTER to continue...'
		READ(*,*)
14	CONTINUE
	END
	PROGRAM D7R13
C	Driver for routine DES
	DIMENSION IN(64),KEY(64),IOUT(64),ICMP(64)
	INTEGER HEX2IN
	CHARACTER HIN(17)*1,HKEY(17)*1,HOUT(17)*1,HCMP(17)*1,
     *		  VERDCT*8,TEXT*60,TEXT2*6,IN2HEX*1
	OPEN(5,FILE='DESTST.DAT',STATUS='OLD')
	READ (5,'(A)') TEXT
	WRITE(*,'(/1X,A)') TEXT
5	READ (5,'(A)') TEXT
	WRITE(*,'(/1X,A,/)') TEXT
	READ (5,*,ERR=99) NCIPHR
	READ (5,'(A)') TEXT2
	IF (TEXT2.EQ.'encode') IDIREC=0
	IF (TEXT2.EQ.'decode') IDIREC=1
10	WRITE(*,'(1X,T9,A,T23,A,T38,A,T56,A)')
     *		'Key','Plaintext','Expected Cipher','Actual Cipher'
	MM=MIN(NCIPHR,16)
	NCIPHR=NCIPHR-16
	DO 16 M=1,MM
		READ(5,'(51A1)') (HKEY(K),K=1,17),(HIN(K),K=1,17),
     *			(HCMP(K),K=1,17)
		DO 12 I=1,16
			J=I+1
			IDUM=HEX2IN(HIN(J))
			JDUM=HEX2IN(HKEY(J))
			DO 11 K=1,4
				L=4*I+1-K
				IN(L)=MOD(IDUM,2)
				IDUM=IDUM/2
				KEY(L)=MOD(JDUM,2)
				JDUM=JDUM/2
11			CONTINUE
12		CONTINUE
		NEWKEY=1
		CALL DES(IN,KEY,NEWKEY,IDIREC,IOUT)
		HOUT(1)=' '
		DO 14 I=1,16
			JDUM=0
			DO 13 J=1,4
				JDUM=JDUM+(2**(4-J))*IOUT(4*(I-1)+J)
13			CONTINUE
			HOUT(I+1)=IN2HEX(JDUM)
14		CONTINUE
		VERDCT='    o.k.'
		DO 15 I=1,17
			IF (HCMP(I).NE.HOUT(I)) VERDCT='   wrong'
15		CONTINUE
		WRITE(*,'(1X,68A1,A)') (HKEY(K),K=1,17),(HIN(K),K=1,17),
     *			(HCMP(K),K=1,17),(HOUT(K),K=1,17),VERDCT
16	CONTINUE
	WRITE(*,'(/1X,A)') 'press RETURN to continue...'
	READ(*,*)
	IF (NCIPHR.LE.0) GOTO 5
	GOTO 10
99	CLOSE(5)
	END
	INTEGER FUNCTION HEX2IN(CH)
C	Converts character representing hexadecimal data to its integer
C	value in a machine-independent way
	CHARACTER CH*1
	IF(CH.GE.'0'.AND.CH.LE.'9')THEN
		HEX2IN=ICHAR(CH)-ICHAR('0')
	ELSE
		HEX2IN=ICHAR(CH)-ICHAR('A')+10
	ENDIF
	RETURN
	END
	CHARACTER*1 FUNCTION IN2HEX(I)
C	Inverse of HEX2IN
	IF(I.LE.9)THEN
		IN2HEX=CHAR(I+ICHAR('0'))
	ELSE
		IN2HEX=CHAR(I-10+ICHAR('A'))
	ENDIF
	RETURN
	END
	PROGRAM D7R14
C	Driver for routines KS and CYFUN in file DESKS.FOR
	DIMENSION KEY(64),KN(48),IR(32),IOUT(32)
	CHARACTER TEXT(64)*1
C	First test routine KS
	DO 11 I=1,64
		KEY(I)=MOD(I,2)
		IF (KEY(I).EQ.0) TEXT(I)='-'
		IF (KEY(I).EQ.1) TEXT(I)='*'
11	CONTINUE
	WRITE(*,*) 'Legend:'
	WRITE(*,'(1X,T12,A)') '-=0   *=1 '
	WRITE(*,*) 'Master key:'
	WRITE(*,'(1X,T12,56A1)') (TEXT(I),I=1,56)
	WRITE(*,*) 'Sub-master keys:'
	DO 13 I=1,16
		CALL KS(KEY,I,KN)
		DO 12 K=1,48
			IF (KN(K).EQ.0) TEXT(K)='-'
			IF (KN(K).EQ.1) TEXT(K)='*'
12		CONTINUE
		WRITE(*,'(1X,I6,T12,48A1)') I,(TEXT(J),J=1,48)
13	CONTINUE
	WRITE(*,*) 'press RETURN to continue...'
	READ(*,*)
C	Now test routine CYFUN
	DO 14 I=1,32
		IR(I)=MOD(I,3)
		IR(I)=MOD(IR(I),2)
		IF (IR(I).EQ.0) TEXT(I)='-'
		IF (IR(I).EQ.1) TEXT(I)='*'
14	CONTINUE
	WRITE(*,*) 'Legend:'
	WRITE(*,'(1X,T12,A)') '-=0   *=1 '
	WRITE(*,*) 'Input to cipher function:'
	WRITE(*,'(1X,T12,32A1)') (TEXT(I),I=1,32)
	WRITE(*,*) 'Ciphered output:'
	DO 16 I=1,16
		CALL KS(KEY,I,KN)
		CALL CYFUN(IR,KN,IOUT)
		DO 15 K=1,32
			IF (IOUT(K).EQ.0) TEXT(K)='-'
			IF (IOUT(K).EQ.1) TEXT(K)='*'
15		CONTINUE
		WRITE(*,'(1X,I6,T12,32A1)') I,(TEXT(J),J=1,32)
16	CONTINUE
	END
	PROGRAM D10R10
C	Driver for routine DF1DIM
	PARAMETER(NDIM=3,NMAX=50)
	EXTERNAL DF1DIM
	COMMON /F1COM/ NCOM,PCOM(NMAX),XICOM(NMAX)
	DIMENSION P(NDIM),XI(NDIM)
	DATA P/0.0,0.0,0.0/
	NCOM=NDIM
	WRITE(*,'(/1X,A)') 'Enter vector direction along which to'
	WRITE(*,'(1X,A)') 'plot the function. Minimum is in the'
	WRITE(*,'(1X,A)') 'direction 1.0,1.0,1.0 - Enter X,Y,Z:'
	READ (*,*) (XI(I),I=1,3)
	DO 11 J=1,NDIM
		PCOM(J)=P(J)
		XICOM(J)=XI(J)
11	CONTINUE
	CALL SCRSHO(DF1DIM)
	END

	SUBROUTINE DFUNC(X,DF)
	DIMENSION X(3),DF(3)
	DO 11 I=1,3
		DF(I)=(X(I)-1.0)**2
11	CONTINUE
	RETURN
	END
	PROGRAM D10R11
C	Driver for routine DFPMIN
	PARAMETER(NDIM=3,FTOL=1.0E-6,PIO2=1.5707963)
	DIMENSION P(NDIM)
	WRITE(*,'(/1X,A)') 'Program finds the minimum of a function'
	WRITE(*,'(1X,A)') 'with different trial starting vectors.'
	WRITE(*,'(1X,A)') 'True minimum is (0.5,0.5,0.5)'
	DO 11 K=0,4
		ANGL=PIO2*K/4.0
		P(1)=2.0*COS(ANGL)
		P(2)=2.0*SIN(ANGL)
		P(3)=0.0
		WRITE(*,'(/1X,A,3(F6.4,A))') 'Starting vector: (',
     *			P(1),',',P(2),',',P(3),')'
		CALL DFPMIN(P,NDIM,FTOL,ITER,FRET)
		WRITE(*,'(1X,A,I3)') 'Iterations:',ITER
		WRITE(*,'(1X,A,3(F6.4,A))') 'Solution vector: (',
     *			P(1),',',P(2),',',P(3),')'
		WRITE(*,'(1X,A,E14.6)') 'Func. value at solution',FRET
11	CONTINUE
	END

	FUNCTION FUNC(X)
	DIMENSION X(3)
	FUNC=1.0-BESSJ0(X(1)-0.5)*BESSJ0(X(2)-0.5)*BESSJ0(X(3)-0.5)
	END

	SUBROUTINE DFUNC(X,DF)
	PARAMETER (NMAX=50)
	DIMENSION X(3),DF(NMAX)
	DF(1)=BESSJ1(X(1)-0.5)*BESSJ0(X(2)-0.5)*BESSJ0(X(3)-0.5)
	DF(2)=BESSJ0(X(1)-0.5)*BESSJ1(X(2)-0.5)*BESSJ0(X(3)-0.5)
	DF(3)=BESSJ0(X(1)-0.5)*BESSJ0(X(2)-0.5)*BESSJ1(X(3)-0.5)
	RETURN
	END
