function P2()
% Checks out BballDesign, Ave,Var

% Benchmark Ave and Var...
h = 10; r = 10; 
clc
disp('  M     N    p      Ave         Var        Ave         Var')
disp('                 Benchmark   Benchmark    Error       Error')
disp('---------------------------------------------------------------')
for M =  [100 ]
    for N = [25 50 ]
        for p = [10 11 200 ]
          tic
          AYou = Ave(h,r,M,N,p);
          tYouAve = toc;
          tic
          ACVL = AveCVL(h,r,M,N,p);
          tCVLAve = toc;
          tic
          VYou = Var(h,r,M,N,p);
          tYouVar = toc;
          tic
          VCVL = VarCVL(h,r,M,N,p);
          tCVLVar = toc;
          errAve = abs(AYou-ACVL)/ACVL;
          errVar = abs(VYou-VCVL)/VCVL;
          fprintf(' %3d  %3d  %3d    %6.3f      %6.3f    %10.3e   %10.3e\n',...
                      M,N,p,tYouAve/tCVLAve,tYouVar/tCVLVar,errAve,errVar)
        end
    end
end

% Check BballDesign....

tic
disp(' ')
disp(' ')
disp('  M     N    p       h_opt      r_opt           Ave             Var')
disp('-----------------------------------------------------------------------')

M = 45; N = 25; p = 5;
[h,r] = BballDesign(M,N,p);
A = AveCVL(h,r,M,N,p);
V = VarCVL(h,r,M,N,p);
fprintf('%3d   %3d  %3d     %8.4f   %8.4f    %12.8f    %12.8f\n',M,N,p,h,r,A,V)

M = 45; N = 25; p = 20;
[h,r] = BballDesign(M,N,p);
A = AveCVL(h,r,M,N,p);
V = VarCVL(h,r,M,N,p);
fprintf('%3d   %3d  %3d     %8.4f   %8.4f    %12.8f    %12.8f\n',M,N,p,h,r,A,V)

M = 45; N = 25; p = 40;
[h,r] = BballDesign(M,N,p);
A = AveCVL(h,r,M,N,p);
V = VarCVL(h,r,M,N,p);
fprintf('%3d   %3d  %3d     %8.4f   %8.4f    %12.8f    %12.8f\n',M,N,p,h,r,A,V)

M = 25; N = 45; p = 40;
[h,r] = BballDesign(M,N,p);
A = AveCVL(h,r,M,N,p);
V = VarCVL(h,r,M,N,p);
fprintf('%3d   %3d  %3d     %8.4f   %8.4f    %12.8f    %12.8f\n',M,N,p,h,r,A,V)

totalTime = toc;
fprintf('\nTotal Time for Four Designs = %6.2f\n',totalTime)


function [h,r] = BballDesignCVL(M,N,p)

% Initial guess based on playing around...
x0 = [(M+N)/2,(M+N)/2];

z = fminsearch(@(x) GoodLightingCVL(x,M,N,p),x0);
h = z(1); r = z(2);

function z = GoodLightingCVL(x,M,N,p)
h = x(1); r = x(2);
lambda = 1/2; % Parameter to play with. 
z = lambda*VarCVL(h,r,M,N,p) - (1-lambda)*AveCVL(h,r,M,N,p);

        
function mu = AveCVL(h,r,M,N,p)
%Useful x and y values for the tile centers...
xTile = (-M:M-1)'+.5;  eM = ones(2*M,1);
yTile = (-N:N-1)'+.5;  eN = ones(2*N,1);
L = zeros(2*N,2*M);
for k=0:floor(p/2)
     % Location of the k-th light...
     xk = r*cos(2*pi*k/p); yk = r*sin(2*pi*k/p);
     % Matrix of illumination values due to the kth light..
     Lk = 1./(( (eN*xTile' - xk).^2 + (yTile*eM' - yk).^2 + h^2).^(3/2));
     % Exploit symmetry across the x-axis...
     if (rem(p,2)==0 && (k==0 || k==p/2)) || (rem(p,2)==1 && k==0)
        L = L + Lk;
     else
        L = L + Lk + Lk(2*N:-1:1,:);
     end
end
mu = sum(sum(L))*h/(4*M*N);

function sigma = VarCVL(h,r,M,N,p)
xTile = (-M:M-1)'+.5;  eM = ones(2*M,1);
yTile = (-N:N-1)'+.5;  eN = ones(2*N,1);
L = zeros(2*N,2*M);
for k=0:floor(p/2)
     % Location of the k-th light...
     xk = r*cos(2*pi*k/p); yk = r*sin(2*pi*k/p);
     % Matrix of illumination values due to the kth light..
     Lk = 1./(( (eN*xTile' - xk).^2 + (yTile*eM' - yk).^2 + h^2).^(3/2));
     % Exploit symmetry...
     if (rem(p,2)==0 && (k==0 || k==p/2)) || (rem(p,2)==1 && k==0)
        L = L + Lk;
     else
        L = L + Lk + Lk(2*N:-1:1,:);
     end
end
L = h*L;
mu = sum(sum(L))/(4*M*N);
% Exploit symmetry..
sigma = sqrt(sum(sum((L(1:N,:)-mu).^2))/(2*M*N));


% -1 if do not exploit symmetry
% -1 if do not vectorize
% -2, -3 for very unreasonable h, r or if code did not run