You need to use spherical trigonometry, part of spherical geometry for full accuracy. However, since you are dealing with only a small piece of the sphere, euclidian geometry will do if you remember one thing.
As latitude increases, the lines of longitude get closer together. For example, near the North Pole, the latitude lines are almost touching. So condition your latitude differences, diminishing them by mulitlying by a factor of cos(latitude). That will give you good enough accuracy for your app.
$n = 60.406505416667;
$s = 60.400570555556;
$e = 5.3351572222222;
$w = 5.3190577777778;
$rotn = 3.7088732260919;
$a = ($e + $w) / 2.0;
$b = ($n + $s) / 2.0;
$squish = cos(deg2rad($b));
$x = $squish * ($e - $w) / 2.0;
$y = ($n - $s) / 2.0;
$ne = array(
$a + ($x * cos(deg2rad($rotn)) - $y * sin(deg2rad($rotn))) /$squish,
$b + $x * sin(deg2rad($rotn)) + $y *cos(deg2rad($rotn))
);
$nw = array(
$a - ($x * cos(deg2rad($rotn)) + $y * sin(deg2rad($rotn))) /$squish,
$b - $x * sin(deg2rad($rotn)) + $y *cos(deg2rad($rotn))
);
$sw = array(
$a - ($x * cos(deg2rad($rotn)) - $y * sin(deg2rad($rotn))) /$squish,
$b - $x * sin(deg2rad($rotn)) - $y *cos(deg2rad($rotn))
);
$se = array(
$a + ($x * cos(deg2rad($rotn)) + $y * sin(deg2rad($rotn))) /$squish,
$b + $x * sin(deg2rad($rotn)) - $y *cos(deg2rad($rotn))
);
print_r(array(
'sw'=>$sw,
'se'=>$se,
'ne'=>$ne,
'nw'=>$nw,
));
My $squish
variable is the cos(lat) I mentioned. There is de-squishing for the relative part of horizontal lengths. The sine table looks like this:
NE: (a + x cos A - y sin A, b + x sin A + y cos A)
NW: (a - x cos A - y sin A, b - x sin A + y cos A)
SW: (a - x cos A + y sin A, b - x sin A - y cos A)
SE: (a + x cos A + y sin A, b + x sin A - y cos A)
Perhaps tttppp could account for the differences from tttppp's table.