How to fix these codes so that I can have a more systematic data for vertices, edges, and faces?

This picture below is just an example. The data for vertices, edges, and faces are arranged accordingly. It can easily put them in the Surface Evolver Software to generate a sqaure shape and get the minimal surface. Notice the color.

Sample: https://www.youtube.com/watch?v=6bXOXxd7hAE&t=457s

NOW, This is the problem

These codes will generate a sinusoidal wave curve with a triangular mesh. The data for vertices, edges and faces are not systematic. I want to arranged them properly. Also, the face data here shows that one of its edges is the same all throughout. It must not be the same. Here are the codes:

clear
MeshSize        = 0.04;
LensStr = 'osci';
radius      = 1;
%%
n = ceil(radius/MeshSize); % radial resolution
[X,F] = CreateDiskMesh([0 0 0], radius, n);
[X,F] = CompactMesh(X,F,1e-9,false);
x = X(:,1);
y = X(:,2);
f = 10; %define modulation frequency
A = .2; % modulation amplitude
theta = atan2(y,x);
Z   = A*cos(f*theta);
X(:,3) = Z;
fig = figure();
set(fig, 'Name', LensStr);
ax = axes('Parent',fig);
hold(ax, 'on');
trisurf(F, X(:,1), X(:,2), X(:,3), 'Parent', ax);
view(ax,3);
xlabel('x')
ylabel('y')
zlabel('z')
axis(ax, 'equal');

edges= [F(:,[1,2]);
    F(:,[2,3]);
    F(:,[1,3])];
edges=unique(sort(edges,2),'rows');
%%
function [X,F] = AddMesh(X,F,X1,F1)
F = [F; F1+size(X,1)];
X = [X; X1];
end
%%
function [X,F] = CreateDiskMesh(xyzc, a, n)
% radius of the inner/outer spherical parts
W = allVL1(3, n); % FEX file
% Connectivity
XY0 = W*[0 sqrt(3)/2 sqrt(3)/2;
    0  -0.5      +0.5].' *(1/n);
F0 = delaunay(XY0);
phi = XY0(:,2)./XY0(:,1)*(pi/(2*sqrt(3)));
phi(XY0(:,1)==0) = 0;
r = XY0(:,1)*(2/sqrt(3));
XY0 = r.*[cos(phi),sin(phi)];
XY0 = XY0*a;
xyzc = xyzc(:).';
F=[];
X=[];
for k=0:5
    theta = k*pi/3;
    R = [cos(theta), -sin(theta);
        sin(theta),  cos(theta)];
    XYZ = XY0*R';
    XYZ(:,3) = 0;
    XYZ = XYZ + xyzc;
    Fk = F0;
    [X,F] = AddMesh(X, F, XYZ, Fk);
end
[X,F] = CompactMesh(X,F,1e-9,true);
end
%%
function [V,F] = CompactMesh(V,F,mergetol,orientface)
if mergetol == 0
    [V,I,J] = unique(V, 'rows');
else
    [V,I,J] = uniquetol(V, mergetol, 'DataScale', 1, 'ByRows', true);
end
F = J(F);
if length(I)<length(J)
    remove = F(:,1)==F(:,2) | F(:,2)==F(:,3) | F(:,3)==F(:,1);
    F(remove,:) = [];
end
if orientface
    T = permute(reshape(V(F,:),[],3,3),[3 2 1]);
    N = cross_dim1(T(:,2,:)-T(:,1,:),T(:,3,:)-T(:,2,:));
    N = reshape(N,3,[]);
    Nz = N(3,:);
    keep = Nz > 0.1*max(Nz(:));
    F = F(keep,:);
end
end
function v = allVL1(n, L1, L1ops, MaxNbSol)
% All integer permutations with sum criteria
%
% function v=allVL1(n, L1); OR
% v=allVL1(n, L1, L1opt);
% v=allVL1(n, L1, L1opt, MaxNbSol);
% 
% INPUT
%    n: length of the vector
%    L1: target L1 norm
%    L1ops: optional string ('==' or '<=' or '<')
%           default value is '=='
%    MaxNbSol: integer, returns at most MaxNbSol permutations.
%    When MaxNbSol is NaN, allVL1 returns the total number of all possible
%    permutations, which is useful to check the feasibility before getting
%    the permutations.
% OUTPUT:
%    v: (m x n) array such as: sum(v,2) == L1,
%       (or <= or < depending on L1ops)                            
%       all elements of v is naturel numbers {0,1,...}
%       v contains all (=m) possible combinations
%       v is sorted by sum (L1 norm), then by dictionnary sorting criteria
%    class(v) is same as class(L1) 
% Algorithm:
%    Recursive
% Remark:
%    allVL1(n,L1-n)+1 for natural numbers defined as {1,2,...}
% Example:
%    This function can be used to generate all orders of all
%    multivariable polynomials of degree p in R^n:
%         Order = allVL1(n, p)
% Author: Bruno Luong
% Original, 30/nov/2007
% Version 1.1, 30/apr/2008: Add H1 line as suggested by John D'Errico
%         1.2, 17/may/2009: Possibility to get the number of permutations
%                           alone (set fourth parameter MaxNbSol to NaN)
%         1.3, 16/Sep/2009: Correct bug for number of solution
%         1.4, 18/Dec/2010: + non-recursive engine
%         1.5: 01/Aug/2020: fix bug of AllVL1(1,1) returns wrong result
global MaxCounter;
if nargin<3 || isempty(L1ops)
    L1ops = '==';
