PROGRAM SPHASE
C     DIAGONAL CORRECTION VIA PHASE INTEGRAL                                    
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
      IMPLICIT INTEGER*4 (I-N)                                                  
      PARAMETER (NT0=12,NT3=3*NT0,NM0=4000)
      DIMENSION EIT(2)
      COMMON/INTEGRAL/G,A,C,D,O(2),IT1,IT2,IT3
c
C     1 a.u.= 2.1947E5 cm^-1                                                    
      CM=219470.0d0
C     PI=3.142592653589                                                         
      PI=REAL(4)*ATAN(REAL(1))                                                  
      DLIMIT=0.01/CM                                                            
      AITGLIM=10000000.0-1.0                                                    
      WRITE(*,*)' Phase integrals'                                              
      WRITE(*,*)' E=G*q+(A/2)*q^2+(C/6)*q^3+(D/24)*q^4'
      WRITE(*,*)' Input G, A, C and D in cm-1:'
      READ(*,*)G,A,C,D
      IF(A.LT.0.0)THEN                                                          
       WRITE(*,*)' Absolute value of frequency used for definition of q'        
      ENDIF                                           
c     
c     get atomic units  
      G=G/CM
      A=A/CM
      C=C/CM
      D=D/CM
      ALFAC=SQRT(ABS(A)*0.5)*PI                                                 
      DO 6 NI=0,1                                                               
      ALIMIT=(REAL(NI)+0.5)*ALFAC                                               
C     int (O1..O2)SQRT(E-V(q))dq = SQRT(Hbar*Wharm/2)*(ni+1/2)*pi               
      DE=1.0/CM                                                                 
      IT1=0                                                                     
      EN1=ABS(A)*(REAL(NI)+0.40)                                                
      EN2=EN1+DE                                                                
7     AI1=AITG(EN1,QI,QA)
      AI2=AITG(EN2,QI,QA)
      IF(AI1.GT.AITGLIM.OR.AI2.GT.AITGLIM)THEN                                  
       WRITE(*,*)'Integral too big ...'
       GOTO 888
      ENDIF                                                                     
      IF(ABS(AI2-AI1).LT.0.0000001)THEN                                         
       WRITE(*,*)'Difference too small'
       GOTO 6
      ENDIF                                                                     
      EN3=EN2+(EN2-EN1)*(ALIMIT-AI2)/(AI2-AI1)                                  
      IT1=IT1+1                                                                 
      IF(IT1.GT.100)THEN                                                        
       WRITE(*,*)' Stable state has not been found after 100 iterations'        
       GOTO 888
      ENDIF                                                                     
      IF(ABS(EN3-EN2).LT.DLIMIT)THEN
       IF(NI.EQ.0)THEN
        Q0I=QI
        Q0A=QA
       ELSE
        Q1I=QI
        Q1A=QA
       ENDIF  
       GOTO 6                                          
      ENDIF
      EN1=EN2                                                                   
      EN2=EN3                                                                   
      GOTO 7                                                                    
6     EIT(NI+1)=EN3                                                             
888   WRITE(6,6006)EIT(2)*CM,Q1I,Q1A,EIT(1)*CM,Q0I,Q0A,
     1 EIT(2)*CM-EIT(1)*CM,2.0D0*EIT(1)*CM,IT3
6006  FORMAT('Upper level ',F15.3,'   q = ',F7.3,' - ',F6.3,/,
     1       'Lower level ',F15.3,'   q = ',F7.3,' - ',F6.3,/,
     2       'Difference  ',F15.3,/,
     3       'Lower * 2   ',F15.3,/,/,
     4       'Number of iterations :',I5,/)
      STOP
      END                                                                       
c
      FUNCTION AITG(E,QI,QA)
      IMPLICIT REAL*8 (A-H,O-Z)                                                 
      IMPLICIT INTEGER*4 (I-N)                                                  
      COMMON/INTEGRAL/G,A,C,D,O(2),IT1,IT2,IT3
C     int(o1,o2)sqrt(E-V(q))dq                                                  
      CM=219470.0                                                               
      O(1)=0.0                                                                  
      O(2)=0.0                                                                  
      OLIMIT=0.000001                                                           
      ELIMIT=0.001/CM                                                           
      SLIMIT=0.000001                                                           
      IT2=0                                                                     
      DO 1 I=1,2                                                                
      DEL=0.025*REAL(I*2-3)                                                     
      ON=0.0                                                                    
2     ON=ON+DEL                                                                 
      IT2=IT2+1                                                                 
      IF(IT2.GT.1000)THEN                                                       
       WRITE(*,*)I,' Cross-section of energy surface not found '                
       VN=ON*(G+ON*(0.5*A+ON*(C/6.0+D*ON/24.0)))
       WRITE(*,3000)A*CM,C*CM,D*CM,ON,DEL,E*CM,VN*CM                            
3000  FORMAT(' A = ',F10.1,' B = ',F10.1,' C = ',F10.1,/,                       
     1       ' ON, DEL, E, VN: ',4F15.6)                                        
       AITG=10000000.0                                                          
       RETURN                                                                   
      ENDIF                                                                     
      VN=ON*(G+ON*(0.5*A+ON*(C/6.0+D*ON/24.0)))
      IF(VN.GE.E)THEN                                                           
       DEL=0.5*DEL                                                              
       ON=ON-2.0*DEL                                                            
       IF(ABS(2.0*DEL).LT.OLIMIT.AND.ABS(VN-E).LT.ELIMIT)GOTO 1                 
       GOTO 2                                                                   
      ENDIF                                                                     
      GOTO 2                                                                    
1     O(I)=ON                                                                   
      IF(O(1).LT.-100.OR.O(2).GT.100)THEN                                       
       WRITE(*,*)' No solution for this energy'                                 
       AITG=10000000.0                                                          
       RETURN                                                                   
      ENDIF                                                                     
      NS=1000                                                                   
      SUMO=0.0                                                                  
      IT3=0                                                                     
4     NS=NS+100                                                                 
      IT3=IT3+1                                                                 
      IF(IT3.GT.100)THEN                                                        
       WRITE(*,*)'Impossible to do numerical integration.'                      
       AITG=10000000.0                                                          
       RETURN                                                                   
      ENDIF                                                                     
      QI=O(1)
      QA=O(2)
      H=(O(2)-O(1))/REAL(NS)                                                    
      SUM=0.0                                                                   
      X=O(1)                                                                    
      DO 3 I=1,NS-1                                                             
      X=X+H                                                                     
      COEF=2.0+2.0*MOD(I+1,2)                                                   
3     SUM=SUM+COEF*SQRT(E-X*(G+X*(A/2.0+X*(C/6.0+X*D/24.0))))
      SUM=SUM*H/3.0                                                             
      IF(ABS(SUM-SUMO).GT.SLIMIT)THEN                                           
       SUMO=SUM                                                                 
       GOTO 4                                                                   
      ENDIF                                                                     
      AITG=SUM                                                                  
      RETURN                                                                    
      END