To simplify the formulas, I will work in terms of latitudes a1 and a2 and longitudes b1 and b2 in degrees and assume the arguments to the sin and cos functions are in degrees (to eliminate that pesky conversion to radians by dividing by 57.2958). Also let d = 24900/360, which is the factor for converting degrees to miles.

Then the great circle distance formula in miles and degrees is

R = d * arccos[sin(a1)sin(a2) + cos(a1)cos(a2)cos(b2 - b1)].

Assume (a1,b1) is given for the lower right corner of the square. I assume the goal is to find (a2,b2) for the upper left corner such that the distance along the bottom side between (a1,b1) and (a1,b2) is 2r miles and the distance along the right side between (a2,b1) and (a1,b1) is also 2r miles.

To find the distance along the bottom side between (a1,b1) and (a1,b2), set a2 = a1 and R = 2r in the great circle formula and solve for b2:

2r = d * arccos[sin(a1)^2 + cos(a1)^2 * cos(b2 - b1)],

cos(2r/d) = sin(a1)^2 + cos(a1)^2 * cos(b2 - b1),

cos(b2 - b1) = (cos(2r/d) - sin(a1)^2)/cos(a1)^2,

and finally

b2 = b1 + arccos[(cos(2r/d) - sin(a1)^2)/cos(a1)^2].[**]

To find the distance along the right side between (a1,b1) and (a2,b1), set b2 = b1 and R = 2r in the great circle formula and solve for a2:

2r = d * arccos[sin(a1)sin(a2) + cos(a1)cos(a2)cos(b2 - b1)],

and since b2 - b1 = 0, cos(b2 - b1) = cos(0) = 1, and

cos(2r/d) = sin(a1)sin(a2) + cos(a1)cos(a2) = cos(a2 - a1)

using an angle difference identity, which yields

2r/d = a2 - a1 and

a2 = a1 + 2r/d.[***]

Comparing equations [***] and [**] shows converting differences in latitudes to miles can be done directly without reference to the longitude, but converting differences in longitudes to miles depends on the latitude. These are known facts but we've proved them here.

You said you wanted to find a square (pin cushion) where the distance from the center point to side is r. Equation [**] finds a pin cushion where the center of the bottom side to the right side is r. To use the center point, use (a1+a2)/2 instead of a1 in the formula:

b2 = b1 + arccos[(cos(2r/d) - sin((a1+a2)/2)^2)/cos((a1+a2)/2)^2].[**]