  function CenterOfPop()
% Center of population for states based on 2010 Census
% Acquire the data...
A = importdata('Counties.txt','\t');
clc
disp('State  Population  #Counties   theta_f    phi_f      theta_g    phi_g        d')
disp('---------------------------------------------------------------------------------')
% n is the total number of counties...
[n,m] = size(A.data);
% As we scan the data set, currentState houses the 2-letter abbreviation
% of the current state and istart is the row in A.data where the data for the 
% the current state begins
istart = 1;
currentState = A.textdata{2,1};
for k=1:n
    State = A.textdata{k+1,1};
    % Check to see if we have are done stepping through the rows associated
    % with the current state..
    if k==n || ~strcmp(State,currentState)
        % istop is the last row associated with the current state
        if k==n
            istop = n;
        else
            istop = k-1;
        end
        % Get the current state's vector of populations, latitudes and
        % longitudes...
        StatePops = A.data(istart:istop,1);
        StateLats = A.data(istart:istop,7);
        StateLongs= A.data(istart:istop,8);
        % Compute the population..
        TotalStatePop = sum(StatePops);
        % The Flat model...
        theta_f = sum(StatePops.*StateLats)/TotalStatePop;
        phi_f   = sum(StatePops.*StateLongs)/TotalStatePop;
        fprintf('  %s    %8d      %3d     %8.3f  %8.3f',currentState,TotalStatePop,(istop-istart+1),theta_f,phi_f)
        % The Globe model...
        [x,y,z] = SphereToCart(StateLats,StateLongs);
        xBar = sum(x.*StatePops)/TotalStatePop;
        yBar = sum(y.*StatePops)/TotalStatePop;
        zBar = sum(z.*StatePops)/TotalStatePop;
        [theta_g,phi_g] = CartToSphere(xBar,yBar,zBar);
        % The separation of the two centers
        d = dist(theta_g,phi_g,theta_f,phi_f);
        fprintf('    %8.3f  %8.3f   %8d\n',theta_g,phi_g,round(d*5280))
        istart = k;
        currentState = State;
    end
end

  function [x,y,z] = SphereToCart(theta,phi)
% A point on the unit sphere has latitude theta and longitude phi in
% radians.
% (x,y,z) is its location in Cartesian coordinates
theta = theta*pi/180;
phi = phi*pi/180;
x = cos(theta).*cos(phi);
y = cos(theta).*sin(phi);
z = sin(theta);

  function [theta,phi] = CartToSphere(x,y,z)
% A point on the unit sphere has Cartesian coordinates (x,y,z)
% theta is its latitude (radians) and phi is its longitude (radians)
theta = asin(z);
phi   = -acos(x./cos(theta));
theta = 180*theta/pi;
phi   = 180*phi/pi;

  function d = dist(theta1,phi1,theta2,phi2)
% Point P1 on the surface of the Earth has latitude theta1 and longitude
% phi1. Point P2 on the surface of the Earth has latitude theta2 and longitude
% phi2. d is the great circle distance between P1 and P2 in miles
[x1,y1,z1] = SphereToCart(theta1,phi1);
[x2,y2,z2] = SphereToCart(theta2,phi2);
d = 2*asin(sqrt((x1-x2)^2+(y1-y2)^2+(z1-z2)^2)/2)*3950;

