Derivation of intermediate points (great circle geometry) in the geosphere package

I have started using the geosphere package kindly provided by Robert  
Hijmans [http://cran.r-project.org/web/packages/geosphere/]. My  
intention is to derive a global map of the air routes density using  
the open access flight routes [http://openflights.org/data.html].

I have discovered some strange results when trying to derive the  
intermediate points on a path:

gcIntermediate(p1, p2, n)
*this uses the formulas described in  
[http://williams.best.vwh.net/avform.htm#Intermediate]

I discovered that, depending on the location of points considering  
quadrants, the current formula will (systematically?) shift the  
coordinates to some new quadrant, which is obviously a bug:

library(geosphere)
gcIntermediate(c(0,0), c(10,10), 5)  # this one is OK

lon      lat
[1,] 1.803343 1.830237
[2,] 3.422141 3.468623
[3,] 4.961631 5.019001
[4,] 6.508449 6.565756
[5,] 8.150886 8.192447

gcIntermediate(c(-10,10), c(10,10), 5)  # this one needs to be  
fixed: ifelse(lon>0, lon+180, 180-lon); -1*lat

lon       lat
[1,]  173.6232 -10.08957
[2,]  176.8920 -10.13646
[3,]  180.0000 -10.15108
[4,] -176.8920 -10.13646
[5,] -173.6232 -10.08957

gcIntermediate(c(-20,-10), c(5,10), 5)  # this is OK

lon       lat
[1,] -15.4984322 -6.467336
[2,] -11.4054190 -3.175910
[3,]  -7.5000000  0.000000
[4,]  -3.5945810  3.175910
[5,]   0.4984322  6.467336

gcIntermediate(c(-20,-25), c(10,5), 5) # this needs to be fixed:   
for(lon2<lon1) {ifelse(lon>0, 180-lon, 180+lon)}; -1*lat

lon        lat
[1,]  165.7968 20.0356468
[2,]  170.9455 15.1760517
[3,]  175.7253 10.3444590
[4,] -179.6410  5.4446120
[5,] -174.9524  0.3684525


Is somebody really good in geosphere geometry?

BR & thanx

T. Hengl
http://home.medewerker.uva.nl/t.hengl/