Télécharger elpdbe.eso

Retour à la liste

Numérotation des lignes :

  1. C ELPDBE SOURCE KK2000 14/04/09 21:15:24 8027
  2. SUBROUTINE ELPDBE(CZ,N,NP,CH10,CH11,XM)
  3. C
  4. C =====================================================================
  5. C
  6. C CALCUL DES FONCTIONS DE HANKEL PAR SERIE
  7. C
  8. C
  9. C =====================================================================
  10. C
  11. IMPLICIT INTEGER(I-N)
  12. IMPLICIT REAL*8(A-B,D-H,O-Z)
  13. IMPLICIT COMPLEX*16(C)
  14. C Include contenant quelques constantes dont XPI :
  15. -INC CCREEL
  16. C
  17. C
  18. DIMENSION CZ(*)
  19. DIMENSION CH10(*)
  20. DIMENSION CH11(*)
  21. C
  22. CI = (0d0,1d0)
  23. C Constante d'Euler
  24. CGAM = .5772156649015328606d0
  25. C
  26. C
  27. XERR = 0D0
  28. DO 1000 I=1,N
  29.  
  30. CJ0 = 0d0
  31. CJ1 = 0d0
  32.  
  33. C
  34.  
  35. IF ( ABS ( CZ(I)) .LT. XM) THEN
  36.  
  37. CT = CZ(I) **CMPLX(2) / CMPLX(4.d0)
  38. CAL0 = CMPLX(1d0)
  39. CAL1 = CZ(I) /CMPLX(2.d0)
  40.  
  41. DO 1100 IK=1,NP
  42. K= IK - 1
  43. CJ0 = CJ0 + CAL0
  44. CAL0 = CMPLX(-1.d0)*CT* CAL0/ ((K+CMPLX(1))**CMPLX(2))
  45. CJ1 = CJ1 + CAL1
  46. CAL1 = CMPLX(-1.d0)*CT* CAL1/ ((K+CMPLX(1))*(K+CMPLX(2)))
  47. 1100 CONTINUE
  48. C
  49. XERR = MAX (XERR,ABS(CAL0/CJ0))
  50. XERR = MAX (XERR,ABS(CAL1/CJ1))
  51.  
  52. CT = CZ(I) **CMPLX(2) / CMPLX(4.d0)
  53. CY0 = LOG ( CZ(I) / CMPLX(2.d0))
  54. CY0 = CJ0 * ( CY0 + CGAM)
  55. CS0 = CMPLX(1d0)
  56. CS0P = CMPLX(0d0)
  57. CBET0 = CT
  58.  
  59. CY1 = LOG ( CZ(I) / CMPLX(2.d0))
  60. CY1 = CJ1 * CY1 - CMPLX(1.d0)/CZ(I)
  61. CS1 = CMPLX(-2.d0)*CGAM + CMPLX(1.d0)
  62. CS1P = CMPLX(0d0)
  63. CBET1 = CMPLX(-1.d0)*CZ(I) *CS1/CMPLX(4.d0)
  64.  
  65.  
  66. DO 1110 IK=1,NP
  67.  
  68. CY0 = CY0 + CBET0
  69. CS0P = CS0 + CMPLX(1.d0)/( IK + 1)
  70. CBET0= CMPLX(-1.d0)*CT* CBET0/ ( (IK+CMPLX(1)) *
  71. & (IK + CMPLX(1))) * ( CS0P/CS0)
  72. CS0 = CS0P
  73.  
  74. K= IK - 1
  75. CY1 = CY1 + CBET1
  76. CS1P = CS1 +( CMPLX(1.D0)/( K + CMPLX(1))) +
  77. & ( CMPLX(1.D0)/( K + CMPLX(2)))
  78. CBET1= CMPLX(-1.d0)*CT* CBET1/
  79. & ( ( K+CMPLX(1)) *( K+ CMPLX(2))) * ( CS1P/CS1)
  80. CS1 = CS1P
  81.  
  82. 1110 CONTINUE
  83. XERR = MAX (XERR,ABS(CBET0/CY0))
  84. XERR = MAX (XERR,ABS(CBET1/CY1))
  85. CY0 = CY0 *CMPLX(2.d0) / XPI
  86. CY1 = CY1 *CMPLX(2.d0) / XPI
  87. CH10(I) = CJ0 + CI*CY0
  88. CH11(I) = CJ1 + CI*CY1
  89. ELSE
  90. C8 = CZ(I) *CMPLX(8d0)
  91. CAL0 = CMPLX(1d0)
  92. CAL1 = CMPLX(1d0)
  93. CP0= CAL0
  94. CP1= CAL1
  95. CQ0=CMPLX(0.D0)
  96. CQ1=CMPLX(0.D0)
  97.  
  98. DO 1200 K=1,NP
  99. I1= 2*K -1
  100. I2= 2*K
  101. CBET0= CMPLX(-1.D0)*CAL0*(CMPLX(2)*I1 -
  102. & CMPLX(1))**CMPLX(2)/ (I1 *C8)
  103. CBET1= CAL1*(CMPLX(4) - (CMPLX(2)*I1 -
  104. & CMPLX(1))**CMPLX(2))/ (I1 *C8)
  105. CAL0= CBET0*(CMPLX(2)*I2 - CMPLX(1))**CMPLX(2)/ (I2 *C8)
  106. CAL1= CMPLX(-1.D0)*CBET1*(CMPLX(4) -
  107. & (CMPLX(2)*I2 - CMPLX(1))**CMPLX(2))/ (I2 *C8)
  108. CP0 = CP0 + CAL0
  109. CQ0 = CQ0 + CBET0
  110. CP1 = CP1 + CAL1
  111. CQ1 = CQ1 + CBET1
  112. 1200 CONTINUE
  113. XERR = MAX (XERR,ABS(CBET0/CQ0))
  114. XERR = MAX (XERR,ABS(CBET1/CQ1))
  115. XERR = MAX (XERR,ABS(CAL0/CP0))
  116. XERR = MAX (XERR,ABS(CAL1/CP1))
  117.  
  118. CCOEF=(CMPLX(2) / ( XPI * CZ(I))) ** CMPLX(.5D0)
  119. CKI0= CZ(I) - XPI /CMPLX(4.d0)
  120. CH10(I)=CCOEF* (CP0 + CI * CQ0) * EXP ( CI * CKI0)
  121. CKI1= CZ(I) - (XPI*CMPLX(3.d0) /CMPLX(4.d0))
  122. CH11(I)=CCOEF* (CP1 + CI * CQ1) * EXP ( CI * CKI1)
  123. ENDIF
  124.  
  125. 1000 CONTINUE
  126.  
  127. C
  128. RETURN
  129. END
  130.  
  131.  
  132.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales