clear all

% Compare MATALB Bessel function to a customly created Bessel function

% Inputs:
x = 0.1:0.1:2; % input to bessel function
m = 0; %order of bessel function
n = 6;      % Newton-Cotes rule
N = 840;  % Newton-Cotes N  

% Compare with built in MATLAB
J = besselj(m,x);

% Create custom Bessel function
lwr = 0;    %lower limit of integration
uppr= pi;   %upper limit of integration
theta = linspace(0,pi,N+1); 

for ii = 1:length(x)    
    % Create function to be integrated
    A = cos(x(ii)*sin(theta)-m*theta)/pi;    

    %Integrate with newton-cotes method
    J_m(ii) = nc_method(A,lwr,uppr,n);
end

% Difference between MATLAB and custom Bessel
compare = J-J_m;

% eps is the Spacing of floating point numbers in MATLAB.
% Any number smaller than in magnitude than eps is essentially zero
% Try 1 + eps, this is equal to 1.
v_eps = [eps, eps]; 
x_eps = [x(1),x(end)];

figure(1),
    subplot(211),plot(x,J,x,J_m,':')
        title('Comparison between MATLAB and Custom Bessel function')
        legend('MATLAB','Custom')
    subplot(212),plot(x,compare,'-o',x_eps,v_eps,x_eps,-v_eps)
        title('Difference between MATLAB and Custom Bessel')
        legend('difference','+eps','-eps')

%figure(2),plot(theta,A)