%!
% Author : Christophe JORSSEN
% ('libre' is the french word for 'free' if you want to contact me ;-))
% Created the : Sat 20 March 2004
% Version : 1.0 $Revision: 1.1 $
% Modified : Dimanche 25 décembre 2005
% par Jean-Paul VIGNAULT
% pour utilisation par jps2ps
% fichier original : http://ftp.jussieu.fr/ftp/pub/tex-archive/graphics/pstricks/contrib/pst-math
/Gammaln {dup dup dup 5.5 add dup ln 3 -1 roll .5 add mul sub neg 1.000000000190015
0 1 5 {
[76.18009172947146 -86.50532032941677 24.0140982483091 -1.231739572450155
.1208650973866179E-2 -.5395239384953E-5 2.5066282746310005] exch get
4 -1 roll 1 add dup 5 1 roll div add} for
4 -1 roll div 2.5066282746310005 mul ln add exch pop} bind def
/BETA {2 copy add Gammaln neg exch Gammaln 3 -1 roll Gammaln Exp} bind def
/Horner {aload length
dup 2 add -1 roll
exch 1 sub {
dup 4 1 roll
mul add exch
} repeat
pop
} bind def
/Bessel_j0 {dup abs 8 lt {
dup mul dup [57568490574 -13362590354 651619640.7 -11214424.18 77392.33017 -184.9052456] Horner
exch [57568490411 1029532985 9494680.718 59272.64853 267.8532712 1] Horner
div}
{abs dup .636619772 exch div sqrt exch dup .785398164 sub exch 8 exch div dup dup mul dup
[1 -1.098628627E-2 .2734510407E-4 -.2073370639E-5 .2093887211E-6] Horner
3 index Cos mul
exch [-.1562499995E-1 .1430488765E-3 -.6911147651E-5 .7621095161E-6 -.934945152E-7] Horner
4 -1 roll Sin mul 3 -1 roll mul neg add mul}
ifelse} bind def
/Bessel_y0 {dup 8 lt {
dup dup mul dup [-2957821389 7062834065 -512359803.6 10879881.29 -86327.92757 228.4622733] Horner
exch [40076544269 745249964.8 7189466.438 47447.26470 226.1030244 1] Horner
div exch dup ln exch Bessel_j0 .636619772 mul mul add}
{dup .636619772 exch div sqrt exch dup .785398164 sub exch 8 exch div dup dup mul dup
[1 -.1098628627E-2 .2734510407E-4 -.2073370639E-5 .2093887211E-6] Horner
3 index Sin mul
exch [-.1562499995E-1 .1430488765E-3 -.6911147651E-5 .7621095161E-6 -.934945152E-7] Horner
4 -1 roll Cos mul 3 -1 roll mul add mul}
ifelse} bind def
/Bessel_j1 {dup abs 8 lt {
dup dup mul dup 3 -2 roll [72362614232 -7895059235 242396853.1 -2972611.439 15704.48260 -30.16036606] Horner mul
exch [144725228442 2300535178 18583304.74 99447.43394 376.9991397 1] Horner
div}
{dup abs dup .636619772 exch div sqrt exch dup 2.356194491 sub exch 8 exch div dup dup mul dup
[1 .183105E-2 -.3516396496E-4 .2457520174E-5 -.240337019E-6] Horner
3 index Cos mul
exch [.04687499995 6.2002690873E-3 .8449199096E-5 -.88228987E-6 .105787412E-6] Horner
4 -1 roll Sin mul 3 -1 roll mul neg add mul exch dup abs div mul}
ifelse} bind def
/Bessel_y1 {dup 8 lt {
dup dup dup mul dup [-.4900604943E13 .1275274390E13 -.5153428139E11 .7349264551E9 -.4237922726E7 .8511937935E4] Horner
exch [.2499580570E14 .4244419664E12 .3733650367E10 .2245904002E8 .1020426050E6 .3549632885E3 1] Horner
div mul exch dup dup ln exch Bessel_j1 mul exch 1 exch div sub .636619772 mul add}
{dup .636619772 exch div sqrt exch dup 2.356194491 sub exch 8 exch div dup dup mul dup
[1 .183105E-2 -.3516396496E-4 .2457520174E-5 -.240337019E-6] Horner
3 index Sin mul
exch [.04687499995 -.2002690873E-3 .8449199096E-5 6.88228987E-6 .105787412E-6] Horner
4 -1 roll Cos mul 3 -1 roll mul add mul}
ifelse} bind def
% En cours...
/Bessel_yn {dup 0 eq {pop Bessel_y0}{dup 1 eq {pop Bessel_y1}{
exch dup Bessel_y0 exch dup Bessel_y1 exch 2 exch div {
mul 3 -1 roll mul 2 index sub pstack} for
} ifelse } ifelse } bind def
|