PROGRAM RAD1
	DIMENSION T(8)
	WRITE(6,1003)
        open(22,file='radi.out')
        PRINT*,'Type t(i),i=1,8'
	READ(5,*)(T(I),I=1,8)
1001    FORMAT(8A8)
	WRITE(6,1001)(T(I),I=1,8)
999     CONTINUE
        PRINT*,'Type N,L,IZ'
	READ(5,*)N,L,IZ
	IF(N.EQ.0)STOP
2001    FORMAT(3I10)
	WRITE(6,2002)N,L,IZ
	PI=3.1416
1003    FORMAT(1H1,1X,'RADIAL WAVEFUNCTIONS AND RADIAL DISTRIBUTION', 
     1              /,'FUNCTIONS FOR HYDROGEN LIKE ORBITALS',///)
2002    FORMAT(1X,'N=',I3,2X,'L=',I3,2X,'IZ=',I3,//)
	WRITE(6,1005)
1005    FORMAT(8X,'R/CM',11X,'RNL',12X,'RSQD',11X,'RDF',/)
	A0=0.529E-8
	DO 10 J=1,50
	R=(J-1)*0.25E-8
	F=4*PI*R**2
	RHO=2*IZ*R/(N*A0)
	FNL=IFAC(N+L)
	FNORM=-SQRT(((2*IZ/(N*A0))**3*IFAC(N-L-1)/(2*FNL**3*N)))
	CALL LAGUERRE(PLANL ,N,L,RHO)
	RHOL=1
	IF(L.GT.0)RHOL=RHO**L
	RNL=FNORM*RHOL*PLANL*EXP(-RHO/2)
	RSQ=RNL*RNL
	RDF=F*RSQ
	WRITE(22,106)R,RNL,RSQ,RDF
	WRITE(6,1006)R,RNL,RSQ,RDF
106     FORMAT(1X,4E15.4)
1006    FORMAT(1X,4E15.4,/)
10      CONTINUE
	GO TO 999
	END



	SUBROUTINE LAGUERRE(PLANL,N,L,RHO)
	LIM=N-L-1
	IF(LIM.GE.0)GO TO 1
	WRITE(6,1001)N,L
1001    FORMAT(1X,'FAULT IN Q, NOS N=',I3,2X,'L=',I3,/)
	STOP
1       CONTINUE
	SUM=0.0
	LUP=LIM+1
	DO 3 K1=1,LUP
	K=K1-1
	IF1=IFAC(N+L)
	RHOK=1.0
	IF(K.GT.0)RHOK=RHO**K
3       SUM=SUM+(-1)**(K+1)*IF1*IF1*RHOK/(IFAC(N-L-1-K)*IFAC(2*L+1+K)
     1  *IFAC(K))
	PLANL=SUM
	RETURN
	END
	FUNCTION IFAC(N)
	IPROD=1
	IF(N.EQ.0)GO TO 40
20      DO 30 I=1,N
30      IPROD=IPROD*I
40      IFAC=IPROD
	RETURN
	END