@@ -13,70 +13,78 @@ arc3d <- function(from, to, center, radius, n, circle = 50, base = 0, plot = TRU
1313 from <- fixarg(from )
1414 to <- fixarg(to )
1515 center <- fixarg(center )
16-
16+
1717 m <- max(nrow(from ), nrow(to ), nrow(center ), length(base ))
1818 base <- rep_len(base , m )
1919
2020 result <- matrix (NA_real_ , nrow = 1 , ncol = 3 )
2121
2222 for (j in seq_len(m )) {
23- from1 <- getrow(from , j )
24- to1 <- getrow(to , j )
25- center1 <- getrow(center , j )
26- base1 <- base [j ]
27- logr1 <- log(veclen(from1 - center1 ))
28- logr2 <- log(veclen(to1 - center1 ))
29- A <- normalize(from1 - center1 )
30- B <- normalize(to1 - center1 )
31- steps <- if (base1 < = 0 ) 4 * abs(base1 ) + 1 else 4 * base1 - 1
32- for (k in seq_len(steps )) {
33- if (k %% 2 ) {
34- A1 <- A * (- 1 )^ (k %/% 2 )
35- B1 <- B * (- 1 )^ (k %/% 2 + (base1 > 0 ))
36- } else {
37- A1 <- B * (- 1 )^ (k %/% 2 + (base1 < = 0 ))
38- B1 <- A * (- 1 )^ (k %/% 2 )
39- }
40- theta <- acos(sum(A1 * B1 ))
41- if (isTRUE(all.equal(theta , pi )))
42- warning(" Arc " j , " points are opposite each other! Arc is not well defined."
43- if (missing(n ))
44- n1 <- ceiling(circle * theta / (2 * pi ))
45- else
46- n1 <- n
47-
48- if (missing(radius )) {
49- pretheta <- (k %/% 2 )* pi - (k %% 2 == 0 )* theta
50- if (k == 1 )
51- totaltheta <- (steps %/% 2 )* pi - (steps %% 2 == 0 )* theta + theta
52- p1 <- pretheta / totaltheta
53- p2 <- (pretheta + theta )/ totaltheta
54- radius1 <- exp(seq(from = (1 - p1 )* logr1 + p1 * logr2 ,
55- to = (1 - p2 )* logr1 + p2 * logr2 ,
56- length.out = n1 + 1 ))
57- } else
58- radius1 <- rep_len(radius , n1 )
59- arc <- matrix (NA_real_ , nrow = n1 + 1 , ncol = 3 )
60- p <- seq(0 , 1 , length.out = n1 + 1 )
61- arc [1 ,] <- center1 + radius1 [1 ]* A1
62- arc [n1 + 1 ,] <- center1 + radius1 [n1 + 1 ]* B1
63- AB <- veclen(A1 - B1 )
64- for (i in seq_len(n1 )[- 1 ]) {
65- ptheta <- p [i ]* theta
66- phi <- pi / 2 + (0.5 - p [i ])* theta
67- q <- (sin(ptheta ) / sin(phi ))/ AB
68- D <- (1 - q )* A1 + q * B1
69- arc [i ,] <- center1 + radius1 [i ] * normalize(D )
70- }
71- if (k == 1 )
72- result <- rbind(result , arc )
73- else
74- result <- rbind(result [- nrow(result ), ,drop = FALSE ], arc )
75- }
76- result <- rbind(result , result [1 ,])
23+ from1 <- getrow(from , j )
24+ to1 <- getrow(to , j )
25+ center1 <- getrow(center , j )
26+ # The "arc" might be a straight line
27+ if (isTRUE(all.equal(from1 , center1 )) ||
28+ isTRUE(all.equal(to1 , center1 )) ||
29+ isTRUE(all.equal(normalize(from1 - center1 ),
30+ normalize(to1 - center1 )))) {
31+ result <- rbind(result , from1 , to1 )
32+ } else {
33+ base1 <- base [j ]
34+ logr1 <- log(veclen(from1 - center1 ))
35+ logr2 <- log(veclen(to1 - center1 ))
36+ A <- normalize(from1 - center1 )
37+ B <- normalize(to1 - center1 )
38+ steps <- if (base1 < = 0 ) 4 * abs(base1 ) + 1 else 4 * base1 - 1
39+ for (k in seq_len(steps )) {
40+ if (k %% 2 ) {
41+ A1 <- A * (- 1 )^ (k %/% 2 )
42+ B1 <- B * (- 1 )^ (k %/% 2 + (base1 > 0 ))
43+ } else {
44+ A1 <- B * (- 1 )^ (k %/% 2 + (base1 < = 0 ))
45+ B1 <- A * (- 1 )^ (k %/% 2 )
46+ }
47+ theta <- acos(sum(A1 * B1 ))
48+ if (isTRUE(all.equal(theta , pi )))
49+ warning(" Arc " j , " points are opposite each other! Arc is not well defined."
50+ if (missing(n ))
51+ n1 <- ceiling(circle * theta / (2 * pi ))
52+ else
53+ n1 <- n
54+
55+ if (missing(radius )) {
56+ pretheta <- (k %/% 2 )* pi - (k %% 2 == 0 )* theta
57+ if (k == 1 )
58+ totaltheta <- (steps %/% 2 )* pi - (steps %% 2 == 0 )* theta + theta
59+ p1 <- pretheta / totaltheta
60+ p2 <- (pretheta + theta )/ totaltheta
61+ radius1 <- exp(seq(from = (1 - p1 )* logr1 + p1 * logr2 ,
62+ to = (1 - p2 )* logr1 + p2 * logr2 ,
63+ length.out = n1 + 1 ))
64+ } else
65+ radius1 <- rep_len(radius , n1 + 1 )
66+ arc <- matrix (NA_real_ , nrow = n1 + 1 , ncol = 3 )
67+ p <- seq(0 , 1 , length.out = n1 + 1 )
68+ arc [1 ,] <- center1 + radius1 [1 ]* A1
69+ arc [n1 + 1 ,] <- center1 + radius1 [n1 + 1 ]* B1
70+ AB <- veclen(A1 - B1 )
71+ for (i in seq_len(n1 )[- 1 ]) {
72+ ptheta <- p [i ]* theta
73+ phi <- pi / 2 + (0.5 - p [i ])* theta
74+ q <- (sin(ptheta ) / sin(phi ))/ AB
75+ D <- (1 - q )* A1 + q * B1
76+ arc [i ,] <- center1 + radius1 [i ] * normalize(D )
77+ }
78+ if (k == 1 )
79+ result <- rbind(result , arc )
80+ else
81+ result <- rbind(result [- nrow(result ), ,drop = FALSE ], arc )
82+ }
83+ }
84+ result <- rbind(result , result [1 ,])
7785 }
7886 if (plot )
79- lines3d(result [c(- 1 , - nrow(result )), , drop = FALSE ], ... )
87+ lines3d(result [c(- 1 , - nrow(result )), , drop = FALSE ], ... )
8088 else
81- result [c(- 1 , - nrow(result )), , drop = FALSE ]
89+ result [c(- 1 , - nrow(result )), , drop = FALSE ]
8290}
0 commit comments