--- a/ChangeLog	Wed May 14 18:09:56 2008 +0200
+++ b/ChangeLog	Sun Apr 27 22:34:17 2008 +0200
@@ -1,3 +1,8 @@
+2008-05-20  David Bateman  <dbateman@free.fr>
+
+	* configure.in (AC_CHECK_FUNCS): Add expm1, lgammaf, lgammaf_r,
+	log1pf and tgammaf. Also check for libfftw3f.
+
 2008-04-09  Rafael Laboissiere  <rafael@debian.org>
 
 	* example/octave.desktop.in: Drop the Encoding key, which is

--- a/PROJECTS	Wed May 14 18:09:56 2008 +0200
+++ b/PROJECTS	Sun Apr 27 22:34:17 2008 +0200
@@ -117,7 +117,6 @@
       - minres
       - qmr
       - symmlq 
-      - spaugment
 
 -------
 Strings:
@@ -143,6 +142,9 @@
 
   * Template functions for mixed-type ops.
 
+  * Convert other functions for use with the floating point type
+  including quad, dasrt, daspk, etc.
+
 ------------
 Input/Output:
 ------------

--- a/configure.in	Wed May 14 18:09:56 2008 +0200
+++ b/configure.in	Sun Apr 27 22:34:17 2008 +0200
@@ -562,7 +562,8 @@
   with_fftw3=no
   AC_CHECK_HEADER(fftw3.h, [have_fftw3_header=yes])
   if test "$have_fftw3_header" = yes; then
-    AC_CHECK_LIB(fftw3, fftw_plan_dft_1d, [FFTW_LIBS="-lfftw3"; with_fftw3=yes])
+    AC_CHECK_LIB(fftw3, fftw_plan_dft_1d, [
+      AC_CHECK_LIB(fftw3f, fftwf_plan_dft_1d, [FFTW_LIBS="-lfftw3 -lfftw3f"; with_fftw3=yes])])
   fi
 fi
 
@@ -1319,15 +1320,15 @@
 ### Checks for functions and variables.
 
 AC_CHECK_FUNCS(atexit basename bcopy bzero canonicalize_file_name \
-  chmod dup2 endgrent endpwent execvp expm1 fcntl fork getcwd \
+  chmod dup2 endgrent endpwent execvp expm1 expm1f fcntl fork getcwd \
   getegid geteuid getgid getgrent getgrgid getgrnam getpgrp getpid \
   getppid getpwent getpwuid gettimeofday getuid getwd _kbhit kill \
-  lgamma lgamma_r link localtime_r log1p lstat memmove mkdir mkfifo \
-  mkstemp on_exit pipe poll putenv raise readlink realpath rename \
-  resolvepath rindex rmdir round select setgrent setlocale setpwent \
-  setvbuf sigaction siglongjmp sigpending sigprocmask sigsuspend \
+  lgamma lgammaf lgamma_r lgammaf_r link localtime_r log1p log1pf lstat \
+  memmove mkdir mkfifo mkstemp on_exit pipe poll putenv raise readlink \
+  realpath rename resolvepath rindex rmdir round select setgrent setlocale \
+  setpwent setvbuf sigaction siglongjmp sigpending sigprocmask sigsuspend \
   snprintf stat strcasecmp strdup strerror stricmp strncasecmp \
-  strnicmp strptime strsignal symlink tempnam tgamma trunc umask \
+  strnicmp strptime strsignal symlink tempnam tgamma tgammaf trunc umask \
   uname unlink usleep utime vfprintf vsprintf vsnprintf waitpid \
   _chmod _snprintf x_utime _utime32)
 
@@ -1539,7 +1540,7 @@
 
 ### Check for nonstandard but common math functions that we need.
 
-AC_CHECK_FUNCS(acosh asinh atanh erf erfc exp2 log2)
+AC_CHECK_FUNCS(acosh acoshf asinh asinhf atanh atanhf erf erff erfc erfcf exp2 exp2f log2 log2f)
 
 ### Checks for OS specific cruft.
 

