Sources Esope | Cast3M

elpdbe
C ELPDBE    SOURCE    KK2000    14/04/09    21:15:24     8027
      SUBROUTINE ELPDBE(CZ,N,NP,CH10,CH11,XM)
C
C =====================================================================
C
C     CALCUL DES FONCTIONS DE HANKEL PAR SERIE
C
C
C =====================================================================
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-B,D-H,O-Z)
      IMPLICIT COMPLEX*16(C)
C    Include contenant quelques constantes dont XPI :
-INC CCREEL
C
C
      DIMENSION CZ(*)
      DIMENSION CH10(*)
      DIMENSION CH11(*)
C
      CI = (0d0,1d0)
C     Constante d'Euler
      CGAM = .5772156649015328606d0
C
C
      XERR = 0D0
      DO  1000  I=1,N
 
         CJ0 = 0d0
         CJ1 = 0d0
 
C
 
         IF ( ABS ( CZ(I)) .LT. XM) THEN
 
            CT = CZ(I) **CMPLX(2) / CMPLX(4.d0)
            CAL0 = CMPLX(1d0)
            CAL1 = CZ(I) /CMPLX(2.d0)
 
            DO  1100  IK=1,NP
               K= IK - 1
               CJ0 = CJ0 + CAL0
               CAL0 = CMPLX(-1.d0)*CT* CAL0/ ((K+CMPLX(1))**CMPLX(2))
               CJ1    = CJ1    + CAL1
               CAL1 = CMPLX(-1.d0)*CT* CAL1/ ((K+CMPLX(1))*(K+CMPLX(2)))
 1100       CONTINUE
C
            XERR = MAX (XERR,ABS(CAL0/CJ0))
            XERR = MAX (XERR,ABS(CAL1/CJ1))
 
            CT = CZ(I) **CMPLX(2) / CMPLX(4.d0)
            CY0 = LOG ( CZ(I) / CMPLX(2.d0))
            CY0 = CJ0  * ( CY0  + CGAM)
            CS0   = CMPLX(1d0)
            CS0P  = CMPLX(0d0)
            CBET0 = CT
 
            CY1 = LOG ( CZ(I) / CMPLX(2.d0))
            CY1 = CJ1 *  CY1  - CMPLX(1.d0)/CZ(I)
            CS1   = CMPLX(-2.d0)*CGAM + CMPLX(1.d0)
            CS1P  = CMPLX(0d0)
            CBET1 = CMPLX(-1.d0)*CZ(I) *CS1/CMPLX(4.d0)
 
 
            DO  1110  IK=1,NP
 
               CY0 = CY0 + CBET0
               CS0P =  CS0 + CMPLX(1.d0)/( IK + 1)
               CBET0= CMPLX(-1.d0)*CT* CBET0/ ( (IK+CMPLX(1)) *
     &                 (IK + CMPLX(1))) * ( CS0P/CS0)
               CS0  =  CS0P
 
               K= IK - 1
               CY1 = CY1 + CBET1
               CS1P =  CS1 +( CMPLX(1.D0)/( K + CMPLX(1))) +
     &                      ( CMPLX(1.D0)/( K + CMPLX(2)))
               CBET1= CMPLX(-1.d0)*CT* CBET1/
     &                    ( ( K+CMPLX(1)) *( K+ CMPLX(2))) * ( CS1P/CS1)
               CS1  =  CS1P
 
 1110       CONTINUE
            XERR = MAX (XERR,ABS(CBET0/CY0))
            XERR = MAX (XERR,ABS(CBET1/CY1))
            CY0 = CY0 *CMPLX(2.d0) / XPI
            CY1 = CY1 *CMPLX(2.d0) / XPI
            CH10(I) = CJ0 + CI*CY0
            CH11(I) = CJ1 + CI*CY1
         ELSE
            C8 = CZ(I) *CMPLX(8d0)
            CAL0 = CMPLX(1d0)
            CAL1 = CMPLX(1d0)
            CP0= CAL0
            CP1= CAL1
            CQ0=CMPLX(0.D0)
            CQ1=CMPLX(0.D0)
 
            DO  1200  K=1,NP
               I1= 2*K -1
               I2= 2*K
               CBET0= CMPLX(-1.D0)*CAL0*(CMPLX(2)*I1 -
     &                CMPLX(1))**CMPLX(2)/ (I1 *C8)
               CBET1=     CAL1*(CMPLX(4) - (CMPLX(2)*I1 -
     &                          CMPLX(1))**CMPLX(2))/ (I1 *C8)
               CAL0=  CBET0*(CMPLX(2)*I2 - CMPLX(1))**CMPLX(2)/ (I2 *C8)
               CAL1= CMPLX(-1.D0)*CBET1*(CMPLX(4) -
     &              (CMPLX(2)*I2 - CMPLX(1))**CMPLX(2))/ (I2 *C8)
               CP0    = CP0 +  CAL0
               CQ0    = CQ0 + CBET0
               CP1    = CP1 +  CAL1
               CQ1    = CQ1 + CBET1
 1200       CONTINUE
            XERR = MAX (XERR,ABS(CBET0/CQ0))
            XERR = MAX (XERR,ABS(CBET1/CQ1))
            XERR = MAX (XERR,ABS(CAL0/CP0))
            XERR = MAX (XERR,ABS(CAL1/CP1))
 
            CCOEF=(CMPLX(2) / ( XPI * CZ(I))) ** CMPLX(.5D0)
            CKI0= CZ(I) - XPI /CMPLX(4.d0)
            CH10(I)=CCOEF* (CP0 + CI * CQ0) * EXP ( CI * CKI0)
            CKI1= CZ(I) - (XPI*CMPLX(3.d0) /CMPLX(4.d0))
            CH11(I)=CCOEF* (CP1 + CI * CQ1) * EXP ( CI * CKI1)
         ENDIF
 
 1000 CONTINUE
 
C
      RETURN
      END