program hydrogenic

c.......Solves Hydrogenic functions 
c.......by Dario Mitnik
 
c        implicit real*8 (a-h,o-z)
	parameter (mxpoints=300,mxpntsnod=100)
        character*30  outputfile
        common/radialcoeff/rhotor,radcof
        parameter(ifact=100)
        common /facts / rlf(0:ifact)
	common/blkfuncdat/enerhyd,iz,nq,lq,iaunits
	common/blkwave/radius(2*mxpoints),psi(2*mxpoints),rmax,nmpts
	data one,two/1.0,2.0/


c------- factorial subroutine generates the log of factorials
        call factt

	ifile=9
	print*,'give Z:'
	read*,iz
5	print*,'give the principal quantum number n:'
	read*,nq
	print*,'give the orbital quantum number l:'
	read*,lq
	if (lq.ge.nq) then
		print*,' Error!! try again '
		go to 5
	endif
        outputfile(1:5) = 'wave.'
50      format(a25)
        print*,'give the name of the output file (extension only)'
        read(*,50) outputfile(6:30)
	print*,' '
	ifile = ifile + 1
	open(unit=ifile,file=outputfile,status='unknown')
	
	iaunits=0	! if "1"=a.u.  else=Rydberg
	
c-------hydrogenic energy
	if (iaunits.eq.1) then
		enerhyd = -((one*iz/nq)**2)/two
 		print*,'energy = : ',enerhyd,' (a.u.)'
 	else
		enerhyd = -(one*iz/nq)**2 
		print*,'energy = : ',enerhyd,' (Rydbergs)'
	endif
	print*,'enerhyd = ::  ',enerhyd	
	print*,'change? ("0"=No) else give the new energy'
	read*,rener
	if (rener.ne.0) enerhyd = rener
	print*,' '
c-------maximum radius maxrad
	if (iaunits.eq.1) then
c....... exp(-85) is the minimum exponent allowed in single precision
c....... exp(-708) is the minimum exponent allowed in double precision
		rmax1 = 85./(sqrt(2*abs(enerhyd)))
	else
		rmax1 = 85./sqrt(abs(enerhyd))
	endif
c	rprom = (3.*nq*nq - lq*(lq+1.) ) / (2.*iz)
c	rmax2 = rprom*40./3.
c	print*,' 10 * 4/3<n> = :',rmax2,
c     +         ' - maximum radius allowed is :',rmax1
	rmax2 = 2.*nq*nq/iz
	print*,' 2/z * n^2 = ',rmax2,
     +          ' - maximum radius allowed is :',rmax1
	rmax = min(rmax1,(rmax2+rmax1)/2.)
	print*,' maximum radius =: ',rmax
	print*,'change? ("0"=No)  else give the new maximum radius'
	read*,rchange
	if (rchange.ne.0) rmax = rchange
	print*,' '
	
c-------number of points nmpts
	nmpts = ( nq - lq )*mxpntsnod
	print*,'number of points / 2. =:',nmpts
	print*,'change? ("0"=No)  else give the new number of points'
	read*,ichange
	if (ichange.ne.0) nmpts = ichange
	if (nmpts.gt.mxpoints) then
		print*,'nmpts = :',nmpts,'  reduced to ',mxpoints
		nmpts=mxpoints
	endif
	print*,' '
	
c-------matching ratio rmatch
      	rmatch1 = rmax2/10.
        rmatch2 = rmax/10.
        rmatch = max(rmatch1,rmatch2)
c        print*,'4/3 * <n>= :',rmatch1,' -  rmax/10=:',rmatch2
        print*,'2/z * n^2= :',rmatch1,' -  rmax/10=:',rmatch2
        print*,'the position of the matching point is:',rmatch
        print*,'change? ("0"=No)  else give the new matching point'
        read*,rchange
        if (rchange.ne.0) rmatch = rchange
        print*,' '