--- a/libcruft/ChangeLog	Wed May 14 18:09:56 2008 +0200
+++ b/libcruft/ChangeLog	Sun Apr 27 22:34:17 2008 +0200
@@ -1,3 +1,191 @@
+2008-05-20  Jaroslav Hajek <highegg@gmail.com>
+
+	* qrupdate/cch1dn.f, qrupdate/cchinx.f, qrupdate/cqhqr.f, 
+	qrupdate/cqrinc.f, qrupdate/cqrinr.f, qrupdate/cqrqhu.f, 
+	qrupdate/cqrqhv.f, qrupdate/sch1dn.f, qrupdate/schinx.f, 
+	qrupdate/sqhqr.f, qrupdate/sqrinc.f, qrupdate/sqrinr.f, 
+	qrupdate/sqrqhu.f: Convert DOUBLE PRECISION constants to REAL.
+	* qrupdate/cqrinr.f, qrupdate/sqrinr.f: Correct EXTERNAL
+	declarations.
+	* qrupdate/sqrinr.f: Convert DOUBLE PRECISION calls to
+	REAL counterparts.
+
+2008-05-20  David Bateman  <dbateman@free.fr>
+
+	* Makefile.in (MISC_OBJ): Add misc/smachar.o
+	* Makerules.in (CRUFT_CSRC, CRUFT_CPICOBJ): Add CEXTRA, allowing
+	objects files with no corresponding source file in the
+	distribution.
+
+	* amos/cacai.f, amos/cacon.f, amos/cbesh.f, amos/cbesi.f,
+	amos/cbesj.f, amos/cbesk.f, amos/cbesy.f, amos/cbinu.f,
+	amos/cbuni.f, amos/cbunk.f, amos/cunk1.f amos/cunk2.f,
+	amos/crati.f, amos/cshch.f, amos/cuni1.f, amos/cuoik.f,
+	amos/cairy.f, amos/cbiry.f, amos/ckscl.f, amos/cs1s2.f,
+	amos/cuchk.f, amos/cuni2.f, amos/cwrsk.f, amos/casyi.f,
+	amos/cbknu.f, amos/cmlri.f, amos/cseri.f, amos/cunhj.f,
+	amos/cunik.f: New files.
+	* amos/Makefile.in (FSRC): Add them.
+
+	* blas-xtra/xsdot.f, blas-xtra/xsnrm2.f, blas-xtra/xscnrm2.f,
+	blas-xtra/xcdotc.f, blas-xtra/xcdotu.f: New files
+	* blas-xtra/Makefile.in (FSRC): Add them.
+
+	* blas/sasum.f, blas/saxpy.f, blas/scabs1.f, blas/scopy.f,
+	blas/sger.f, blas/smach.f, blas/snrm2.f, blas/srot.f,
+	blas/sswap.f, blas/ssymv.f, blas/ssyr.f, blas/ssyr2.f,
+	blas/ssyr2k.f, blas/stbsv.f, blas/strmm.f, blas/strmv.f,
+	blas/strsv.f, blas/scasum.f, blas/scnrm2.f, blas/caxpy.f,
+	blas/ccopy.f, blas/cdotc.f, blas/cdotu.f, blas/, blas/csrot.f,
+	blas/csscal.f, blas/cgemm.f, blas/cgemv.f, blas/cgerc.f,
+	blas/cgeru.f, blas/chemv.f, blas/cher.f, blas/cher2.f,
+	blas/cher2k.f, blas/cherk.f, blas/cscal.f, blas/cswap.f,
+	blas/ctbsv.f, blas/ctrmm.f, blas/ctrmv.f, blas/, blas/ctrsm.f,
+	blas/ctrsv.f: New files
+	* blas/Makefile.in (FSRC): Add them.
+
+	* fftpack/zfftb.f, zfftb1.f, fftpack/zfftf.f, fftpack/zfftf1.f,
+	fftpack/zffti.f, fftpack/zffti1.f, fftpack/zpassb.f,
+	fftpack/zpassb2.f, fftpack/zpassb3.f, fftpack/zpassb4.f,
+	fftpack/zpassb5.f, fftpack/zpassf.f, fftpack/zpassf2.f,
+	fftpack/zpassf3.f, fftpack/zpassf4.f, fftpack/zpassf5.f: Rename
+	function (c -> z | add z).
+	* fftpack/cfftb.f, cfftb1.f, fftpack/cfftf.f, fftpack/cfftf1.f,
+	fftpack/cffti.f, fftpack/cffti1.f, fftpack/passb.f,
+	fftpack/passb2.f, fftpack/passb3.f, fftpack/passb4.f,
+	fftpack/passb5.f, fftpack/passf.f, fftpack/passf2.f,
+	fftpack/passf3.f, fftpack/passf4.f, fftpack/passf5.f: New files
+	for single precision.
+	* fftpack/Makefile.in (FSRC): Add new files.
+
+	* lapack-xtra/xclange.f, lapack-xtra/xslamch.f,
+	lapack-xtra/xslange.f: New files.
+	* lapack-xtra/Makefile.in (FSRC): Add them.
+
+	* lapack/cbdsqr.f, lapack/csrscl.f, lapack/cgbcon.f,
+	lapack/cgbtf2.f, lapack/cgbtrf.f, lapack/cgbtrs.f,
+	lapack/cgebak.f, lapack/cgebal.f, lapack/cgebd2.f,
+	lapack/cgebrd.f, lapack/cgecon.f, lapack/cgeesx.f, lapack/cgeev.f,
+	lapack/cgehd2.f, lapack/cgehrd.f, lapack/cgelq2.f,
+	lapack/cgelqf.f, lapack/cgelsd.f, lapack/cgelss.f,
+	lapack/cgelsy.f, lapack/cgeqp3.f, lapack/cgeqpf.f,
+	lapack/cgeqr2.f, lapack/cgeqrf.f, lapack/cgesv.f, lapack/cgesvd.f,
+	lapack/cgetf2.f, lapack/cgetrf.f, lapack/cgetri.f,
+	lapack/cgetrs.f, lapack/cggbal.f, lapack/cgtsv.f, lapack/cgttrf.f,
+	lapack/cgttrs.f, lapack/cgtts2.f, lapack/cheev.f, lapack/chetd2.f,
+	lapack/chetrd.f, lapack/chseqr.f, lapack/clabrd.f,
+	lapack/clacgv.f, lapack/clacn2.f, lapack/clacon.f,
+	lapack/clacpy.f, lapack/cladiv.f, lapack/clahqr.f,
+	lapack/clahr2.f, lapack/clahrd.f, lapack/claic1.f,
+	lapack/clals0.f, lapack/clalsa.f, lapack/clalsd.f,
+	lapack/clange.f, lapack/clanhe.f, lapack/clanhs.f,
+	lapack/clantr.f, lapack/claqp2.f, lapack/claqps.f,
+	lapack/claqr0.f, lapack/claqr1.f, lapack/claqr2.f,
+	lapack/claqr3.f, lapack/claqr4.f, lapack/claqr5.f, lapack/clarf.f,
+	lapack/clarfb.f, lapack/clarfg.f, lapack/clarft.f,
+	lapack/clarfx.f, lapack/clartg.f, lapack/clarz.f, lapack/clarzb.f,
+	lapack/clarzt.f, lapack/clascl.f, lapack/claset.f, lapack/clasr.f,
+	lapack/classq.f, lapack/claswp.f, lapack/clatbs.f,
+	lapack/clatrd.f, lapack/clatrs.f, lapack/clatrz.f,
+	lapack/clauu2.f, lapack/clauum.f, lapack/cpbcon.f,
+	lapack/cpbtf2.f, lapack/cpbtrf.f, lapack/cpbtrs.f,
+	lapack/cpocon.f, lapack/cpotf2.f, lapack/cpotrf.f,
+	lapack/cpotri.f, lapack/cpotrs.f, lapack/cptsv.f, lapack/cpttrf.f,
+	lapack/cpttrs.f, lapack/cptts2.f, lapack/crot.f, lapack/csteqr.f,
+	lapack/ctrcon.f, lapack/ctrevc.f, lapack/ctrexc.f,
+	lapack/ctrsen.f, lapack/ctrsyl.f, lapack/ctrti2.f,
+	lapack/ctrtri.f, lapack/ctrtrs.f, lapack/ctzrzf.f,
+	lapack/cung2l.f, lapack/cung2r.f, lapack/cungbr.f,
+	lapack/cunghr.f, lapack/cungl2.f, lapack/cunglq.f,
+	lapack/cungql.f, lapack/cungqr.f, lapack/cungtr.f,
+	lapack/cunm2r.f, lapack/cunmbr.f, lapack/cunml2.f,
+	lapack/cunmlq.f, lapack/cunmqr.f, lapack/cunmr3.f,
+	lapack/cunmrz.f, lapack/sbdsqr.f, lapack/sgbcon.f,
+	lapack/sgbtf2.f, lapack/sgbtrf.f, lapack/sgbtrs.f,
+	lapack/sgebak.f, lapack/sgebal.f, lapack/sgebd2.f,
+	lapack/sgebrd.f, lapack/sgecon.f, lapack/sgeesx.f, lapack/sgeev.f,
+	lapack/sgehd2.f, lapack/sgehrd.f, lapack/sgelq2.f,
+	lapack/sgelqf.f, lapack/sgelsd.f, lapack/sgelss.f,
+	lapack/sgelsy.f, lapack/sgeqp3.f, lapack/sgeqpf.f,
+	lapack/sgeqr2.f, lapack/sgeqrf.f, lapack/sgesv.f, lapack/sgesvd.f,
+	lapack/sgetf2.f, lapack/sgetrf.f, lapack/sgetri.f,
+	lapack/sgetrs.f, lapack/sggbak.f, lapack/sggbal.f,
+	lapack/sgghrd.f, lapack/sgtsv.f, lapack/sgttrf.f, lapack/sgttrs.f,
+	lapack/sgtts2.f, lapack/shgeqz.f, lapack/shseqr.f,
+	lapack/slabad.f, lapack/slabrd.f, lapack/slacn2.f,
+	lapack/slacon.f, lapack/slacpy.f, lapack/sladiv.f, lapack/slae2.f,
+	lapack/slaed6.f, lapack/slaev2.f, lapack/slaexc.f, lapack/slag2.f,
+	lapack/slahqr.f, lapack/slahr2.f, lapack/slahrd.f,
+	lapack/slaic1.f, lapack/slaln2.f, lapack/slals0.f,
+	lapack/slalsa.f, lapack/slalsd.f, lapack/slamc1.f,
+	lapack/slamc2.f, lapack/slamc3.f, lapack/slamc4.f,
+	lapack/slamc5.f, lapack/slamch.f, lapack/slamrg.f,
+	lapack/slange.f, lapack/slanhs.f, lapack/slanst.f,
+	lapack/slansy.f, lapack/slantr.f, lapack/slanv2.f,
+	lapack/slapy2.f, lapack/slapy3.f, lapack/slaqp2.f,
+	lapack/slaqps.f, lapack/slaqr0.f, lapack/slaqr1.f,
+	lapack/slaqr2.f, lapack/slaqr3.f, lapack/slaqr4.f,
+	lapack/slaqr5.f, lapack/slarf.f, lapack/slarfb.f, lapack/slarfg.f,
+	lapack/slarft.f, lapack/slarfx.f, lapack/slartg.f, lapack/slarz.f,
+	lapack/slarzb.f, lapack/slarzt.f, lapack/slas2.f, lapack/slascl.f,
+	lapack/slasd0.f, lapack/slasd1.f, lapack/slasd2.f,
+	lapack/slasd3.f, lapack/slasd4.f, lapack/slasd5.f,
+	lapack/slasd6.f, lapack/slasd7.f, lapack/slasd8.f,
+	lapack/slasda.f, lapack/slasdq.f, lapack/slasdt.f,
+	lapack/slaset.f, lapack/slasq1.f, lapack/slasq2.f,
+	lapack/slasq3.f, lapack/slasq4.f, lapack/slasq5.f,
+	lapack/slasq6.f, lapack/slasr.f, lapack/slasrt.f, lapack/slassq.f,
+	lapack/slasv2.f, lapack/slaswp.f, lapack/slasy2.f,
+	lapack/slatbs.f, lapack/slatrd.f, lapack/slatrs.f,
+	lapack/slatrz.f, lapack/slauu2.f, lapack/slauum.f,
+	lapack/slazq3.f, lapack/slazq4.f, lapack/sorg2l.f,
+	lapack/sorg2r.f, lapack/sorgbr.f, lapack/sorghr.f,
+	lapack/sorgl2.f, lapack/sorglq.f, lapack/sorgql.f,
+	lapack/sorgqr.f, lapack/sorgtr.f, lapack/sorm2r.f,
+	lapack/sormbr.f, lapack/sorml2.f, lapack/sormlq.f,
+	lapack/sormqr.f, lapack/sormr3.f, lapack/sormrz.f,
+	lapack/spbcon.f, lapack/spbtf2.f, lapack/spbtrf.f,
+	lapack/spbtrs.f, lapack/spocon.f, lapack/spotri.f,
+	lapack/spotrs.f, lapack/sptsv.f, lapack/spttrf.f, lapack/spttrs.f,
+	lapack/sptts2.f, lapack/srscl.f, lapack/ssteqr.f, lapack/ssterf.f,
+	lapack/ssyev.f, lapack/ssytd2.f, lapack/ssytrd.f, lapack/stgevc.f,
+	lapack/strcon.f, lapack/strevc.f, lapack/strexc.f,
+	lapack/strsen.f, lapack/strsyl.f, lapack/strti2.f,
+	lapack/strtri.f, lapack/strtrs.f, lapack/stzrzf.f,
+	lapack/scsum1.f: New files
+	* lapack/Makefile.in (FSRC): Add them.
+
+	* misc/r1mach.f: New file
+	* misc/machar.cc: Modify to allow to be build twice, once for
+	double precision and once for single precision.
+	* misc/Makefile.in (FSRC): Add it.
+	(CEXTRA): Add smachar.c, and target for smachar.o
+	(MAKEDEPS): Include CEXTRA.
+
+	* qrupdate/sch1up.f, qrupdate/cch1up.f, qrupdate/sqrinc.f,
+	qrupdate/cqrinc.f, qrupdate/sqrdec.f, qrupdate/cqrdec.f,
+	qrupdate/sqrinr.f, qrupdate/cqrinr.f, qrupdate/sqrder.f,
+	qrupdate/cqrder.f, qrupdate/sqrshc.f, qrupdate/cqrshc.f,
+	qrupdate/sqr1up.f, qrupdate/cqr1up.f, qrupdate/sch1dn.f,
+	qrupdate/cch1dn.f, qrupdate/schinx.f, qrupdate/cchinx.f,
+	qrupdate/schdex.f, qrupdate/cchdex.f, qrupdate/sqrqhu.f,
+	qrupdate/cqrqhu.f, qrupdate/sqrqhv.f, qrupdate/cqrqhv.f,
+	qrupdate/sqhqr.f, qrupdate/cqhqr.f: New files.
+	* qrupdate/Makefile.in (FSRC): Add them.
+
+	* slatec-fn/acosh.f, slatec-fn/albeta.f, slatec-fn/algams.f,
+	slatec-fn/alngam.f, slatec-fn/alnrel.f, slatec-fn/asinh.f,
+	slatec-fn/atanh.f, slatec-fn/betai.f, slatec-fn/csevl.f,
+	slatec-fn/erf.f, slatec-fn/erfc.f, slatec-fn/gami.f,
+	slatec-fn/gamit.f, slatec-fn/gamlim.f, slatec-fn/gamma.f,
+	slatec-fn/gamr.f, slatec-fn/inits.f, slatec-fn/pchim.f,
+	slatec-fn/pchst.f, slatec-fn/r9gmit.f, slatec-fn/r9lgic.f,
+	slatec-fn/r9lgit.f, slatec-fn/r9lgmc.f, slatec-fn/xacosh.f,
+	slatec-fn/xasinh.f, slatec-fn/xatanh.f, slatec-fn/xbetai.f,
+	slatec-fn/xerf.f, slatec-fn/xerfc.f, slatec-fn/xgamma.f,
+	slatec-fn/xsgmainc.f: New files.
+	* slatec-fn/Makefile.in (FSRC): Add them.	
+
 2008-04-20  Jaroslav Hajek <highegg@gmail.com>
 
 	* qrupdate/dch1dn.f, qrupdate/dchdex.f, qrupdate/dchinx.f,

