C FOBEY0 SOURCE CHAT 05/01/13 00:05:37 5004 SUBROUTINE FOBEY0(XR1,NXR1,XR2) IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C Include contenant quelques constantes dont XPI : -INC CCREEL C C C XR1( NXR1) : TABLEAU DE REELS TOUS > 0 C XR2( NXR1) : TABLEAU DE REELS C C XR2(I) = Y0 ( XR1(I) ) C Y0 FONCTION DE BESSEL DE DEUXIEME ESPECE D'ORDRE 0 C C APPROXIMATION POLYNOMIALE PAR SECTEUR C REFERENCE : ABRAMOWITZ HANDBOOK OF MATHEMATICAL FONCTIONS C PRECISION E = 1.D-8 C DIMENSION XR1(NXR1) DIMENSION XR2(NXR1) DIMENSION XB(1) DIMENSION XBB(1) C A0= .36746691D0 A2= .60559366D0 A4= -.74350384D0 A6= .25300117D0 A8= -.04261214D0 A10= .00427916D0 A12=-.00024846D0 C B0= .79788456D0 B1= -.00000077D0 B2= -.00552740D0 B3= -.00009512D0 B4= .00137237D0 B5= -.00072805D0 B6= .00014476D0 C C0= -.78539816D0 C1= -.04166397D0 C2= -.00003954D0 C3= .00262573D0 C4= -.00054125D0 C5= -.00029333D0 C6= .00013558D0 C DO 100 I=1,NXR1 IF (XR1(I) .LE.3D0) THEN XB(1)= XR1(I) CALL FOBEJ0(XB,1,XBB) A00 = (2.D0 / XPI)* (LOG( XB(1) / 2.D0))* XBB(1) Y2= (XR1(I) / 3D0)**2 XR2(I)= A00+A0+Y2*(A2+Y2*(A4+Y2*(A6+Y2*(A8+Y2*(A10+Y2*A12))) $ )) ELSE Y= 3D0/XR1(I) F0= B0+Y*(B1+Y*(B2+Y*(B3+Y*(B4+Y*(B5+Y*B6))))) T0 = XR1(I)+C0+Y*(C1+Y*(C2+Y*(C3+Y*(C4+Y*(C5+Y*C6))))) XR2(I) = (1.D0 / ( XR1(I) ** .5D0))* F0 * SIN (T0) ENDIF 100 CONTINUE C RETURN END