Sources Esope | Cast3M

fobej1
C FOBEJ1    SOURCE    CHAT      05/01/13    00:05:31     5004
      SUBROUTINE  FOBEJ1(XR1,NXR1,XR2)
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C
C
C  XR1( NXR1) : TABLEAU DE REELS TOUS >-3
C  XR2( NXR1) : TABLEAU DE REELS
C
C           XR2(I) = J1 (  XR1(I) )
C      J1 FONCTION DE BESSEL DE PREMIERE ESPECE D'ORDRE 1
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=         .5D0
      A2= -.56249985D0
      A4=  .21093573D0
      A6= -.03954289D0
      A8=  .00443319D0
      A10=-.00031761D0
      A12= .00001109D0
C
      B0=  .79788456D0
      B1=  .00000156D0
      B2=  .01659667D0
      B3=  .00017105D0
      B4= -.00249511D0
      B5=  .00113653D0
      B6= -.00020033D0
C
      C0=-2.35619816D0
      C1=  .12499612D0
      C2=  .00005650D0
      C3= -.00637879D0
      C4=  .00074348D0
      C5=  .00079824D0
      C6= -.00029166D0
C
      DO 100  I=1,NXR1
        IF (XR1(I) .LE. 3D0) THEN
          Y= XR1(I) / 3D0
          BB     = A0 + (A2 * (Y**2))+ ( A4 * (Y**4 )) + (A6 * (Y**6 ))
     &                + (A8 * (Y**8))+ ( A10* (Y**10)) + (A12* (Y**12))
          XR2(I)= XR1(I) *  BB
        ELSE
          Y= 3D0/XR1(I)
          F1 = B0 + (B1 * Y)+ ( B2 * (Y**2 )) + (B3 * (Y**3 ))
     &            + (B4 * (Y**4))+ ( B5* (Y**5)) + (B6* (Y**6))
          T1 =  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)) * F1 * COS(T1)
        ENDIF
100   CONTINUE
C
      RETURN
      END