function profiles = profiles_from_AFS2D(AFS,V,bndpts,hx,hy)
% computes the feasible profiles from a 2D-AFS for a three component system
% 
% input:
% AFS ...... area of feasible solutions computed by facpack
% V ........ matrix of the right singular vectors
% bndpts ... number of nearly equidistant boundary points for each AFS-subset
% hx ....... stepsize for the discretization
% hy ....... stepsize for the discretization
% 
% example:
% profiles = profiles_from_AFS2D(AFS, V, 20, 0.1, 0.1)

col = [0     0.4470     0.7410;
    0.8500    0.3250    0.0980;
    0.9290    0.6940    0.1250;
    0.4940    0.1840    0.5560;
    0.4660    0.6740    0.1880;
    0.3010    0.7450    0.9330];

profiles = [];

parts = 0;
for i=1:length(AFS)
    if ~isempty(AFS{i})
        dimP(i) = length(AFS{i});
        parts = parts+1;
    end
end

L = zeros(parts,max(dimP));
for i=1:parts
    l(i) = 0;
    for j=1:dimP(i)-1
        L(i,j) = norm(AFS{i}(j,:)-AFS{i}(j+1,:),2);
        l(i) = l(i)+L(i,j);
    end
end
lenP = l;
LENP = L;

if parts == 2
    parts = 1;
end

for seg = 1:parts
    if parts == 1
        SolT = [];
        for j = 1:2
            % the extremal points
            [m,id(1)] = min(AFS{j}(:,1));
            [m,id(2)] = max(AFS{j}(:,1));
            [m,id(3)] = min(AFS{j}(:,2));
            [m,id(4)] = max(AFS{j}(:,2));

            SolT = [SolT; AFS{j}(id,:)];

            % the boundary points
            SolT = [SolT; AFS{j}(1,:)];
            lb = size(SolT,1);
            lt = 1;
            hr = lenP(j)/bndpts;
            i = 1;
            while i < bndpts
                len = sum(LENP(j,1:lt));
                if hr*i < len
                    lb = lb+1;
                    al = (len-hr*i)/LENP(j,lt);
                    SolT(lb,:) = al*AFS{j}(lt,:)+(1-al)*AFS{j}(lt+1,:);
                    i = i+1;
                else
                    lt = lt+1;
                end
            end
        end

        % the inner points
        for j=1:2
            lockup{j} = zeros(dimP(j)-1,2);
            for i=1:dimP(j)-1
                lockup{j}(i,:) = [min(AFS{j}(i:i+1,2)) max(AFS{j}(i:i+1,2))];
            end
        end
        ny = ceil((max(AFS{1}(:,2))-min(AFS{1}(:,2)))/hy)+1;
        gridy = linspace(min(AFS{1}(:,2)),max(AFS{1}(:,2)),ny);

        for i=2:ny-1
            if (gridy(i) > max(AFS{2}(:,2))) || (gridy(i) < min(AFS{2}(:,2)))
                cuts = [];
                cutp = [];
                for j=1:dimP(1)-1
                    if (gridy(i) > lockup{1}(j,1)) && (gridy(i) < lockup{1}(j,2))
                        al = (AFS{1}(j+1,2)-gridy(i))/(AFS{1}(j+1,2)-AFS{1}(j,2));

                        cuts = [cuts; j];
                        cutp = [cutp; al*AFS{1}(j,:)+(1-al)*AFS{1}(j+1,:)];
                    end
                end
                nc = length(cuts);

                nx = ceil(abs(cutp(1,1)-cutp(2,1))/hx)+1;
                gridx = linspace(cutp(1,1),cutp(2,1),nx);
                for j=2:nx-1
                    lb = lb+1;
                    SolT(lb,:) = [gridx(j), gridy(i)];
                end
            else
                cuts = [];
                cutp = [];
                for ii = 1:2
                    cuts{ii} = [];
                    cutp{ii} = [];
                    for j=1:dimP(ii)-1
                        if (gridy(i) > lockup{ii}(j,1)) && (gridy(i) < lockup{ii}(j,2))
                            al = (AFS{ii}(j+1,2)-gridy(i))/(AFS{ii}(j+1,2)-AFS{ii}(j,2));

                            cuts{ii} = [cuts{ii}; j];
                            cutp{ii} = [cutp{ii}; al*AFS{ii}(j,:)+(1-al)*AFS{ii}(j+1,:)];
                        end
                    end
                end
                nc = length(cuts);

                nx = ceil(abs(cutp{1}(1,1)-cutp{1}(2,1))/hx)+1;
                gridx = linspace(cutp{1}(1,1),cutp{1}(2,1),nx);
                for j=2:nx-1
                    if (gridx(j) < min(cutp{2}(:,1))) || (gridx(j) > max(cutp{2}(:,1)))
                        lb = lb+1;
                        SolT(lb,:) = [gridx(j), gridy(i)];
                    end
                end
            end 
        end
        profiles = [ones(lb,1) SolT]*V(:,1:3)';
        
    elseif parts >= 3
        afs = AFS{seg};
        dP = dimP(seg);
        lP = lenP(seg);
        LP = LENP(seg,:);

        % the extremal points
        [m,id(1)] = min(afs(:,1));
        [m,id(2)] = max(afs(:,1));
        [m,id(3)] = min(afs(:,2));
        [m,id(4)] = max(afs(:,2));

        SolT = afs(id,:);

        % the boundary points
        SolT(5,:) = afs(1,:);
        lb = 5;
        lt = 1;
        hr = lP/bndpts;
        i = 1;

        while (i < bndpts) && (lt < size(afs,1))
            len = sum(LP(1:lt));
            if hr*i <= len
                lb = lb+1;
                al = (len-hr*i)/LP(lt);
                SolT(lb,:) = al*afs(lt,:)+(1-al)*afs(lt+1,:);
                i = i+1;
            else
                lt = lt+1;
            end
        end

        % the inner points
        lockup = zeros(dP-1,2);
        for i=1:dP-1
            lockup(i,:) = [min(afs(i:i+1,2)) max(afs(i:i+1,2))];
        end
        ny = ceil((max(afs(:,2))-min(afs(:,2)))/hy)+1;
        gridy = linspace(min(afs(:,2)),max(afs(:,2)),ny);

        for i=2:ny-1
            cuts = [];
            cutp = [];
            for j=1:dP-1
                if (gridy(i) > lockup(j,1)) && (gridy(i) < lockup(j,2))
                    al = (afs(j+1,2)-gridy(i))/(afs(j+1,2)-afs(j,2));

                    cuts = [cuts; j];
                    cutp = [cutp; al*afs(j,:)+(1-al)*afs(j+1,:)];
                end
            end
            nc = length(cuts);

            nx = ceil(abs(cutp(1,1)-cutp(2,1))/hx)+1;
            gridx = linspace(cutp(1,1),cutp(2,1),nx);
            for j=2:nx-1
                lb = lb+1;
                SolT(lb,:) = [gridx(j), gridy(i)];
            end

        end
        profiles{seg} = [ones(lb,1) SolT]*V(:,1:3)';
    end
    
    figure(seg)
    plot(1:size(V,1), profiles{seg}, 'color', col(seg,:))
    axis tight
    yl = get(gca, 'ylim');
    set(gca, 'ylim',[min(0, yl(1)), yl(2)]);
end
