% Script C1
% Average distance to the Sun.

P = input('Enter P: ');
A = input('Enter A: ');
N = 100000;
x_old = P;
y_old = 0;
d_old = P;
d = 0;
a = pi*sqrt(P*A)*(P+A)/2;
% Compue average distance ...
for k=1:N
    % Get the k-th orbit point...
    tau_k = k*(2*pi/N);
    xk = (P-A)/2 + ((P+A)/2)*cos(tau_k);
    yk = sqrt(P*A)*sin(tau_k);
    % Distance to Sun associated with the kth sector...
    dk = sqrt(xk^2+yk^2);
    % Area of the kth sector...
    s1 = dk;
    s2 = d_old;
    s3 = sqrt((xk-x_old)^2+(yk-y_old)^2);
    s = (s1+s2+s3)/2;
    Area_k = sqrt(s*(s-s1)*(s-s2)*(s-s3));
    % Update the running sum...
    Area_fract = Area_k/a;
    d = d + Area_fract*((dk+d_old)/2);  
    % For the next sector, (xk,yk) is (x_old,y_old)...
    x_old = xk;
    y_old = yk;
    d_old = dk;
end
% Output the results...
fprintf('P = %12.3f\n',P)
fprintf('A = %12.3f\n',A)
fprintf('N = %1d\n\n',N)
fprintf('      d = %12.3f\n',d)
fprintf('(P+A)/2 = %12.3f\n',(P+A)/2)

