The Tissot's indicatrix of an Azimuthal Equidistant projection shows ellipses being stretched along latitudes.

How can I calculate the degree (amount) of this distortion? That is, what is the ratio of major axis to minor axis of a Tissot ellipse as a function of latitude (assuming the center of projection is the north pole)?

I can interpolate this well and get what I need. But would still be good to know what the analytical formula is. Feels like a function of cos(latitude) since an Azimuthal Equidistant projects latitudes with linearly increasing distance from center instead of cosine.