[R] Plotted circle does not go through desired points - very long email with code

[Monica Pisica]

Hi,

 

I am trying to solve a problem related to Poincare
circles ( for more info http://www.ms.uky.edu/~droyster/courses/spring08/math6118/Classnotes/Chapter09.pdf).
In a nutshell i am trying to replicate the method in the above pdf section
9.2.1. that explains in broad terms how to draw the arc inside a circle that
goes through 2 previously set points on the first circle.

 

I think i came up with some code to do that in R but even
if i think the theory is sound enough and i cannot pinpoint my mistake, the
graphic does not go through my 2 points – although i think it should.

 

So i start with:

1.     
A matrix called mat 10 x 10 that has only 0 and
1 as values, 0 means no link between points, 1 means a link …. So if i plot all
these points on a circle of center O (0,0) (called "base circle") and
radius r = 20, i will know which points i want to link with an arc inside the
base circle.

 

mat <- matrix(0, nrow = 10, ncol = 10, byrow = TRUE,
dimnames = list(c("A", "B", "C", "D",
"E", "F", "G", "H", "K",
"L"), c("A", "B", "C", "D",
"E", "F", "G", "H", "K",
"L")))

 

mat[1,] <- c(0,1,0,0,1,0,1,0,0,1)

mat[2,] <- c(1,0,1,0,1,0,0,1,0,0)

mat[3,] <- c(0,1,0,1,0,1,0,0,0,0)

mat[4,] <- c(0,0,1,0,0,1,0,1,1,0)

mat[5,] <- c(1,1,0,0,0,0,1,0,1,1)

mat[6,] <- c(0,0,1,1,0,0,0,1,0,1)

mat[7,] <- c(1,0,0,0,1,0,0,0,1,0)

mat[8,] <- c(0,1,0,1,0,1,0,0,1,0)

mat[9,] <- c(0,0,0,1,1,0,1,1,0,0)

mat[10,] <- c(1,0,0,0,1,1,0,0,0,0)

 

2.     
The rest of the code – i will try to comment it
as much as i could ;-)

 

library(plotrix)

library(CircStats)

library(spatstat)

 

plot.new()  #
draw.circle will give an error without a plot open – although for now i don't
want to plot it.

r=20

cr <- draw.circle(0, 0, r, nv=200)

n <- dim(mat)[1]

p <- round(seq(1, 200, 200/n))

xy <- data.frame(x = cr$x[p], y = cr$y[p])  # get coordinates for the 10 points in mat.

row.names(xy) <- dimnames(mat)[[2]]

 

# decide to draw the arc (circle) between the 8th
point in mat (point "H") and its first dependency (point "B")

 

i=8

j=1

 

m <- xy$y[i]/xy$x[i]

a.m <- atan(m) # angle in radians

b.m <-( pi/2)+a.m

m1 <- tan(b.m) # slope of tangent line on
"H"

b1 <- xy$y[i] -(m1*xy$x[i])

 

# decide on 2 random points on the tangent line

x1a <- -100

x1b <- 100

y1a <- m1*x1a + b1

y1b <- m1*x1b + b1

 

# Get the point dependency for "H"

mati <- mat[i,][mat[i,] ==1]

cmati <- names(mati)

n1 <- length(cmati)

 

# Get the first dependency: j=1

pt2 <- xy[cmati[j],]

 

# Get coordinates of the middle of segment
"H"-"B"

xm <- min(xy$x[i], pt2$x) + abs((min(xy$x[i], pt2$x)-max(xy$x[i],
pt2$x))/2)

ym <- min(xy$y[i], pt2$y) + abs((min(xy$y[i],
pt2$y)-max(xy$y[i], pt2$y))/2)

 

# find the perpendicular line equation that goes through
(xm, ym) 

m2 <- (xy$y[i]-pt2$y)/(xy$x[i]-pt2$x)

a.m2 <- atan(m2) # angle in radians

b.m2 <- (pi/2) +a.m2

m3 <- tan(b.m2) # slope of perpendicular line

b3 <- ym -(m3*xm) 

 

# Choose again some random points on this line:

x3a <- -100

x3b <- 100

y3a <- m3*x3a + b3

y3b <- m3*x3b + b3

 

# Build the spatial lines (the 2 perpendicular lines) and
find their intersection:

win <- owin(range(-100,x1a, x1b, x3a, x3b, xm, 100),
range(-100,y1a, y1b, y3a, y3b, ym, 100))

ln1 <- psp(x1a, y1a, x1b, y1b, window = win)

ln3 <- psp(x3a, y3a, x3b, y3b, window = win)

 

ppt <- crossing.psp(ln1, ln3)

 

# find the radius of the second circle that has center
at (ppt$x, ppt$y) and goes through B and H


r.crp <- sqrt((pt2$x-ppt$x)^2 + (pt2$y-ppt$y)^2)

r1.crp <- sqrt((xy$x[i]-ppt$x)^2 + (xy$y[i]-ppt$y)^2)

 

# r.crp should be equal with r1.crp – but there are some
very little difference  but not enough in
my opinion that the second circle not go through points H and B.

> r.crp==r1.crp

[1] FALSE

> r.crp-r1.crp

[1] 1.421085e-14

 

# now let's draw the plot: you will see that the second
circle goes through H but not through point B – although it should:

 

plot(cr$x, cr$y, type = "l", col =
"red", xlim = c(-70, 70), ylim=c(-70, 70))

points(xy$x, xy$y, cex=2, col = "blue")

identify(xy$x, xy$y, labels = row.names(xy))

crp <- draw.circle(ppt$x, ppt$y, r.crp, lwd=2, border
= "green")

lines(c(x1a, x1b), c(y1a, y1b), col = "blue")

lines(c(x3a, x3b), c(y3a, y3b), col = "cyan")

lines(c(xy$x[i], pt2$x), c(xy$y[i], pt2$y), col =
"violet", lty=2)

 

# and to test that the difference between r1.crp and
r.crp does not matter for the graphic:

crp <- draw.circle(ppt$x, ppt$y, r1.crp, lwd=2,
border = "orange", lty=2)

 

##############################

 

So my question is: Where is the mistake that makes the
second circle (green) not to go through my points H and B?

 

Thanks so much for your insights and i am sorry this is
such a long email – but it is quite a bit of code as well.

 

Monica