Monthly Archives: April 2014

What If You Dig A Hole Through The Earth?

It suddenly struck me that that tiny pea, pretty and blue, was the Earth. I put up my thumb and shut one eye, and my thumb blotted out the planet Earth. I didn’t feel like a giant. I felt very, very small (Neil Armstrong)

Where would you come out if you dig a hole straight downward from where you live through the center of the Earth? Supposing you survive to the extremely high pressure and temperature of the nucleus, would you find water or land at the other side? It maybe sound a ridiculous question (I will not refute that) but in this post I will estimate how many people would die drowned and how many will find land in the antipode of where they live. At least, knowledge does not take up any space.

I found a database of the United Nations with very useful information for my experiment: longitude, latitude and population of all capital cities of the world in 2011(1). I assumed that capital cities are a good sample of where people live. Maybe is not the best one since some very big countries are represented by only a city but is a good way to obtain a quick estimation. On the other hand, capital cities represent approximately 7% of the world population so in this sense is a very good sample.

Process is simple: loading the xls file, calculating the antipode of each point and checking where it is. Google provides information about country where a coordinate belongs. For coordinates on the sea, no information is returned. Once you have this, is easy to calculate proportion of people that will find water. My estimation is around 77% of people will find water on the other side. Taking into account that all people leave from land and approximately 70% of the Earth’s surface is water, this figure seems to be small but since both poles are symmetrical and are uninhabited, the estimation makes sense. Here you have a map with the result of the experimet. Points are capital cities and size is related with population. In blue, capitals with antipode on the sea and in brown, capitals with antipode in land:
WorldMapR
By the way, I am one of the 23% of lucky people that would find land in the other side. I live in Madrid, Spain:
madrid
An if some rainy afternoon having little to do I dig a hole through the Earth, I will appear in a place called Weber, in New Zealand:
weber
This estimation can be very silly but the physics involved in the experiment are very interesting as you can see here. By the way, there is a film called Total Recall (2012) where the only way to travel between last two cities in the world is using an elevator through the Earth. Here you have the code:

library(xlsx)
library(ggmap)
library(mapdata)
library(ggplot2)
#The xls file is in http://esa.un.org/unpd/wup/CD-ROM/WUP2011-F13-Capital_Cities.xls
CapitalCities <- read.xlsx("WUP2011-F13-Capital_Cities.xls", sheetName="Capital_Cities", startRow=13,header=TRUE)
names(CapitalCities) = gsub("\\.", "", names(CapitalCities))
#Obtain symmetric coordinates for each capital
CapitalCities$LatitudeSym <- -CapitalCities$Latitude
CapitalCities$LongitudeSym <- -sign(CapitalCities$Longitude)*(180-abs(CapitalCities$Longitude))
CapitalCities$DigResult <- apply(CapitalCities, 1, function(x) {unlist(revgeocode(c(as.numeric(x[11]),as.numeric(x[10]))))})
CapitalCities$Drowned <- is.na(CapitalCities$DigResult)*1
#Percentage of population saved
sum(CapitalCities$Drowned*CapitalCities$Populationthousands)/sum(CapitalCities$Populationthousands)
world <- map_data("world")
opt <- theme(legend.position="none",
axis.ticks=element_blank(),
axis.title=element_blank(),
axis.text =element_blank(),
plot.title = element_text(size = 35),
panel.background = element_rect(fill="turquoise1"))
p <- ggplot()
p <- p + geom_polygon(data=world, aes(x=long, y=lat, group = group),colour="white", fill="lightgoldenrod2" )
p <- p + geom_point(data=CapitalCities, aes(x=Longitude, y=Latitude, color=Drowned, size = Populationthousands)) + scale_size(range = c(2, 20), name="Population (thousands)")
p <- p + labs(title = "What if you dig a hole through the Earth?")
p <- p + scale_colour_gradient(low = "brown", high = "blue")
p <- p + annotate("rect", xmin = -135, xmax = -105, ymin = -70, ymax = -45, fill = "white")
p <- p + annotate("text", label = "Drowned", x = -120, y = -60, size = 6, colour = "blue")
p <- p + annotate("text", label = "Saved", x = -120, y = -50, size = 6, colour = "brown")
p <- p + geom_point(aes(x = -120, y = -65), size=8, colour="blue")
p <- p + geom_point(aes(x = -120, y = -55), size=8, colour = "brown")
p + opt
# Get a map of Spain, centered and signed in Madrid
madrid <- geocode('Madrid, Spain')
map.madrid <- get_map( location = as.numeric(madrid), color = "color", maptype = "roadmap", scale = 2, zoom = 6)
ggmap(map.madrid) + geom_point(aes(x = lon, y = lat), data = madrid, colour = 'red', size = 4)
# Get a map of New Zealand, centered and signed in Weber (the antipode of Madrid)
weber <- geocode('Weber, New Zealand')
map.weber <- get_map( location = as.numeric(weber), color = "color", maptype = "roadmap", scale = 2, zoom = 6)
ggmap(map.weber) + geom_point(aes(x = lon, y = lat), data = weber, colour = 'red', size = 4)

