function [R, M, C, Np, f, g, Vmax] = mecrow_model()
  mecrow_model_data;
  f = @(F, theta, phase)  interp2_periodic(F_vec, theta_vec, psi_mat, Np, F, theta, phase);
  g = @(F, theta, phase)  interp2_periodic(F_vec, theta_vec, tau_mat, Np, F, theta, phase);
end

function z = interp2_periodic(X, Theta, Z, Np, x, theta, phase);
  if phase > 1
    z = interp2_periodic(X, Theta, Z, Np, x, theta + (phase-1)*2*pi/(3*Np), 1);
  else
    if theta >= 0 && theta <= 2*pi/Np
      z = interp2(X, Theta, Z, x, theta);
    elseif theta > 2*pi
      z = interp2_periodic(X, Theta, Z, Np, x, theta - 2*pi/Np, 1);
    elseif theta < 0
      z = interp2_periodic(X, Theta, Z, Np, x, theta + 2*pi/Np, 1);
    end
  end
end