function W_ind = upoclust_newdist2(eps)
%UPOCLUST Summary of this function goes here
%   Inputs: eps - Connectivity threshold.

load('ks22f90h25t120_sortedupos.mat'); %Load UPO data.
%distfunction=sprintf('upo_%s',funHandle);
%paramfile=sprintf('param_%s.mat',funHandle);
%load(paramfile); %Load parameters, such as Fourier Modes, how many points to interpolate, etc.

%%%%%TEMP PARAMS%%%%
nph=120; %Interpolate 120 phase points between -L/2 and L/2.
phase=gpuArray(complex(exp((-2i*pi/nph).*(1:nph))));
ExtraPoints=5; %Interpolate 5 points between each trajectory point.
%%%%%%%%%%%%%%%%%%%%
nupo=length(upo);

%eval(sprintf('W=%s(upo(1).a,upo(1).T,A1,A2,ExtraPoints,L,h);',interpfunction));
[~,aa]=ksfmetd2(upo(1).a, L, h, upo(1).T, 1);
W_t=upo_pchip(aa,ExtraPoints);
W=gpuArray(complex(W_t(1:2:end,:)+1i*W_t(2:2:end,:))); %Combine real and imaginary parts.
clear W_t;

W_ind=1;
dW=size(W,1);
sW=size(W,2);

for i=2:20
    [~,aa]=ksfmetd2(upo(i).a, L, h, upo(i).T, 1);
    W_t=upo_pchip(aa,ExtraPoints);
    W_t=gpuArray(complex(W_t(1:2:end,:)+1i*W_t(2:2:end,:))); %After interpolating in all dimensions, combine real and imaginary parts.
    sWt=size(W_t,2);

    for nWt=1:sWt %For all points in period of W_t.

        minD=gpuArray(Inf);
        Wt=W_t(:,nWt);
        normWi=sqrt(sum(abs(W).^2));

        dotWiWt=sum(bsxfun(@times,conj(Wt),W),1); %Dot product along every element of W, with Wt.
        ph_ang=abs(acos(bsxfun(@rdivide,bsxfun(@times,dotWiWt,phase'),normWi)));
        [minColVal,minRowIdx]=min(ph_ang); %minColIdx(minRowIdx) gives which point in W.
        [~,minColIdx]=min(minColVal);

        Wi=W(:,minColIdx);
        Wtt=Wt.*phase(minRowIdx(minColIdx));
        minVal=norm(Wi-Wtt);

        rdotWiWt=sum(bsxfun(@times,-Wt,W),1); %Just -Wt in here, because -conj(conj(Wt))=-Wt.
        rph_ang=abs(acos(bsxfun(@rdivide,bsxfun(@times,rdotWiWt,phase'),normWi)));
        [minColVal,minRowIdx]=min(rph_ang); %minColIdx(minRowIdx) gives which point in W.
        [~,minColIdx]=min(minColVal);

        rminVal=norm(W(:,minColIdx)+conj(Wt).*exp(-2i*pi*minRowIdx(minColIdx)/nph));

        minD=min(minD,min(minVal,rminVal));

        minD;

        if minD>eps;
            W=[W, W_t];
            W_ind=[W_ind, i];
            sW=sW+sWt;
            break; %If any min of nWt are larger than eps, then break and move onto next orbit.
        end
        %If you get to the end of this loop without breaking out,
        %then it means that the orbit can be represented by a combination of other orbits in the set.

    end




    %%%%%

end

%save('W_new.mat','eps','funHandle','W_ind')

end