c-------printing the radial hydrogenic wavefunction from r0 to rmatch
	rhotor = (2.*iz/nq)
	radcof = rhotor**1.5
	r0 = 0	 	
 	rdel = (rmatch - r0) / nmpts
 	r = r0 - rdel
 	do 100 i=1,nmpts
 		r = r + rdel
 		hy1 = hydradial(r)
 		write(ifile,25) r,hy1
100	continue
c	printing the radial hydrogenic wavefunction from rmatch to rmax
	rdel = (rmax - r) / nmpts
	do 200 i=2,nmpts
		r = r + rdel
		hy1 = hydradial(r)
		write(ifile,25) r,hy1
200	continue
25	format(3(2x,1pg15.7))
	write(ifile,*) ' '
	
c-------solving the same problem using the Numerov Method
	call hydnumerov(r0,rmatch)

c.......normalization
        call normalization

c.......begin the function with positive slop
	if (psi(2).gt.0) then 
        	do 400 i=1,2*nmpts
400        	write(ifile,40) radius(i),psi(i)
	else
        	do 500 i=1,2*nmpts
        	psi(i) = -1*psi(i)
500        	write(ifile,40) radius(i),psi(i)
	endif
40	format(2x,f12.5,2x,1pg15.6)
	write(ifile,*) ' '

	close(unit=ifile)
	print*,' continue? (no = "0") '
	read*,icont
	if (icont.ne.0) go to 5
	stop
	end	
	

c--------------------------------------------------------------------
	subroutine hydnumerov(r0,rmatch)
c        implicit real*8 (a-h,o-z)
	
c.......driver subroutine for the Numerov method

	parameter (mxpoints=300)
	common/blkfuncdat/enerhyd,iz,nq,lq,iaunits
	common/blkwave/radius(2*mxpoints),psi(2*mxpoints),rmax,nmpts
	
	data yepsilon/2.e-19/ ! value of the wavefunction at rmax
	
	h = ( rmatch - r0 ) / nmpts
	
c-------outgoing integration r0 --> rmatch
	ioutwards = 1	
	x0 = r0
	y0 = 0
	y1 = h**(lq+1)
	print*, ' First step, outgoing integration:'
	call numerov(x0,y0,x1,y1,h,rmatch,ioutwards)
	youtw = psi(nmpts)
	print*,'function at rmatch=:',radius(nmpts),' = :',youtw
		
c--------inward integration x0 = rmax --> rmatch
	ioutwards=0
	x0 =  rmax
	y0 = yepsilon
	h = (rmatch-x0)/nmpts
	x1 = x0 + h
c........ initial value of the function at x=rmax-|h|
	if (iaunits.eq.1) then
c...............y(x) = const * exp(-sqrt(-2*E)*x)   (a.u.)
 		eners = sqrt(-2*enerhyd)
 	else
c---------------y(x) = const * exp(-sqrt(-E)*x)    (Ryd)
		eners = sqrt(-1*enerhyd)
	endif
	expo = (x0-x1)*eners
	expo = exp(expo)
	y1 = y0 * expo
	print*, ' Second step, inward integration:'
	radius(2*nmpts) = x0
	psi(2*nmpts) = y0
	call numerov(x0,y0,x1,y1,h,rmatch,ioutwards)
	yinwd = psi(nmpts)
	print*,'function at rmatch=:',radius(nmpts),' = :',yinwd

c.......matchinf factor at r(nmpts)
	wvfactor = yinwd/youtw
	print*,'  matching ratio = :',wvfactor
	
c.......multiplying the inward wavefunction by the matching factor
	isgfactor = 1
	if (wvfactor.lt.0) isgfactor = -1
	do 300 i=1,nmpts-1
300	psi(i) = psi(i)*wvfactor*isgfactor
c.......multiplying the outward wavefunction by the sign of the matching factor
	if (wvfactor.lt.0) then
	do 305 i=nmpts,2*nmpts
305	psi(i) = psi(i)*-1
	endif

	return
	end 

