This works:
x = 500;
y = 700;
% a = p(1), b = p(2)
fn = @(p) norm(x - 110972.*p(1).*cos(p(2)/2)) + norm(y - 110972.*p(2).*cos(p(1)/2));
[ab, fval] = fminsearch(fn, [1 1]);
fprintf(1,'\n\ta = %f\n\tb = %f\n\n',ab)
x2 = 110972.*ab(1).*cos(ab(2)/2)
y2 = 110972.*ab(2).*cos(ab(1)/2)
producing:
a = 0.004506
b = 0.006308
x2 =
500.0003e+000
y2 =
700.0000e+000
