Esope source files

   1 : C FOBEJ1    SOURCE    CHAT      05/01/13    00:05:31     5004           
   2 :       SUBROUTINE  FOBEJ1(XR1,NXR1,XR2)                                          
   3 :       IMPLICIT INTEGER(I-N)
   4 :       IMPLICIT REAL*8 (A-H,O-Z)                                                 
   5 : C                                                                               
   6 : C                                                                               
   7 : C  XR1( NXR1) : TABLEAU DE REELS TOUS >-3                                       
   8 : C  XR2( NXR1) : TABLEAU DE REELS                                                
   9 : C                                                                               
  10 : C           XR2(I) = J1 (  XR1(I) )                                             
  11 : C      J1 FONCTION DE BESSEL DE PREMIERE ESPECE D'ORDRE 1                       
  12 : C                                                                               
  13 : C       APPROXIMATION POLYNOMIALE PAR SECTEUR                                   
  14 : C       REFERENCE : ABRAMOWITZ HANDBOOK OF MATHEMATICAL FONCTIONS               
  15 : C       PRECISION  E = 1.E-8                                                    
  16 : C                                                                               
  17 :       DIMENSION XR1(NXR1)                                                       
  18 :       DIMENSION XR2(NXR1)                                                       
  19 : C                                                                               
  20 :       A0=         .5D0                                                          
  21 :       A2= -.56249985D0                                                          
  22 :       A4=  .21093573D0                                                          
  23 :       A6= -.03954289D0                                                          
  24 :       A8=  .00443319D0                                                          
  25 :       A10=-.00031761D0                                                          
  26 :       A12= .00001109D0                                                          
  27 : C                                                                               
  28 :       B0=  .79788456D0                                                          
  29 :       B1=  .00000156D0                                                          
  30 :       B2=  .01659667D0                                                          
  31 :       B3=  .00017105D0                                                          
  32 :       B4= -.00249511D0                                                          
  33 :       B5=  .00113653D0                                                          
  34 :       B6= -.00020033D0                                                          
  35 : C                                                                               
  36 :       C0=-2.35619816D0                                                          
  37 :       C1=  .12499612D0                                                          
  38 :       C2=  .00005650D0                                                          
  39 :       C3= -.00637879D0                                                          
  40 :       C4=  .00074348D0                                                          
  41 :       C5=  .00079824D0                                                          
  42 :       C6= -.00029166D0                                                          
  43 : C                                                                               
  44 :       DO 100  I=1,NXR1                                                          
  45 :         IF (XR1(I) .LE. 3D0) THEN                                               
  46 :           Y= XR1(I) / 3D0                                                       
  47 :           BB     = A0 + (A2 * (Y**2))+ ( A4 * (Y**4 )) + (A6 * (Y**6 ))         
  48 :      &                + (A8 * (Y**8))+ ( A10* (Y**10)) + (A12* (Y**12))         
  49 :           XR2(I)= XR1(I) *  BB                                                  
  50 :         ELSE                                                                    
  51 :           Y= 3D0/XR1(I)                                                         
  52 :           F1 = B0 + (B1 * Y)+ ( B2 * (Y**2 )) + (B3 * (Y**3 ))                  
  53 :      &            + (B4 * (Y**4))+ ( B5* (Y**5)) + (B6* (Y**6))                 
  54 :           T1 =  XR1(I)+C0 + (C1 * Y)+ ( C2 * (Y**2 )) + (C3 * (Y**3 ))          
  55 :      &            + (C4 * (Y**4))+ ( C5* (Y**5)) + (C6* (Y**6))                 
  56 :           XR2(I) = (1. / ( XR1(I) ** .5)) * F1 * COS(T1)                        
  57 :         ENDIF                                                                   
  58 : 100   CONTINUE                                                                  
  59 : C                                                                               
  60 :       RETURN                                                                    
  61 :       END                                                                       
  62 :