c-----------------------------------------------------------------
	subroutine numerov(x0,y0,x1,y1,h,rmatch,ioutwards)
c        implicit real*8 (a-h,o-z)
	
c.......Numerov Formula for Y''(j) = f(j) * Y(j)

c	R(j+1)*Y(j+1) = 2*R(j)*Y(j) - R(j-1)*Y(j-1) + h^2*f(j)*Y(j)
c	where R(j) = 1 - h^2/12 * f(j)

c	for the hydrogenic problem 
c	f(j) = 2 (Veff(j) - E) (in a.u.) 
c	f(j) =   (Veff(j) - E) (in Rydbergs) 
c	and
c	Veff (j) = V(X(j)) + l(l+1)/x^2 = -Z/X(j)  + l(l+1)/x^2 (a.u.)
c	Veff (j) = V(X(j)) + l(l+1)/x^2 = -2Z/X(j) + l(l+1)/x^2 (Ryd)
c
c	f is calculated in function diffeq

	parameter (mxpoints=300)
	common/blkfuncdat/enerhyd,iz,nq,lq,iaunits
	common/blkwave/radius(2*mxpoints),psi(2*mxpoints),rmax,nmpts
	data one,two,twelve/1.0,2.0,12.0/

	h2 = h*h
	iway = 2*ioutwards - 1
	indxinit = (2*nmpts - 1)  - (2*nmpts - 2)*ioutwards 
c....... Loop in x
	do 100 indxwave=indxinit,nmpts,iway
	x1 = h + x0
	x2 = h + x1
	f0 =  diffeq(x0)
	f1 =  diffeq(x1)
	f2 =  diffeq(x2)
	r0 = one - h2*f0/twelve
	r1 = one - h2*f1/twelve
	r2 = one - h2*f2/twelve
	r2y2 = two*r1*y1 - r0*y0 + h2*f1*y1
	y2 = r2y2 / r2
	radius(indxwave) = x1
	psi(indxwave) = y1
c.......storing the new values
	x0 = x1
	y0 = y1
	y1 = y2	
100	continue
	return
	end

c-----------------------------------------------------------------
	function diffeq(x)
c        implicit real*8 (a-h,o-z)
	common/blkfuncdat/enerhyd,iz,nq,lq,iaunits
	
c	for the hydrogenic problem 
c	f(j) = 2 (Veff(j) - E) (in a.u.) 
c	f(j) =   (Veff(j) - E) (in Rydbergs) 
c	and
c	Veff (j) = V(X(j)) + l(l+1)/x^2 = -Z/X(j)  + l(l+1)/x^2 (a.u.)
c	Veff (j) = V(X(j)) + l(l+1)/x^2 = -2Z/X(j) + l(l+1)/x^2 (Ryd)

	if (x.eq.0) then
		diffeq=0.
		go to 200
	endif
	if (iaunits.eq.1) then
c---------------(a.u.)
		veff = -2.*iz/x + 2*lq*(lq+1.)/(x*x) 
		diffeq = veff - 2*enerhyd
	else
c---------------(Ryd)
		veff = -2.*iz/x +  lq*(lq+1.)/(x*x) 
		diffeq = veff -  enerhyd 
	endif
200	return
	end
	
c-----------------------------------------------------------------
	subroutine normalization
c        implicit real*8 (a-h,o-z)
	parameter (mxpoints=300)
	common/blkwave/radius(2*mxpoints),psi(2*mxpoints),rmax,nmpts
	data yepsilon/2.e-19/ ! value of the wavefunction at rmax
	
	sum = 0.
	h1 = radius(2) - radius(1)
	do 100 i=2,nmpts-1
	if (abs(psi(i)).lt.yepsilon/h1) go to 100
	sum = sum + psi(i)*psi(i)
