BTW, I just wanted to add that to get the elliptical path itself, you don't have to use
p.Y -= (c_radius1 - c_radius2) * (float)Math.Sin(a)

you can simply multiply the Y-component with a scale factor:
p.Y = p.Y * c_radius2/c_radius1 //---(where c_radius2/c_radius1 < 1)

But this only works to get the location on the path. in your case, since you're transforming all the points, this will squash the square figure itself, as well as the circle. What you can do is to define the square separately, in so-called object-space, rotate it there by the appropriate amount, and than translate it's "pivot point" to the location on the ellipse, in world-space. Another potential problem is that the square figure orients itself as if it was still on a circular path, which is the most apparent when the eccentricity is high. You can get the normal vector to a point on an ellipse made from a vertically squashed circle by stretching the same circle, along with the appropriate radius vector by the same amount (c_radius1/c_radius2), and then normalizing the radius vector.
Once normalized, than sin(figureAngle) = normal.Y, and cos(figureAngle) = normal.X.

So if:
circle:
x = cos(a);
y = sin(a);

ellipse:
x = circle.x;
y = circle.y * c_radius2/c_radius1;

Then:
normal:
x = circle.x;
y = circle.y * c_radius1/c_radius2;

or:

x = circle.x * c_radius2/c_radius1;;
y = circle.y * c_radius1/c_radius2 * c_radius2/c_radius1;

x = circle.x * c_radius2/c_radius1;
y = circle.y;

normal = normalize(normal);