%% YOURNAME (DUEDATE)
%
% *DO NOT CHANGE ANYTHING IN THIS SCRIPT AT ALL!*
%
% (Besides your name and due date above as needed.)
%
% Instead create a new function file SiWhile.m that
% contains the requested function SiWhile.  Then use this
% script to test whether your function works correctly.
%
%%
% <latex>
% % for Octave, or for Matlab published to LaTeX for proclatex
% \input{../setup.tex}
% </latex>
%
%% FOR, IF, WHILE, Exercise 11

if ~exist('___code___','var') ; clear ; end
format compact
more off


%% EXERCISE DESCRIPTION:
%
% Create a Matlab function SiWhile(x) that evaluates the
% sine integral of x.
%
% This function should sum the Taylor series to the
% maximum possible accuracy.  It should also be careful
% that the round-off errors in Matlab numbers do not
% produce unacceptably bad values.  If necessary, it
% should instead return the asymptotic expression
% |sign(x)| |pi/2| |-| |cos(x)/x|
%
% Unlike Si, SiWhile should use a |while| loop to do the
% sum.


%% Additional file: SiWhile.m
% <include>SiWhile.m</include>


%% TESTS:


%%  Try x = 1

% set the given x
x=1;
fprintf('x: %.2f\n',x)
% https://keisan.casio.com/exec/system/1180573427
SiExact=0.9460830703671830149414;
% run the script
SiValue=SiWhile(x);
% find the signed error
SiError=SiValue-SiExact
% print out the results
fprintf('Value:%19.16f\n',SiValue)
fprintf('Exact:%19.16f\n',SiExact)
fprintf('Error:%19.16f\n',SiError)
if abs(SiError) > 0
    fprintf('*** The accuracy is less than expected!')
end


%% Try x = 5 pi

% set the given x
x=5*pi;
fprintf('x: %.2f\n',x)
% https://keisan.casio.com/exec/system/1180573427
SiExact=1.63396484610283521331;
% run the script
SiValue=SiWhile(x);
% find the signed error
SiError=SiValue-SiExact
% print out the results
fprintf('Value:%19.16f\n',SiValue)
fprintf('Exact:%19.16f\n',SiExact)
fprintf('Error:%19.16f\n',SiError)
disp('Error in the largest term x^15/(15!*15)*eps:')
fprintf('      %19.16f\n',x^15/(factorial(15)*15)*eps)
if abs(SiError) > 7e-12
    fprintf('*** The accuracy is less than expected!')
end


%% Try x = 34

% set the given x
x=34;
fprintf('x: %.2f\n',x)
% https://keisan.casio.com/exec/system/1180573427
SiExact=1.595256185182468624967;
% run the script
SiValue=SiWhile(x);
% find the signed error
SiError=SiValue-SiExact
% print out the results
fprintf('Value:%19.16f\n',SiValue)
fprintf('Exact:%19.16f\n',SiExact)
fprintf('Error:%19.16f\n',SiError)
disp('Error in the largest term x^33/(33!*33)*eps:')
fprintf('      %19.16f\n',x^33/(factorial(33)*33)*eps)
if abs(SiError) > 1e-04
    fprintf('*** The accuracy is less than expected!')
end


%% Try x = 34.6

% set the given x
x=34.6;
fprintf('x: %.2f\n',x)
% https://keisan.casio.com/exec/system/1180573427
SiExact=1.599659517262426255204;
% run the script
SiValue=SiWhile(x);
% find the signed error
SiError=SiValue-SiExact
% print out the results
fprintf('Value:%19.16f\n',SiValue)
fprintf('Exact:%19.16f\n',SiExact)
fprintf('Error:%19.16f\n',SiError)
disp('Error in the largest term x^33/(33!*33)*eps:')
fprintf('      %19.16f\n',x^33/(factorial(33)*33)*eps)
if abs(SiError) > 3e-06
    fprintf('*** The accuracy is less than expected!')
end


%% Try x = 34.7 (after the cut-off)

% set the given x
x=34.7;
fprintf('x: %.2f\n',x)
% https://keisan.casio.com/exec/system/1180573427
SiExact=1.599393177611930938636;
% run the script
SiValue=SiWhile(x);
% find the signed error
SiError=SiValue-SiExact
% print out the results
fprintf('Value:%19.16f\n',SiValue)
fprintf('Exact:%19.16f\n',SiExact)
fprintf('Error:%19.16f\n',SiError)
if SiValue ~= pi/2-cos(x)/x
    fprintf('*** The value is not the asymptotic one!\n')
end
if abs(SiError) > 8e-5
    fprintf('*** The accuracy is less than expected!')
end


%% Try x = 36 (after the cut-off)

% set the given x
x=36;
fprintf('x: %.2f\n',x)
% https://keisan.casio.com/exec/system/1180573427
SiExact=1.575107209556777000758;
% run the script
SiValue=SiWhile(x);
% find the signed error
SiError=SiValue-SiExact
% print out the results
fprintf('Value:%19.16f\n',SiValue)
fprintf('Exact:%19.16f\n',SiExact)
fprintf('Error:%19.16f\n',SiError)
if SiValue ~= pi/2-cos(x)/x
    fprintf('*** The value is not the asymptotic one!\n')
end
if abs(SiError) > 8e-04
    fprintf('*** The accuracy is less than expected!')
end


%% Try x = 100 (after the cut-off)

% set the given x
x=100;
fprintf('x: %.2f\n',x)
% https://keisan.casio.com/exec/system/1180573427
SiExact=1.562225466889056293352;
% run the script
SiValue=SiWhile(x);
% find the signed error
SiError=SiValue-SiExact
% print out the results
fprintf('Value:%19.16f\n',SiValue)
fprintf('Exact:%19.16f\n',SiExact)
fprintf('Error:%19.16f\n',SiError)
if SiValue ~= pi/2-cos(x)/x
    fprintf('*** The value is not the asymptotic one!\n')
end
if abs(SiError) > 6e-05
    fprintf('*** The accuracy is less than expected!')
end


%% Test negative x values

for x = [1 10 100]
    SiValue=SiWhile(x);
    SiValue2=SiWhile(-x);
    if SiValue2 ~= -SiValue
        fprintf(...
            '*** Si is not antisymmetric at x = %.0f:',x)
        fprintf('    -Si(x):f19.16',-SiValue)
        fprintf('    Si(-x):f19.16',SiValue2)
        fprintf('     Error:f19.16',SiValue2+SiValue)
    end
end
