Sources Esope | Cast3M

fobej0
C FOBEJ0    SOURCE    CHAT      05/01/13    00:05:28     5004
      SUBROUTINE  FOBEJ0(XR1,NXR1,XR2)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C
C  XR1( NXR1) : TABLEAU DE REELS TOUS >-3
C  XR2( NXR1) : TABLEAU DE REELS
C
C           XR2(I) = J0 (  XR1(I) )
C      J0 FONCTION DE BESSEL DE PREMIERE ESPECE D'ORDRE 0
C
C       APPROXIMATION POLYNOMIALE PAR SECTEUR
C       REFERENCE : ABRAMOWITZ HANDBOOK OF MATHEMATICAL FONCTIONS
C       PRECISION  E = 1.E-8
C
      DIMENSION XR1(NXR1)
      DIMENSION XR2(NXR1)
C
      A0=         1.D0
      A2= -2.2499997D0
      A4=  1.2656208D0
      A6=  -.3163866D0
      A8=   .0444479D0
      A10= -.0039444D0
      A12=  .0002100D0
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
          Y= XR1(I) / 3D0
          XR2(I) = A0 + (A2 * (Y**2))+ ( A4 * (Y**4 )) + (A6 * (Y**6 ))
     &                + (A8 * (Y**8))+ ( A10* (Y**10)) + (A12* (Y**12))
        ELSE
          Y= 3D0/XR1(I)
          F0 = B0 + (B1 * Y)+ ( B2 * (Y**2 )) + (B3 * (Y**3 ))
     &            + (B4 * (Y**4))+ ( B5* (Y**5)) + (B6* (Y**6))
          T0 =  XR1(I)+C0 + (C1 * Y)+ ( C2 * (Y**2 )) + (C3 * (Y**3 ))
     &            + (C4 * (Y**4))+ ( C5* (Y**5)) + (C6* (Y**6))
          XR2(I) = (1. / ( XR1(I) ** .5)) * F0 * COS(T0)
        ENDIF
100   CONTINUE
C
      RETURN
      END