100	continue
	sum = h1*(sum + (psi(1)**2)/2.)
	if (abs(psi(nmpts)).gt.yepsilon/h1) 
     +      sum = sum + h1*(psi(nmpts)**2)/2.
	
	sum2 = 0.
	h2 = radius(nmpts+2) - radius(nmpts+1)
	do 200 i=nmpts+2,2*nmpts-1
	if (abs(psi(i)).lt.yepsilon/h2) go to 200
 	sum2 = sum2 + psi(i)*psi(i)
c	print*,'i=:',i,'  radius:',radius(i),'   psi^2=:',psi(i)*psi(i)
200	continue
	sum2 = h2*sum2
	if (abs(psi(nmpts+1)).gt.yepsilon/h2) 
     +         sum2 = sum2 + h2*(psi(nmpts+1)**2)/2.
	if (abs(psi(2*nmpts)).gt.yepsilon/h2) 
     +         sum2 = sum2 + h2*(psi(2*nmpts)**2)/2.
	
	cnorm = sqrt(sum + sum2)
	print*,' sqrt(Int(psi*psi)) = :: ',cnorm

	do 300 i=1,2*nmpts
	psi(i) = psi(i)/cnorm
300	continue
	
	return
	end

c-------------------------------------------------------------------------
	function hydradial(r)
c        implicit real*8 (a-h,o-z)
	
c	This function returns the value of the radial 
c	hydrogenic wavefunction Rnl(r)
c	Normalized as Int(0,inf){ |Rnl(r)|^2 r^2 dr } = 1

       	common/radialcoeff/rhotor,radcof
	parameter(ifact=100)
        common /facts / rlf(0:ifact)
	common/blkfuncdat/enerhyd,iz,nq,lq,iaunits
		
c-------  	[  (n-l-1)!	]^1/2
c-------	|  ---------	|
c-------	[  2n[(n+l)!]^3	]

	if (nq-lq-1.gt.ifact.or.nq+lq.gt.ifact)pause 'increase ifact' 
	rf = exp(rlf(nq-lq-1)/2. - (3./2.)*rlf(nq+lq) )
	twon = 2.*nq
	rf = rf/sqrt(twon)
	
c-------  			  [  (n-l-1)!		]^1/2
c-------	- {2Z/(n*a0)}^3/2 |  ---------		|
c-------			  [  2n[(n+l)!]^3	]
c------- 	 
	rf = -1*radcof*rf
	
c------- rho = 2Z/(n*a0) * r
	 rho = rhotor * r

c-------  R(nl) = rf * exp(-rho/2) * rho**l * Laguerre(n+l,2l+1)(rho)
	rexpo = -1*rho/2.
	hydradial = rf * exp(rexpo) * (rho**lq) * aLaguerre(nq,lq,rho)
	return
	end
	
c--------------------------------------------------------------------
	function aLaguerre(nq,lq,rho)
c        implicit real*8 (a-h,o-z)
        parameter(ifact=100)
        common /facts / rlf(0:ifact)

c	Associated Laguerre polynom L(n+l,2l+1)

	nr = nq - lq - 1
	ALaguerre = 0.
	do 100 k=0,nr
	if (nr-k.gt.ifact.or.2*lq+1+k.gt.ifact)pause 'increase ifact' 
	rf = exp(2*rlf(nq+lq) - rlf(nr-k) - rlf(2*lq+1+k) - rlf(k) )
	aLaguerre = aLaguerre + (-1)**(k+1) * rf * (rho**k) 
100	continue
	return
	end
	
c--------------------------------------------------------------------
	subroutine factt
c        implicit real*8 (a-h,o-z)

c       ***  calculates the logs of factorials required by the racah
c       ***  coefficient routine dracah

        parameter(ifact=100)
        common /facts / rlf(0:ifact)

        rlf(0)=1.
        rlf(1)=1.
        x=2.
        do 10 i=2,30
        rlf(i)=rlf(i-1)*x
        x=x+1.
10      continue
        do 20 i=0,30
        rlf(i)=log(rlf(i))
20      continue
        x=30.
        do 30 i=31,ifact
        rlf(i)=rlf(i-1)+log(x)
        x=x+1.
30      continue
	return
	end