      PROGRAM HORBS
C	Program to calculate the wavefunction of Hydrogen like orbitals
C 	and to write out a 2D slice of the 3D cube onto a file in format
C	suitable for reading by UNIMAP. 
C	Henry Rzepa.  October, 1989.
      COMMON /GRID/ G(25,25,25)
C	The dimension of the grid is is a maximum of 25, since anything
C	larger would take too much computer time!
      CHARACTER IONN*4, IOFF*4
      DATA PI /3.14159265/
      IONN = CHAR(27)//'[7m'
      IOFF = CHAR(27)//'[0m'
C	The above shows how too turn inverse video on and off on a
C	terminal, to highlight an error! CHAR(27) is the Escape 
C	character, 7=inverse, 0=attribute off, 1=bold, 5=blink etc
C
      PI = SQRT(1./PI)
C       Read size of grid (max 25)
      WRITE(*,*) 'Enter dimension of grid (Odd number, Max 25)'
      READ(5,*) IG
      IF(IG.GT.25) IG=25
      IGH = (IG/2) 
C     READ THE QUANTUM NUMBERS OF THE ORBITALS
   5  WRITE(*,4)IGH,IGH
   4  FORMAT(' Enter N,L, M, Z Plane (-',I2,' to +',I2') and scale')
      READ(5,*,END=10) N,L,M,IZ,SCALE
      IF(N.EQ.0) GO TO 10
      IF(N.GT.4) WRITE(*,6)IONN,IOFF
    6 FORMAT(1x,A4,'N Too large for program',A4)
      IF(L.GT.(N-1)) WRITE(*,*)'VALUE OF L INCORRECT'
      IF(IABS(M).GT.L)WRITE(*,*)'VALUE OF M INCORRECT'
      IF(N.GT.4 .OR. L.GT.(N-1) .OR. IABS(M).GT.L) GOTO 5
      IF(N.EQ.4.AND.L.LT.3) WRITE(*,*)'4SPD ORBITALS NOT AVAILABLE'
      IF(SCALE.LT.1.E-6)SCALE=1.0
C	Set the origin at the centre of the cube
      XO = (IGH+1)*SCALE
      YO = (IGH+1)*SCALE
      ZO = (IGH+1)*SCALE
C       Build up a grid of points through a slice of the Z plane.
      DO 130 L1 = 1,IG
      DO 130 I = 1,IG
      DO 130 J = 1,IG
      XI=FLOAT(I)*SCALE
      YI=FLOAT(J)*SCALE
      ZI=FLOAT(L1)*SCALE
C	Calculate the distance from the origin
      R=SQRT((XI-XO)**2+(YI-YO)**2+(ZI-ZO)**2)
C 	Make sure division by zero is not attempted!
      IF(R.LT.1.E-6) GOTO 65
      TEMP=SQRT((YI-YO)**2+(ZI-ZO)**2)
C	Define SIN PHI (SP), COS PHI, SIN PSI (SS) COS PSI (CS)
      SP= (TEMP/R)
      CP= (XI-XO)/R
      IF(TEMP.LT.1.E-6)GO TO 70
      SS= (ZI-ZO)/TEMP
      CS= (YI-YO)/TEMP
      GOTO 80
   65 SP=0
      CP=0
      SS=0
      CS=0.
      GOTO 80
   70 SS = 0
      CS = 0.
   80 CONTINUE