(1) United Nations, Department of Economic and Social Affairs, Population Division (2012). World Urbanization Prospects: The 2011 Revision, CD-ROM Edition.

Advertisements

The Pythagorean Tree Is In Bloom

There is geometry in the humming of the strings, there is music in the spacing of the spheres (Pythagoras)

Spring is here and I will be on holiday next week. I cannot be more happy! It is time to celebrate so I have drawn another fractal. It is called the Pythagorean Tree:

PythagoreanTree

Here you have the code. See you soon:

library("grid")
l=0.15 #Length of the square
grid.newpage()
gr <- rectGrob(width=l, height=l, name="gr") #Basic Square
pts <- data.frame(level=1, x=0.5, y=0.1, alfa=0) #Centers of the squares
for (i in 2:10) #10=Deep of the fractal. Feel free to change it
{
  df<-pts[pts$level==i-1,]
  for (j in 1:nrow(df))
  {
    pts <- rbind(pts, 
                 c(i, 
                   df[j,]$x-2*l*((1/sqrt(2))^(i-1))*sin(df[j,]$alfa+pi/4)-0.5*l*((1/sqrt(2))^(i-2))*sin(df[j,]$alfa+pi/4-3*pi/4), 
                   df[j,]$y+2*l*((1/sqrt(2))^(i-1))*cos(df[j,]$alfa+pi/4)+0.5*l*((1/sqrt(2))^(i-2))*cos(df[j,]$alfa+pi/4-3*pi/4),                  
                   df[j,]$alfa+pi/4))
    pts <- rbind(pts, 
                 c(i, 
                   df[j,]$x-2*l*((1/sqrt(2))^(i-1))*sin(df[j,]$alfa-pi/4)-0.5*l*((1/sqrt(2))^(i-2))*sin(df[j,]$alfa-pi/4+3*pi/4), 
                   df[j,]$y+2*l*((1/sqrt(2))^(i-1))*cos(df[j,]$alfa-pi/4)+0.5*l*((1/sqrt(2))^(i-2))*cos(df[j,]$alfa-pi/4+3*pi/4),                  
                   df[j,]$alfa-pi/4))
  }
}
for (i in 1:nrow(pts))
{
  grid.draw(editGrob(gr, vp=viewport(x=pts[i,]$x, y=pts[i,]$y, w=((1/sqrt(2))^(pts[i,]$level-1)), h=((1/sqrt(2))^(pts[i,]$level-1)), angle=pts[i,]$alfa*180/pi), 
            gp=gpar(col=0, lty="solid", fill=rgb(139*(nrow(pts)-i)/(nrow(pts)-1), 
                                                 (186*i+69*nrow(pts)-255)/(nrow(pts)-1), 
                                                  19*(nrow(pts)-i)/(nrow(pts)-1), 
                                                  alpha= (-110*i+200*nrow(pts)-90)/(nrow(pts)-1), max=255))))
}

The Collatz Fractal

It seems to me that the poet has only to perceive that which others do not perceive, to look deeper than others look. And the mathematician must do the same thing (Sofia Kovalevskaya)

