empty_space_function.m~

%% F(r) integration

% SobolPoints
N = 2^11; % Number of points to sum over to estimate inner integral
P = sobolset(2); % set dimension

%% Integrate
r = 0:0.1:2.5;
%r = (0.1:0.1:2).*1000;
F = zeros(length(r),1);

Fc = zeros(length(r),1);
Fin = zeros(length(r),1);

% The smaller num_integrals, more larger integral_result
% num_integrals>12  -> integral_result = 0
% increasing 
num_integrals = 128;

% calculate F(r)
parfor i =1:length(r)
    X_full = net(P,N);
    % Outer most summation
    for n=1:num_integrals
        % Calculate front constant
        c = (-1)^(n-1)*(2*r(i))^(2*n)/factorial(n);
        if isnan(c)
            disp('NAN');
        end
        % Integrals summation
        integral_result = InnerIntegral(r(i), N, X_full, n);
        %integral_result = 1;
        F(i) = F(i) + c*integral_result;
        Fc(i) = Fc(i) + c;
        Fin(i) = Fin(i) + integral_result;
        
        fprintf('n = %d\n',n);
    end
end

figure(1);plot(r,F);
%axis([0 2.5 0 2]);
figure(2);plot(r,Fc);
figure(3);plot(r,Fin);