C 	Fill grid
C	The Radial part is a function RADXX
C	The angular part is a function ANGXX
      IF(N.EQ.1) G(J,I,L1)=RAD1S(R)*ANGS(SP,CP,SS,CS,PI)
      IF(N.EQ.2.AND.L.EQ.0.AND.M.EQ.0)
     +         G(J,I,L1)=RAD2S(R)*ANGS(SP,CP,SS,CS,PI)
      IF(N.EQ.2.AND.L.EQ.1.AND.M.EQ.0)
     +         G(J,I,L1)=RAD2P(R)*ANGPX(SP,CP,SS,CS,PI)
      IF(N.EQ.2.AND.L.EQ.1.AND.M.EQ.1)
     +         G(J,I,L1)=RAD2P(R)*ANGPY(SP,CP,SS,CS,PI)
      IF(N.EQ.2.AND.L.EQ.1.AND.M.EQ.-1)
     +         G(J,I,L1)=RAD2P(R)*ANGPZ(SP,CP,SS,CS,PI)
      IF(N.EQ.3.AND.L.EQ.0.AND.M.EQ.0)
     +         G(J,I,L1)=RAD3S(R)*ANGS(SP,CP,SS,CS,PI)
      IF(N.EQ.3.AND.L.EQ.1.AND.M.EQ.0)
     +         G(J,I,L1)=RAD3P(R)*ANGPX(SP,CP,SS,CS,PI)
      IF(N.EQ.3.AND.L.EQ.1.AND.M.EQ.1)
     +         G(J,I,L1)=RAD3P(R)*ANGPY(SP,CP,SS,CS,PI)
      IF(N.EQ.3.AND.L.EQ.1.AND.M.EQ.-1)
     +         G(J,I,L1)=RAD3P(R)*ANGPZ(SP,CP,SS,CS,PI)
      IF(N.EQ.3.AND.L.EQ.2.AND.M.EQ.0)
     +         G(J,I,L1)=RAD3D(R)*ANGDZR(SP,CP,SS,CS,PI)
      IF(N.EQ.3.AND.L.EQ.2.AND.M.EQ.1)
     +         G(J,I,L1)=RAD3D(R)*ANGDXZ(SP,CP,SS,CS,PI)
      IF(N.EQ.3.AND.L.EQ.2.AND.M.EQ.-1)
     +         G(J,I,L1)=RAD3D(R)*ANGDYZ(SP,CP,SS,CS,PI)
      IF(N.EQ.3.AND.L.EQ.2.AND.M.EQ.2)
     +         G(J,I,L1)=RAD3D(R)*ANGDXY2(SP,CP,SS,CS,PI)
      IF(N.EQ.3.AND.L.EQ.2.AND.M.EQ.-2)
     +         G(J,I,L1)=RAD3D(R)*ANGDXY(SP,CP,SS,CS,PI)
      IF(N.EQ.4.AND.L.EQ.3.AND.M.EQ.0)
     +   G(J,I,L1)=RAD4F(R)*ANGFZR(SP,CP,SS,CS,PI)
      IF(N.EQ.4.AND.L.EQ.3.AND.M.EQ.1)
     +   G(J,I,L1)=RAD4F(R)*ANGFXZY(SP,CP,SS,CS,PI)
      IF(N.EQ.4.AND.L.EQ.3.AND.M.EQ.-1)
     +   G(J,I,L1)=RAD4F(R)*ANGFYZX(SP,CP,SS,CS,PI)
      IF(N.EQ.4.AND.L.EQ.3.AND.M.EQ.2)
     +   G(J,I,L1)=RAD4F(R)*ANGFZXY(SP,CP,SS,CS,PI)
      IF(N.EQ.4.AND.L.EQ.3.AND.M.EQ.-2)
     +   G(J,I,L1)=RAD4F(R)*ANGFXYZ(SP,CP,SS,CS,PI)
      IF(N.EQ.4.AND.L.EQ.3.AND.M.EQ.3)
     +   G(J,I,L1)=RAD4F(R)*ANGFXR(SP,CP,SS,CS,PI)
      IF(N.EQ.4.AND.L.EQ.3.AND.M.EQ.-3)
     +   G(J,I,L1)=RAD4F(R)*ANGFYR(SP,CP,SS,CS,PI)
  130 CONTINUE
