How is the cieplot 1931 chromaticity diagram calculated in matlab?

Question

I'm completely baffled by the chromaticity diagram calculation in matlab.

It is found here code found here.

The xyz to srgb conversion is standard.. though it seems to be missing an illuminant adaptation.

function [rgb] = xyz2srgb(xyz)
    M = [ 3.2406 -1.5372 -0.4986; -0.9689 1.8758 0.0415; 0.0557 -0.2040 1.0570 ];
    [rows cols ] = size(xyz);
    rgb = M*xyz;
    for c = 1:cols
        for ch = 1:3
            if rgb(ch,c) <= 0.0031308
                rgb(ch,c) = 12.92*rgb(ch,c);
            else
                rgb(ch,c) = 1.055*(rgb(ch,c)^(1.0/2.4)) - 0.055;
            end
            % clip RGB
            if rgb(ch,c) < 0
                rgb(ch,c) = 0;
            elseif rgb(ch,c) > 1
                rgb(ch,c) = 1;
            end
        end
    end
end

What makes no sense to me is the xyY to xyz conversion. It's seen below and steps through the visible spectrum. I haven't been able to find any other implementations like this one.. which seems unique to the patch function in matlab:

w2 = mod(w,N)+1;
a1 = atan2(y(w)-e,x(w)-e); % start angle
a2 = atan2(y(w2)-e,x(w2)-e); % end angle
r1 = ((x(w)-e) ˆ 2 + (y(w)- e) ˆ 2) ˆ 0.5;% start radius
r2 = ((x(w2)-e) ˆ 2 + (y(w2)-e) ˆ 2) ˆ 0.5; % end radius
for c = 1:steps % colourfulness
% patch polygon
xyz(1,1) = e+r1*cos(a1)*c/steps;
xyz(1,2) = e+r1*sin(a1)*c/steps;
xyz(1,3) = 1 - xyz(1,1) - xyz(1,2);
xyz(2,1) = e+r1*cos(a1)*(c-1)/steps;
xyz(2,2) = e+r1*sin(a1)*(c-1)/steps;
xyz(2,3) = 1 - xyz(2,1) - xyz(2,2);
xyz(3,1) = e+r2*cos(a2)*(c-1)/steps;
xyz(3,2) = e+r2*sin(a2)*(c-1)/steps;
xyz(3,3) = 1 - xyz(3,1) - xyz(3,2);
xyz(4,1) = e+r2*cos(a2)*c/steps;
xyz(4,2) = e+r2*sin(a2)*c/steps;
xyz(4,3) = 1 - xyz(4,1) - xyz(4,2);