--- a/libcruft/Makefile.in	Wed May 14 18:09:56 2008 +0200
+++ b/libcruft/Makefile.in	Sun Apr 27 22:34:17 2008 +0200
@@ -75,7 +75,7 @@
 # FIXME -- this should build the shared library directly from
 # a normal archive file (created from PIC code, though).
 
-MISC_OBJ := misc/machar.o misc/f77-extern.o \
+MISC_OBJ := misc/machar.o misc/smachar.o misc/f77-extern.o \
 	misc/f77-fcn.o misc/lo-error.o misc/quit.o misc/cquit.o
 
 CRUFT_FSRC := $(foreach dir, $(SUBDIRS), $(wildcard $(srcdir)/$(dir)/*.f))

--- a/libcruft/Makerules.in	Wed May 14 18:09:56 2008 +0200
+++ b/libcruft/Makerules.in	Sun Apr 27 22:34:17 2008 +0200
@@ -29,7 +29,7 @@
 DLL_CXXDEFS = @CRUFT_DLL_DEFS@
 
 CRUFT_FSRC = $(addprefix $(srcdir)/, $(FSRC))
-CRUFT_CSRC = $(addprefix $(srcdir)/, $(CSRC))
+CRUFT_CSRC = $(addprefix $(srcdir)/, $(CSRC) $(CEXTRA))
 CRUFT_CXXSRC = $(addprefix $(srcdir)/, $(CXXSRC))
 
 CRUFT_SRC = $(CRUFT_FSRC) $(CRUFT_CSRC) $(CRUFT_CXXSRC)
@@ -61,9 +61,9 @@
     CRUFT_FPICOBJ := $(CRUFT_FOBJ)
   endif
   ifdef CPICFLAG
-    CRUFT_CPICOBJ := $(addprefix pic/, $(CRUFT_COBJ))
+    CRUFT_CPICOBJ := $(addprefix pic/, $(CRUFT_COBJ) $(CEXTRA))
   else
-    CRUFT_CPICOBJ := $(CRUFT_COBJ)
+    CRUFT_CPICOBJ := $(CRUFT_COBJ) $(CEXTRA)
   endif
   ifdef CXXPICFLAG
     CRUFT_CXXPICOBJ := $(addprefix pic/, $(CRUFT_CXXOBJ))

--- a/libcruft/amos/Makefile.in	Wed May 14 18:09:56 2008 +0200
+++ b/libcruft/amos/Makefile.in	Sun Apr 27 22:34:17 2008 +0200
@@ -26,7 +26,11 @@
 
 EXTERNAL_DISTFILES = $(DISTFILES)
 
-FSRC = dgamln.f xzabs.f xzexp.f xzlog.f xzsqrt.f zacai.f zacon.f \
+FSRC = cacai.f cacon.f cbesh.f cbesi.f cbesj.f cbesk.f cbesy.f cbinu.f \
+  cbuni.f cbunk.f cunk1.f cunk2.f crati.f cshch.f cuni1.f \
+  cuoik.f cairy.f cbiry.f ckscl.f cs1s2.f cuchk.f cuni2.f cwrsk.f \
+  casyi.f cbknu.f cmlri.f cseri.f cunhj.f cunik.f dgamln.f gamln.f \
+  xzabs.f xzexp.f xzlog.f xzsqrt.f zacai.f zacon.f \
   zairy.f zasyi.f zbesh.f zbesi.f zbesj.f zbesk.f zbesy.f zbinu.f \
   zbiry.f zbknu.f zbuni.f zbunk.f zdiv.f zkscl.f zmlri.f zmlt.f \
   zrati.f zs1s2.f zseri.f zshch.f zuchk.f zunhj.f zuni1.f zuni2.f \

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libcruft/amos/cacai.f	Sun Apr 27 22:34:17 2008 +0200
@@ -0,0 +1,90 @@
+      SUBROUTINE CACAI(Z, FNU, KODE, MR, N, Y, NZ, RL, TOL, ELIM, ALIM)
+C***BEGIN PROLOGUE  CACAI
+C***REFER TO  CAIRY
+C
+C     CACAI APPLIES THE ANALYTIC CONTINUATION FORMULA
+C
+C         K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN)
+C                 MP=PI*MR*CMPLX(0.0,1.0)
+C
+C     TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT
+C     HALF Z PLANE FOR USE WITH CAIRY WHERE FNU=1/3 OR 2/3 AND N=1.
+C     CACAI IS THE SAME AS CACON WITH THE PARTS FOR LARGER ORDERS AND
+C     RECURRENCE REMOVED. A RECURSIVE CALL TO CACON CAN RESULT IF CACON
+C     IS CALLED FROM CAIRY.
+C
+C***ROUTINES CALLED  CASYI,CBKNU,CMLRI,CSERI,CS1S2,R1MACH
+C***END PROLOGUE  CACAI
+      COMPLEX CSGN, CSPN, C1, C2, Y, Z, ZN, CY
+      REAL ALIM, ARG, ASCLE, AZ, CPN, DFNU, ELIM, FMR, FNU, PI, RL,
+     * SGN, SPN, TOL, YY, R1MACH
+      INTEGER INU, IUF, KODE, MR, N, NN, NW, NZ
+      DIMENSION Y(N), CY(2)
+      DATA PI / 3.14159265358979324E0 /
+      NZ = 0
+      ZN = -Z
+      AZ = CABS(Z)
+      NN = N
+      DFNU = FNU + FLOAT(N-1)
+      IF (AZ.LE.2.0E0) GO TO 10
+      IF (AZ*AZ*0.25E0.GT.DFNU+1.0E0) GO TO 20
+   10 CONTINUE
+C-----------------------------------------------------------------------
+C     POWER SERIES FOR THE I FUNCTION
+C-----------------------------------------------------------------------
+      CALL CSERI(ZN, FNU, KODE, NN, Y, NW, TOL, ELIM, ALIM)
+      GO TO 40
+   20 CONTINUE
+      IF (AZ.LT.RL) GO TO 30
+C-----------------------------------------------------------------------
+C     ASYMPTOTIC EXPANSION FOR LARGE Z FOR THE I FUNCTION
+C-----------------------------------------------------------------------
+      CALL CASYI(ZN, FNU, KODE, NN, Y, NW, RL, TOL, ELIM, ALIM)
+      IF (NW.LT.0) GO TO 70
+      GO TO 40
+   30 CONTINUE
+C-----------------------------------------------------------------------
+C     MILLER ALGORITHM NORMALIZED BY THE SERIES FOR THE I FUNCTION
+C-----------------------------------------------------------------------
+      CALL CMLRI(ZN, FNU, KODE, NN, Y, NW, TOL)
+      IF(NW.LT.0) GO TO 70
+   40 CONTINUE
+C-----------------------------------------------------------------------
+C     ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
+C-----------------------------------------------------------------------
+      CALL CBKNU(ZN, FNU, KODE, 1, CY, NW, TOL, ELIM, ALIM)
+      IF (NW.NE.0) GO TO 70
+      FMR = FLOAT(MR)
+      SGN = -SIGN(PI,FMR)
+      CSGN = CMPLX(0.0E0,SGN)
+      IF (KODE.EQ.1) GO TO 50
+      YY = -AIMAG(ZN)
+      CPN = COS(YY)
+      SPN = SIN(YY)
+      CSGN = CSGN*CMPLX(CPN,SPN)
+   50 CONTINUE
+C-----------------------------------------------------------------------
+C     CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
+C     WHEN FNU IS LARGE
+C-----------------------------------------------------------------------
+      INU = INT(FNU)
+      ARG = (FNU-FLOAT(INU))*SGN
+      CPN = COS(ARG)
+      SPN = SIN(ARG)
+      CSPN = CMPLX(CPN,SPN)
+      IF (MOD(INU,2).EQ.1) CSPN = -CSPN
+      C1 = CY(1)
+      C2 = Y(1)
+      IF (KODE.EQ.1) GO TO 60
+      IUF = 0
+      ASCLE = 1.0E+3*R1MACH(1)/TOL
+      CALL CS1S2(ZN, C1, C2, NW, ASCLE, ALIM, IUF)
+      NZ = NZ + NW
+   60 CONTINUE
+      Y(1) = CSPN*C1 + CSGN*C2
+      RETURN
+   70 CONTINUE
+      NZ = -1
+      IF(NW.EQ.(-2)) NZ=-2
+      RETURN
+      END

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libcruft/amos/cacon.f	Sun Apr 27 22:34:17 2008 +0200
@@ -0,0 +1,149 @@
+      SUBROUTINE CACON(Z, FNU, KODE, MR, N, Y, NZ, RL, FNUL, TOL, ELIM,
+     * ALIM)
+C***BEGIN PROLOGUE  CACON
+C***REFER TO  CBESK,CBESH
+C
+C     CACON APPLIES THE ANALYTIC CONTINUATION FORMULA
+C
+C         K(FNU,ZN*EXP(MP))=K(FNU,ZN)*EXP(-MP*FNU) - MP*I(FNU,ZN)
+C                 MP=PI*MR*CMPLX(0.0,1.0)
+C
+C     TO CONTINUE THE K FUNCTION FROM THE RIGHT HALF TO THE LEFT
+C     HALF Z PLANE
+C
+C***ROUTINES CALLED  CBINU,CBKNU,CS1S2,R1MACH
+C***END PROLOGUE  CACON
+      COMPLEX CK, CONE, CS, CSCL, CSCR, CSGN, CSPN, CSS, CSR, C1, C2,
+     * RZ, SC1, SC2, ST, S1, S2, Y, Z, ZN, CY
+      REAL ALIM, ARG, ASCLE, AS2, BSCLE, BRY, CPN, C1I, C1M, C1R, ELIM,
+     * FMR, FNU, FNUL, PI, RL, SGN, SPN, TOL, YY, R1MACH
+      INTEGER I, INU, IUF, KFLAG, KODE, MR, N, NN, NW, NZ
+      DIMENSION Y(N), CY(2), CSS(3), CSR(3), BRY(3)
+      DATA PI / 3.14159265358979324E0 /
+      DATA CONE / (1.0E0,0.0E0) /
+      NZ = 0
+      ZN = -Z
+      NN = N
+      CALL CBINU(ZN, FNU, KODE, NN, Y, NW, RL, FNUL, TOL, ELIM, ALIM)
+      IF (NW.LT.0) GO TO 80
+C-----------------------------------------------------------------------
+C     ANALYTIC CONTINUATION TO THE LEFT HALF PLANE FOR THE K FUNCTION
+C-----------------------------------------------------------------------
+      NN = MIN0(2,N)
+      CALL CBKNU(ZN, FNU, KODE, NN, CY, NW, TOL, ELIM, ALIM)
+      IF (NW.NE.0) GO TO 80
+      S1 = CY(1)
+      FMR = FLOAT(MR)
+      SGN = -SIGN(PI,FMR)
+      CSGN = CMPLX(0.0E0,SGN)
+      IF (KODE.EQ.1) GO TO 10
+      YY = -AIMAG(ZN)
+      CPN = COS(YY)
+      SPN = SIN(YY)
+      CSGN = CSGN*CMPLX(CPN,SPN)
+   10 CONTINUE
+C-----------------------------------------------------------------------
+C     CALCULATE CSPN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
+C     WHEN FNU IS LARGE
+C-----------------------------------------------------------------------
+      INU = INT(FNU)
+      ARG = (FNU-FLOAT(INU))*SGN
+      CPN = COS(ARG)
+      SPN = SIN(ARG)
+      CSPN = CMPLX(CPN,SPN)
+      IF (MOD(INU,2).EQ.1) CSPN = -CSPN
+      IUF = 0
+      C1 = S1
+      C2 = Y(1)
+      ASCLE = 1.0E+3*R1MACH(1)/TOL
+      IF (KODE.EQ.1) GO TO 20
+      CALL CS1S2(ZN, C1, C2, NW, ASCLE, ALIM, IUF)
+      NZ = NZ + NW
+      SC1 = C1
+   20 CONTINUE
+      Y(1) = CSPN*C1 + CSGN*C2
+      IF (N.EQ.1) RETURN
+      CSPN = -CSPN
+      S2 = CY(2)
+      C1 = S2
+      C2 = Y(2)
+      IF (KODE.EQ.1) GO TO 30
+      CALL CS1S2(ZN, C1, C2, NW, ASCLE, ALIM, IUF)
+      NZ = NZ + NW
+      SC2 = C1
+   30 CONTINUE
+      Y(2) = CSPN*C1 + CSGN*C2
+      IF (N.EQ.2) RETURN
+      CSPN = -CSPN
+      RZ = CMPLX(2.0E0,0.0E0)/ZN
+      CK = CMPLX(FNU+1.0E0,0.0E0)*RZ
+C-----------------------------------------------------------------------
+C     SCALE NEAR EXPONENT EXTREMES DURING RECURRENCE ON K FUNCTIONS
+C-----------------------------------------------------------------------
+      CSCL = CMPLX(1.0E0/TOL,0.0E0)
+      CSCR = CMPLX(TOL,0.0E0)
+      CSS(1) = CSCL
+      CSS(2) = CONE
+      CSS(3) = CSCR
+      CSR(1) = CSCR
+      CSR(2) = CONE
+      CSR(3) = CSCL
+      BRY(1) = ASCLE
+      BRY(2) = 1.0E0/ASCLE
+      BRY(3) = R1MACH(2)
+      AS2 = CABS(S2)
+      KFLAG = 2
+      IF (AS2.GT.BRY(1)) GO TO 40
+      KFLAG = 1
+      GO TO 50
+   40 CONTINUE
+      IF (AS2.LT.BRY(2)) GO TO 50
+      KFLAG = 3
+   50 CONTINUE
+      BSCLE = BRY(KFLAG)
+      S1 = S1*CSS(KFLAG)
+      S2 = S2*CSS(KFLAG)
+      CS = CSR(KFLAG)
+      DO 70 I=3,N
+        ST = S2
+        S2 = CK*S2 + S1
+        S1 = ST
+        C1 = S2*CS
+        ST = C1
+        C2 = Y(I)
+        IF (KODE.EQ.1) GO TO 60
+        IF (IUF.LT.0) GO TO 60
+        CALL CS1S2(ZN, C1, C2, NW, ASCLE, ALIM, IUF)
+        NZ = NZ + NW
+        SC1 = SC2
+        SC2 = C1
+        IF (IUF.NE.3) GO TO 60
+        IUF = -4
+        S1 = SC1*CSS(KFLAG)
+        S2 = SC2*CSS(KFLAG)
+        ST = SC2
+   60   CONTINUE
+        Y(I) = CSPN*C1 + CSGN*C2
+        CK = CK + RZ
+        CSPN = -CSPN
+        IF (KFLAG.GE.3) GO TO 70
+        C1R = REAL(C1)
+        C1I = AIMAG(C1)
+        C1R = ABS(C1R)
+        C1I = ABS(C1I)
+        C1M = AMAX1(C1R,C1I)
+        IF (C1M.LE.BSCLE) GO TO 70
+        KFLAG = KFLAG + 1
+        BSCLE = BRY(KFLAG)
+        S1 = S1*CS
+        S2 = ST
+        S1 = S1*CSS(KFLAG)
+        S2 = S2*CSS(KFLAG)
+        CS = CSR(KFLAG)
+   70 CONTINUE
+      RETURN
+   80 CONTINUE
+      NZ = -1
+      IF(NW.EQ.(-2)) NZ=-2
+      RETURN
+      END

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libcruft/amos/cairy.f	Sun Apr 27 22:34:17 2008 +0200
@@ -0,0 +1,336 @@
+      SUBROUTINE CAIRY(Z, ID, KODE, AI, NZ, IERR)
+C***BEGIN PROLOGUE  CAIRY
+C***DATE WRITTEN   830501   (YYMMDD)
+C***REVISION DATE  890801   (YYMMDD)
+C***CATEGORY NO.  B5K
+C***KEYWORDS  AIRY FUNCTION,BESSEL FUNCTIONS OF ORDER ONE THIRD
+C***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
+C***PURPOSE  TO COMPUTE AIRY FUNCTIONS AI(Z) AND DAI(Z) FOR COMPLEX Z
+C***DESCRIPTION
+C
+C         ON KODE=1, CAIRY COMPUTES THE COMPLEX AIRY FUNCTION AI(Z) OR
+C         ITS DERIVATIVE DAI(Z)/DZ ON ID=0 OR ID=1 RESPECTIVELY. ON
+C         KODE=2, A SCALING OPTION CEXP(ZTA)*AI(Z) OR CEXP(ZTA)*
+C         DAI(Z)/DZ IS PROVIDED TO REMOVE THE EXPONENTIAL DECAY IN
+C         -PI/3.LT.ARG(Z).LT.PI/3 AND THE EXPONENTIAL GROWTH IN
+C         PI/3.LT.ABS(ARG(Z)).LT.PI WHERE ZTA=(2/3)*Z*CSQRT(Z)
+C
+C         WHILE THE AIRY FUNCTIONS AI(Z) AND DAI(Z)/DZ ARE ANALYTIC IN
+C         THE WHOLE Z PLANE, THE CORRESPONDING SCALED FUNCTIONS DEFINED
+C         FOR KODE=2 HAVE A CUT ALONG THE NEGATIVE REAL AXIS.
+C         DEFINITIONS AND NOTATION ARE FOUND IN THE NBS HANDBOOK OF
+C         MATHEMATICAL FUNCTIONS (REF. 1).
+C
+C         INPUT
+C           Z      - Z=CMPLX(X,Y)
+C           ID     - ORDER OF DERIVATIVE, ID=0 OR ID=1
+C           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
+C                    KODE= 1  RETURNS
+C                             AI=AI(Z)                ON ID=0 OR
+C                             AI=DAI(Z)/DZ            ON ID=1
+C                        = 2  RETURNS
+C                             AI=CEXP(ZTA)*AI(Z)       ON ID=0 OR
+C                             AI=CEXP(ZTA)*DAI(Z)/DZ   ON ID=1 WHERE
+C                             ZTA=(2/3)*Z*CSQRT(Z)
+C
+C         OUTPUT
+C           AI     - COMPLEX ANSWER DEPENDING ON THE CHOICES FOR ID AND
+C                    KODE
+C           NZ     - UNDERFLOW INDICATOR
+C                    NZ= 0   , NORMAL RETURN
+C                    NZ= 1   , AI=CMPLX(0.0,0.0) DUE TO UNDERFLOW IN
+C                              -PI/3.LT.ARG(Z).LT.PI/3 ON KODE=1
+C           IERR   - ERROR FLAG
+C                    IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
+C                    IERR=1, INPUT ERROR   - NO COMPUTATION
+C                    IERR=2, OVERFLOW      - NO COMPUTATION, REAL(ZTA)
+C                            TOO LARGE WITH KODE=1.
+C                    IERR=3, CABS(Z) LARGE      - COMPUTATION COMPLETED
+C                            LOSSES OF SIGNIFCANCE BY ARGUMENT REDUCTION
+C                            PRODUCE LESS THAN HALF OF MACHINE ACCURACY
+C                    IERR=4, CABS(Z) TOO LARGE  - NO COMPUTATION
+C                            COMPLETE LOSS OF ACCURACY BY ARGUMENT
+C                            REDUCTION
+C                    IERR=5, ERROR              - NO COMPUTATION,
+C                            ALGORITHM TERMINATION CONDITION NOT MET
+C
+C
+C***LONG DESCRIPTION
+C
+C         AI AND DAI ARE COMPUTED FOR CABS(Z).GT.1.0 FROM THE K BESSEL
+C         FUNCTIONS BY
+C
+C            AI(Z)=C*SQRT(Z)*K(1/3,ZTA) , DAI(Z)=-C*Z*K(2/3,ZTA)
+C                           C=1.0/(PI*SQRT(3.0))
+C                           ZTA=(2/3)*Z**(3/2)
+C
+C         WITH THE POWER SERIES FOR CABS(Z).LE.1.0.
+C
+C         IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
+C         MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z IS LARGE, LOSSES
+C         OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. CONSEQUENTLY, IF
+C         THE MAGNITUDE OF ZETA=(2/3)*Z**1.5 EXCEEDS U1=SQRT(0.5/UR),
+C         THEN LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR
+C         FLAG IERR=3 IS TRIGGERED WHERE UR=R1MACH(4)=UNIT ROUNDOFF.
+C         ALSO, IF THE MAGNITUDE OF ZETA IS LARGER THAN U2=0.5/UR, THEN
+C         ALL SIGNIFICANCE IS LOST AND IERR=4. IN ORDER TO USE THE INT
+C         FUNCTION, ZETA MUST BE FURTHER RESTRICTED NOT TO EXCEED THE
+C         LARGEST INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF ZETA
+C         MUST BE RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2,
+C         AND U3 ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE
+C         PRECISION ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE
+C         PRECISION ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMIT-
+C         ING IN THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT THE MAG-
+C         NITUDE OF Z CANNOT EXCEED 3.1E+4 IN SINGLE AND 2.1E+6 IN
+C         DOUBLE PRECISION ARITHMETIC. THIS ALSO MEANS THAT ONE CAN
+C         EXPECT TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES,
+C         NO DIGITS IN SINGLE PRECISION AND ONLY 7 DIGITS IN DOUBLE
+C         PRECISION ARITHMETIC. SIMILAR CONSIDERATIONS HOLD FOR OTHER
+C         MACHINES.
+C
+C         THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
+C         BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
+C         ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
+C         SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
+C         ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
+C         ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
+C         CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
+C         HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
+C         ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
+C         SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
+C         THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
+C         0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
+C         THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
+C         COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
+C         BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
+C         COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
+C         MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
+C         THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
+C         OR -PI/2+P.
+C
+C***REFERENCES  HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
+C                 AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
+C                 COMMERCE, 1955.
+C
+C               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
+C                 AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
+C
+C               A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
+C                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
+C                 1018, MAY, 1985
+C
+C               A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
+C                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
+C                 MATH. SOFTWARE, 1986
+C
+C***ROUTINES CALLED  CACAI,CBKNU,I1MACH,R1MACH
+C***END PROLOGUE  CAIRY
+      COMPLEX AI, CONE, CSQ, CY, S1, S2, TRM1, TRM2, Z, ZTA, Z3
+      REAL AA, AD, AK, ALIM, ATRM, AZ, AZ3, BK, CK, COEF, C1, C2, DIG,
+     * DK, D1, D2, ELIM, FID, FNU, RL, R1M5, SFAC, TOL, TTH, ZI, ZR,
+     * Z3I, Z3R, R1MACH, BB, ALAZ
+      INTEGER ID, IERR, IFLAG, K, KODE, K1, K2, MR, NN, NZ, I1MACH
+      DIMENSION CY(1)
+      DATA TTH, C1, C2, COEF /6.66666666666666667E-01,
+     * 3.55028053887817240E-01,2.58819403792806799E-01,
+     * 1.83776298473930683E-01/
+      DATA  CONE / (1.0E0,0.0E0) /
+C***FIRST EXECUTABLE STATEMENT  CAIRY
+      IERR = 0
+      NZ=0
+      IF (ID.LT.0 .OR. ID.GT.1) IERR=1
+      IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1
+      IF (IERR.NE.0) RETURN
+      AZ = CABS(Z)
+      TOL = AMAX1(R1MACH(4),1.0E-18)
+      FID = FLOAT(ID)
+      IF (AZ.GT.1.0E0) GO TO 60
+C-----------------------------------------------------------------------
+C     POWER SERIES FOR CABS(Z).LE.1.
+C-----------------------------------------------------------------------
+      S1 = CONE
+      S2 = CONE
+      IF (AZ.LT.TOL) GO TO 160
+      AA = AZ*AZ
+      IF (AA.LT.TOL/AZ) GO TO 40
+      TRM1 = CONE
+      TRM2 = CONE
+      ATRM = 1.0E0
+      Z3 = Z*Z*Z
+      AZ3 = AZ*AA
+      AK = 2.0E0 + FID
+      BK = 3.0E0 - FID - FID
+      CK = 4.0E0 - FID
+      DK = 3.0E0 + FID + FID
+      D1 = AK*DK
+      D2 = BK*CK
+      AD = AMIN1(D1,D2)
+      AK = 24.0E0 + 9.0E0*FID
+      BK = 30.0E0 - 9.0E0*FID
+      Z3R = REAL(Z3)
+      Z3I = AIMAG(Z3)
+      DO 30 K=1,25
+        TRM1 = TRM1*CMPLX(Z3R/D1,Z3I/D1)
+        S1 = S1 + TRM1
+        TRM2 = TRM2*CMPLX(Z3R/D2,Z3I/D2)
+        S2 = S2 + TRM2
+        ATRM = ATRM*AZ3/AD
+        D1 = D1 + AK
+        D2 = D2 + BK
+        AD = AMIN1(D1,D2)
+        IF (ATRM.LT.TOL*AD) GO TO 40
+        AK = AK + 18.0E0
+        BK = BK + 18.0E0
+   30 CONTINUE
+   40 CONTINUE
+      IF (ID.EQ.1) GO TO 50
+      AI = S1*CMPLX(C1,0.0E0) - Z*S2*CMPLX(C2,0.0E0)
+      IF (KODE.EQ.1) RETURN
+      ZTA = Z*CSQRT(Z)*CMPLX(TTH,0.0E0)
+      AI = AI*CEXP(ZTA)
+      RETURN
+   50 CONTINUE
+      AI = -S2*CMPLX(C2,0.0E0)
+      IF (AZ.GT.TOL) AI = AI + Z*Z*S1*CMPLX(C1/(1.0E0+FID),0.0E0)
+      IF (KODE.EQ.1) RETURN
+      ZTA = Z*CSQRT(Z)*CMPLX(TTH,0.0E0)
+      AI = AI*CEXP(ZTA)
+      RETURN
+C-----------------------------------------------------------------------
+C     CASE FOR CABS(Z).GT.1.0
+C-----------------------------------------------------------------------
+   60 CONTINUE
+      FNU = (1.0E0+FID)/3.0E0
+C-----------------------------------------------------------------------
+C     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
+C     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
+C     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
+C     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
+C     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
+C     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
+C     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
+C     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
+C-----------------------------------------------------------------------
+      K1 = I1MACH(12)
+      K2 = I1MACH(13)
+      R1M5 = R1MACH(5)
+      K = MIN0(IABS(K1),IABS(K2))
+      ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0)
+      K1 = I1MACH(11) - 1
+      AA = R1M5*FLOAT(K1)
+      DIG = AMIN1(AA,18.0E0)
+      AA = AA*2.303E0
+      ALIM = ELIM + AMAX1(-AA,-41.45E0)
+      RL = 1.2E0*DIG + 3.0E0
+      ALAZ=ALOG(AZ)
+C-----------------------------------------------------------------------
+C     TEST FOR RANGE
+C-----------------------------------------------------------------------
+      AA=0.5E0/TOL
+      BB=FLOAT(I1MACH(9))*0.5E0
+      AA=AMIN1(AA,BB)
+      AA=AA**TTH
+      IF (AZ.GT.AA) GO TO 260
+      AA=SQRT(AA)
+      IF (AZ.GT.AA) IERR=3
+      CSQ=CSQRT(Z)
+      ZTA=Z*CSQ*CMPLX(TTH,0.0E0)
+C-----------------------------------------------------------------------
+C     RE(ZTA).LE.0 WHEN RE(Z).LT.0, ESPECIALLY WHEN IM(Z) IS SMALL
+C-----------------------------------------------------------------------
+      IFLAG = 0
+      SFAC = 1.0E0
+      ZI = AIMAG(Z)
+      ZR = REAL(Z)
+      AK = AIMAG(ZTA)
+      IF (ZR.GE.0.0E0) GO TO 70
+      BK = REAL(ZTA)
+      CK = -ABS(BK)
+      ZTA = CMPLX(CK,AK)
+   70 CONTINUE
+      IF (ZI.NE.0.0E0) GO TO 80
+      IF (ZR.GT.0.0E0) GO TO 80
+      ZTA = CMPLX(0.0E0,AK)
+   80 CONTINUE
+      AA = REAL(ZTA)
+      IF (AA.GE.0.0E0 .AND. ZR.GT.0.0E0) GO TO 100
+      IF (KODE.EQ.2) GO TO 90
+C-----------------------------------------------------------------------
+C     OVERFLOW TEST
+C-----------------------------------------------------------------------
+      IF (AA.GT.(-ALIM)) GO TO 90
+      AA = -AA + 0.25E0*ALAZ
+      IFLAG = 1
+      SFAC = TOL
+      IF (AA.GT.ELIM) GO TO 240
+   90 CONTINUE
+C-----------------------------------------------------------------------
+C     CBKNU AND CACAI RETURN EXP(ZTA)*K(FNU,ZTA) ON KODE=2
+C-----------------------------------------------------------------------
+      MR = 1
+      IF (ZI.LT.0.0E0) MR = -1
+      CALL CACAI(ZTA, FNU, KODE, MR, 1, CY, NN, RL, TOL, ELIM, ALIM)
+      IF (NN.LT.0) GO TO 250
+      NZ = NZ + NN
+      GO TO 120
+  100 CONTINUE
+      IF (KODE.EQ.2) GO TO 110
+C-----------------------------------------------------------------------
+C     UNDERFLOW TEST
+C-----------------------------------------------------------------------
+      IF (AA.LT.ALIM) GO TO 110
+      AA = -AA - 0.25E0*ALAZ
+      IFLAG = 2
+      SFAC = 1.0E0/TOL
+      IF (AA.LT.(-ELIM)) GO TO 180
+  110 CONTINUE
+      CALL CBKNU(ZTA, FNU, KODE, 1, CY, NZ, TOL, ELIM, ALIM)
+  120 CONTINUE
+      S1 = CY(1)*CMPLX(COEF,0.0E0)
+      IF (IFLAG.NE.0) GO TO 140
+      IF (ID.EQ.1) GO TO 130
+      AI = CSQ*S1
+      RETURN
+  130 AI = -Z*S1
+      RETURN
+  140 CONTINUE
+      S1 = S1*CMPLX(SFAC,0.0E0)
+      IF (ID.EQ.1) GO TO 150
+      S1 = S1*CSQ
+      AI = S1*CMPLX(1.0E0/SFAC,0.0E0)
+      RETURN
+  150 CONTINUE
+      S1 = -S1*Z
+      AI = S1*CMPLX(1.0E0/SFAC,0.0E0)
+      RETURN
+  160 CONTINUE
+      AA = 1.0E+3*R1MACH(1)
+      S1 = CMPLX(0.0E0,0.0E0)
+      IF (ID.EQ.1) GO TO 170
+      IF (AZ.GT.AA) S1 = CMPLX(C2,0.0E0)*Z
+      AI = CMPLX(C1,0.0E0) - S1
+      RETURN
+  170 CONTINUE
+      AI = -CMPLX(C2,0.0E0)
+      AA = SQRT(AA)
+      IF (AZ.GT.AA) S1 = Z*Z*CMPLX(0.5E0,0.0E0)
+      AI = AI + S1*CMPLX(C1,0.0E0)
+      RETURN
+  180 CONTINUE
+      NZ = 1
+      AI = CMPLX(0.0E0,0.0E0)
+      RETURN
+  240 CONTINUE
+      NZ = 0
+      IERR=2
+      RETURN
+  250 CONTINUE
+      IF(NN.EQ.(-1)) GO TO 240
+      NZ=0
+      IERR=5
+      RETURN
+  260 CONTINUE
+      IERR=4
+      NZ=0
+      RETURN
+      END

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libcruft/amos/casyi.f	Sun Apr 27 22:34:17 2008 +0200
@@ -0,0 +1,126 @@
+      SUBROUTINE CASYI(Z, FNU, KODE, N, Y, NZ, RL, TOL, ELIM, ALIM)
+C***BEGIN PROLOGUE  CASYI
+C***REFER TO  CBESI,CBESK
+C
+C     CASYI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY
+C     MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z) IN THE
+C     REGION CABS(Z).GT.MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL RETURN.
+C     NZ.LT.0 INDICATES AN OVERFLOW ON KODE=1.
+C
+C***ROUTINES CALLED  R1MACH
+C***END PROLOGUE  CASYI
+      COMPLEX AK1, CK, CONE, CS1, CS2, CZ, CZERO, DK, EZ, P1, RZ, S2,
+     * Y, Z
+      REAL AA, ACZ, AEZ, AK, ALIM, ARG, ARM, ATOL, AZ, BB, BK, DFNU,
+     * DNU2, ELIM, FDN, FNU, PI, RL, RTPI, RTR1, S, SGN, SQK, TOL, X,
+     * YY, R1MACH
+      INTEGER I, IB, IL, INU, J, JL, K, KODE, KODED, M, N, NN, NZ
+      DIMENSION Y(N)
+      DATA PI, RTPI  /3.14159265358979324E0 , 0.159154943091895336E0 /
+      DATA CZERO, CONE / (0.0E0,0.0E0), (1.0E0,0.0E0) /
+C
+      NZ = 0
+      AZ = CABS(Z)
+      X = REAL(Z)
+      ARM = 1.0E+3*R1MACH(1)
+      RTR1 = SQRT(ARM)
+      IL = MIN0(2,N)
+      DFNU = FNU + FLOAT(N-IL)
+C-----------------------------------------------------------------------
+C     OVERFLOW TEST
+C-----------------------------------------------------------------------
+      AK1 = CMPLX(RTPI,0.0E0)/Z
+      AK1 = CSQRT(AK1)
+      CZ = Z
+      IF (KODE.EQ.2) CZ = Z - CMPLX(X,0.0E0)
+      ACZ = REAL(CZ)
+      IF (ABS(ACZ).GT.ELIM) GO TO 80
+      DNU2 = DFNU + DFNU
+      KODED = 1
+      IF ((ABS(ACZ).GT.ALIM) .AND. (N.GT.2)) GO TO 10
+      KODED = 0
+      AK1 = AK1*CEXP(CZ)
+   10 CONTINUE
+      FDN = 0.0E0
+      IF (DNU2.GT.RTR1) FDN = DNU2*DNU2
+      EZ = Z*CMPLX(8.0E0,0.0E0)
+C-----------------------------------------------------------------------
+C     WHEN Z IS IMAGINARY, THE ERROR TEST MUST BE MADE RELATIVE TO THE
+C     FIRST RECIPROCAL POWER SINCE THIS IS THE LEADING TERM OF THE
+C     EXPANSION FOR THE IMAGINARY PART.
+C-----------------------------------------------------------------------
+      AEZ = 8.0E0*AZ
+      S = TOL/AEZ
+      JL = INT(RL+RL) + 2
+      YY = AIMAG(Z)
+      P1 = CZERO
+      IF (YY.EQ.0.0E0) GO TO 20
+C-----------------------------------------------------------------------
+C     CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF
+C     SIGNIFICANCE WHEN FNU OR N IS LARGE
+C-----------------------------------------------------------------------
+      INU = INT(FNU)
+      ARG = (FNU-FLOAT(INU))*PI
+      INU = INU + N - IL
+      AK = -SIN(ARG)
+      BK = COS(ARG)
+      IF (YY.LT.0.0E0) BK = -BK
+      P1 = CMPLX(AK,BK)
+      IF (MOD(INU,2).EQ.1) P1 = -P1
+   20 CONTINUE
+      DO 50 K=1,IL
+        SQK = FDN - 1.0E0
+        ATOL = S*ABS(SQK)
+        SGN = 1.0E0
+        CS1 = CONE
+        CS2 = CONE
+        CK = CONE
+        AK = 0.0E0
+        AA = 1.0E0
+        BB = AEZ
+        DK = EZ
+        DO 30 J=1,JL
+          CK = CK*CMPLX(SQK,0.0E0)/DK
+          CS2 = CS2 + CK
+          SGN = -SGN
+          CS1 = CS1 + CK*CMPLX(SGN,0.0E0)
+          DK = DK + EZ
+          AA = AA*ABS(SQK)/BB
+          BB = BB + AEZ
+          AK = AK + 8.0E0
+          SQK = SQK - AK
+          IF (AA.LE.ATOL) GO TO 40
+   30   CONTINUE
+        GO TO 90
+   40   CONTINUE
+        S2 = CS1
+        IF (X+X.LT.ELIM) S2 = S2 + P1*CS2*CEXP(-Z-Z)
+        FDN = FDN + 8.0E0*DFNU + 4.0E0
+        P1 = -P1
+        M = N - IL + K
+        Y(M) = S2*AK1
+   50 CONTINUE
+      IF (N.LE.2) RETURN
+      NN = N
+      K = NN - 2
+      AK = FLOAT(K)
+      RZ = (CONE+CONE)/Z
+      IB = 3
+      DO 60 I=IB,NN
+        Y(K) = CMPLX(AK+FNU,0.0E0)*RZ*Y(K+1) + Y(K+2)
+        AK = AK - 1.0E0
+        K = K - 1
+   60 CONTINUE
+      IF (KODED.EQ.0) RETURN
+      CK = CEXP(CZ)
+      DO 70 I=1,NN
+        Y(I) = Y(I)*CK
+   70 CONTINUE
+      RETURN
+   80 CONTINUE
+      NZ = -1
+      RETURN
+   90 CONTINUE
+      NZ=-2
+      RETURN
+      END

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/libcruft/amos/cbesh.f	Sun Apr 27 22:34:17 2008 +0200
@@ -0,0 +1,331 @@
+      SUBROUTINE CBESH(Z, FNU, KODE, M, N, CY, NZ, IERR)
+C***BEGIN PROLOGUE  CBESH
+C***DATE WRITTEN   830501   (YYMMDD)
+C***REVISION DATE  890801   (YYMMDD)
+C***CATEGORY NO.  B5K
+C***KEYWORDS  H-BESSEL FUNCTIONS,BESSEL FUNCTIONS OF COMPLEX ARGUMENT,
+C             BESSEL FUNCTIONS OF THIRD KIND,HANKEL FUNCTIONS
+C***AUTHOR  AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
+C***PURPOSE  TO COMPUTE THE H-BESSEL FUNCTIONS OF A COMPLEX ARGUMENT
+C***DESCRIPTION
+C
+C         ON KODE=1, CBESH COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
+C         HANKEL (BESSEL) FUNCTIONS CY(J)=H(M,FNU+J-1,Z) FOR KINDS M=1
+C         OR 2, REAL, NONNEGATIVE ORDERS FNU+J-1, J=1,...,N, AND COMPLEX
+C         Z.NE.CMPLX(0.0E0,0.0E0) IN THE CUT PLANE -PI.LT.ARG(Z).LE.PI.
+C         ON KODE=2, CBESH COMPUTES THE SCALED HANKEL FUNCTIONS
+C
+C         CY(I)=H(M,FNU+J-1,Z)*EXP(-MM*Z*I)       MM=3-2M,      I**2=-1.
+C
+C         WHICH REMOVES THE EXPONENTIAL BEHAVIOR IN BOTH THE UPPER
+C         AND LOWER HALF PLANES. DEFINITIONS AND NOTATION ARE FOUND IN
+C         THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS (REF. 1).
+C
+C         INPUT
+C           Z      - Z=CMPLX(X,Y), Z.NE.CMPLX(0.,0.),-PI.LT.ARG(Z).LE.PI
+C           FNU    - ORDER OF INITIAL H FUNCTION, FNU.GE.0.0E0
+C           KODE   - A PARAMETER TO INDICATE THE SCALING OPTION
+C                    KODE= 1  RETURNS
+C                             CY(J)=H(M,FNU+J-1,Z),      J=1,...,N
+C                        = 2  RETURNS
+C                             CY(J)=H(M,FNU+J-1,Z)*EXP(-I*Z*(3-2M))
+C                                  J=1,...,N  ,  I**2=-1
+C           M      - KIND OF HANKEL FUNCTION, M=1 OR 2
+C           N      - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1
+C
+C         OUTPUT
+C           CY     - A COMPLEX VECTOR WHOSE FIRST N COMPONENTS CONTAIN
+C                    VALUES FOR THE SEQUENCE
+C                    CY(J)=H(M,FNU+J-1,Z)  OR
+C                    CY(J)=H(M,FNU+J-1,Z)*EXP(-I*Z*(3-2M))  J=1,...,N
+C                    DEPENDING ON KODE, I**2=-1.
+C           NZ     - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW,
+C                    NZ= 0   , NORMAL RETURN
+C                    NZ.GT.0 , FIRST NZ COMPONENTS OF CY SET TO ZERO
+C                              DUE TO UNDERFLOW, CY(J)=CMPLX(0.0,0.0)
+C                              J=1,...,NZ WHEN Y.GT.0.0 AND M=1 OR
+C                              Y.LT.0.0 AND M=2. FOR THE COMPLMENTARY
+C                              HALF PLANES, NZ STATES ONLY THE NUMBER
+C                              OF UNDERFLOWS.
+C           IERR    -ERROR FLAG
+C                    IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
+C                    IERR=1, INPUT ERROR   - NO COMPUTATION
+C                    IERR=2, OVERFLOW      - NO COMPUTATION, FNU+N-1 TOO
+C                            LARGE OR CABS(Z) TOO SMALL OR BOTH
+C                    IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
+C                            BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
+C                            REDUCTION PRODUCE LESS THAN HALF OF MACHINE
+C                            ACCURACY
+C                    IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
+C                            TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
+C                            CANCE BY ARGUMENT REDUCTION
+C                    IERR=5, ERROR              - NO COMPUTATION,
+C                            ALGORITHM TERMINATION CONDITION NOT MET
+C
+C***LONG DESCRIPTION
+C
+C         THE COMPUTATION IS CARRIED OUT BY THE RELATION
+C
+C         H(M,FNU,Z)=(1/MP)*EXP(-MP*FNU)*K(FNU,Z*EXP(-MP))
+C             MP=MM*HPI*I,  MM=3-2*M,  HPI=PI/2,  I**2=-1
+C
+C         FOR M=1 OR 2 WHERE THE K BESSEL FUNCTION IS COMPUTED FOR THE
+C         RIGHT HALF PLANE RE(Z).GE.0.0. THE K FUNCTION IS CONTINUED
+C         TO THE LEFT HALF PLANE BY THE RELATION
+C
+C         K(FNU,Z*EXP(MP)) = EXP(-MP*FNU)*K(FNU,Z)-MP*I(FNU,Z)
+C         MP=MR*PI*I, MR=+1 OR -1, RE(Z).GT.0, I**2=-1
+C
+C         WHERE I(FNU,Z) IS THE I BESSEL FUNCTION.
+C
+C         EXPONENTIAL DECAY OF H(M,FNU,Z) OCCURS IN THE UPPER HALF Z
+C         PLANE FOR M=1 AND THE LOWER HALF Z PLANE FOR M=2.  EXPONENTIAL
+C         GROWTH OCCURS IN THE COMPLEMENTARY HALF PLANES.  SCALING
+C         BY EXP(-MM*Z*I) REMOVES THE EXPONENTIAL BEHAVIOR IN THE
+C         WHOLE Z PLANE FOR Z TO INFINITY.
+C
+C         FOR NEGATIVE ORDERS,THE FORMULAE
+C
+C               H(1,-FNU,Z) = H(1,FNU,Z)*CEXP( PI*FNU*I)
+C               H(2,-FNU,Z) = H(2,FNU,Z)*CEXP(-PI*FNU*I)
+C                         I**2=-1
+C
+C         CAN BE USED.
+C
+C         IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
+C         MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
+C         LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
+C         CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
+C         LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
+C         IERR=3 IS TRIGGERED WHERE UR=R1MACH(4)=UNIT ROUNDOFF. ALSO
+C         IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
+C         LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
+C         MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
+C         INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
+C         RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
+C         ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
+C         ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
+C         ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
+C         THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
+C         TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
+C         IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
+C         SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
+C
+C         THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
+C         BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
+C         ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
+C         SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
+C         ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
+C         ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
+C         CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
+C         HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
+C         ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
+C         SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
+C         THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
+C         0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
+C         THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
+C         COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
+C         BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
+C         COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
+C         MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
+C         THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
+C         OR -PI/2+P.
+C
+C***REFERENCES  HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
+C                 AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
+C                 COMMERCE, 1955.
+C
+C               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
+C                 BY D. E. AMOS, SAND83-0083, MAY, 1983.
+C
+C               COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
+C                 AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
+C
+C               A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
+C                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
+C                 1018, MAY, 1985
+C
+C               A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
+C                 ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
+C                 MATH. SOFTWARE, 1986
+C
+C***ROUTINES CALLED  CACON,CBKNU,CBUNK,CUOIK,I1MACH,R1MACH
+C***END PROLOGUE  CBESH
+C
+      COMPLEX CY, Z, ZN, ZT, CSGN
+      REAL AA, ALIM, ALN, ARG, AZ, CPN, DIG, ELIM, FMM, FN, FNU, FNUL,
+     * HPI, RHPI, RL, R1M5, SGN, SPN, TOL, UFL, XN, XX, YN, YY, R1MACH,
+     * BB, ASCLE, RTOL, ATOL
+      INTEGER I, IERR, INU, INUH, IR, K, KODE, K1, K2, M,
+     * MM, MR, N, NN, NUF, NW, NZ, I1MACH
+      DIMENSION CY(N)
+C
+      DATA HPI /1.57079632679489662E0/
+C
+C***FIRST EXECUTABLE STATEMENT  CBESH
+      NZ=0
+      XX = REAL(Z)
+      YY = AIMAG(Z)
+      IERR = 0
+      IF (XX.EQ.0.0E0 .AND. YY.EQ.0.0E0) IERR=1
+      IF (FNU.LT.0.0E0) IERR=1
+      IF (M.LT.1 .OR. M.GT.2) IERR=1
+      IF (KODE.LT.1 .OR. KODE.GT.2) IERR=1
+      IF (N.LT.1) IERR=1
+      IF (IERR.NE.0) RETURN
+      NN = N
+C-----------------------------------------------------------------------
+C     SET PARAMETERS RELATED TO MACHINE CONSTANTS.
+C     TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
+C     ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
+C     EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL    AND
+C     EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL       ARE INTERVALS NEAR
+C     UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
+C     RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
+C     DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
+C     FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU
+C-----------------------------------------------------------------------
+      TOL = AMAX1(R1MACH(4),1.0E-18)
+      K1 = I1MACH(12)
+      K2 = I1MACH(13)
+      R1M5 = R1MACH(5)
+      K = MIN0(IABS(K1),IABS(K2))
+      ELIM = 2.303E0*(FLOAT(K)*R1M5-3.0E0)
+      K1 = I1MACH(11) - 1
+      AA = R1M5*FLOAT(K1)
+      DIG = AMIN1(AA,18.0E0)
+      AA = AA*2.303E0
+      ALIM = ELIM + AMAX1(-AA,-41.45E0)
+      FNUL = 10.0E0 + 6.0E0*(DIG-3.0E0)
+      RL = 1.2E0*DIG + 3.0E0
+      FN = FNU + FLOAT(NN-1)
+      MM = 3 - M - M
+      FMM = FLOAT(MM)
+      ZN = Z*CMPLX(0.0E0,-FMM)
+      XN = REAL(ZN)
+      YN = AIMAG(ZN)
+      AZ = CABS(Z)
+C-----------------------------------------------------------------------
+C     TEST FOR RANGE
+C-----------------------------------------------------------------------
+      AA = 0.5E0/TOL
+      BB=FLOAT(I1MACH(9))*0.5E0
+      AA=AMIN1(AA,BB)
+      IF(AZ.GT.AA) GO TO 240
+      IF(FN.GT.AA) GO TO 240
+      AA=SQRT(AA)
+      IF(AZ.GT.AA) IERR=3
+      IF(FN.GT.AA) IERR=3
+C-----------------------------------------------------------------------
+C     OVERFLOW TEST ON THE LAST MEMBER OF THE SEQUENCE
+C-----------------------------------------------------------------------
+      UFL = R1MACH(1)*1.0E+3
+      IF (AZ.LT.UFL) GO TO 220
+      IF (FNU.GT.FNUL) GO TO 90
+      IF (FN.LE.1.0E0) GO TO 70
+      IF (FN.GT.2.0E0) GO TO 60
+      IF (AZ.GT.TOL) GO TO 70
+      ARG = 0.5E0*AZ
+      ALN = -FN*ALOG(ARG)
+      IF (ALN.GT.ELIM) GO TO 220
+      GO TO 70
+   60 CONTINUE
+      CALL CUOIK(ZN, FNU, KODE, 2, NN, CY, NUF, TOL, ELIM, ALIM)
+      IF (NUF.LT.0) GO TO 220
+      NZ = NZ + NUF
+      NN = NN - NUF
+C-----------------------------------------------------------------------
+C     HERE NN=N OR NN=0 SINCE NUF=0,NN, OR -1 ON RETURN FROM CUOIK
+C     IF NUF=NN, THEN CY(I)=CZERO FOR ALL I
+C-----------------------------------------------------------------------
+      IF (NN.EQ.0) GO TO 130
+   70 CONTINUE
+      IF ((XN.LT.0.0E0) .OR. (XN.EQ.0.0E0 .AND. YN.LT.0.0E0 .AND.
+     * M.EQ.2)) GO TO 80
+C-----------------------------------------------------------------------
+C     RIGHT HALF PLANE COMPUTATION, XN.GE.0. .AND. (XN.NE.0. .OR.
+C     YN.GE.0. .OR. M=1)
+C-----------------------------------------------------------------------
+      CALL CBKNU(ZN, FNU, KODE, NN, CY, NZ, TOL, ELIM, ALIM)
+      GO TO 110
+C-----------------------------------------------------------------------
+C     LEFT HALF PLANE COMPUTATION
+C-----------------------------------------------------------------------
+   80 CONTINUE
+      MR = -MM
+      CALL CACON(ZN, FNU, KODE, MR, NN, CY, NW, RL, FNUL, TOL, ELIM,
+     * ALIM)
+      IF (NW.LT.0) GO TO 230
+      NZ=NW
+      GO TO 110
+   90 CONTINUE
+C-----------------------------------------------------------------------
+C     UNIFORM ASYMPTOTIC EXPANSIONS FOR FNU.GT.FNUL
+C-----------------------------------------------------------------------
+      MR = 0
+      IF ((XN.GE.0.0E0) .AND. (XN.NE.0.0E0 .OR. YN.GE.0.0E0 .OR.
+     * M.NE.2)) GO TO 100
+      MR = -MM
+      IF (XN.EQ.0.0E0 .AND. YN.LT.0.0E0) ZN = -ZN
+  100 CONTINUE
+      CALL CBUNK(ZN, FNU, KODE, MR, NN, CY, NW, TOL, ELIM, ALIM)
+      IF (NW.LT.0) GO TO 230
+      NZ = NZ + NW
+  110 CONTINUE
+C-----------------------------------------------------------------------
+C     H(M,FNU,Z) = -FMM*(I/HPI)*(ZT**FNU)*K(FNU,-Z*ZT)
+C
+C     ZT=EXP(-FMM*HPI*I) = CMPLX(0.0,-FMM), FMM=3-2*M, M=1,2
+C-----------------------------------------------------------------------
+      SGN = SIGN(HPI,-FMM)
+C-----------------------------------------------------------------------
+C     CALCULATE EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
+C     WHEN FNU IS LARGE
+C-----------------------------------------------------------------------
+      INU = INT(FNU)
+      INUH = INU/2
+      IR = INU - 2*INUH
+      ARG = (FNU-FLOAT(INU-IR))*SGN
+      RHPI = 1.0E0/SGN
+      CPN = RHPI*COS(ARG)
+      SPN = RHPI*SIN(ARG)
+C     ZN = CMPLX(-SPN,CPN)
+      CSGN = CMPLX(-SPN,CPN)
+C     IF (MOD(INUH,2).EQ.1) ZN = -ZN
+      IF (MOD(INUH,2).EQ.1) CSGN = -CSGN
+      ZT = CMPLX(0.0E0,-FMM)
+      RTOL = 1.0E0/TOL
+      ASCLE = UFL*RTOL
+      DO 120 I=1,NN
+C       CY(I) = CY(I)*ZN
+C       ZN = ZN*ZT
+        ZN=CY(I)
+        AA=REAL(ZN)
+        BB=AIMAG(ZN)
+        ATOL=1.0E0
+        IF (AMAX1(ABS(AA),ABS(BB)).GT.ASCLE) GO TO 125
+          ZN = ZN*CMPLX(RTOL,0.0E0)
+          ATOL = TOL
+  125   CONTINUE
+        ZN = ZN*CSGN
+        CY(I) = ZN*CMPLX(ATOL,0.0E0)
+        CSGN = CSGN*ZT
+  120 CONTINUE
+      RETURN
+  130 CONTINUE
+      IF (XN.LT.0.0E0) GO TO 220
+      RETURN
+  220 CONTINUE
+      IERR=2
+      NZ=0
+      RETURN
+  230 CONTINUE
+      IF(NW.EQ.(-1)) GO TO 220
+      NZ=0
+      IERR=5
+      RETURN
+  240 CONTINUE
+      NZ=0
+      IERR=4
+      RETURN
+      END