Sources Esope | Cast3M

fobey0
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