function [fi,K5] = intsec(nvr1,nvr2,G,fi,dfi)
%
% Function intsec calculates the intersection of two vector spaces
% stored in fi and dfi with dimension nvr2 and nvr1, respectively.
% The resulting space is stored in fi.
%
ntry2 = zeros(1,G);
ntry3 = zeros(1,G);
%
% The normals to the spaces spanned by fi and dfi in the sum of both
% spaces, are calculated and stored in dfin and fin.
%
nvr12 = nvr1 + nvr2;
K6 = 1;
for K5 = nvr2:nvr12
  while K6 <= nvr1
    for I4 = 1:K5
      if I4 <= nvr2
        vec(I4) = dfi(K6,1:G)*conj(fi(1:G,I4));
      else
        vec(I4) = dfi(K6,1:G)*dfin(I4-nvr2,1:G)';
      end
    end
    I7 = K5 - nvr2;
    for I4 = 1:G
      vec2(I4) = dfi(K6,I4);
      for I5 = 1:nvr2
        vec2(I4) = vec2(I4) - vec(I5)*fi(I4,I5);
      end
      if I7 ~= 0
        for I6 = 1:I7
          vec2(I4) = vec2(I4) - vec(I6+nvr2)*dfin(I6,I4);
        end
      end
    end
    I4 = 1;
    while I4 <= G
      if (abs(vec2(I4)))^2 < 0.001
        I4 = I4 + 1;
      else
        break;
      end
    end
    if I4 <= G
      rnorm = sum(abs(vec2(1:G).').^2);
      rnorm = 1/sqrt(rnorm);
       dfin(I7+1,1:G) = vec2(1:G)*rnorm;
      ntry2(K6) = 1;
      K6 = K6 + 1;
      break;
    end
    K6 = K6 + 1;
  end
  if K6 > nvr1
    break;
  end
end
K6 = 1;
for K5 = nvr1:nvr12
  while K6 <= nvr2
    for I4 = 1:K5
      if I4 <= nvr1
        vec(I4) = conj(dfi(I4,1:G))*fi(1:G,K6);
      else
        vec(I4) = fin(1:G,I4-nvr1)'*fi(1:G,K6);
      end
    end
    I7 = K5 - nvr1;
    for I4 = 1:G
      vec2(I4) = fi(I4,K6);
      for I5 = 1:nvr1
        vec2(I4) = vec2(I4) - vec(I5)*dfi(I5,I4);
      end
      if I7 ~= 0
        for I6 = 1:I7
          vec2(I4) = vec2(I4) - vec(I6+nvr1)*fin(I4,I6);
        end
      end
    end
    I4 = 1;
    while I4 <= G
      if (abs(vec2(I4)))^2 < 0.001
        I4 = I4 + 1;
      else
        break;
      end
    end
    if I4 <= G
      rnorm = sum(abs(vec2(1:G).').^2);
      rnorm = 1/sqrt(rnorm);
      fin(1:G,I7+1) = vec2(1:G).'*rnorm;
      ntry3(K6) = 1;
      K6 = K6 + 1;
      break;
    end
    K6 = K6 + 1;
  end
  if K6 > nvr2
    break;
  end
end
K5 = 0;
K6 = 0;
for I4 = 1:G
  K5 = K5 + ntry2(I4);
  K6 = K6 + ntry3(I4);
end
%
% The K5 normals to fi and the K6 normals to dfi are stored in
% dfin and fin. The wanted intersection of spaces fi and dfi is
% to be simultaneously orthogonal to the space spanned by fin and dfin.
%
I7 = K6;
for I4 = 1:K5
  vec(1:I7) = dfin(I4,1:G)*conj(fin(1:G,1:I7));
  for I5 = 1:G
    vec2(I5) = dfin(I4,I5);
    for I6 = 1:I7
      vec2(I5) = vec2(I5) - vec(I6)*fin(I5,I6);
    end
  end
  I5 = 1;
  while I5 <= G
    if abs(vec2(I5)) < 0.001
      I5 = I5 + 1;
    else
      break;
    end
  end
  if I5 <= G
    I7 = I7 + 1;
    rnorm = 0;
    rnorm = sum(abs(vec2(1:G).').^2);
    rnorm = 1/sqrt(rnorm);
    fin(1:G,I7) = vec2(1:G).'*rnorm;
  end
end
%
% An orthonormal basis for the sum of the spaces spanned by fin
% and dfin, is now stored in fin with I7 basis vectors. By 
% orthogonalizing the basis vectors of dfi to those of fin, we
% find the intersection of fi and dfi.
%
K5 = 0;
K6 = I7;
for I4 = 1:nvr1
  for I5 = 1:I7
    if I5 <= K6
      vec(I5) = dfi(I4,1:G)*conj(fin(1:G,I5));
    else
      vec(I5) = dfi(I4,1:G)*dfin(I5-K6,1:G)';
    end
  end
  for I5 = 1:G
    vec2(I5) = dfi(I4,I5);
    for I6 = 1:I7
      if I6 <= K6
        vec2(I5) = vec2(I5) - vec(I6)*fin(I5,I6);
      else
        vec2(I5) = vec2(I5) - vec(I6)*dfin(I6-K6,I5); 
      end
    end
  end
  I5 = 1;
  while I5 <= G
    if abs(vec2(I5)) < 0.001
      I5 = I5 + 1;
    else
      break;
    end
  end
  if I5 <= G     
    I7 = I7 + 1;
    K5 = K5 + 1;
    rnorm = sum(abs(vec2(1:G).').^2);
    rnorm = 1/sqrt(rnorm);
    dfin(I7-K6,1:G) = vec2(1:G)*rnorm;
  end
end
if K5 ~= 0
%
% The intersection of fi and dfi is stored in fi.
%
  fi(1:G,1:K5) = dfin(1:K5,1:G).'; 
end
return             