You are an expert in writing code in many different languages. Your goal is to perform code completion for the following code, keeping in mind the rest of the code and the file meta data.
    Metadata details:
    filename: refine_p_surface.m
    path of file: ./repos/Isogeometric-Analysis/routines_for_iga
    
    the code of the file until where you have to start completion: %--------------------------------------------------------------------------
%  Author: Isaac J. Lee (ijlee2@ices.utexas.edu)
%  
%  
%  Summary:
%  
%  This routine does a p-refinement (degree elevation) of the given knot
%  knot vector in 2D.
%  
%  
%  Warning:
%  
%  The knot vector is assumed to be open. The nodes are assumed to be
%  ordered along the direction of 1, then along that of 2.
%  
%  
%  Instructions:
%  
%  Type the following onto Matlab's command window or in a code,
%  
%      [knots1_new, knots2_new, nodes_new, p1_new, p2_new] = refine_p_surface(knots1, knots2, nodes, p1, p2, q1, q2);
%  
%  where,
%  
%      knots1, knots2 are the knot vectors (column vectors)
%      nodes is an array of nodes
%      p1, p2 are the degree of the B-splines
%      q1, q2 are the degree increments
%  
%  
%  Outputs:
%  
%  1. New knot vectors
%  
%  2. New nodes array
%  
%  3. New degrees
%--------------------------------------------------------------------------
function [knots1_new, knots2_new, nodes_new, p1_new, p2_new] = refine_p_surface(knots1, knots2, nodes, p1, p2, q1, q2)
    % Some 
    
function [knots1_new, knots2_new, nodes_new, p1_new, p2_new] = refine_p_surface(knots1, knots2, nodes, p1, p2, q1, q2)
    % Some useful constants
    constant_p1p1 = p1 + 1;
    constant_p2p1 = p2 + 1;
    constant_p1pq1 = p1 + q1;
    constant_p2pq2 = p2 + q2;
    constant_p1pq1_p1 = constant_p1pq1 + 1;
    constant_p2pq2_p1 = constant_p2pq2 + 1;
    
    numKnots1 = size(knots1, 1);
    numKnots2 = size(knots2, 1);
    numNodes1 = numKnots1 - constant_p1p1;
    numNodes2 = numKnots2 - constant_p2p1;
    numNodes  = numNodes1 * numNodes2;
    numDimensions = size(nodes, 2);
    
    numKnots1_unique = size(unique(knots1), 1);
    numKnots2_unique = size(unique(knots2), 1);
    numKnots1_new = numKnots1 + q1 * numKnots1_unique;
    numKnots2_new = numKnots2 + q2 * numKnots2_unique;
    numNodes1_new = numNodes1 + q1 * (numKnots1_unique - 1);
    numNodes2_new = numNodes2 + q2 * (numKnots2_unique - 1);
    
    % Check for errors
    if (numNodes ~= size(nodes, 1))
        fprintf('Error: For the given knot vector, there should be %d nodes.\n\n', numNodes);
        
        knots1_new = knots1;
        knots2_new = knots2;
        nodes_new = nodes;
        p1_new = p1;
        p2_new = p2;
        
        return;
        
    elseif (q1 == 0 && q2 == 0)
        fprintf('Warning: Please specify a positive number for q1 or q2 for degree elevation.\n\nNo operation has been performed.\n\n');
        
        knots1_new = knots1;
        knots2_new = knots2;
        nodes_new = nodes;
        p1_new = p1;
        p2_new = p2;
        
        return;
        
    end