How beautiful is this fractal! In previous posts I colored plots using module of complex numbers generated after some iterations. In this occasion I have used the escape-time algorithm, a very well known coloring algorithm which is very easy to implement in R.
Collatz07
Those who want to know more about this fractal can go here. For coloring, I chose a simple scale from red to yellow resulting a fractal interpretation of my country’s flag. You can choose another scale or use a RColorBrewer palette as I did in this previous post. Choosing another x or y ranges you can zoom particular areas of the fractal.

Try yourself and send me your pictures!

library(ggplot2)
xrange <- seq(-8, 8, by = 0.01)
yrange <- seq(-3, 3, by = 0.01)
f  <- function (z) {1/4*(2+7*z-(2+5*z)*cos(pi*z))}
z <- outer(xrange, 1i*yrange,'+')
t <- mat.or.vec(nrow(z), ncol(z))
for (k in 1:10)
{
  z <- f(z)
  t <- t + (is.finite(z)+0)
}
## Supressing texts, titles, ticks, background and legend.
opt <- theme(legend.position="none",
             panel.background = element_blank(),
             axis.ticks=element_blank(),
             axis.title=element_blank(),
             axis.text =element_blank())
z <- data.frame(expand.grid(x=xrange, y=yrange), z=as.vector(t))
ggplot(z, aes(x=x, y=y, color=z)) + geom_tile() + scale_colour_gradient(low="red", high="yellow") + opt

Do Not Play With Mr. Penney

Facts do not speak (Henry Poincare)

Mr. Penney is my best friend. He is maths teacher and loves playing. Yesterday we were in his office at the university when he suggested me a game:

When you toss a coin three times, you can obtain eight different sequences of tails and heads: TTT, TTH, THT, HTT, THH, HTH, HHT and HHH. Using a fair coin, all sequences have the same chances to appear. Choose one sequence and I will then choose another one. I will toss a coin until either your or my sequence appears as a consecutive subsequence of the coin toss outcomes. The player whose sequence appears first wins. I will repeat this procedure 100 times. The one with more games won is the winner of the game.  Don’t worry: I will not toss the coin manually. I will simulate using my computer. What’s your bet?

Ok, my bet is THT, I said. After some seconds, Mr. Penney said: My bet is TTH.

This was the result of the first round:

Round1Another chance? told me Mr. Penney. Of course! Now my bet is TTH! I said. In fact, I was thinking Take that! Now I chose your previous bet. Do you think I am foolish?. After some seconds, Mr. Penney said: My bet now is HTT.

This was the result of the second round:

Round2

Another chance? told me Mr. Penney. At this point, I was very suspicious but I wanted the last chance so I told him Of course! Now my bet is HTT! I wanted to try my strategy one more time. After some seconds, Mr. Penney said: My bet now is HHT.

This was the result of the third round:

Round3Ok, I give it up! What’s the trick? I said. And Mr. Penney explained it to me. You can find the explanation here. This is the last time I play with you! I told him once he finished the explanation.

Here you have the code. Feel free to play:

library(gridExtra)
Me     <- "TTH"
Penney <- "HTT"
results <- data.frame(play= numeric(0), Penney = integer(0), Me = character(0))
for (i in 1:100) {
play <- c()
repeat {play <- do.call(paste, c(play, as.list(sample(c("H","T"), 1)), sep=""))
  if (grepl(Penney, play)|grepl(Me, play)) {
    results <- rbind(results, data.frame(play= i, Penney = as.numeric(grepl(Penney, play)), Me = as.numeric(grepl(Me, play))))
    break}}}
grid.newpage()
table <- rbind(
  c("Me", Me, sum(results$Me), if(sum(results$Penney) > sum(results$Me)) "Loser" else "Winner"),
  c("Penney", Penney, sum(results$Penney), if(sum(results$Penney) > sum(results$Me)) "Winner" else "Loser"))
grid.table(table,
           cols = c("Player", "Bet", "Games Won", "Result"),
           gpar.colfill = gpar(fill="palegreen3",col="White"),
           gpar.corefill =  gpar(fill="palegreen",col="White"),
           gpar.rowfill = gpar(fill=NA, col=NA))