R Under development (unstable) (2024-07-28 r86931 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(Bessel)
> 
> ### Replicate some testing "ideas" from  zqcbh.f (TOMS 644 test program)
> 
> ##  zqcbh generates sequences of H Bessel functions for kind 2 from
> ##  zbesh and checks them against analytic continuation formulas
> ##  in the (z,fnu) space:
> ##
> ##  kode = 1 tests (analytic continuation formulae, i**2 = -1):
> ##
> ##  H(fnu,2,z)= -exp(i*pi*fnu) * H(fnu,1,-z),    -pi < arg(z) <= 0
> ##
> ##            = 2*cos(pi*fnu) * H(fnu,2,-z) + exp(i*pi*fnu)*H(fnu,1,-z),
> ##
> ##                                                0 < arg(z) <= pi
> ##  kode = 2 tests for kinds 1 and 2:
> ##
> ##         exp(-i*z) * H(fnu,1,z) = [exp(-i*z)* H(fnu,1,z)]
> ##
> ##         exp( i*z) * H(fnu,2,z) = [exp( i*z)* H(fnu,2,z)]
> ##
> ##  where the left side is computed with kode = 1 and the right side
> ##  with kode = 2.
> 
> N <- 100 # speed ...
> set.seed(1)
> for(nu in round(sort(rlnorm(20)), 2)) {
+     cat("nu =", format(nu),"  ")
+     for(n in 1:5) {
+         cat(".")
+         z <- complex(re = rnorm(N),
+                      im = rnorm(N))
+         b1  <- BesselH(1, z, nu=nu)
+         b1. <- BesselH(1, z, nu=nu, expon.scaled = TRUE)
+         b2  <- BesselH(2, z, nu=nu)
+         b2. <- BesselH(2, z, nu=nu, expon.scaled = TRUE)
+         exp.i.pn <- exp(1i*pi*nu)
+         b2.alt <-
+             ifelse(0 < Arg(z),
+                    2*cos(pi*nu) * BesselH(2,-z,nu) + exp.i.pn * BesselH(1,-z,nu),
+                    ## else Arg(z) <= 0:
+                    -exp.i.pn * BesselH(1,-z,nu))
+         stopifnot(all.equal(b2,  b2.alt,            tol = 30e-16),
+                   all.equal(exp(-1i * z) * b1, b1., tol = 80e-16),
+                   all.equal(exp( 1i * z) * b2, b2., tol = 80e-16))
+     }
+     cat("\n")
+ }
nu = 0.11   .....
nu = 0.43   .....
nu = 0.44   .....
nu = 0.53   .....
nu = 0.54   .....
nu = 0.74   .....
nu = 0.96   .....
nu = 0.98   .....
nu = 1.2   .....
nu = 1.39   .....
nu = 1.48   .....
nu = 1.63   .....
nu = 1.78   .....
nu = 1.81   .....
nu = 2.09   .....
nu = 2.27   .....
nu = 2.57   .....
nu = 3.08   .....
nu = 4.53   .....
nu = 4.93   .....
> 
> 
> cat('Time elapsed: ', proc.time(),'\n') # "stats"
Time elapsed:  1.53 0.17 1.68 NA NA 
> 
> if(!interactive()) warnings()
> 
> proc.time()
   user  system elapsed 
   1.53    0.17    1.68