Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

How to match concentric ellipses' edges with some vertical line in R plot

Tags:

plot

r

Background: My ultimate goal in asking the following question is to turn the following into a function. In short, I have a curve placed above multiple concentric ellipses (ovals). The core issue is how to match the curve above these ellipses so if the curve changes the ellipses change accordingly.

Details:

More specifically, I want to have 11 concentric ellipses (see below image). At the outer edge of these 11 ellipses, I need 13 line segments to go up to attach to the curve line (see image below). The curve shown below is a Cauchy distribution (see last line of R code below the image). But this curve can be any other symmetric curve (e.g., normal distribution, t-distribution).

Question:

How can I determine (match) the coordinates of the ellipses with the line segment that must exactly connect to the outer edge of the ellipses?

Note: Most probably, the problem has to do with the math behind ellipses.

(I've tried something and have provided annotated R codes further below).

enter image description here

Here is my R code:

if(!require(library(plotrix))){install.packages('plotrix') }
library(plotrix)                ## A package for drawing ellipses ##

# In the "plotrix":

# "x" and "y" are the coordinates of the center.
# "a" and "b" give the radii of the ovals.


plot(1, ty='n', ann = F, xlim = c(-4, 6), ylim = c(-3, 1.5) ) ## A platform for ellipses


draw.ellipse(x = rep(1, 11), y = rep(-1.2, 11), 
         a = seq(1, 6, by = .4), b = seq(1/4.5, 6/4.5 , by = .4/4.5 ), 
         lty = 2, border = 'gray60' )       ## Draw multiple Concentric ellipses ##


AA <- seq(-4, 6, len = 13)     ## A range of values on the x-xis just like "xlim" ##
BB <- dcauchy( AA, 1, .95)*5   ## The Height for the AA according to a distribution ##

segments(AA, rep(-1.2, length(AA) ), AA, BB, lty = 3, lwd = 2, col= 'green4' )

curve(dcauchy(x, 1, .95)*5, -4, 6, add = T, col ='magenta', lwd = 3)
like image 780
rnorouzian Avatar asked Mar 29 '17 18:03

rnorouzian


1 Answers

library(plotrix)

#AVAILABE DATA
el_x = 1 #Centre of ellipse
el_y = -1.2 #Centre of ellipse
n = 13 #The number of ellipses (13 in this example)
el_a = seq(1, 6, length.out = n) #'a' values
#You likely have 'b' too, but I am computing for now
el_b = seq(1/4.5, 6/4.5, length.out = n)

#Retain only every other 'a' value
if (n %%2 !=0){ #Odd 'n'
    el_a_2 = el_a[seq_along(el_a) %% 2 != 0]
} else { #Even 'n'
    el_a_2 = el_a[seq_along(el_a) %% 2 == 0] #Modify if necessary
}

#Store curve in a variable for now
cc = curve(dcauchy(x, 1, .95)*5, min(el_x - el_a), max(el_x + el_a))

#Calculate x-values (subtract and add 'a' to 'el_x') of ellipse extremes
AA <- c(el_x - el_a_2, el_x, el_x + el_a_2) #The x-values you need for segments
AA = sort(AA) #To plot lines and points as you want

#Calculate y-values from AA
BB <- dcauchy(AA, 1, .95)*5  #The y-values you need for segments

#Raise the curve to have 0.2 distance (OPTIONAL)
Raise_curve = max(el_b + el_y) + 0.2 - (min(BB) - max(el_b + el_y))

BB = BB + Raise_curve #Raise BB

plot(1, type = 'n', ann = F,
    xlim = c(min(el_x - el_a), max(el_x + el_a)),
    ylim = c(min(el_y - el_b), max(BB)))

draw.ellipse(x = rep(el_x, n), y = rep(el_y, n), 
         a = el_a, b = el_b, 
         lty = 2, border = 'gray60')

segments(x0 = AA, x1 = AA, y0 = el_y, y1 = BB, lty = 3, lwd = 2,
    col = c(rep('green3', floor(length(AA)-3)/2),
        rep('navy', 3 ),
        rep('green3', length(AA) - 3 - floor(length(AA)-3)/2)))

lines(cc$x, cc$y+ Raise_curve, lwd = 2, col = "red")

#ADDITIONAL EXAMPLE
points(AA, rep(-1.2, length(AA) ), pch = c(rep(22, 6), 21, rep(22, 6)),
    col = c(rep('blue', 6), 'red', rep('blue', 6)),
    bg = c(rep('blue', 6), 'red', rep('blue', 6)), cex = 3)

xx <- seq(AA[7], AA[9], length.out = 1000)
yy = dcauchy(xx, 1, .95)*5
yy = yy+Raise_curve
polygon(c(AA[7], xx, AA[9]), c(el_y, yy, el_y), border = NA,
                       col= adjustcolor('green', alpha.f = .2) )

enter image description here

like image 71
d.b Avatar answered Oct 11 '22 09:10

d.b