cutting a sphere with 3 planes and calculating 8 sub-volumes using R-language -
i want calculate 8 sub-volumes (v) resulting intersection of sphere center coordinate (10.5,10.5,10.5) , r=1 3 palnes (x=10, y=10 , z=10). using r-language , code. can calculate volume of 2 parts of sphere correctly. don't understand why not calculate correct value other 6 sub-volumes.
rm(list=ls(all=true)) # center of sphere x <- 10.5 y <- 10.5 z <- 10.5 r <- 1 # box size , number of walls in each direction u <- 3 lbx <- 30 lby <- 30 lbz <- 30 v <- u^2 w <- u^3 ui <- u - 1 # initiation of walls lx <- c() ly <- c() lz <- c() # position of walls for(i in 1:ui){ lx[i] <- i*lbx/u ly[i] <- i*lby/u lz[i] <- i*lbz/u } # calculating absolute distances between center of sphere , wall # calculation of xmin , xmid , xmax calculation of integrals dx <- sapply(x,function(x) lx - x) abs.dx <- apply(abs(dx),2,min) dx <- abs.dx*(2*(apply(replace(dx,which(dx<0),inf),2,min)==abs.dx)-1) dy <- sapply(y,function(y) ly - y) abs.dy <- apply(abs(dy),2,min) dy <- abs.dy*(2*(apply(replace(dy,which(dy<0),inf),2,min)==abs.dy)-1) dz <- sapply(z,function(z) lz - z) abs.dz <- apply(abs(dz),2,min) dz <- abs.dz*(2*(apply(replace(dz,which(dz<0),inf),2,min)==abs.dz)-1) #distance between radius center intersection wall d <- (dx)^2+(dy)^2+(dz)^2 dxy <- (dx)^2+(dy)^2 dyz <- (dy)^2+(dz)^2 dxz <- (dx)^2+(dz)^2 d1 <- apply(matrix(rbind(dxy,dyz,dxz),nrow=3,ncol=length(x)),2,min) d2 <- matrix(rbind(dx,dy,dz),nrow=3,ncol=length(x)) abs.d2 <- matrix(rbind(abs.dx,abs.dy,abs.dz),nrow=3,ncol=length(x)) d11 <- apply(abs.d2,2,min) x <- c() y <- c() z <- c() vol <- c() library('pracma') for(j in 1:length(x)) { if (d[j] < r[j]^2) { f1 <- function(x, y) (sqrt(r[j]^2 -x^2 - y^2) - abs.dz[j]) f2 <- function(x, y) (2*(sqrt(r[j]^2 - x^2 - y^2))) f3 <- function(x, y) (sqrt(r[j]^2 -x^2 - y^2)) f4 <- function(x, y) (sqrt(r[j]^2 -x^2 - y^2)- dz[j]) xmin <- (-sqrt(r[j]^2-dy[j]^2-dz[j]^2)) xmid <- dx[j] xmax <- (sqrt(r[j]^2-dy[j]^2-dz[j]^2)) ymin <- function(x) (-(sqrt(r[j]^2-dz[j]^2-x^2))) ymid <- dy[j] ymax <- function(x) (sqrt(r[j]^2-dz[j]^2-x^2)) vol1 <- integral2(f1, xmin, xmid, ymin, ymid) vol2 <- integral2(f1, xmid, xmax, ymin, ymid) vol3 <- integral2(f1, xmin, xmid, ymid, ymax) vol34 <- integral2(f2, xmin, xmid, ymid, ymax) vol5 <- integral2(f1, xmid, xmax, ymin, ymid) vol56 <- integral2(f2, xmid, xmax, ymin, ymid) vol7 <- integral2(f1, xmid, xmax, ymid, ymax) vol78 <- integral2(f2, xmid, xmax, ymid, ymax) # finding position of calculated volumes vol8.1 <- vol1$q vol[[length(vol)+1]]=vol8.1 vol8.2 <- vol2$q vol[[length(vol)+1]]=vol8.2 vol8.3 <- vol3$q vol[[length(vol)+1]]=vol8.3 vol8.4 <- vol34$q - vol8.3 vol[[length(vol)+1]]=vol8.4 vol8.5 <- vol5$q vol[[length(vol)+1]]=vol8.5 vol8.6 <- vol56$q - vol8.5 vol[[length(vol)+1]]=vol8.6 vol8.7 <- vol7$q vol[[length(vol)+1]]=vol8.7 vol8.8 <- vol78$q - vol8.7 vol[[length(vol)+1]]=vol8.8 } } newdata <- cbind(x,y,z,vol) newdata <- as.data.frame(newdata) newdata #autocad check shows different values #volume1:0.0019 - bounding box:x: 9.7929--10.0000/y:9.7929--10.0000/z:9.7929--10.0000 #volume2:0.0609 - bounding box:x: 9.6340--10.0000/y:9.6340--10.0000/z:10.000--11.2071 #volume3:0.5308 - bounding box:x: 9.5000--10.0000/y:10.000--11.36607z:10.000--11.3660 #volume4:0.5308 - bounding box:x: 10.000--11.3660/y:9.5000--10.0000/z:10.000--11.3660 #volume5:0.0609 - bounding box:x: 10.000--11.2071/y:9.6340--10.0000/z:9.6340--10.0000 #volume6:0.0609 - bounding box:x: 9.6340--10.0000/y:10.000--11.2071/z:9.6340--10.0000 #volume7:2.4117 - bounding box:x: 10.000--11.5000/y:10.000--11.5000/z:10.000--11.5000 #volume8:0.5308 - bounding box:x: 10.000--11.3660/y:10.000--11.3660/z:9.5000--10.0000 the results of r following:
vol 1 0.001900686 2 0.060883156 3 0.051049453 4 0.282041684 5 0.060883156 6 0.399931778 7 0.520997769 8 2.067153172 does have idea? why result of these 2 integrals correct
vol1 <- integral2(f1, xmin, xmid, ymin, ymid) vol2 <- integral2(f1, xmid, xmax, ymin, ymid) and rest, volumes don't know how r calculates them.
Comments
Post a Comment