PROGRAM ORB

*************************************************************************
*   An Orbital calculator.      
*************************************************************************
  
      CHARACTER *2 ANS
      CHARACTER *1 X
      DIMENSION RDF(50)

     
      PRINT *, ' Welcome to orbital calculator.'
1     PRINT *, ' For which orbital do you want to do calculations?'
      

      READ (*,800) ANS
800   FORMAT (A2)
      
      IF (ANS.NE.'1s'.AND.ANS.NE.'2s'.AND.ANS.NE.'2p') GOTO 1

      PRINT *, ' Right, you want info on the ',ANS,' orbital.'
      PRINT *, ' Here it is.'
      PRINT *

      PI=3.1415926

*************************************************************************
* Since we are going to store the RDFs in an array, we need to have
*  a variable which will count for us as we place successive values
*  in the array. This we call ICOUNT, and it must first be initialised
*  (set to zero), then incremented on each pass through the loop.
*************************************************************************


      ICOUNT=0
           
      DO 100 R=0, 2.0, 0.05
      ICOUNT=ICOUNT+1

*************************************************************************
*   We shall use FUNCTIONS for each type of orbital
*************************************************************************

      IF (ANS.EQ.'1s') PSI=ONES(R)
      IF (ANS.EQ.'2s') PSI=TWOS(R)
      IF (ANS.EQ.'2p') PSI=TWOP(R)

      RDF(ICOUNT)=(PSI*PSI)*4.0*PI*R*R
100   CONTINUE

*************************************************************************
*  To scale the array, we call a SUBROUTINE. This is much like a
*  FUNCTION, but does not return a value to the line from which it
*  was called.
*************************************************************************

      CALL SCALE(RDF,ICOUNT)

      ICOUNT=0
      
      DO 101 R=0.0, 4.0, 0.1
      ICOUNT=ICOUNT+1

*************************************************************************       
* Instead of printing out the numercial value of the rdf, we shall
*  print some Xs on the screen, in a row proportional to the size
*  of the rdf. To do this, we need to store the character we shall
*  print out in the variable X
*************************************************************************                
      X='X'   

*************************************************************************
* Now we just print the variable X a number of times equal to rdf      
* The WRITE command rather than PRINT allows us more control 
* of how the data is output. The FORMAT line does the formating 
* of the data. 
*************************************************************************

      WRITE(*,801) R, (X, I=1,50*RDF(ICOUNT))
801   FORMAT(F5.2,5x,100A1)
101   CONTINUE

      STOP
      END


*************************************************************************
*   The 1s calculator.      
*************************************************************************   

      FUNCTION ONES(R)
      
      ONES=EXP(-R/2.0)
      
      RETURN 
      END
      

*************************************************************************
*   The 2s calculator.      
*************************************************************************

      FUNCTION TWOS(R)
      
      FAC=(1.0/SQRT(32.0))*(2.0-R)
      TWOS=FAC*EXP(-R/2.0)

      RETURN
      END

*************************************************************************
*   The 2p calculator.      
*************************************************************************
      
      FUNCTION TWOP(R)
      
      FAC=(1.0/SQRT(972.0))*(6-6*R+R*R)
      TWOP=FAC*EXP(-R/2.0)
      
      RETURN
      END

*************************************************************************
*   The orbital scaler.      
*************************************************************************

      SUBROUTINE SCALE(RDF,ICOUNT)
      
      DIMENSION RDF(50)
      
      BIGGEST=0.0

************************************    
*  First find the largest value in the array....
************************************

       DO 100 I=1,ICOUNT
       IF (BIGGEST.LT.RDF(I)) BIGGEST=RDF(I)
100   CONTINUE

************************************      
*  ...and then divide all values 
*  by that value so that the greatest
*  value in the array is now 1.0
************************************

      DO 101 I=1,ICOUNT
      RDF(I)=RDF(I)/BIGGEST
101   CONTINUE
      
      RETURN
      END