23/07/26 20:51:37.21 0pDrMthqH.net
表面積計算と３D画像作成機能を追加。
abcd2ABCD=\(AB,BC,CA,DA,DB,DC,print=FALSE){
c=AB
a=BC
b=CA
p=DA
q=DB
r=DC
abc=sort(c(a,b,c))
if(sum(abc[1:2])<=abc[3]) return(NA)
st=c*exp(1i*acos((c^2+a^2-b^2)/(2*c*a)))
s=Re(st)
t=Im(st)
x0=(a^2+q^2-r^2)/(2*a)
P= p^2-(x0-s)^2
Q=q^2-x0^2
y0=(-P+Q+t^2)/(2*t)
if(Q-(-P+Q+t^2)^2/(4*t^2)<0) return(NA)
z0=sqrt(Q-(-P+Q+t^2)^2/(4*t^2))
D=c(x0,y0,z0)
A=c(s,t,0)
B=c(0,0,0)
C=c(a,0,0)
v=rbind(A,B,C,D)
Vol=abs(det(rbind(v[1,]-v[4,],v[2,]-v[4,],v[3,]-v[4,])))/6
ABC2S3d=\(A,B,C){
cross<-function(x,y)c(x[2]*y[3]-x[3]*y[2],x[3]*y[1]-x[1]*y[3],x[1]*y[2]-x[2]*y[1])
sqrt(sum(cross(A-C,B-C)^2))/2
}
S=ABC2S3d(A,B,C)+ABC2S3d(B,C,D)+ABC2S3d(C,D,A)+ABC2S3d(D,A,B)
if(print){
library(rgl)
open3d()
plot3d(v,type='n',col=2,size=1,xlab='x',ylab='y',zlab='z')
for(i in 1:4) points3d(v[i,],col=i,size=8)
connect <- function(x) segments3d(rbind(v[x[1],],v[x[2],]),col=8,lwd=2)
combn(4,2,connect)
}
return(list(D=D,Vol=Vol,S=S))
}