# v_besselratioi

## PURPOSE

V_BESSELRATIOI calculate the inverse Bessel function ratio

## SYNOPSIS

function s=v_besselratioi(r,v,p)

## DESCRIPTION

    V_BESSELRATIOI calculate the inverse Bessel function ratio
      Inputs:
          r    value of the Bessel function ratio besseli(v+1,s)./besseli(v,s) (scalar or matrix)
          v    denominator Bessel function order [0] [for now v=0 always]
          p    digits precision <=14 [5]

     Outputs:
          s    value of s such that r=besseli(v+1,s)./besseli(v,s)

     This function implements function VKAPPA from [1] for which a FORTRAN
     implementation is available in [2].

     Refs:
     [1] G. W. Hill, Evaluation and Inversion of the Ratios of Modified Bessel
         Functions, I 1 (x)/I 0 (x) and I 1.5 (x)/I 0.5 (x),
         ACM Transactions on Mathematical Software (TOMS), Vol 7, pp199-208, 1981
     [2] ACM Collected Algorithms (CALGO), http://calgo.acm.org/

## CROSS-REFERENCE INFORMATION

This function calls:

- v_besselratio V_BESSELRATIO calculate the Bessel function ratio besseli(v+1,x)./besseli(v,x)

This function is called by:

## SOURCE CODE

    0001 function s=v_besselratioi(r,v,p)
    0002 %V_BESSELRATIOI calculate the inverse Bessel function ratio
    0003 %  Inputs:
    0004 %      r    value of the Bessel function ratio besseli(v+1,s)./besseli(v,s) (scalar or matrix)
    0005 %      v    denominator Bessel function order [0] [for now v=0 always]
    0006 %      p    digits precision <=14 [5]
    0007 %
    0008 % Outputs:
    0009 %      s    value of s such that r=besseli(v+1,s)./besseli(v,s)
    0010 %
    0011 % This function implements function VKAPPA from [1] for which a FORTRAN
    0012 % implementation is available in [2].
    0013 %
    0014 % Refs:
    0015 % [1] G. W. Hill, Evaluation and Inversion of the Ratios of Modified Bessel
    0016 %     Functions, I 1 (x)/I 0 (x) and I 1.5 (x)/I 0.5 (x),
    0017 %     ACM Transactions on Mathematical Software (TOMS), Vol 7, pp199-208, 1981
    0018 % [2] ACM Collected Algorithms (CALGO), http://calgo.acm.org/
    0019
    0020 %      Copyright (C) Mike Brookes 2017-2018
    0021 %      Version: $Id: v_besselratioi.m 10865 2018-09-21 17:22:45Z dmb $
    0022 %
    0023 %   VOICEBOX is a MATLAB toolbox for speech processing.
    0024 %   Home page: http://www.ee.ic.ac.uk/hp/staff/dmb/voicebox/voicebox.html
    0025 %
    0026 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0027 %   This program is free software; you can redistribute it and/or modify
    0028 %   it under the terms of the GNU Lesser General Public License as published by
    0029 %   the Free Software Foundation; either version 3 of the License, or
    0030 %   (at your option) any later version.
    0031 %
    0032 %   This program is distributed in the hope that it will be useful,
    0033 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
    0034 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    0035 %   GNU Lesser General Public License for more details.
    0036 %
    0037 %   You can obtain a copy of the GNU Lesser General Public License from
    0038 %   https://www.gnu.org/licenses/ .
    0039 %
    0040 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    0041 if nargin<3 || isempty(p)
    0042     p=5;
    0043 else
    0044     p=min(max(ceil(p),1),13);       % digits precision required
    0045 end
    0046 if nargin<2 || isempty(v)
    0047     v=0;                            % default to I1/I0
    0048 end
    0049 sr=size(r);
    0050 r=r(:); % convert to a column vector
    0051 nr=prod(sr);
    0052 if v~=0
    0053     error('v must be zero (for now)');
    0054 end
    0055 s=repmat(-1,nr,1);      % default is -1 for r<0 or r>=1
    0056 y=2./(1-r);
    0057 m1=r>=0 & r<1;          % m1 is mask for valid range
    0058 mn=m1;                  % mn is mask for using Newton iteration
    0059 m=m1 & r<=0.85;         % use inverse taylor series for 0<=r<=0.85
    0060 if any(m)
    0061     rm=r(m);
    0062     mn(m)= rm>=0.642;   % we will need a Newton iteration for 0.642<=r<=0.85
    0063     xm=rm.^2;
    0064     sm=(((rm-5.6076).*rm+5.0797).*rm-4.6494).*y(m).*xm-1;
    0065     s(m)=((((sm.*xm+15.0).*xm+60.0).*xm/360.0+1.0).*xm-2.0).*rm./(xm-1);
    0066 end
    0067 m=m1 & ~m;              % use continued fraction for 0.85<r<1
    0068 if any(m)
    0069     rm=r(m);
    0070     mn(m)=rm<0.95;      % we will need a Newton iteration for 0.85<r<0.95
    0071     ym=y(m);
    0072     mc=rm<0.95;
    0073     if any(mc)
    0074         rm(mc)=(-2326.0.*rm(mc)+4317.5526).*rm(mc) - 2001.035224;
    0075     end
    0076     if any(~mc)
    0077         rm(~mc)=32.0./(120.0*rm(~mc)-131.5+ym(~mc));
    0078     end
    0079     s(m)=(ym+1.0+3.0./(ym-5.0-12.0./(ym-10.0-rm)))*0.25;
    0080 end
    0081 if any(mn)              % we need to do one or two Newton iterations
    0082     rmn=r(mn);
    0083     smn=s(mn);
    0084     ymn=y(mn);
    0085     ymn= ((0.00048*ymn-0.1589).*ymn+0.744).*ymn - 4.2932;   % Gradient approximation ds/dr
    0086     smn=smn+(v_besselratio(smn,0,p+1)-rmn).*ymn;              % do a Newton iteration
    0087     mr=(rmn>=0.75) & (rmn<=0.875);                          % need a second Newton iteration for 0.75<=r<=0.875
    0088     if any(mr)
    0089         smn(mr)=smn(mr)+(v_besselratio(smn(mr),0,p+1)-rmn(mr)).*ymn(mr); % do a second Newton iteration
    0090     end
    0091     s(mn)=smn;                                              % update the main s() array
    0092 end
    0093 s=reshape(s,sr);                % put back into the same shape is the input r
    0094 if ~nargout
    0095     plot(r,s);
    0096     ylabel('x');
    0097     xlabel(sprintf('I_{%d}(x)/I_{%d}(x)',v+1,v));
    0098 end