function GS = qc3d_var_interp_Plot(Q, alpha, Val,fcc, C)

% need to run as GS = qc3d_var_interp_Plot(Qp05, alpha, Valp05/3,Val6p05, C);
% with data taken from 3d_T=0.05and0.1.mat

% val(NQ, Na) values of 4th order vertex. alpha runs from 0 to pi/2 (latter
% corresponsd to angle pi between vectors
% (1/N)sum_{j\ne 0} u(alpha_{0j})+u_0/2N
% non-electroni repulsion term gives C(1 - 1/(2N)) contribution to
% denominator

 gm = (sqrt(5) + 1)/2;
 alphaicosa = asin(1/sqrt(gm^2 + 1));
d = 1e-6; % "infinitesimal" for numerical derivatives

NQ = length(Q);
Na = length(alpha);

D = zeros(NQ, 11);

%C = 10:80;
NC = length(C);

GS = zeros(NQ, NC);


% original coloring
for iC = 1:NC
    for iQ = 1:NQ
        % cubic; N=3. 2 neighbors with alpha=pi/2;
        N = 3;
        D(iQ,1) = (2* interp1(alpha,Val(iQ,:), pi/4) + Val(iQ, 1)/2)/3 + (1-1/2/N)*C(iC);
        
        % FCC lattice (BCC reciprocal), N=4. 3 neighbors with alpha=acos(1/3)
        N = 4;
        al = acos(1/3);
        
        D(iQ,2) = (3* interp1(alpha,Val(iQ,:), al/2) + Val(iQ,1)/2)/4 + (1-1/2/N)*C(iC) -4/16*abs(fcc(iQ) + C(iC));
        
        % BCC lattice (FCC reciprocal), N=6. 4 neighbors with alpha=pi/3, 1 with alpha=pi/2;
        N = 6;
        
        D(iQ,3) = (4* interp1(alpha,Val(iQ,:), pi/6) + interp1(alpha,Val(iQ,:), pi/4) + Val(iQ,1)/2)/6 + (1-1/2/N)*C(iC);
        
        % icosahedral 1 in momentum space, N=6 (vertices). 5 neighbors with alppha = 63.4o
        N = 6;
        
        D(iQ,4) = (5* interp1(alpha,Val(iQ,:), alphaicosa) + Val(iQ,1)/2)/6 + (1-1/2/N)*C(iC);
        
        % dodecahedral in momentum space, N=10, 3 neighbors with alpha 41.8o, 6 neighbors with alpha 70.5o
        N = 10;
        
        D(iQ,5) = (3*interp1(alpha,Val(iQ,:), 41.8/180*pi/2) + 6*interp1(alpha,Val(iQ,:), 70.5/180*pi/2)+ Val(iQ,1)/2)/10 + (1-1/2/N)*C(iC);
        
        % stripe - lamellar; N = 1
        N = 1;
        D(iQ,6) =  (0* min(Val(iQ,:)) + Val(iQ,1)/2)/N + (1-1/2/N)*C(iC);
        
        % columnar-A; N = 2
        N = 2;
        
        D(iQ,7) = (min(Val(iQ,:)) + Val(iQ,1)/2)/2 + (1-1/2/N)*C(iC);
        
        % icosahedral 2 in momentum space, N=15 (edges). 4 neighbors with
        % alppha = 60o, 4 neighbors with alpha = 72o, 4 neighbors with alpha =
        % 36o, and 2 neighbors with alpha = 90.
        N = 15;
        
        D(iQ,8) = (4* interp1(alpha,Val(iQ,:), pi/6) + 4*interp1(alpha,Val(iQ,:), pi/5) ...
            + 4*interp1(alpha,Val(iQ,:), pi/10) + 2*interp1(alpha,Val(iQ,:), pi/4) + Val(iQ,1)/2)/N + (1-1/2/N)*C(iC);
        
        % rombohedral; N = 3
        N = 3;
        
        D(iQ,9) = (2* min(Val(iQ,:)) + Val(iQ,1)/2)/N + (1-1/2/N)*C(iC);
        
        % icosahedral -1 in momentum space, N=5 (dropped 1 vertex). 4 neighbors with alppha = 63.4o
        N = 5;
        
        D(iQ,10) = (4* interp1(alpha,Val(iQ,:), 63.4/180*pi/2) + Val(iQ,1)/2)/N + (1-1/2/N)*C(iC);
        
        
        % columnar-B; N = 3
        N = 3;
        
        D(iQ,11) = (2* interp1(alpha,Val(iQ,:), pi/6) + Val(iQ,1)/2)/N + (1-1/2/N)*C(iC);

        
        if min(D(iQ, :) > 0)
            [~,ind] = sort(D(iQ,:));
            GS(iQ,iC) = ind(1);
        end
        
        if GS(iQ,iC) == 4 % test if icosahedron is distorted
            f = interp1(alpha,Val(iQ,:), [alphaicosa - d,alphaicosa, alphaicosa + d], 'spline');
            der1 = (f(3) - f(1))/2/d;
            der2 = (f(3) - 2*f(2) + f(1))/d^2;
            if  der2 < 0
                GS(iQ, iC) = -1;
            elseif (der2 > 0 && der1 < -2/3*der2)
                GS(iQ, iC) = 4.9;
                
            end
            
        end

% better plotting for production        
%        figure(17); plot3(Q(iQ), C(iC), GS(iQ, iC), 's', 'MarkerSize', GS(iQ, iC)+.01, 'MarkerFaceColor', [GS(iQ, iC)/10,0, 1 - GS(iQ, iC)/10], 'MarkerEdgeColor', 'none' )
%        hold on
                
    end
end
hold off;
