Understanding the output of `areaPolygon()` from p

From this Question I obtained the great function areaPolygon() which gives me the area within a polygon of coordinates. however, when I experiment with the function the calculations seem odd:

I start by creating a set of two points

 require(fields)
 coords <- c(11.3697193956209, 47.233380520521, 11.3723606043791, 
 47.235179479479)
 coords <- matrix(coords, nrow=2, ncol=2, byrow=TRUE)

then I check the distance between these two:

 rdist.earth(coords,coords,miles=FALSE)[1,2]

obtaining: 0.2827821 kilometres (which will be the diagonal of the rectangle)

I go on with creating a rectangle

 polygon <- matrix(coords, nrow=2, ncol=2)
 polygon <- rbind(polygon, polygon)
 polygon[4,2] <- polygon[1,2]
 polygon[4,1] <- polygon[2,1]
 polygon[3,2] <- polygon[2,2]
 polygon[3,1] <- polygon[1,1]
 polygon <- rbind(polygon, polygon[1,])

see if this looks good: plot(polygon)

Fourth step: I calculate the area within the polygon.

geosphere::areaPolygon(polygon)
[1] 31.99288 #from the help file I know this ought to be square metres.

however, I would have expected 200*200=40000 m² since the sides of my recangle are 200 by 200 metres. this can be checked via

rdist.earth(polygon,coords,miles=FALSE)

           [,1]         [,2]
 [1,] 0.0000000 2.827821e-01
 [2,] 0.2827821 9.504539e-05
 [3,] 0.2002671 1.996434e-01
 [4,] 0.1996501 2.002671e-01

so coming now to my question (finally) what am I doing wrong? thank you very much for your help!

You've made an invalid polygon! If you plot it with type="l" you'll see a bow-tie:

> plot(polygon,type="l")

Half the bow-tie will have a negative area, the other half positive, so your result is the difference in the sizes of the halves. They won't be exactly the same because of the earth being spherical...

You just need to reorder your points in polygon.