blob: e0d0595a587eb17cf7c9e4f974828e06a2088029 (
plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
|
! "bessj0" function that returns the Bessel function J0(x) for any real x
!
!
! REFERENCE: "Numerical recipes in FORTRAN,
! The Art of Scientific Computing, Second Edition"
! Press, Teukolsky, Vetterling, Flannery
! Cambridge University Press
! ISBN 0-521-43064-X
! PAGES: 225 ("FUNCTION bessj0")
!
! INPUT: "x" the input argument (REAL)
! OUTPUT: "bessj0" is the value of J0(x) (REAL)
REAL bessj0 x ;
REAL ax xx z ;
DOUBLE p1 p2 p3 p4 p5 := 1.D0 -.1098628627D-2 .2734510407D-4
-.2073370639D-5 .2093887211D-6 ;
DOUBLE q1 q2 q3 q4 q5 := -.1562499995D-1 .1430488765D-3 -.6911147651D-5
.7621095161D-6 -.934945152D-7 ;
DOUBLE r1 r2 r3 r4 r5 r6 := 57568490574.D0 -13362590354.D0
651619640.7D0 -11214424.18D0 77392.33017D0 -184.9052456D0 ;
DOUBLE s1 s2 s3 s4 s5 s6 := 57568490411.D0 1029532985.D0
9494680.718D0 59272.64853D0 267.8532712D0 1.D0 ;
DOUBLE y ;
:: >>x<< ;
IF x ABS 8. < THEN
EVALUATE y := x x * R_TO_D ;
EVALUATE bessj0 := r6 y * r5 + y * r4 + y * r3 + y * r2 + y * r1 +
s6 y * s5 + y * s4 + y * s3 + y * s2 + y * s1 + /
D_TO_R ;
ELSE
EVALUATE ax := x ABS ;
EVALUATE z xx := 8. ax / ax .785398164 - ;
EVALUATE y := z z * R_TO_D ;
EVALUATE bessj0 := p5 y * p4 + y * p3 + y * p2 + y * p1 +
xx COS R_TO_D *
q5 y * q4 + y * q3 + y * q2 + y * q1 +
xx SIN z * R_TO_D * -
.636619772 ax / SQRT R_TO_D *
D_TO_R ;
ENDIF ;
:: <<bessj0>> ;
QUIT " Function *bessj0* XREF " .
|
