# R/crownsurface.R In Maeswrap: Wrapper Functions for MAESTRA/MAESPA

#### Defines functions CWfunCW.boxCW.cylCW.coneCW.ellipsoidCW.paraboloidSAfunS.cylS.boxS.coneS.ellipsoidS.paraboloid

#----------------------------------------------------------------------------------------------------#

# Find SA with CW and CL

# Paraboloid
S.paraboloid <- function(CW,CL){

a <- CW/2
h <- CL

pi*a/(6*h^2) * ((a^2 + 4*h^2)^1.5 - a^3) + pi*a^2
}

# Ellipsoid surface area (actually prolate/oblate spheroid)
# From Maximilian Posch.
S.ellipsoid <- function(CW,CL){

r1 <- CW/2  # rotational axis radius
r2 <- CL/2

omega <- function(x){
if(identical(x,1))x<-1.0001
if(x < 1)res <- 1+asin(sqrt(1-x^2))/(x*sqrt(1-x^2))
if(x > 1)res <- 1+log(x + sqrt(x^2-1))/(x*sqrt(x^2-1))
res
}
omega_V <- Vectorize(omega)

S <- 2*pi*r1^2*omega_V(r1/r2)

return(S)
}

# Cone, includes bottom circle.
S.cone <- function(CW,CL)pi*(CW/2)*sqrt((CW/2)^2 + CL^2) + pi*(CW/2)^2

# Box
S.box <- function(CW,CL)2*CW^2 + 4*CW*CL

# Cylinder
S.cyl <- function(CW,CL)2*pi*(CW/2)*CL + 2*pi*(CW/2)^2

SAfun <- function(CW,CL, crownform){

firstlet <- strsplit(crownform,"")[[1]][1]

res <- switch(firstlet,
c = S.cone(CW,CL),
e = S.ellipsoid(CW,CL),
p = S.paraboloid(CW,CL))

return(res)
}

#-----------------------------------------------------------------------------------------------#

# find CW if CL and SA given

# Radius of paraboloid given surface area and ratio of height to radius.
# Includes lower circle.
# b is h/r1 (height divided by radius)
CW.paraboloid <- function(SA, CWCL){

b <- 2/CWCL

2*sqrt(SA/ ( pi/(6*b^2) * ((1+4*b^2)^1.5 - 1) + pi )  )

}

CW.ellipsoid <- function(SA, CWCL){

# From Posch:
omega <- function(x){

res <- rep(NA, length(x))

res[x<1] <- 1+asin(sqrt(1-x[x<1]^2))/(x[x<1]*sqrt(1-x[x<1]^2))
res[x>1] <- 1+log(x[x>1]+sqrt(x[x>1]^2-1))/(x[x>1]*sqrt(x[x>1]^2-1))
res
}

CW <- 2*sqrt(SA/(2*pi*omega(CWCL)))

return(CW)
}

# Cone
CW.cone <- function(SA, CWCL){

b <- 2/CWCL

2*sqrt(SA / (pi*(1+sqrt(1+b^2))))

}

# Cylinder
CW.cyl <- function(SA, CWCL){

gam <- 2/CWCL

r <- sqrt(SA/(2*pi*(1+gam)))
2*r
}

# Box
CW.box <- function(SA, CWCL){

b <- 2/CWCL

2*sqrt(SA/(8+8*b))

}

CWfun <- function(shape=c("BOX","CONE","PARA","ELIP","CYL"),
SA, CWCL){

shape <- match.arg(shape)
cw <- switch(shape,
BOX = CW.box(SA, CWCL),
CONE = CW.cone(SA, CWCL),
PARA = CW.paraboloid(SA, CWCL),
ELIP = CW.ellipsoid(SA, CWCL),
CYL = CW.cyl(SA, CWCL)
)

return(cw)
}

#
# ## Test:
# CW.cone(S.cone(CW=2,CL=10), CWCL=2/10)
# CW.ellipsoid(S.ellipsoid(CW=2,CL=10), CWCL=2/10)
# CW.box(S.box(CW=2,CL=10), CWCL=2/10)
# CW.paraboloid(S.paraboloid(CW=2,CL=10), CWCL=2/10)
# CW.cyl(S.cyl(CW=2,CL=10), CWCL=2/10)

## Try the Maeswrap package in your browser

Any scripts or data that you put into this service are public.

Maeswrap documentation built on May 29, 2017, 9:35 p.m.