Add additional functon for calculating zeros of Hankel function

See #57 (comment)
The result seems to be more or less identical to the one obtained with
sphbesselh_zeros() in terms of stability -- the obtained values for z, p differ
quite a lot.
One advantage of sphhankel_laplace() is that it will not lead to Inf or Nan, but
only to large numerical errors in the reproduced sound field for orders >85.

Maybe it would be a good idea to try to calculate z directly via an formula
without using the zeros() function.

SFS_general/sphhankel_laplace.m

@@ -0,0 +1,35 @@
+function [z,p] = sphhankel_laplace(order)
+
+hankel = struct(...
+  'zeros',[],...
+  'poles',[],...
+  'scale',[],...
+  'delay',[],...
+  'nominator',[],...
+  'denominator',[]);
+
+hankel = repmat(hankel,order+1,1);
+for n=0:order
+  % nominator
+  hankel(n+1).nominator = zeros(1, n+1);
+  for k=0:n-1
+    hankel(n+1).nominator(k+1) = ...
+      ((2*n-k-1)*(2*n - k))/(2*(n-k))*hankel(n).nominator(k+1);
+  end
+  hankel(n+1).nominator(n+1) = 1;
+end
+
+for n=0:order
+  % flip nominator polynoms
+  hankel(n+1).nominator = hankel(n+1).nominator(end:-1:1);
+  % zeros
+  hankel(n+1).zeros = roots(hankel(n+1).nominator);
+  % denominator
+  hankel(n+1).denominator = zeros(1, n+2);
+  hankel(n+1).denominator(1) = 1;
+  % poles
+  hankel(n+1).poles = roots(hankel(n+1).denominator);
+end
+
+z = hankel(order).zeros;
+p = hankel(order).poles;

SFS_time_domain/driving_functions_imp/driving_function_imp_nfchoa_pw.m

@@ -72,7 +72,8 @@
 
 %% ===== Computation =====================================================
 % Find spherical hankel function zeros
-[z,p] = sphbesselh_zeros(N);
+%[z,p] = sphbesselh_zeros(N);
+[z,p] = sphhankel_laplace(N);
 
 % Get the delay and weighting factors
 if strcmp('2D',dimension)