end
n = floor(n); % make sure n is integer
if n<1
    v = [];
    return
end
if nargin<4  || isempty(MaxNbSol)
    MaxCounter = Inf;
else
    MaxCounter = MaxNbSol;
end
Counter(0);
switch L1ops
    case {'==' '='}
        if isnan(MaxCounter)
            % return the number of solutions
            v = nchoosek(n+L1-1,L1); % nchoosek(n+L1-1,n-1)
        else
            v = allVL1eq(n, L1);
        end
    case '<=' % call allVL1eq for various sum targets
        if isnan(MaxCounter)
            % return the number of solutions
            %v = nchoosek(n+L1,L1)*factorial(n-L1); BUG <- 16/Sep/2009: 
            v = 0;
            for j=0:L1
                v = v + nchoosek(n+j-1,j);
            end
            % See pascal's 11th identity, the sum doesn't seem to
            % simplify to a fix formula
        else
            v = cell2mat(arrayfun(@(j) allVL1eq(n, j), (0:L1)', ...
                'UniformOutput', false));
        end
    case '<'
        v = allVL1(n, L1-1, '<=', MaxCounter);
    otherwise
        error('allVL1: unknown L1ops')
end
end % allVL1
%%
function v = allVL1eq(n, L1)
global MaxCounter;
n = feval(class(L1),n);
s = n+L1;
sd = double(n)+double(L1);
notoverflowed = double(s)==sd;
if isinf(MaxCounter) && notoverflowed
    v = allVL1nonrecurs(n, L1);
else
    v = allVL1recurs(n, L1);
end
end % allVL1eq
%% Recursive engine
function v = allVL1recurs(n, L1, head)
% function v=allVL1eq(n, L1);
% INPUT
%    n: length of the vector
%    L1: desired L1 norm
%    head: optional parameter to by concatenate in the first column
%          of the output
% OUTPUT:
%    if head is not defined
%      v: (m x n) array such as sum(v,2)==L1
%         all elements of v is naturel numbers {0,1,...}
%         v contains all (=m) possible combinations
%         v is (dictionnary) sorted
% Algorithm:
%    Recursive
global MaxCounter;
if n==1
    if Counter < MaxCounter
        v = L1;
    else
        v = zeros(0,1,class(L1));
    end
else % recursive call
    v = cell2mat(arrayfun(@(j) allVL1recurs(n-1, L1-j, j), (0:L1)', ...
        'UniformOutput', false));
end
if nargin>=3 % add a head column
    v = [head+zeros(size(v,1),1,class(head)) v];
end
end % allVL1recurs
%%
function res=Counter(newval)
persistent counter;
if nargin>=1
    counter = newval;
    res = counter;
else
    res = counter;
    counter = counter+1;
end
end % Counter
%% Non-recursive engine
function v = allVL1nonrecurs(n, L1)
% function v=allVL1eq(n, L1);
% INPUT
%    n: length of the vector
%    L1: desired L1 norm
% OUTPUT:
%    if head is not defined
%      v: (m x n) array such as sum(v,2)==L1
%         all elements of v is naturel numbers {0,1,...}
%         v contains all (=m) possible combinations
%         v is (dictionnary) sorted
% Algorithm:
%    NonRecursive
% Chose (n-1) the splitting points of the array [0:(n+L1)]
p = n+L1-1;
if p ~= 1 % bug occurs for allVL1nonrecurs(1,1) since
    % nchoosek take 1 vector length and not vector it self
    s = nchoosek(1:p,n-1);
else
    if n == 1
        s = zeros(1,0);
    elseif n == 2
        s = 1;
    end
end
m = size(s,1);
s1 = zeros(m,1,class(L1));
s2 = (n+L1)+s1;
v = diff([s1 s s2],1,2); % m x n
v = v-1;
end % allVL1nonrecurs
function c = cross_dim1(a,b)
% c = cross_dim1(a,b)
% Calculate cross product along the first dimension
% NOTE: auto expansion allowed
c = zeros(max(size(a),size(b)));
c(1,:) = a(2,:).*b(3,:)-a(3,:).*b(2,:);
c(2,:) = a(3,:).*b(1,:)-a(1,:).*b(3,:);
c(3,:) = a(1,:).*b(2,:)-a(2,:).*b(1,:);
end % cross_dim1