C	Write out the grid to a file
      CALL SCPLOT(IG)
      GO TO 5
  10  STOP
      END
      FUNCTION RAD1S(R)
      RAD1S=2*EXP(-R)
      RETURN
      END
      FUNCTION RAD2S(R)
      RAD2S=SQRT(1./8.)*(2.-R)*EXP(-R/2.)
      RETURN
      END
      FUNCTION RAD2P(R)
      RAD2P=SQRT(1./24.)*R*EXP(-R/2.)
      RETURN
      END
      FUNCTION RAD3S(R)
      RAD3S=(2./81.)*SQRT(1./3.)*(27.-18.*R+ 2.*R*R)*EXP(-R/3.)
      RETURN
      END
      FUNCTION RAD3P(R)
      RAD3P=(1./81.)*SQRT(8./3.)*(6.-R)*R*EXP(-R/3.)
      RETURN
      END
      FUNCTION RAD3D(R)
      RAD3D=(1./81.)*SQRT(8./15.)*R*R*EXP(-R/3.)
      RETURN
      END
      FUNCTION RAD4F(R)
      RAD4F=(1./(768.*SQRT(35.)))*R*R*R*EXP(-R/4.)
      RETURN
      END
      FUNCTION ANGS(SP,CP,SS,CS,PI)
      ANGS= PI*0.5
      RETURN
      END
      FUNCTION ANGPX(SP,CP,SS,CS,PI)
      ANGPX=PI*SQRT(3./4.)*CP
      RETURN
      END
      FUNCTION ANGPY(SP,CP,SS,CS,PI)
      ANGPY=PI*SQRT(3./4.)*SP*CS
      RETURN
      END
      FUNCTION ANGPZ(SP,CP,SS,CS,PI)
      ANGPZ=PI*SQRT(3./4.)*SP*SS
      RETURN
      END
      FUNCTION ANGDZR(SP,CP,SS,CS,PI)
      ANGDZR=PI*SQRT(5./16.)*(3*CP*CP-1)
      RETURN
      END
      FUNCTION ANGDXZ(SP,CP,SS,CS,PI)
      ANGDXZ=PI*SQRT(15./4.)*SP*CP*CS
      RETURN
      END
      FUNCTION ANGDYZ(SP,CP,SS,CS,PI)
      ANGDYZ=PI*(15./4.)*SP*CP*SS
      RETURN
      END
      FUNCTION ANGDXY2(SP,CP,SS,CS,PI)
      ANGDXY2=PI*SQRT(15./16.)*SP*SP*(2*CS*CS-1.)
      RETURN
      END
      FUNCTION ANGDXY(SP,CP,SS,CS,PI)
      ANGDXY=PI*SQRT(16./15.)*SP*SP*(2.*CS*SS)
      RETURN
      END
      FUNCTION ANGFZR(SP,CP,SS,CS,PI)
      ANGFZR=PI*SQRT(7./4.)*(5.*CP*CP*CP-3.*CP)
      RETURN
      END
      FUNCTION ANGFXZY(SP,CP,SS,CS,PI)
      ANGFXZY=PI*SQRT(105./16.)*SP*CS*(CP*CP-SP*SP*SS*SS)
      RETURN
      END
      FUNCTION ANGFYZX(SP,CP,SS,CS,PI)
      ANGFYZX=PI*SQRT(105./16.)*SP*SS*(CP*CP-SP*SP*CS*CS)
      RETURN
      END
      FUNCTION ANGFZXY(SP,CP,SS,CS,PI)
      ANGFZXY=PI*SQRT(105./16.)*SP*SP*CP*(2.*CS*CS-1.)
      RETURN
      END
      FUNCTION ANGFXYZ(SP,CP,SS,CS,PI)
      ANGFXYZ=PI*SQRT(105./16.)*SP*SP*CP*(2.*SS*CS)
      RETURN
      END
      FUNCTION ANGFXR(SP,CP,SS,CS,PI)
      ANGFXR=PI*SQRT(7./16.)*SP*CS*(5.*SP*SP*CS*CS-3.)
      RETURN
      END
      FUNCTION ANGFYR(SP,CP,SS,CS,PI)
      ANGFYR=PI*SQRT(7./16.)*SP*SS*(5.*SP*SP*SS*SS-3.)
      RETURN
      END
      SUBROUTINE SCPLOT(IG)
      COMMON /GRID/ G(25,25,25)
      OPEN(UNIT=20,FILE='DATA.DAT',FORM='FORMATTED',STATUS='UNKNOWN')
C	Normalise so that the maximum value is 1.0
      AMAX=0.0
C      GOTO 4
      DO 1 J=1,IG
      DO 1 K=1,IG 
      DO 1 L=1,IG
    1  IF(G(K,J,L).GT.AMAX) AMAX= G(K,J,L)
      DO 2 J=1,IG
      DO 2 K=1,IG
      DO 2 L=1,IG
    2 G(K,J,L) = G(K,J,L)/AMAX
    4 WRITE(20,3) (((G(K,J,L1),K=1,IG),J=1,IG),L1=1,IG)
      RETURN
    3 FORMAT(F10.5)
      END
