function [E0, Th, Ph, Rho2] = qcsimSinterp(N, Iter, alpha, Val)
% throw at random points on a unit sphere, each point representing a
% fourier component of density. The points have amplitude and location.
% Evolution: 1. keep positins fixed, vary amplitude; 2. keep amplitude
% fixed vary position. Always move in the direction of decreasing energy,
% thugh can also explore drift-diffusion, to get out of minima. For
% evolution pick a random point.

r = -1; % GL coeff in front of second order (nonint) term

d = 1e-6; % "infinitesimal" for numerical derivatives

Cr = .0001; %Gradient step in the rho direction
Ct = .3; %Gradient step in the theta direction
Cp = .3; %Gradient step in the phi direction

clf(1);
indx = 0;
for i = 1:N
    theta = rand*pi;
    if rand < sin(theta)
        indx = indx + 1;
        Th(indx) = theta;
        Ph(indx) = 2*pi*rand;
        Rho2(indx) = .01 + .00001*rand;
        plot3(sin(Th(indx))*cos(Ph(indx)), sin(Th(indx))*sin(Ph(indx)), cos(Th(indx)), '*')
        hold on
    end
end
hold off

M = indx; % initial number of Bragg spots 

% typical distance between the "particles" is 2/sqrt(M)


Rho2old = Rho2;
i = 0; err = 1;
figure(1);
while M > 0 && i < Iter && err > 1e-10;
    i
    indx = ceil(rand*M);
    th = Th(indx);
    ph = Ph(indx);
    rho2 = Rho2(indx);
    Thp = [Th(1:indx-1), Th(indx+1:M)];
    Php = [Ph(1:indx-1), Ph(indx+1:M)];
    Rho2p = [Rho2(1:indx-1), Rho2(indx+1:M)];

    
    alphaij = acos(sin(th)*cos(ph)*sin(Thp).*cos(Php) + sin(th)*sin(ph)*sin(Thp).*sin(Php) + cos(th)*cos(Thp));
    
    E0 = r*rho2 + 2*sum(interp1(alpha,Val,alphaij/2).*Rho2p)*rho2 + Val(1)*rho2^2/2;
    E1r = r*(rho2+d) + 2*sum(interp1(alpha,Val,alphaij/2).*Rho2p)*(rho2+d) + Val(1)*(rho2+d)^2/2; % change in rho direction
    
    rho2new = rho2 - (E1r - E0)/d*Cr;
    
    if rho2new > 0 
        Rho2(indx) = rho2new; %updade Rho
        alphaijt = acos(sin(th+d)*cos(ph)*sin(Thp).*cos(Php) + sin(th+d)*sin(ph)*sin(Thp).*sin(Php) + cos(th+d)*cos(Thp));
        E1t = r*rho2 + 2*sum(interp1(alpha,Val,alphaijt/2).*Rho2p)*rho2 + Val(1)*rho2^2/2; % change in theta direction
        alphaijp = acos(sin(th)*cos(ph+d)*sin(Thp).*cos(Php) + sin(th)*sin(ph+d)*sin(Thp).*sin(Php) + cos(th)*cos(Thp));
        E1p = r*rho2 + 2*sum(interp1(alpha,Val,alphaijp/2).*Rho2p)*rho2 + Val(1)*rho2^2/2; % change in theta direction
       
        Gt = (E1t - E0)/d; % gradients
        Gp = (E1p - E0)/d;
        
        Th(indx) = Th(indx) - Gt*Ct + .0*rand; % update Theta
        Ph(indx) = Ph(indx) - Gp*Cp + .0*rand; % update Phi
               
        
    else  % a fourier component indx has collapes to nothing
        Rho2 = Rho2p;
        Th = Thp;
        Ph = Php;
        M = M-1; % shrink system size
    end
    
    if mod(i,100) == 0 % display every 100 iterations
        figure(1); clf(1); hold on
        rhomax = max(Rho2);
        for m = 1:M
            plot3([0, sin(Th(m)).*cos(Ph(m))], [0, sin(Th(m)).*sin(Ph(m))], [0, cos(Th(m))], ':*', 'Color', [Rho2(m)/rhomax 0 1-Rho2(m)/rhomax])
            plot3([0, sin(pi - Th(m)).*cos(Ph(m)+pi)],[0, sin(pi - Th(m)).*sin(Ph(m)+pi)], [0, cos(pi - Th(m))], ':*', 'Color', [Rho2(m)/rhomax 0 1-Rho2(m)/rhomax])
        end
        axis([-1, 1,-1, 1,-1, 1])
        
        pause(.001)
    end
    i = i+1;

    if length(Rho2old) ~= length(Rho2)
        err = 1;
    else
        err = norm(Rho2old - Rho2), M
    end
    
    Rho2old = Rho2;
end


% caclulate final enerrgy
E0= 0;
for indx = 1:M

    th = Th(indx);
    ph = Ph(indx);
    rho2 = Rho2(indx);
    Thp = [Th(1:indx-1), Th(indx+1:M)];
    Php = [Ph(1:indx-1), Ph(indx+1:M)];
    Rho2p = [Rho2(1:indx-1), Rho2(indx+1:M)];

    indx
    alphaij = acos(sin(th)*cos(ph)*sin(Thp).*cos(Php) + sin(th)*sin(ph)*sin(Thp).*sin(Php) + cos(th)*cos(Thp));
    alphaij*180/pi
    
    E0 = E0+ r*rho2 + sum(interp1(alpha,Val,alphaij/2).*Rho2p)*rho2 + Val(1)*rho2^2/2;
end
E0
Rho2

disp('perfect icosa')
1/4/((5* interp1(alpha,Val, 63.4349/180*pi/2) + Val(1)/2)/6)

end


function u_ = U(alpha)
    

A = 10;
B = 15;
    
u_ = A + B*cos(2*alpha); % interaction at 0 and pi is the same. alpha is angle between q's, NOT half-angle
    
end
