Esope source files

   1 : C FOBEY1    SOURCE    CHAT      05/01/13    00:05:40     5004           
   2 :       SUBROUTINE  FOBEY1(XR1,NXR1,XR2)
   3 :       IMPLICIT INTEGER(I-N)
   4 :       IMPLICIT REAL*8 (A-H,O-Z)
   5 : C    Include contenant quelques constantes dont XPI :
   6 : -INC CCREEL
   7 : C     
   8 : C     
   9 : C     
  10 : C  XR1( NXR1) : TABLEAU DE REELS TOUS > 0
  11 : C  XR2( NXR1) : TABLEAU DE REELS
  12 : C     
  13 : C           XR2(I) = Y1 (  XR1(I) )
  14 : C      Y1 FONCTION DE BESSEL DE DEUXIEME ESPECE D'ORDRE 1
  15 : C     
  16 : C       APPROXIMATION POLYNOMIALE PAR SECTEUR
  17 : C       REFERENCE : ABRAMOWITZ HANDBOOK OF MATHEMATICAL FONCTIONS
  18 : C      PRECISION E = 1.D-8
  19 : C     
  20 :       DIMENSION XR1(NXR1)
  21 :       DIMENSION XR2(NXR1)
  22 :       DIMENSION XB(1)
  23 :       DIMENSION XBB(1)
  24 : C     
  25 :       A0= -.6366198D0
  26 :       A2=  .2212091D0
  27 :       A4= 2.1682709D0
  28 :       A6=-1.3164827D0
  29 :       A8=  .3123951D0
  30 :       A10=-.0400976D0
  31 :       A12= .0027873D0
  32 : C     
  33 :       B0=  .79788456D0
  34 :       B1=  .00000156D0
  35 :       B2=  .01659667D0
  36 :       B3=  .00017105D0
  37 :       B4= -.00249511D0
  38 :       B5=  .00113653D0
  39 :       B6= -.00020033D0
  40 : C     
  41 :       C0=-2.35619816D0
  42 :       C1=  .12499612D0
  43 :       C2=  .00005650D0
  44 :       C3= -.00637879D0
  45 :       C4=  .00074348D0
  46 :       C5=  .00079824D0
  47 :       C6= -.00029166D0
  48 : C     
  49 :       DO 100  I=1,NXR1
  50 :          IF (XR1(I) .LE. 3.D0) THEN
  51 :             XB(1)= XR1(I)
  52 :             Y2= (XB(1) / 3D0)**2
  52 :             CALL FOBEJ1(XB,1,XBB)
  53 :             A00=(2.D0 /XPI)*XB(1)*(LOG(XB(1)/2.D0))* XBB(1)
  54 :             BB= A00+A0+Y2*(A2+Y2*(A4+Y2*(A6+Y2*(A8+Y2*(A10+Y2*A12)))))
  55 :             XR2(I)= BB / XR1(I)
  56 :          ELSE
  57 :             Y= 3D0/XR1(I)
  58 :             F1= B0+Y*(B1+Y*(B2+Y*(B3+Y*(B4+Y*(B5+Y*B6)))))
  59 :             T1 = XR1(I)+C0+Y*(C1+Y*(C2+Y*(C3+Y*(C4+Y*(C5+Y*C6)))))
  60 :             XR2(I) = (1.D0 / ( XR1(I) ** .5D0)) * F1 * SIN(T1)
  61 :          ENDIF
  62 :  100  CONTINUE
  63 : C     
  64 :       RETURN
  65 :       END
  66 : 
  67 : 
  68 :  
  69 :  
  70 :