Calculates the Ricatti-Bessel function psi_n(rho)=j_n(rho)*rho and its derivative for n=1:nNmax Parameters: - nNmax: scalar integer number of n in series - rho: column vector [R x 1] (with possibly zero components) arguments of the Ricatti-Bessel function Returns: structure with 2 fields each a matrix [R x nNmax] with values for each rho and n=1..nNmax stRBpsi2.psi: psi_n(rho) stRBpsi2.Dpsi: psi'_n(rho)} This file is part of the SPlaC v1.0 package (copyright 2008) Check the README file for further information
0001 n=1:nNmax; 0002 nm1=0:nNmax; 0003 nu=nm1+0.5; 0004 [f ierr] = besselj(nu,rho); 0005 % f is matrix [R x nNmax+1] of cylindrical Bessel 0006 % J_{n+0.5}(rho), n=0..nNmax 0007 0008 if (max(max(ierr)) > 0) 0009 disp 'Error in besselj in GenRBpsi2' 0010 end 0011 0012 % find and exclude indices where rho=0 0013 ind0 = find(rho==0); % ind0 is possibly empty vector 0014 % replace with 1 to avoid error messages 0015 rho(ind0) = 1; 0016 0017 sq=sqrt((pi/2).*rho); % [R x 1] 0018 sqmat=repmat(sq,1,nNmax+1); % [R x nNmax+1] 0019 f=f.*sqmat; 0020 % f is now matrix of spherical Ricati-Bessel 0021 % psi_n(rho), n=0..nNmax or equivalently psi_{n-1}(rho), n=1..nNmax+1 0022 0023 stRBpsi2.psi=f(:,n+1); 0024 0025 rhomat=repmat(rho,1,nNmax); % [R x nNmax] 0026 nmat=repmat(n,length(rho),1); % [R x nNmax] 0027 0028 % Computes: psi_n'=psi_{n-1} - n psi_n/rho 0029 % and check for loss of precision in sum 0030 stRBpsi2.Dpsi=GenCheckSum2Mat(f(:,n), - nmat.*f(:,n+1)./rhomat,'Dpsi','GenRBpsi2'); 0031 0032 % Replace proper value for cases where rho=0, i.e. psi_n'(0)=0 0033 if (~isempty(ind0)) 0034 stRBpsi2.Dpsi(ind0,:)=zeros(length(ind0),nNmax); 0035 end 0036