function y = mysinintegral(x);

% mysinintegral.m - double-precision Sine-integral special function
%
% Input:
%       x       3-D array, with entries of the form pi*sqrt(integer)
%               if x does not have this form, mysintintegral(x) will NOT return Si(x)
% Output:
%       y       3-D array, double-precision computation of Si(x) component-by-component
%
% Remark:       make sure matlab is correctly linked to Maple's mfun.m
%
% By Laurent Demanet, 2005

% precomputation of all the matrix elements
maxx = ceil(max(max(max(x)))/pi);
newx = sqrt(0:(maxx^2))*pi;
newy = mfun('Si',newx);

% lookup assignment
y = zeros(size(x));
for m = 1:size(x,1)
    for n = 1:size(x,2)
        for p =1:size(x,3)
            index = 1 + round((x(m,n,p)/pi)^2);
            y(m,n,p) = newy(index);
        end
    end
end

