1 : C ELPDBE SOURCE CHAT 05/01/12 23:36:17 5004 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 : CXPI = XPI 24 : C Constante d'Euler 25 : CGAM = .5772156649015328606d0 26 : C 27 : C 28 : XERR = 0D0 29 : DO 1000 I=1,N 30 : 31 : CJ0 = 0d0 32 : CJ1 = 0d0 33 : 34 : C 35 : 36 : IF ( ABS ( CZ(I)) .LT. XM) THEN 37 : 38 : CT = CZ(I) **2 / 4.d0 39 : CAL0 = 1d0 40 : CAL1 = CZ(I) /2.d0 41 : 42 : DO 1100 IK=1,NP 43 : K= IK - 1 44 : CJ0 = CJ0 + CAL0 45 : CAL0 = -1.d0*CT* CAL0/ ( ( K+1) **2) 46 : CJ1 = CJ1 + CAL1 47 : CAL1 = -1.d0*CT* CAL1/ ( ( K+1) *(K + 2)) 48 : 1100 CONTINUE 49 : C 50 : XERR = MAX (XERR,ABS(CAL0/CJ0)) 51 : XERR = MAX (XERR,ABS(CAL1/CJ1)) 52 : 53 : CT = CZ(I) **2 / 4.d0 54 : CY0 = LOG ( CZ(I) / 2.d0) 55 : CY0 = CJ0 * ( CY0 + CGAM) 56 : CS0 = 1d0 57 : CS0P = 0d0 58 : CBET0 = CT 59 : 60 : CY1 = LOG ( CZ(I) / 2.d0) 61 : CY1 = CJ1 * CY1 - 1.d0/CZ(I) 62 : CS1 = -2.d0*CGAM + 1.d0 63 : CS1P = 0d0 64 : CBET1 =-1.d0*CZ(I) *CS1/4.d0 65 : 66 : 67 : DO 1110 IK=1,NP 68 : 69 : CY0 = CY0 + CBET0 70 : CS0P = CS0 + 1.d0/( IK + 1) 71 : CBET0= -1.d0*CT* CBET0/ ( (IK+1) *(IK + 1)) * ( CS0P/CS0) 72 : CS0 = CS0P 73 : 74 : K= IK - 1 75 : CY1 = CY1 + CBET1 76 : CS1P = CS1 +( 1.D0/( K + 1)) + (1.D0/( K + 2 )) 77 : CBET1= -1.d0*CT* CBET1/ ( ( K+1) *( K+ 2)) * ( CS1P/CS1) 78 : CS1 = CS1P 79 : 80 : 1110 CONTINUE 81 : XERR = MAX (XERR,ABS(CBET0/CY0)) 82 : XERR = MAX (XERR,ABS(CBET1/CY1)) 83 : CY0 = CY0 *2.d0 / CXPI 84 : CY1 = CY1 *2.D0 / CXPI 85 : CH10(I) = CJ0 + CI*CY0 86 : CH11(I) = CJ1 + CI*CY1 87 : ELSE 88 : C8 = CZ(I) *8d0 89 : CAL0 = 1d0 90 : CAL1 = 1d0 91 : CP0= CAL0 92 : CP1= CAL1 93 : CQ0=0.D0 94 : CQ1=0.D0 95 : 96 : DO 1200 K=1,NP 97 : I1= 2*K -1 98 : I2= 2*K 99 : CBET0= -1.D0*CAL0*(2*I1 - 1)**2/ (I1 *C8) 100 : CBET1= CAL1*(4 - (2*I1 - 1)**2)/ (I1 *C8) 101 : CAL0= CBET0*(2*I2 - 1)**2/ (I2 *C8) 102 : CAL1= -1.D0*CBET1*(4 - (2*I2 - 1)**2)/ (I2 *C8) 103 : CP0 = CP0 + CAL0 104 : CQ0 = CQ0 + CBET0 105 : CP1 = CP1 + CAL1 106 : CQ1 = CQ1 + CBET1 107 : 1200 CONTINUE 108 : XERR = MAX (XERR,ABS(CBET0/CQ0)) 109 : XERR = MAX (XERR,ABS(CBET1/CQ1)) 110 : XERR = MAX (XERR,ABS(CAL0/CP0)) 111 : XERR = MAX (XERR,ABS(CAL1/CP1)) 112 : 113 : CCOEF=(2 / ( CXPI * CZ(I))) ** .5D0 114 : CKI0= CZ(I) - CXPI /4.d0 115 : CH10(I)=CCOEF* (CP0 + CI * CQ0) * EXP ( CI * CKI0) 116 : CKI1= CZ(I) - (CXPI*3.d0 /4.d0) 117 : CH11(I)=CCOEF* (CP1 + CI * CQ1) * EXP ( CI * CKI1) 118 : ENDIF 119 : 120 : 1000 CONTINUE 121 : 122 : C 123 : RETURN 124 : END 125 : 126 : 127 : 128 : 129 : 130 : 131 : 132 : 133 : 134 : 135 :
© Cast3M 2003 - All rights reserved.
Disclaimer