% EX1p9abcd.m
% (a) Using 4-digit decimal arithmetic to compute
% the surface area of earth based on formula:
% A = 4 pi r^2.
% (b) Recalculate for r = r+1, & find difference in the areas.
% (c) Calculate diff. in areas using dA = 8*pi*r*h from calculus.
% (d) Redo using 6-digit decimal precision.
% Function required: rd(x,p).m to round numbers to p-digit precision.
% K. Ming Leung

clear all;
r0 = 6370; exact = 8*pi*r0,
p = 4;      % set 4-digit decimal precision
myPi = rd(pi,p),
p4 = rd(4*myPi,p), % 4 is a machine number

r0 = 6370;  % is a machine number
r0s = rd(r0*r0,p),
A0 = rd(p4*r0s,p),
B0 = rd( rd(p4*r0,p)*r0, p ),

r1 = 6370+1;    % are all machine numbers
r1s = rd(r1*r1,p),
A1 = rd(p4*r1s,p),
B1 = rd( rd(p4*r1,p)*r1, p ),

% Changes in the areas
dA = rd(A1-A0,p), 
dB = rd(B1-B0,p),

% Compute using formula from calculus
h = 1;
diffA = rd(rd(rd(8*myPi,p)*r0,p)*h,p),

p = 10;      % set 6-digit decimal precision
myPi = rd(pi,p),
p4 = rd(4*myPi,p), 

r0 = 6370;  % is a machine number
r0s = rd(r0*r0,p);
A0 = rd(p4*r0s,p),
B0 = rd( rd(p4*r0,p)*r0, p ),

r1 = 6370+1;    % are all machine numbers
r1s = rd(r1*r1,p);
A1 = rd(p4*r1s,p),
B1 = rd( rd(p4*r1,p)*r1, p ),

% Changes in the areas
dA = rd(A1-A0,p), 
dB = rd(B1-B0,p),

% Compute using formula from calculus
h = 1;
diffA = rd(rd(rd(8*myPi,p)*r0,p)*h,p),

exactdiffA = 8*pi*r0*h,     % Double-precision result for dA = 8*pi*r0*h
exactdA = 4*pi*(r1^2-r0^2), % Double-precision result for the difference in the areas.

