% fit_Bessel.m

%

% This MATLAB program fits a function f(x) that satisfies the

% boundary conditions f(R_bc) = 0 to a Bessel function expansion.

% The user enters vectors containing the values of r and

% f(r) at grid points, and then computes the expansion in

% terms of the Bessel functions.  This program uses an

% expansion in the first five Bessel functions of the first

% kind of order zero.

%

% K. Beers. MIT ChE. 7/31/2002

function [iflag,coeff] = fit_Bessel(r,f,R_bc,iplot,fun_name);

iflag = 0;

J_roots_guess = [2.4048, 5.5201, 8.6537, 11.7915, 14.9309];

% First, we need to find the first leading zeros of

% the Bessel function.

i_truncate = length(J_roots_guess);

J_roots = zeros(i_truncate,1);

for k=1:i_truncate

    J_roots(k) = fzero('besselj(0,x)',J_roots_guess(k));

end

% Now, we evaluate the coefficients by integration.

coeff = zeros(i_truncate,1);

for k=1:i_truncate

    var1 = 2 / (R_bc * besselj(1,J_roots(k)))^2;

    J_func = besselj(0,J_roots(k)*r/R_bc);

    integrand = J_func.*f.*r;

    coeff(k) = var1 * trapz(r,integrand);

end

% We now calculate the representation of f(r) as a linear

% combination of Bessel functions.

f_approx = 0*r;

for k = 1:i_truncate

    f_approx = f_approx + ...

        coeff(k)*besselj(0,J_roots(k)*r/R_bc);

end

% Now, we plot the function f(r)and its Bessel function representation.

if(iplot)

    figure;

    subplot(2,1,1);

    plot(r,f);

    hold on;

    plot(r,f_approx,'-.');

    legend('exact','Bessel expansion');

    title('f(r) and its expansion in J_0(r) - five terms');

    xlabel('r');

    ylabel('f(r)');

    %

    subplot(2,1,2);

    f_error = f - f_approx;

    plot(r,f_error);

    xlabel('r');

    ylabel('Error = exact f(r) - approx f(r)');

    if(nargin==5)

        gtext(['f(r) = ', fun_name]);

    end

end

iflag = 1;

return;