Can we use the distance function to calculate lengths of arbitrary paths on ellipsoids of user-defined dimension?

I wanted to use the distance function to calculate path lengths along an ellipsoid of arbitrary semimajor and semiminor axes. When I deviate substantially from the wgs84 ellipsoid, the solution is no longer the same for the full pole to pole ellipse perimeter as the Ramanujan formula or discrete integral calculus. Is there a reason why that is? Is the reference ellipsoid object not defined correctly? Is there a domain over which the distance function is considered valid that I’m exceeding here? Code:

clc
close all
clear all

% used for integration way
dt = 0.001;
t = 0:dt:360;

refEllipse = referenceEllipsoid('wgs84');

% first check perimeter along the WGS84 globe - should be ellipse perimeter
a = refEllipse.SemimajorAxis;
b = refEllipse.SemiminorAxis;
P1 = 4*distance(0,0,90,0,refEllipse); % along lon=0, 1/4 of the ellipse
p = (b^2/a);
e = sqrt(1 - b^2/a^2);
r = p./(1+e.*cosd(t));
x = r.*cosd(t);
y = r.*sind(t);
P2 = sum(sqrt(diff(x).^2+diff(y).^2));
h = ((a-b)^2)/((a+b)^2);
P3 = pi*(a+b)*(1+3*h/(10+sqrt(4-3*h)));
disp('Verify for WGS84 the perimeter is all within 1 centimeter:')
disp([num2str(P1) ' ' num2str(P2) ' ' num2str(P3) ' m']);              

nEllipses = 12;
Perimeter1 = nan(nEllipses,nEllipses);
Perimeter2 = nan(nEllipses,nEllipses);
Perimeter3 = nan(nEllipses,nEllipses);
for a = 1:nEllipses
    for b = 1:a
        refEllipse.SemimajorAxis = a;
        refEllipse.SemiminorAxis = b;
        Perimeter1(a,b) = 4*distance(0,0,90,0,refEllipse);

          % just using discrete integral ds=sqrt(dx^2 + dy^2)
          p = (b^2/a);
          e = sqrt(1 - b^2/a^2);
          r = p./(1+e.*cosd(t));
          x = r.*cosd(t);
          y = r.*sind(t);
          Perimeter2(a,b) = sum(sqrt(diff(x).^2+diff(y).^2));

          % Ramanujan formula 
          % http://www.mathsisfun.com/geometry/ellipse-perimeter.html
          h = ((a-b)^2)/((a+b)^2);
          Perimeter3(a,b) = pi*(a+b)*(1+3*h/(10+sqrt(4-3*h)));
      end
  end

disp('Compare highly elliptical results:');
disp(['Max difference distance and Ramanujan formula: ' num2str(max(max(abs(Perimeter1-Perimeter3))))]);
disp(['Max difference discrete integral and Ramanujan formula: ' num2str(max(max(abs(Perimeter2-Perimeter3))))]);