I was running into problems when using my fitted terms in another program: the same coefficients did not give the same surface.
I found the cause: the cosine and sine terms are inverted in zernfun. M>0 should be a cosine, while M<0 should be a sine.
I've corrected this, starting at line 193
z(:,idx_pos) = y(:,idx_pos).*cos(theta*m(idx_pos)');
z(:,idx_neg) = y(:,idx_neg).*sin(theta*-m(idx_neg)');
% note the required sign change "-m" in the sine term