Search Site | Contact Details | FAQ

ADAS Subroutine r8fmon1

C
       REAL*8 FUNCTION R8FMON1(E1,E2,L)                                 
       IMPLICIT REAL*8(A-H,O-Z)                                         
C-----------------------------------------------------------------------
C
C  ****************** FORTRAN77 FUNCTION: R8FMON1 **********************
C
C PURPOSE: CALCULATES THE MONOPOLE INTEGRAL |<E1,L|1/r>|E2,L>|^2  
C          
C
C NOTE: CREATED BY  ALAN BURGESS FOR USE IN THE DIPOLE INTEGRAL 
C       I(KAPPA1,L1,KAPPA2,L2,1) EVALUATION AS DEFINED IN PHIL. 
C       TRANS. ROY. SOC. A226,255,1970, WHERE E1=KAPPA1**2 AND 
C       E2=KAPPA2**2. APPLIES TO POSITIVE ELECTRON ENERGIES, .  
c       THAT IS THE FREE-FREE CASE.                                              
C
C CALLING PROGRAMS: R8MON1
C
C INPUT:  (R*8)  E1      = KAPPA1**2 WHERE KAPPA1 IS SCALED INITIAL
C                          ELECTRON WAVE NUMBER 
C INPUT:  (I*4)  L       = ORBITAL ANGULAR OMENTUM OF INITIAL ELECTRON 
C INPUT:  (R*8)  E2      = KAPPA2**2 WHERE KAPPA2 IS SCALED INITIAL
C                          ELECTRON WAVE NUMBER 
C
C OUTPUT: (R*8)  R8MON1 = |<E1,L|1/r>|E2,L>|^2 
C
C ROUTINES:
C          ROUTINE    SOURCE   BRIEF DESCRIPTION
C          -------------------------------------------------------------
C          ARGAM      ADAS     CALCULATES ARGGAMMA(L+1+I*A)
C
C UNIX-IDL PORT:
C
C VERSION: 1.1                          DATE: 17-04-07
C MODIFIED: HUGH SUMMERS
C               - FIRST FULLY COMMENTED RELEASE
C
C-----------------------------------------------------------------------
       IF(E1+E2-1.0D-40)28,28,29                                        
   28  R8FMON1=1.0D50                                                   
       RETURN                                                           
   29  CONTINUE                                                         
       VMAX=200.0D0                                                     
       X1=DSQRT(E1)                                                     
       X2=DSQRT(E2)                                                     
       X3=X1+X2                                                         
       X4=X3*X3                                                         
       X5=X1*X2                                                         
       X6=X2-X1                                                         
       X7=4.0D0/X4                                                      
       PI=3.141592653589793D0                                           
       IF(E1-E2)1,1,2                                                   
    1  ETA=1.0D0/X2                                                     
       GO TO 3                                                          
    2  ETA=1.0D0/X1                                                     
    3  G=0.5D0*PI*DEXP(-PI*ETA)                                         
       A1=1.0D0                                                         
       A2=1.0D0                                                         
       MG=0                                                             
       MA1=0                                                            
       MA2=0                                                            
       M=-1                                                             
    4  M=M+1                                                            
       EM=M                                                             
       T=EM+EM+1.0D0                                                    
       G=G*X7/(T*(T+1.0D0))                                             
       EMM=EM*EM                                                        
       A1=A1*(1.0D0+EMM*E1)                                             
       A2=A2*(1.0D0+EMM*E2)                                             
   30  IF(G-0.015625D0) 31,32,32                                        
   31  G=64.0D0*G                                                       
       MG=MG-1                                                          
       GO TO 30                                                         
   32  IF(G-64.0D0) 34,34,33                                            
   33  G=0.015625D0*G                                                   
       MG=MG+1                                                          
       GO TO 32                                                         
   34  IF(A1-64.0D0) 36,36,35                                           
   35  A1=0.015625D0*A1                                                 
       MA1=MA1+1                                                        
       GO TO 34                                                         
   36  IF(A2-64.0D0) 38,38,37                                           
   37  A2=0.015625D0*A2                                                 
       MA2=MA2+1                                                        
       GO TO 36                                                         
   38  CONTINUE                                                         
       IF(M-L)4,5,5                                                     
    5  G=G*(T+1.0D0)                                                    
       IF(X1-300.0D0)7,6,6                                              
    6  B=PI/X1                                                          
       A1=1.5D0*A1/(B*(3.0D0-B*(3.0D0-B*(2.0D0-B))))                    
       GO TO 9                                                          
    7  IF(X1-0.2D0)9,9,8                                                
    8  B=-PI/X1                                                         
       A1=A1/(1.0D0-DEXP(B+B))                                          
    9  IF(X2-300.0D0)11,10,10                                           
   10  B=PI/X2                                                          
       A2=1.5D0*A2/(B*(3.0D0-B*(3.0D0-B*(2.0D0-B))))                    
       GO TO 13                                                         
   11  IF(X2-0.2)13,13,12                                               
   12  B=-PI/X2                                                         
       A2=A2/(1.0D0-DEXP(B+B))                                          
   13  G=G*DSQRT(A1*A2)*(8.0D0)**(MG+MG+MA1+MA2)                        
       S0=1.0D0                                                         
       S1=0.0D0                                                         
       U=L                                                              
       V=0.0D0                                                          
       W=U+U+1.0D0                                                      
       T0=1.0D0                                                         
       T1=0.0D0                                                         
   14  U=U+1.0D0                                                        
       V=V+1.0D0                                                        
       W=W+1.0D0                                                        
       IF(V-VMAX)21,21,20                                               
   20  R8FMON1=0.0D0                                                    
       RETURN                                                           
   21  CONTINUE                                                         
       U0=U*U*X5+1.0D0                                                  
       U1=U*X6                                                          
       T=T0*U0-T1*U1                                                    
       T1=T0*U1+T1*U0                                                   
       T0=T                                                             
       T=X7/(V*W)                                                       
       T0=T*T0                                                          
       T1=T*T1                                                          
       S0=S0+T0                                                         
       S1=S1+T1                                                         
       S=S0*S0+S1*S1                                                    
       T=T0*T0+T1*T1                                                    
       IF(S-1.0D24*T)14,15,15                                           
   15  R8FMON1=G*DSQRT(S)                                               
       IV=V                                                             
       RETURN                                                           
      END                                                               
      INTEGER             L
      REAL*8              E1,          E2
© Copyright 1995-2018 The ADAS Project
Comments and questions to: adas-at-adas.ac.uk