Nothing
R/trio.power.R
In trio: Testing of SNPs and SNP Interactions in Case-Parent Trio Studies
Defines functions
vt4r
vt3r
vt2r
vt1r
vt4
vt3
vt2
vt1
print.trio.power
trio.power
Documented in
print.trio.power
trio.power
trio.power<-function(maf=0.5,RR=1.5,alpha=5*10^(-8),n=NULL,beta=NULL,model=c("additive","dominant","recessive"),
test=c("gTDT","Score","aTDT")){
p <- maf
if(is.null(n)&&is.null(beta))
stop("Either n or beta has to be speficied.")
if(!is.null(n) && !is.null(beta))
stop("Only one of n or beta is allowed to be specified.")
if(!is.null(n) && any(n<1)){
n<-n[-which(n<1)]
if(length(n)==0)
stop("n must be an positive integer.")
else
warning("Some n were deleted, as they were not positive integers.",call.=FALSE)
}
if(!is.null(beta) && (any(beta >= 1) | any(beta <= 0))){
del<-which((beta >= 1) | (beta <= 0))
beta<-beta[-del]
if(length(beta)==0)
stop("beta must be between 0 and 1.")
else
warning("Some beta were deleted, as they were not between 0 and 1.", call.=FALSE)
}
if(any(p>=1) | any(p<= 0)){
del<-which((p>=1) | (p<= 0))
p<-p[-del]
if(length(p)==0)
stop("maf must be between 0 and 1.")
else
warning("Some maf were deleted, as they were not between 0 and 1.",call.=FALSE)
}
if(any(alpha>=1) | any(alpha<=0)){
del<-which((alpha >= 1) | (alpha <= 0))
alpha<-alpha[-del]
if(length(alpha)==0)
stop("alpha must be between 0 and 1.")
else
warning("Some alpha were deleted, as they were not between 0 and 1.",call.=FALSE)
}
if(any(RR<=0.5)){
del<-which(RR<=0.5)
RR<-RR[-del]
if(length(RR)==0)
stop("RR must be greater than 0.5.")
else
warning("Some RR were deleted, as they were not greater than 0.5.",call.=FALSE)
}
test<-match.arg(test,several.ok=TRUE)
model<-match.arg(model,several.ok=TRUE)
single<-function(p,RR,alpha,n=NULL,beta=NULL,model,test){
outa<-FALSE
if(test=="aTDT"){
outa<-TRUE
model<-"additive"
test<-"Score"
}
if(test=="gTDT"){
if(model=="additive"){
q<-1-p
R<-(RR*2-1)*p^2+2*RR*p*q+q^2
prob<-c(2*p^3*q*(RR*2-1),2*p^3*q*RR,p^2*q^2*(RR*2-1),2*p^2*q^2*RR,p^2*q^2,2*p*q^3*RR,2*p*q^3)
s<-prob/R
u<-c(1,0,2,1,0,1,0)
v<-c(0,1,0,1,2,0,1)
e1<-t(u)%*%s
e2<-t(v)%*%s
e1pe2<-2*p*q/R*(p*(RR*2-1)+RR+q)
e1me2<-2*p*q/R*(p*(RR*2-1)+(1-2*p)*RR-q)
d1<-sum(u^2*s)
d12<-sum(u*v*s)
d2<-sum(v^2*s)
e1e2<-e1*e2
me<-logit(e2/(e1pe2))*sqrt(e1*e2)/sqrt(e1pe2)
mu<-me
vt1<-(sqrt(e2/e1)*((log(e2/e1)/2-1)))*sqrt(e1pe2)-log((e2/e1))*sqrt(e1e2)/(2*sqrt(e1pe2))
vt2<-(sqrt(e1/e2)*(1+(log(e2/e1)/2)))*sqrt(e1pe2)-log((e2/e1))*sqrt(e1e2)/(2*sqrt(e1pe2))
vt<-c(vt1,vt1*vt2,vt2)
V<-vt[1]^2*(d1-e1^2)+vt[2]*2*(d12-e1e2)+vt[3]^2*(d2-e2^2)
va<-V/e1pe2^2
std<-sqrt(va)
}
else{
if(model=="dominant"){
q<-1-p
R<-RR*(p^2+2*p*q)+q^2
prob<-c(2*p*q^3,2*p*q^3*RR,p^2*q^2,3*p^2*q^2*RR)/R
d1<-prob[1]
d2<-prob[2]
d3<-prob[3]
d4<-prob[4]
tmp<- d1/3 - d2 + d3 - d4/3
tmp2<-2*d1 + 2*d3
s<- ((d1 + d2)*(((tmp/tmp2)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
tmp/tmp2))/(((tmp/tmp2)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
tmp/tmp2 + 1) + ((d3 + d4)*(((tmp/tmp2)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - tmp/tmp2))/(((tmp/tmp2)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - tmp/tmp2 + 1/3) -
((d1 + d2)*(((tmp/tmp2)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
tmp/tmp2)^2)/(((tmp/tmp2)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2)-
tmp/tmp2 + 1)^2 - ((d3 + d4)*(((tmp/tmp2)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - tmp/tmp2)^2)/(((tmp/tmp2)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - tmp/tmp2 + 1/3)^2
be <-log(((tmp/tmp2)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) - tmp/tmp2)
chi<-be*sqrt(s)
mu<-chi
vt1<-vt1(d1,d2,d3,d4)
vt2<-vt2(d1,d2,d3,d4)
vt3<-vt3(d1,d2,d3,d4)
vt4<-vt4(d1,d2,d3,d4)
vt<-c(vt1,vt2,vt3,vt4)
}
if(model=="recessive"){
q<-1-p
R<-RR*p^2+2*p*q+q^2
prob<-c(2*p^3*q,2*p^3*q*RR,3*p^2*q^2,p^2*q^2*RR)/R
r1<-prob[1]
r2<-prob[2]
r3<-prob[3]
r4<-prob[4]
tmp<-3*r1 - r2 + r3 - 3*r4
tmp2<-2*r1 + 2*r3
s <- - ((r1 + r2)*(tmp/tmp2 - ((3*r2 + 3*r4)/(r1 + r3) +
1/tmp2^2*tmp^2)^(1/2)))/(((3*r2 + 3*r4)/(r1 + r3) +
(tmp/tmp2)^2)^(1/2) - tmp/tmp2 + 1)^2 - ((r3 + r4)*((3*tmp)/tmp2 -
3*((3*r2 + 3*r4)/(r1 + r3) + (tmp/tmp2)^2)^(1/2)))/(((3*r2 + 3*r4)/(r1 + r3) +
(tmp/tmp2)^2)^(1/2) - tmp/tmp2 + 3)^2
be<- log(((3*r2 + 3*r4)/(r1 + r3) + (tmp/tmp2)^2)^(1/2) - tmp/tmp2)
chi<-be*sqrt(s)
mu<-chi
vt1r<-vt1r(r1,r2,r3,r4)
vt2r<-vt2r(r1,r2,r3,r4)
vt3r<-vt3r(r1,r2,r3,r4)
vt4r<-vt4r(r1,r2,r3,r4)
vt<-c(vt1r,vt2r,vt3r,vt4r)
}
varmat<-prob%*%t((-prob))
diag(varmat)<-prob*(1-prob)
va<-sum((varmat%*%vt)*vt)
std<-sqrt(va)
}
}
if(test=="Score"){
if(model=="additive"){
q<-1-p
R<-(RR*2-1)*p^2+2*RR*p*q+q^2
prob<-c(2*p^3*q*(RR*2-1),2*p^3*q*RR,p^2*q^2*(RR*2-1),2*p^2*q^2*RR,p^2*q^2,2*p*q^3*RR,2*p*q^3)/R
u<-c(1,0,2,1,0,1,0)
v<-c(0,1,0,1,2,0,1)
mu<-(t(u)%*%prob-t(v)%*%prob)/sqrt(t(u)%*%prob+t(v)%*%prob)
e1pe2<-2*p*q/R*(p*(RR*2-1)+RR+q)
e1me2<-2*p*q/R*(p*(RR*2-1)+(1-2*p)*RR-q)
d1<-2*p*q/R*((1-q^2)*(RR*2-1)+(1-2*p*q)*RR+(1-p^2))
d2<-2*p*q/R*((1-q^2)*(RR*2-1)+(1+2*p*q)*RR+(1-p^2))
d12<-2*p*q/R*(-(1-q^2)*(RR*2-1)+(p-q)*RR+(1-p^2))
vari<-(d1-e1me2^2/4)/e1pe2+(e1me2*d12)/e1pe2^2+1/4*(e1me2)^2*d2/e1pe2^3
std<-sqrt(vari)
}
else{
if(model=="dominant"){
q<-1-p
R<-RR*(p^2+2*p*q)+q^2
prob<-c(2*p*q^3,2*p*q^3*RR,p^2*q^2,3*p^2*q^2*RR)/R
d1<-prob[1]
d2<-prob[2]
d3<-prob[3]
d4<-prob[4]
mu<-(2*(d2-d1)+d4-3*d3)/sqrt(4*d2+4*d1+3*d3+3*d4)
div<-sqrt((4*d1+4*d2+3*d3+3*d4)^3)
vt1<- -4*d1-12*d2-8*d4
vt2<-12*d1+4*d2+12*d3+4*d4
vt3<- -9*d1-15*d2-9/2*d3-21/2*d4
vt4<-7*d1+d2+15/2*d3+3/2*d4
vt<-c(vt1,vt2,vt3,vt4)/div
}
if(model=="recessive"){
q<-1-p
R<-RR*p^2+2*p*q+q^2
prob<-c(2*p^3*q,2*p^3*q*RR,3*p^2*q^2,p^2*q^2*RR)/R
r1<-prob[1]
r2<-prob[2]
r3<-prob[3]
r4<-prob[4]
mu<-(2*(r2-r1)+3*r4-r3)/(sqrt(4*(r1+r2)+3*(r3+r4)))
vt1r<- -(4*(r1 + 3*r2 + r3 + 3*r4))/(4*r1 + 4*r2 + 3*r3 + 3*r4)^(3/2)
vt2r<-(4*(3*r1 + r2 + 2*r3))/(4*r1 + 4*r2 + 3*r3 + 3*r4)^(3/2)
vt3r<- -(r1 + 7*r2 + (3*r3)/2 + (15*r4)/2)/(4*r1 + 4*r2 + 3*r3 + 3*r4)^(3/2)
vt4r<-(3*(10*r1 + 6*r2 + 7*r3 + 3*r4))/(2*(4*r1 + 4*r2 + 3*r3 + 3*r4)^(3/2))
vt<-c(vt1r,vt2r,vt3r,vt4r)
}
varmat<-prob%*%t((-prob))
diag(varmat)<-prob*(1-prob)
va<-sum((varmat%*%vt)*vt)
std<-sqrt(va)
}
}
if(!is.null(n)){
z<-qnorm(1-alpha/2)
power<-pnorm(-z,sd=std,mean=sqrt(n)*mu)+(1-pnorm(z,sd=std,mean=sqrt(n)*mu))
if(outa){
model<-"multiplicative"
test="aTDT"
}
out<-c(model=model,power=power,alpha=alpha,n=n,test=test,RR=RR,p=p,calc="power")
return(out)
}
if(!is.null(beta)){
size<-as.vector(((std*qnorm(beta)+qnorm(1-alpha/2))/mu)^2)
if(outa){
model<-"multiplicative"
test<-"aTDT"
}
out<-c(model=model,size=size,beta=beta,alpha=alpha,test=test,RR=RR,p=p,calc="size")
return(out)
}
}
if(!is.null(n)){
vars<-expand.grid(n,p,RR,alpha,model,test,stringsAsFactors=FALSE)
wa<-which(vars[,6]=="aTDT")
le<-length(wa)/length(model)
if(le>=1){
aw<-wa[-c(1:le)]
vars<-vars[-aw,]
}
out<-t(mapply(single,n=vars[,1],p=vars[,2],RR=vars[,3],alpha=vars[,4],model=vars[,5],test=vars[,6]))
out<-as.data.frame(out)
out$n<-as.numeric(as.character(out$n))
out$alpha<-as.numeric(as.character(out$alpha))
out$RR<-as.numeric(as.character(out$RR))
out$p<-as.numeric(as.character(out$p))
out$power<-as.numeric(as.character(out$power))
}
if(!is.null(beta)){
vars<-expand.grid(beta,p,RR,alpha,model,test,stringsAsFactors=FALSE)
wa<-which(vars[,6]=="aTDT")
le<-length(wa)/length(model)
if(le>=1){
aw<-wa[-c(1:le)]
vars<-vars[-aw,]}
out<-t(mapply(single,beta=vars[,1],p=vars[,2],RR=vars[,3],alpha=vars[,4],model=vars[,5],test=vars[,6]))
out<-as.data.frame(out)
out$size<-as.numeric(as.character(out$size))
out$beta<-as.numeric(as.character(out$beta))
out$alpha<-as.numeric(as.character(out$alpha))
out$RR<-as.numeric(as.character(out$RR))
out$p<-as.numeric(as.character(out$p))
}
class(out) <- "trio.power"
out
}
print.trio.power<-function(x,digits=4,...){
if(x$calc[1]=="power"){
power<-format.pval(as.numeric(as.character(x$power)),digits=digits)
out<-data.frame(Test=x$test,Model=x$model, MAF=x$p,alpha=x$alpha,RR=x$RR,Trios=x$n,power=power)
cat(" Trio studies power calculation", "\n\n")
print(format(out,digits=digits))
}
if(x$calc[1]=="size"){
size<-ceiling(as.numeric(as.character(x$size)))
out<-data.frame(Test=x$test,Model=x$model,MAF=x$p,alpha=x$alpha,RR=x$RR,beta=x$beta,Trios=size)
cat(" Trio studies sample size calculation", "\n\n")
print(format(out,digits=digits,trim=FALSE)
)
if(any(x$RR==1)){
warning("RR should differ from one to get a finite sample size.",call.=FALSE)}
}
}
vt1<-function(d1,d2,d3,d4){
h<-(1/3*d1-d2+d3-1/3*d4)/(2*(d1+d3))
be<-log(sqrt((d2+d4)/(3*(d1+d3))+h^2)-h)
summe<-((d1 + d2)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3)))/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1) + ((d3 + d4)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - (d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3)))/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - (d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1/3) - ((d1 + d2)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - (d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))^2)/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - (d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1)^2 - ((d3 + d4)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - (d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))^2)/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - (d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1/3)^2
ablb<- -(1/(3*(2*d1 + 2*d3)) + ((4*(d1/3 - d2 + d3 - d4/3)^2)/(2*d1 + 2*d3)^3 - ((2*d1)/9 - (2*d2)/3 + (2*d3)/3 - (2*d4)/9)/(2*d1 + 2*d3)^2 +
(3*(d2 + d4))/(3*d1 + 3*d3)^2)/(2*((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2)) -
(2*(d1/3 - d2 + d3 - d4/3))/(2*d1 + 2*d3)^2)/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))
abls<-((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1) - ((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1)^2 - (d1 + d2)*(1/3/(2*d1 + 2*d3) +
1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 -
1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) + 3*(d2 + d4)/(3*d1 + 3*d3)^2) -
2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1) - (d3 + d4)*(1/3/(2*d1 + 2*d3) + 1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 - 1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) +
3*(d2 + d4)/(3*d1 + 3*d3)^2) - 2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3) +
3*(d1 + d2)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))*(1/3/(2*d1 + 2*d3) + 1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 - 1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) +
3*(d2 + d4)/(3*d1 + 3*d3)^2) - 2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1)^2 +
3*(d3 + d4)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))*(1/3/(2*d1 + 2*d3) + 1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 - 1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) +
3*(d2 + d4)/(3*d1 + 3*d3)^2) - 2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3)^2 -
2*(d1 + d2)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2*(1/3/(2*d1 + 2*d3) + 1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 - 1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) +
3*(d2 + d4)/(3*d1 + 3*d3)^2) - 2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1)^3 -
2*(d3 + d4)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2*(1/3/(2*d1 + 2*d3) + 1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 - 1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) +
3*(d2 + d4)/(3*d1 + 3*d3)^2) - 2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3)^3
abl<-ablb*sqrt(summe)+be*1/(2*sqrt(summe))*abls
abl
}
vt2<-function(d1,d2,d3,d4){
h<-(1/3*d1-d2+d3-1/3*d4)/(2*(d1+d3))
be<-log(sqrt((d2+d4)/(3*(d1+d3))+h^2)-h)
summe<-((d1 + d2)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3)))/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1) + ((d3 + d4)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3)))/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1/3) - ((d1 + d2)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))^2)/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1)^2 - ((d3 + d4)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))^2)/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1/3)^2
ablb<- -((((2*d1)/3 - 2*d2 + 2*d3 - (2*d4)/3)/(2*d1 + 2*d3)^2 - 1/(3*d1 + 3*d3))/(2*((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)) - 1/(2*d1 + 2*d3))/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))
abls<-((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1) - ((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1)^2 - (1/2*(1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) -
1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3))*(d1 + d2)/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1) - (1/2*(1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) -
1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3))*(d3 + d4)/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3) + 3*(1/2*(1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) -
1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3))*(d1 + d2)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1)^2 + 3*(1/2*(1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) -
1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3))*(d3 + d4)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3)^2 - 2*(1/2*(1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) -
1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3))*(d1 + d2)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1)^3 - 2*(1/2*(1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) -
1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3))*(d3 + d4)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3)^3
abl<-ablb*sqrt(summe)+be*1/(2*sqrt(summe))*abls
abl
}
vt3<-function(d1,d2,d3,d4){
h<-(1/3*d1-d2+d3-1/3*d4)/(2*(d1+d3))
be<-log(sqrt((d2+d4)/(3*(d1+d3))+h^2)-h)
summe<-((d1 + d2)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3)))/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1) + ((d3 + d4)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3)))/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1/3) - ((d1 + d2)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - (d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))^2)/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - (d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1)^2 -
((d3 + d4)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))^2)/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1/3)^2
ablb<- -(1/(2*d1 + 2*d3) + ((4*(d1/3 - d2 + d3 - d4/3)^2)/(2*d1 + 2*d3)^3 - ((2*d1)/3 - 2*d2 + 2*d3 - (2*d4)/3)/(2*d1 + 2*d3)^2 +
(3*(d2 + d4))/(3*d1 + 3*d3)^2)/(2*((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2)) -
(2*(d1/3 - d2 + d3 - d4/3))/(2*d1 + 2*d3)^2)/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))
abls<-((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3) - ((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3)^2 - (d1 + d2)*(1/(2*d1 + 2*d3) +
1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 -
1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) + 3*(d2 + d4)/(3*d1 + 3*d3)^2) -
2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1) - (d3 + d4)*(1/(2*d1 + 2*d3) + 1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 - 1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) +
3*(d2 + d4)/(3*d1 + 3*d3)^2) - 2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3) -
2*(d1 + d2)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2*(1/(2*d1 + 2*d3) + 1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 - 1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) +
3*(d2 + d4)/(3*d1 + 3*d3)^2) - 2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1)^3 -
2*(d3 + d4)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2*(1/(2*d1 + 2*d3) + 1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 - 1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) +
3*(d2 + d4)/(3*d1 + 3*d3)^2) - 2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3)^3 +
3*(d1 + d2)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))*(1/(2*d1 + 2*d3) + 1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 - 1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) +
3*(d2 + d4)/(3*d1 + 3*d3)^2) - 2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1)^2 +
3*(d3 + d4)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))*(1/(2*d1 + 2*d3) + 1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 - 1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) +
3*(d2 + d4)/(3*d1 + 3*d3)^2) - 2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3)^2
abl<-ablb*sqrt(summe)+be*1/(2*sqrt(summe))*abls
abl
}
vt4<-function(d1,d2,d3,d4){
h<-(1/3*d1-d2+d3-1/3*d4)/(2*(d1+d3))
be<-log(sqrt((d2+d4)/(3*(d1+d3))+h^2)-h)
summe<-((d1 + d2)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3)))/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1) + ((d3 + d4)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3)))/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1/3) - ((d1 + d2)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))^2)/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1)^2 - ((d3 + d4)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - (d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))^2)/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - (d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1/3)^2
ablb<- -((((2*d1)/9 - (2*d2)/3 + (2*d3)/3 - (2*d4)/9)/(2*d1 + 2*d3)^2 - 1/(3*d1 + 3*d3))/(2*((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)) - 1/(3*(2*d1 + 2*d3)))/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))
abls<-((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3) - ((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3)^2 -
(1/2*(1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) - 1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/3/(2*d1 + 2*d3))*(d1 + d2)/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1) -
(1/2*(1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) - 1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/3/(2*d1 + 2*d3))*(d3 + d4)/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3) +
3*(1/2*(1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) - 1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/3/(2*d1 + 2*d3))*(d1 + d2)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1)^2 +
3*(1/2*(1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) - 1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/3/(2*d1 + 2*d3))*(d3 + d4)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3)^2 -
2*(1/2*(1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) - 1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/3/(2*d1 + 2*d3))*(d1 + d2)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1)^3 -
2*(1/2*(1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) - 1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/3/(2*d1 + 2*d3))*(d3 + d4)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3)^3
abl<-ablb*sqrt(summe)+be*1/(2*sqrt(summe))*abls
abl
}
vt1r<-function(r1,r2,r3,r4){
h<-(3*r1-r2+r3-3*r4)/(2*(r1+r3))
be<-log(sqrt(3*(r2+r4)/(r1+r3)+h^2)-h)
summe<- -((r1 + r2)*((3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) -
((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)))/(((3*r2 + 3*r4)/(r1 + r3) +
(3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2) - (3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) + 1)^2 -
((r3 + r4)*((3*(3*r1 - r2 + r3 - 3*r4))/(2*r1 + 2*r3) - 3*((3*r2 + 3*r4)/(r1 + r3) +
(3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2)))/(((3*r2 + 3*r4)/(r1 + r3) + (3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2) -
(3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) + 3)^2
ablb<- (3/(2*r1 + 2*r3) - (2*(3*r1 - r2 + r3 - 3*r4))/(2*r1 + 2*r3)^2 + ((3*r2 + 3*r4)/(r1 + r3)^2 -
(18*r1 - 6*r2 + 6*r3 - 18*r4)/(2*r1 + 2*r3)^2 + (4*(3*r1 - r2 + r3 - 3*r4)^2)/(2*r1 + 2*r3)^3)/(2*((3*r2 + 3*r4)/(r1 + r3) +
(3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2)))/((3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) -
((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))
abls<- -(1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) - ((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 1)^2 - (r1 + r2)/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 1)^2*(3/(2*r1 + 2*r3) - 2/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4) + 1/2/((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)*((3*r2 + 3*r4)/(r1 + r3)^2 - 1/(2*r1 + 2*r3)^2*(18*r1 - 6*r2 + 6*r3 - 18*r4) +
4/(2*r1 + 2*r3)^3*(3*r1 - r2 + r3 - 3*r4)^2)) - (r3 + r4)/(((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) - 1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 3)^2*(9/(2*r1 + 2*r3) -
6/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4) + 3/2/((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)*((3*r2 + 3*r4)/(r1 + r3)^2 - 1/(2*r1 + 2*r3)^2*(18*r1 - 6*r2 + 6*r3 - 18*r4) +
4/(2*r1 + 2*r3)^3*(3*r1 - r2 + r3 - 3*r4)^2)) - 2*(r1 + r2)*(1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) - ((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 1)^3*(3/(2*r1 + 2*r3) - 2/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4) +
1/2/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)*((3*r2 + 3*r4)/(r1 + r3)^2 -
1/(2*r1 + 2*r3)^2*(18*r1 - 6*r2 + 6*r3 - 18*r4) + 4/(2*r1 + 2*r3)^3*(3*r1 - r2 + r3 - 3*r4)^2)) -
2*(r3 + r4)*(3/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) - 3*((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 3)^3*(3/(2*r1 + 2*r3) - 2/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4) +
1/2/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)*((3*r2 + 3*r4)/(r1 + r3)^2 -
1/(2*r1 + 2*r3)^2*(18*r1 - 6*r2 + 6*r3 - 18*r4) + 4/(2*r1 + 2*r3)^3*(3*r1 - r2 + r3 - 3*r4)^2))
abl<-ablb*sqrt(summe)+be*1/(2*sqrt(summe))*abls
abl
}
vt2r<-function(r1,r2,r3,r4){
h<-(3*r1-r2+r3-3*r4)/(2*(r1+r3))
be<-log(sqrt(3*(r2+r4)/(r1+r3)+h^2)-h)
summe<- -((r1 + r2)*((3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) - ((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)))/(((3*r2 + 3*r4)/(r1 + r3) + (3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2) -
(3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) + 1)^2 - ((r3 + r4)*((3*(3*r1 - r2 + r3 - 3*r4))/(2*r1 + 2*r3) -
3*((3*r2 + 3*r4)/(r1 + r3) + (3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2)))/(((3*r2 + 3*r4)/(r1 + r3) +
(3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2) - (3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) + 3)^2
ablb<- -(1/(2*r1 + 2*r3) - ((3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3)^2 - 3/(2*(r1 + r3)))/((3*r2 + 3*r4)/(r1 + r3) +
(3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2))/((3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) -
((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))
abls<- -(1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) - ((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 1)^2 - (r1 + r2)*(1/2*(1/(2*r1 + 2*r3)^2*(6*r1 - 2*r2 + 2*r3 - 6*r4) -
3/(r1 + r3))/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 1)^2 - (r3 + r4)*(3/2*(1/(2*r1 + 2*r3)^2*(6*r1 - 2*r2 + 2*r3 - 6*r4) -
3/(r1 + r3))/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
3/(2*r1 + 2*r3))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 3)^2 - 2*(r1 + r2)*(1/2*(1/(2*r1 + 2*r3)^2*(6*r1 - 2*r2 + 2*r3 - 6*r4) -
3/(r1 + r3))/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3))*(1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) - ((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))/(((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) - 1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 1)^3 -
2*(r3 + r4)*(1/2*(1/(2*r1 + 2*r3)^2*(6*r1 - 2*r2 + 2*r3 - 6*r4) - 3/(r1 + r3))/((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) - 1/(2*r1 + 2*r3))*(3/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) -
3*((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))/(((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) - 1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 3)^3
abl<-ablb*sqrt(summe)+be*1/(2*sqrt(summe))*abls
abl
}
vt3r<-function(r1,r2,r3,r4){
h<-(3*r1-r2+r3-3*r4)/(2*(r1+r3))
be<-log(sqrt(3*(r2+r4)/(r1+r3)+h^2)-h)
summe<- -((r1 + r2)*((3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) - ((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)))/(((3*r2 + 3*r4)/(r1 + r3) + (3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2) -
(3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) + 1)^2 - ((r3 + r4)*((3*(3*r1 - r2 + r3 - 3*r4))/(2*r1 + 2*r3) - 3*((3*r2 + 3*r4)/(r1 + r3) +
(3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2)))/(((3*r2 + 3*r4)/(r1 + r3) + (3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2) -
(3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) + 3)^2
ablb <- ((((3*r2)/2 + (3*r4)/2)/(r1 + r3)^2 - (3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3)^2 +
(2*(3*r1 - r2 + r3 - 3*r4)^2)/(2*r1 + 2*r3)^3)/((3*r2 + 3*r4)/(r1 + r3) + (3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2) +
1/(2*r1 + 2*r3) - (6*r1 - 2*r2 + 2*r3 - 6*r4)/(2*r1 + 2*r3)^2)/((3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) -
((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))
abls<- -(3/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) - 3*((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 3)^2 - (r1 + r2)*(1/(2*r1 + 2*r3) - 2/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4) +
1/2/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)*((3*r2 + 3*r4)/(r1 + r3)^2 -
1/(2*r1 + 2*r3)^2*(6*r1 - 2*r2 + 2*r3 - 6*r4) + 4/(2*r1 + 2*r3)^3*(3*r1 - r2 + r3 - 3*r4)^2))/(((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) - 1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 1)^2 -
(r3 + r4)/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 3)^2*(3/(2*r1 + 2*r3) - 6/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4) + 3/2/((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)*((3*r2 + 3*r4)/(r1 + r3)^2 - 1/(2*r1 + 2*r3)^2*(6*r1 - 2*r2 + 2*r3 - 6*r4) +
4/(2*r1 + 2*r3)^3*(3*r1 - r2 + r3 - 3*r4)^2)) - 2*(r1 + r2)*(1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) -
((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))*(1/(2*r1 + 2*r3) -
2/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4) + 1/2/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)*((3*r2 + 3*r4)/(r1 + r3)^2 -
1/(2*r1 + 2*r3)^2*(6*r1 - 2*r2 + 2*r3 - 6*r4) + 4/(2*r1 + 2*r3)^3*(3*r1 - r2 + r3 - 3*r4)^2))/(((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) - 1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 1)^3 -
2*(r3 + r4)*(3/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) - 3*((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))*(1/(2*r1 + 2*r3) - 2/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4) +
1/2/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)*((3*r2 + 3*r4)/(r1 + r3)^2 -
1/(2*r1 + 2*r3)^2*(6*r1 - 2*r2 + 2*r3 - 6*r4) + 4/(2*r1 + 2*r3)^3*(3*r1 - r2 + r3 - 3*r4)^2))/(((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) - 1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 3)^3
abl<-ablb*sqrt(summe)+be*1/(2*sqrt(summe))*abls
abl
}
vt4r<-function(r1,r2,r3,r4){
h<-(3*r1-r2+r3-3*r4)/(2*(r1+r3))
be<-log(sqrt(3*(r2+r4)/(r1+r3)+h^2)-h)
summe<- -((r1 + r2)*((3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) -
((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)))/(((3*r2 + 3*r4)/(r1 + r3) +
(3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2) - (3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) + 1)^2 -
((r3 + r4)*((3*(3*r1 - r2 + r3 - 3*r4))/(2*r1 + 2*r3) - 3*((3*r2 + 3*r4)/(r1 + r3) +
(3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2)))/(((3*r2 + 3*r4)/(r1 + r3) + (3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2) -
(3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) + 3)^2
ablb<-(((9*r1 - 3*r2 + 3*r3 - 9*r4)/(2*r1 + 2*r3)^2 - 3/(2*(r1 + r3)))/((3*r2 + 3*r4)/(r1 + r3) +
(3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2) - 3/(2*r1 + 2*r3))/((3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) -
((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))
abls<- -(3/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) - 3*((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 3)^2 - (r1 + r2)*(1/2*(1/(2*r1 + 2*r3)^2*(18*r1 - 6*r2 + 6*r3 - 18*r4) -
3/(r1 + r3))/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
3/(2*r1 + 2*r3))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 1)^2 - (r3 + r4)*(3/2*(1/(2*r1 + 2*r3)^2*(18*r1 - 6*r2 + 6*r3 - 18*r4) -
3/(r1 + r3))/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
9/(2*r1 + 2*r3))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 3)^2 - 2*(r1 + r2)*(1/2*(1/(2*r1 + 2*r3)^2*(18*r1 - 6*r2 + 6*r3 - 18*r4) -
3/(r1 + r3))/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
3/(2*r1 + 2*r3))*(1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) - ((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 1)^3 - 2*(r3 + r4)*(1/2*(1/(2*r1 + 2*r3)^2*(18*r1 - 6*r2 + 6*r3 - 18*r4) -
3/(r1 + r3))/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
3/(2*r1 + 2*r3))*(3/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) - 3*((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 3)^3
abl<-ablb*sqrt(summe)+be*1/(2*sqrt(summe))*abls
abl
}
Try the trio package in your browser
Any scripts or data that you put into this service are public.
trio documentation built on Nov. 8, 2020, 7:41 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions vt4r vt3r vt2r vt1r vt4 vt3 vt2 vt1 print.trio.power trio.power
Documented in print.trio.power trio.power
trio.power<-function(maf=0.5,RR=1.5,alpha=5*10^(-8),n=NULL,beta=NULL,model=c("additive","dominant","recessive"),
test=c("gTDT","Score","aTDT")){
p <- maf
if(is.null(n)&&is.null(beta))
stop("Either n or beta has to be speficied.")
if(!is.null(n) && !is.null(beta))
stop("Only one of n or beta is allowed to be specified.")
if(!is.null(n) && any(n<1)){
n<-n[-which(n<1)]
if(length(n)==0)
stop("n must be an positive integer.")
else
warning("Some n were deleted, as they were not positive integers.",call.=FALSE)
}
if(!is.null(beta) && (any(beta >= 1) | any(beta <= 0))){
del<-which((beta >= 1) | (beta <= 0))
beta<-beta[-del]
if(length(beta)==0)
stop("beta must be between 0 and 1.")
else
warning("Some beta were deleted, as they were not between 0 and 1.", call.=FALSE)
}
if(any(p>=1) | any(p<= 0)){
del<-which((p>=1) | (p<= 0))
p<-p[-del]
if(length(p)==0)
stop("maf must be between 0 and 1.")
else
warning("Some maf were deleted, as they were not between 0 and 1.",call.=FALSE)
}
if(any(alpha>=1) | any(alpha<=0)){
del<-which((alpha >= 1) | (alpha <= 0))
alpha<-alpha[-del]
if(length(alpha)==0)
stop("alpha must be between 0 and 1.")
else
warning("Some alpha were deleted, as they were not between 0 and 1.",call.=FALSE)
}
if(any(RR<=0.5)){
del<-which(RR<=0.5)
RR<-RR[-del]
if(length(RR)==0)
stop("RR must be greater than 0.5.")
else
warning("Some RR were deleted, as they were not greater than 0.5.",call.=FALSE)
}
test<-match.arg(test,several.ok=TRUE)
model<-match.arg(model,several.ok=TRUE)
single<-function(p,RR,alpha,n=NULL,beta=NULL,model,test){
outa<-FALSE
if(test=="aTDT"){
outa<-TRUE
model<-"additive"
test<-"Score"
}
if(test=="gTDT"){
if(model=="additive"){
q<-1-p
R<-(RR*2-1)*p^2+2*RR*p*q+q^2
prob<-c(2*p^3*q*(RR*2-1),2*p^3*q*RR,p^2*q^2*(RR*2-1),2*p^2*q^2*RR,p^2*q^2,2*p*q^3*RR,2*p*q^3)
s<-prob/R
u<-c(1,0,2,1,0,1,0)
v<-c(0,1,0,1,2,0,1)
e1<-t(u)%*%s
e2<-t(v)%*%s
e1pe2<-2*p*q/R*(p*(RR*2-1)+RR+q)
e1me2<-2*p*q/R*(p*(RR*2-1)+(1-2*p)*RR-q)
d1<-sum(u^2*s)
d12<-sum(u*v*s)
d2<-sum(v^2*s)
e1e2<-e1*e2
me<-logit(e2/(e1pe2))*sqrt(e1*e2)/sqrt(e1pe2)
mu<-me
vt1<-(sqrt(e2/e1)*((log(e2/e1)/2-1)))*sqrt(e1pe2)-log((e2/e1))*sqrt(e1e2)/(2*sqrt(e1pe2))
vt2<-(sqrt(e1/e2)*(1+(log(e2/e1)/2)))*sqrt(e1pe2)-log((e2/e1))*sqrt(e1e2)/(2*sqrt(e1pe2))
vt<-c(vt1,vt1*vt2,vt2)
V<-vt[1]^2*(d1-e1^2)+vt[2]*2*(d12-e1e2)+vt[3]^2*(d2-e2^2)
va<-V/e1pe2^2
std<-sqrt(va)
}
else{
if(model=="dominant"){
q<-1-p
R<-RR*(p^2+2*p*q)+q^2
prob<-c(2*p*q^3,2*p*q^3*RR,p^2*q^2,3*p^2*q^2*RR)/R
d1<-prob[1]
d2<-prob[2]
d3<-prob[3]
d4<-prob[4]
tmp<- d1/3 - d2 + d3 - d4/3
tmp2<-2*d1 + 2*d3
s<- ((d1 + d2)*(((tmp/tmp2)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
tmp/tmp2))/(((tmp/tmp2)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
tmp/tmp2 + 1) + ((d3 + d4)*(((tmp/tmp2)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - tmp/tmp2))/(((tmp/tmp2)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - tmp/tmp2 + 1/3) -
((d1 + d2)*(((tmp/tmp2)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
tmp/tmp2)^2)/(((tmp/tmp2)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2)-
tmp/tmp2 + 1)^2 - ((d3 + d4)*(((tmp/tmp2)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - tmp/tmp2)^2)/(((tmp/tmp2)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - tmp/tmp2 + 1/3)^2
be <-log(((tmp/tmp2)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) - tmp/tmp2)
chi<-be*sqrt(s)
mu<-chi
vt1<-vt1(d1,d2,d3,d4)
vt2<-vt2(d1,d2,d3,d4)
vt3<-vt3(d1,d2,d3,d4)
vt4<-vt4(d1,d2,d3,d4)
vt<-c(vt1,vt2,vt3,vt4)
}
if(model=="recessive"){
q<-1-p
R<-RR*p^2+2*p*q+q^2
prob<-c(2*p^3*q,2*p^3*q*RR,3*p^2*q^2,p^2*q^2*RR)/R
r1<-prob[1]
r2<-prob[2]
r3<-prob[3]
r4<-prob[4]
tmp<-3*r1 - r2 + r3 - 3*r4
tmp2<-2*r1 + 2*r3
s <- - ((r1 + r2)*(tmp/tmp2 - ((3*r2 + 3*r4)/(r1 + r3) +
1/tmp2^2*tmp^2)^(1/2)))/(((3*r2 + 3*r4)/(r1 + r3) +
(tmp/tmp2)^2)^(1/2) - tmp/tmp2 + 1)^2 - ((r3 + r4)*((3*tmp)/tmp2 -
3*((3*r2 + 3*r4)/(r1 + r3) + (tmp/tmp2)^2)^(1/2)))/(((3*r2 + 3*r4)/(r1 + r3) +
(tmp/tmp2)^2)^(1/2) - tmp/tmp2 + 3)^2
be<- log(((3*r2 + 3*r4)/(r1 + r3) + (tmp/tmp2)^2)^(1/2) - tmp/tmp2)
chi<-be*sqrt(s)
mu<-chi
vt1r<-vt1r(r1,r2,r3,r4)
vt2r<-vt2r(r1,r2,r3,r4)
vt3r<-vt3r(r1,r2,r3,r4)
vt4r<-vt4r(r1,r2,r3,r4)
vt<-c(vt1r,vt2r,vt3r,vt4r)
}
varmat<-prob%*%t((-prob))
diag(varmat)<-prob*(1-prob)
va<-sum((varmat%*%vt)*vt)
std<-sqrt(va)
}
}
if(test=="Score"){
if(model=="additive"){
q<-1-p
R<-(RR*2-1)*p^2+2*RR*p*q+q^2
prob<-c(2*p^3*q*(RR*2-1),2*p^3*q*RR,p^2*q^2*(RR*2-1),2*p^2*q^2*RR,p^2*q^2,2*p*q^3*RR,2*p*q^3)/R
u<-c(1,0,2,1,0,1,0)
v<-c(0,1,0,1,2,0,1)
mu<-(t(u)%*%prob-t(v)%*%prob)/sqrt(t(u)%*%prob+t(v)%*%prob)
e1pe2<-2*p*q/R*(p*(RR*2-1)+RR+q)
e1me2<-2*p*q/R*(p*(RR*2-1)+(1-2*p)*RR-q)
d1<-2*p*q/R*((1-q^2)*(RR*2-1)+(1-2*p*q)*RR+(1-p^2))
d2<-2*p*q/R*((1-q^2)*(RR*2-1)+(1+2*p*q)*RR+(1-p^2))
d12<-2*p*q/R*(-(1-q^2)*(RR*2-1)+(p-q)*RR+(1-p^2))
vari<-(d1-e1me2^2/4)/e1pe2+(e1me2*d12)/e1pe2^2+1/4*(e1me2)^2*d2/e1pe2^3
std<-sqrt(vari)
}
else{
if(model=="dominant"){
q<-1-p
R<-RR*(p^2+2*p*q)+q^2
prob<-c(2*p*q^3,2*p*q^3*RR,p^2*q^2,3*p^2*q^2*RR)/R
d1<-prob[1]
d2<-prob[2]
d3<-prob[3]
d4<-prob[4]
mu<-(2*(d2-d1)+d4-3*d3)/sqrt(4*d2+4*d1+3*d3+3*d4)
div<-sqrt((4*d1+4*d2+3*d3+3*d4)^3)
vt1<- -4*d1-12*d2-8*d4
vt2<-12*d1+4*d2+12*d3+4*d4
vt3<- -9*d1-15*d2-9/2*d3-21/2*d4
vt4<-7*d1+d2+15/2*d3+3/2*d4
vt<-c(vt1,vt2,vt3,vt4)/div
}
if(model=="recessive"){
q<-1-p
R<-RR*p^2+2*p*q+q^2
prob<-c(2*p^3*q,2*p^3*q*RR,3*p^2*q^2,p^2*q^2*RR)/R
r1<-prob[1]
r2<-prob[2]
r3<-prob[3]
r4<-prob[4]
mu<-(2*(r2-r1)+3*r4-r3)/(sqrt(4*(r1+r2)+3*(r3+r4)))
vt1r<- -(4*(r1 + 3*r2 + r3 + 3*r4))/(4*r1 + 4*r2 + 3*r3 + 3*r4)^(3/2)
vt2r<-(4*(3*r1 + r2 + 2*r3))/(4*r1 + 4*r2 + 3*r3 + 3*r4)^(3/2)
vt3r<- -(r1 + 7*r2 + (3*r3)/2 + (15*r4)/2)/(4*r1 + 4*r2 + 3*r3 + 3*r4)^(3/2)
vt4r<-(3*(10*r1 + 6*r2 + 7*r3 + 3*r4))/(2*(4*r1 + 4*r2 + 3*r3 + 3*r4)^(3/2))
vt<-c(vt1r,vt2r,vt3r,vt4r)
}
varmat<-prob%*%t((-prob))
diag(varmat)<-prob*(1-prob)
va<-sum((varmat%*%vt)*vt)
std<-sqrt(va)
}
}
if(!is.null(n)){
z<-qnorm(1-alpha/2)
power<-pnorm(-z,sd=std,mean=sqrt(n)*mu)+(1-pnorm(z,sd=std,mean=sqrt(n)*mu))
if(outa){
model<-"multiplicative"
test="aTDT"
}
out<-c(model=model,power=power,alpha=alpha,n=n,test=test,RR=RR,p=p,calc="power")
return(out)
}
if(!is.null(beta)){
size<-as.vector(((std*qnorm(beta)+qnorm(1-alpha/2))/mu)^2)
if(outa){
model<-"multiplicative"
test<-"aTDT"
}
out<-c(model=model,size=size,beta=beta,alpha=alpha,test=test,RR=RR,p=p,calc="size")
return(out)
}
}
if(!is.null(n)){
vars<-expand.grid(n,p,RR,alpha,model,test,stringsAsFactors=FALSE)
wa<-which(vars[,6]=="aTDT")
le<-length(wa)/length(model)
if(le>=1){
aw<-wa[-c(1:le)]
vars<-vars[-aw,]
}
out<-t(mapply(single,n=vars[,1],p=vars[,2],RR=vars[,3],alpha=vars[,4],model=vars[,5],test=vars[,6]))
out<-as.data.frame(out)
out$n<-as.numeric(as.character(out$n))
out$alpha<-as.numeric(as.character(out$alpha))
out$RR<-as.numeric(as.character(out$RR))
out$p<-as.numeric(as.character(out$p))
out$power<-as.numeric(as.character(out$power))
}
if(!is.null(beta)){
vars<-expand.grid(beta,p,RR,alpha,model,test,stringsAsFactors=FALSE)
wa<-which(vars[,6]=="aTDT")
le<-length(wa)/length(model)
if(le>=1){
aw<-wa[-c(1:le)]
vars<-vars[-aw,]}
out<-t(mapply(single,beta=vars[,1],p=vars[,2],RR=vars[,3],alpha=vars[,4],model=vars[,5],test=vars[,6]))
out<-as.data.frame(out)
out$size<-as.numeric(as.character(out$size))
out$beta<-as.numeric(as.character(out$beta))
out$alpha<-as.numeric(as.character(out$alpha))
out$RR<-as.numeric(as.character(out$RR))
out$p<-as.numeric(as.character(out$p))
}
class(out) <- "trio.power"
out
}
print.trio.power<-function(x,digits=4,...){
if(x$calc[1]=="power"){
power<-format.pval(as.numeric(as.character(x$power)),digits=digits)
out<-data.frame(Test=x$test,Model=x$model, MAF=x$p,alpha=x$alpha,RR=x$RR,Trios=x$n,power=power)
cat(" Trio studies power calculation", "\n\n")
print(format(out,digits=digits))
}
if(x$calc[1]=="size"){
size<-ceiling(as.numeric(as.character(x$size)))
out<-data.frame(Test=x$test,Model=x$model,MAF=x$p,alpha=x$alpha,RR=x$RR,beta=x$beta,Trios=size)
cat(" Trio studies sample size calculation", "\n\n")
print(format(out,digits=digits,trim=FALSE)
)
if(any(x$RR==1)){
warning("RR should differ from one to get a finite sample size.",call.=FALSE)}
}
}
vt1<-function(d1,d2,d3,d4){
h<-(1/3*d1-d2+d3-1/3*d4)/(2*(d1+d3))
be<-log(sqrt((d2+d4)/(3*(d1+d3))+h^2)-h)
summe<-((d1 + d2)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3)))/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1) + ((d3 + d4)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - (d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3)))/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - (d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1/3) - ((d1 + d2)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - (d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))^2)/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - (d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1)^2 - ((d3 + d4)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - (d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))^2)/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - (d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1/3)^2
ablb<- -(1/(3*(2*d1 + 2*d3)) + ((4*(d1/3 - d2 + d3 - d4/3)^2)/(2*d1 + 2*d3)^3 - ((2*d1)/9 - (2*d2)/3 + (2*d3)/3 - (2*d4)/9)/(2*d1 + 2*d3)^2 +
(3*(d2 + d4))/(3*d1 + 3*d3)^2)/(2*((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2)) -
(2*(d1/3 - d2 + d3 - d4/3))/(2*d1 + 2*d3)^2)/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))
abls<-((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1) - ((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1)^2 - (d1 + d2)*(1/3/(2*d1 + 2*d3) +
1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 -
1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) + 3*(d2 + d4)/(3*d1 + 3*d3)^2) -
2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1) - (d3 + d4)*(1/3/(2*d1 + 2*d3) + 1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 - 1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) +
3*(d2 + d4)/(3*d1 + 3*d3)^2) - 2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3) +
3*(d1 + d2)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))*(1/3/(2*d1 + 2*d3) + 1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 - 1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) +
3*(d2 + d4)/(3*d1 + 3*d3)^2) - 2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1)^2 +
3*(d3 + d4)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))*(1/3/(2*d1 + 2*d3) + 1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 - 1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) +
3*(d2 + d4)/(3*d1 + 3*d3)^2) - 2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3)^2 -
2*(d1 + d2)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2*(1/3/(2*d1 + 2*d3) + 1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 - 1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) +
3*(d2 + d4)/(3*d1 + 3*d3)^2) - 2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1)^3 -
2*(d3 + d4)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2*(1/3/(2*d1 + 2*d3) + 1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 - 1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) +
3*(d2 + d4)/(3*d1 + 3*d3)^2) - 2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3)^3
abl<-ablb*sqrt(summe)+be*1/(2*sqrt(summe))*abls
abl
}
vt2<-function(d1,d2,d3,d4){
h<-(1/3*d1-d2+d3-1/3*d4)/(2*(d1+d3))
be<-log(sqrt((d2+d4)/(3*(d1+d3))+h^2)-h)
summe<-((d1 + d2)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3)))/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1) + ((d3 + d4)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3)))/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1/3) - ((d1 + d2)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))^2)/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1)^2 - ((d3 + d4)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))^2)/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1/3)^2
ablb<- -((((2*d1)/3 - 2*d2 + 2*d3 - (2*d4)/3)/(2*d1 + 2*d3)^2 - 1/(3*d1 + 3*d3))/(2*((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)) - 1/(2*d1 + 2*d3))/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))
abls<-((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1) - ((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1)^2 - (1/2*(1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) -
1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3))*(d1 + d2)/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1) - (1/2*(1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) -
1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3))*(d3 + d4)/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3) + 3*(1/2*(1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) -
1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3))*(d1 + d2)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1)^2 + 3*(1/2*(1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) -
1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3))*(d3 + d4)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3)^2 - 2*(1/2*(1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) -
1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3))*(d1 + d2)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1)^3 - 2*(1/2*(1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) -
1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3))*(d3 + d4)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3)^3
abl<-ablb*sqrt(summe)+be*1/(2*sqrt(summe))*abls
abl
}
vt3<-function(d1,d2,d3,d4){
h<-(1/3*d1-d2+d3-1/3*d4)/(2*(d1+d3))
be<-log(sqrt((d2+d4)/(3*(d1+d3))+h^2)-h)
summe<-((d1 + d2)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3)))/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1) + ((d3 + d4)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3)))/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1/3) - ((d1 + d2)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - (d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))^2)/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - (d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1)^2 -
((d3 + d4)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))^2)/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1/3)^2
ablb<- -(1/(2*d1 + 2*d3) + ((4*(d1/3 - d2 + d3 - d4/3)^2)/(2*d1 + 2*d3)^3 - ((2*d1)/3 - 2*d2 + 2*d3 - (2*d4)/3)/(2*d1 + 2*d3)^2 +
(3*(d2 + d4))/(3*d1 + 3*d3)^2)/(2*((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2)) -
(2*(d1/3 - d2 + d3 - d4/3))/(2*d1 + 2*d3)^2)/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))
abls<-((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3) - ((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3)^2 - (d1 + d2)*(1/(2*d1 + 2*d3) +
1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 -
1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) + 3*(d2 + d4)/(3*d1 + 3*d3)^2) -
2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1) - (d3 + d4)*(1/(2*d1 + 2*d3) + 1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 - 1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) +
3*(d2 + d4)/(3*d1 + 3*d3)^2) - 2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3) -
2*(d1 + d2)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2*(1/(2*d1 + 2*d3) + 1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 - 1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) +
3*(d2 + d4)/(3*d1 + 3*d3)^2) - 2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1)^3 -
2*(d3 + d4)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2*(1/(2*d1 + 2*d3) + 1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 - 1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) +
3*(d2 + d4)/(3*d1 + 3*d3)^2) - 2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3)^3 +
3*(d1 + d2)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))*(1/(2*d1 + 2*d3) + 1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 - 1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) +
3*(d2 + d4)/(3*d1 + 3*d3)^2) - 2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1)^2 +
3*(d3 + d4)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))*(1/(2*d1 + 2*d3) + 1/2/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)*(4/(2*d1 + 2*d3)^3*(1/3*d1 - d2 + d3 - 1/3*d4)^2 - 1/(2*d1 + 2*d3)^2*(2/3*d1 - 2*d2 + 2*d3 - 2/3*d4) +
3*(d2 + d4)/(3*d1 + 3*d3)^2) - 2/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3)^2
abl<-ablb*sqrt(summe)+be*1/(2*sqrt(summe))*abls
abl
}
vt4<-function(d1,d2,d3,d4){
h<-(1/3*d1-d2+d3-1/3*d4)/(2*(d1+d3))
be<-log(sqrt((d2+d4)/(3*(d1+d3))+h^2)-h)
summe<-((d1 + d2)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3)))/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1) + ((d3 + d4)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3)))/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1/3) - ((d1 + d2)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))^2)/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1)^2 - ((d3 + d4)*(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - (d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))^2)/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - (d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3) + 1/3)^2
ablb<- -((((2*d1)/9 - (2*d2)/3 + (2*d3)/3 - (2*d4)/9)/(2*d1 + 2*d3)^2 - 1/(3*d1 + 3*d3))/(2*((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2)) - 1/(3*(2*d1 + 2*d3)))/(((d1/3 - d2 + d3 - d4/3)^2/(2*d1 + 2*d3)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
(d1/3 - d2 + d3 - d4/3)/(2*d1 + 2*d3))
abls<-((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 + (d2 + d4)/(3*d1 + 3*d3))^(1/2) -
1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3) - ((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3)^2 -
(1/2*(1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) - 1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/3/(2*d1 + 2*d3))*(d1 + d2)/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1) -
(1/2*(1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) - 1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/3/(2*d1 + 2*d3))*(d3 + d4)/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3) +
3*(1/2*(1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) - 1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/3/(2*d1 + 2*d3))*(d1 + d2)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1)^2 +
3*(1/2*(1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) - 1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/3/(2*d1 + 2*d3))*(d3 + d4)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3)^2 -
2*(1/2*(1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) - 1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/3/(2*d1 + 2*d3))*(d1 + d2)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1)^3 -
2*(1/2*(1/(2*d1 + 2*d3)^2*(2/9*d1 - 2/3*d2 + 2/3*d3 - 2/9*d4) - 1/(3*d1 + 3*d3))/(1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/3/(2*d1 + 2*d3))*(d3 + d4)*((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4))^2/((1/(2*d1 + 2*d3)^2*(1/3*d1 - d2 + d3 - 1/3*d4)^2 +
(d2 + d4)/(3*d1 + 3*d3))^(1/2) - 1/(2*d1 + 2*d3)*(1/3*d1 - d2 + d3 - 1/3*d4) + 1/3)^3
abl<-ablb*sqrt(summe)+be*1/(2*sqrt(summe))*abls
abl
}
vt1r<-function(r1,r2,r3,r4){
h<-(3*r1-r2+r3-3*r4)/(2*(r1+r3))
be<-log(sqrt(3*(r2+r4)/(r1+r3)+h^2)-h)
summe<- -((r1 + r2)*((3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) -
((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)))/(((3*r2 + 3*r4)/(r1 + r3) +
(3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2) - (3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) + 1)^2 -
((r3 + r4)*((3*(3*r1 - r2 + r3 - 3*r4))/(2*r1 + 2*r3) - 3*((3*r2 + 3*r4)/(r1 + r3) +
(3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2)))/(((3*r2 + 3*r4)/(r1 + r3) + (3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2) -
(3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) + 3)^2
ablb<- (3/(2*r1 + 2*r3) - (2*(3*r1 - r2 + r3 - 3*r4))/(2*r1 + 2*r3)^2 + ((3*r2 + 3*r4)/(r1 + r3)^2 -
(18*r1 - 6*r2 + 6*r3 - 18*r4)/(2*r1 + 2*r3)^2 + (4*(3*r1 - r2 + r3 - 3*r4)^2)/(2*r1 + 2*r3)^3)/(2*((3*r2 + 3*r4)/(r1 + r3) +
(3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2)))/((3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) -
((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))
abls<- -(1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) - ((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 1)^2 - (r1 + r2)/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 1)^2*(3/(2*r1 + 2*r3) - 2/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4) + 1/2/((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)*((3*r2 + 3*r4)/(r1 + r3)^2 - 1/(2*r1 + 2*r3)^2*(18*r1 - 6*r2 + 6*r3 - 18*r4) +
4/(2*r1 + 2*r3)^3*(3*r1 - r2 + r3 - 3*r4)^2)) - (r3 + r4)/(((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) - 1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 3)^2*(9/(2*r1 + 2*r3) -
6/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4) + 3/2/((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)*((3*r2 + 3*r4)/(r1 + r3)^2 - 1/(2*r1 + 2*r3)^2*(18*r1 - 6*r2 + 6*r3 - 18*r4) +
4/(2*r1 + 2*r3)^3*(3*r1 - r2 + r3 - 3*r4)^2)) - 2*(r1 + r2)*(1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) - ((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 1)^3*(3/(2*r1 + 2*r3) - 2/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4) +
1/2/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)*((3*r2 + 3*r4)/(r1 + r3)^2 -
1/(2*r1 + 2*r3)^2*(18*r1 - 6*r2 + 6*r3 - 18*r4) + 4/(2*r1 + 2*r3)^3*(3*r1 - r2 + r3 - 3*r4)^2)) -
2*(r3 + r4)*(3/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) - 3*((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 3)^3*(3/(2*r1 + 2*r3) - 2/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4) +
1/2/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)*((3*r2 + 3*r4)/(r1 + r3)^2 -
1/(2*r1 + 2*r3)^2*(18*r1 - 6*r2 + 6*r3 - 18*r4) + 4/(2*r1 + 2*r3)^3*(3*r1 - r2 + r3 - 3*r4)^2))
abl<-ablb*sqrt(summe)+be*1/(2*sqrt(summe))*abls
abl
}
vt2r<-function(r1,r2,r3,r4){
h<-(3*r1-r2+r3-3*r4)/(2*(r1+r3))
be<-log(sqrt(3*(r2+r4)/(r1+r3)+h^2)-h)
summe<- -((r1 + r2)*((3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) - ((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)))/(((3*r2 + 3*r4)/(r1 + r3) + (3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2) -
(3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) + 1)^2 - ((r3 + r4)*((3*(3*r1 - r2 + r3 - 3*r4))/(2*r1 + 2*r3) -
3*((3*r2 + 3*r4)/(r1 + r3) + (3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2)))/(((3*r2 + 3*r4)/(r1 + r3) +
(3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2) - (3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) + 3)^2
ablb<- -(1/(2*r1 + 2*r3) - ((3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3)^2 - 3/(2*(r1 + r3)))/((3*r2 + 3*r4)/(r1 + r3) +
(3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2))/((3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) -
((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))
abls<- -(1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) - ((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 1)^2 - (r1 + r2)*(1/2*(1/(2*r1 + 2*r3)^2*(6*r1 - 2*r2 + 2*r3 - 6*r4) -
3/(r1 + r3))/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 1)^2 - (r3 + r4)*(3/2*(1/(2*r1 + 2*r3)^2*(6*r1 - 2*r2 + 2*r3 - 6*r4) -
3/(r1 + r3))/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
3/(2*r1 + 2*r3))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 3)^2 - 2*(r1 + r2)*(1/2*(1/(2*r1 + 2*r3)^2*(6*r1 - 2*r2 + 2*r3 - 6*r4) -
3/(r1 + r3))/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3))*(1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) - ((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))/(((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) - 1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 1)^3 -
2*(r3 + r4)*(1/2*(1/(2*r1 + 2*r3)^2*(6*r1 - 2*r2 + 2*r3 - 6*r4) - 3/(r1 + r3))/((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) - 1/(2*r1 + 2*r3))*(3/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) -
3*((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))/(((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) - 1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 3)^3
abl<-ablb*sqrt(summe)+be*1/(2*sqrt(summe))*abls
abl
}
vt3r<-function(r1,r2,r3,r4){
h<-(3*r1-r2+r3-3*r4)/(2*(r1+r3))
be<-log(sqrt(3*(r2+r4)/(r1+r3)+h^2)-h)
summe<- -((r1 + r2)*((3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) - ((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)))/(((3*r2 + 3*r4)/(r1 + r3) + (3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2) -
(3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) + 1)^2 - ((r3 + r4)*((3*(3*r1 - r2 + r3 - 3*r4))/(2*r1 + 2*r3) - 3*((3*r2 + 3*r4)/(r1 + r3) +
(3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2)))/(((3*r2 + 3*r4)/(r1 + r3) + (3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2) -
(3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) + 3)^2
ablb <- ((((3*r2)/2 + (3*r4)/2)/(r1 + r3)^2 - (3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3)^2 +
(2*(3*r1 - r2 + r3 - 3*r4)^2)/(2*r1 + 2*r3)^3)/((3*r2 + 3*r4)/(r1 + r3) + (3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2) +
1/(2*r1 + 2*r3) - (6*r1 - 2*r2 + 2*r3 - 6*r4)/(2*r1 + 2*r3)^2)/((3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) -
((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))
abls<- -(3/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) - 3*((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 3)^2 - (r1 + r2)*(1/(2*r1 + 2*r3) - 2/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4) +
1/2/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)*((3*r2 + 3*r4)/(r1 + r3)^2 -
1/(2*r1 + 2*r3)^2*(6*r1 - 2*r2 + 2*r3 - 6*r4) + 4/(2*r1 + 2*r3)^3*(3*r1 - r2 + r3 - 3*r4)^2))/(((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) - 1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 1)^2 -
(r3 + r4)/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 3)^2*(3/(2*r1 + 2*r3) - 6/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4) + 3/2/((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)*((3*r2 + 3*r4)/(r1 + r3)^2 - 1/(2*r1 + 2*r3)^2*(6*r1 - 2*r2 + 2*r3 - 6*r4) +
4/(2*r1 + 2*r3)^3*(3*r1 - r2 + r3 - 3*r4)^2)) - 2*(r1 + r2)*(1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) -
((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))*(1/(2*r1 + 2*r3) -
2/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4) + 1/2/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)*((3*r2 + 3*r4)/(r1 + r3)^2 -
1/(2*r1 + 2*r3)^2*(6*r1 - 2*r2 + 2*r3 - 6*r4) + 4/(2*r1 + 2*r3)^3*(3*r1 - r2 + r3 - 3*r4)^2))/(((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) - 1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 1)^3 -
2*(r3 + r4)*(3/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) - 3*((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))*(1/(2*r1 + 2*r3) - 2/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4) +
1/2/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)*((3*r2 + 3*r4)/(r1 + r3)^2 -
1/(2*r1 + 2*r3)^2*(6*r1 - 2*r2 + 2*r3 - 6*r4) + 4/(2*r1 + 2*r3)^3*(3*r1 - r2 + r3 - 3*r4)^2))/(((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) - 1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 3)^3
abl<-ablb*sqrt(summe)+be*1/(2*sqrt(summe))*abls
abl
}
vt4r<-function(r1,r2,r3,r4){
h<-(3*r1-r2+r3-3*r4)/(2*(r1+r3))
be<-log(sqrt(3*(r2+r4)/(r1+r3)+h^2)-h)
summe<- -((r1 + r2)*((3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) -
((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2)))/(((3*r2 + 3*r4)/(r1 + r3) +
(3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2) - (3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) + 1)^2 -
((r3 + r4)*((3*(3*r1 - r2 + r3 - 3*r4))/(2*r1 + 2*r3) - 3*((3*r2 + 3*r4)/(r1 + r3) +
(3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2)))/(((3*r2 + 3*r4)/(r1 + r3) + (3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2) -
(3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) + 3)^2
ablb<-(((9*r1 - 3*r2 + 3*r3 - 9*r4)/(2*r1 + 2*r3)^2 - 3/(2*(r1 + r3)))/((3*r2 + 3*r4)/(r1 + r3) +
(3*r1 - r2 + r3 - 3*r4)^2/(2*r1 + 2*r3)^2)^(1/2) - 3/(2*r1 + 2*r3))/((3*r1 - r2 + r3 - 3*r4)/(2*r1 + 2*r3) -
((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))
abls<- -(3/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) - 3*((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 3)^2 - (r1 + r2)*(1/2*(1/(2*r1 + 2*r3)^2*(18*r1 - 6*r2 + 6*r3 - 18*r4) -
3/(r1 + r3))/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
3/(2*r1 + 2*r3))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 1)^2 - (r3 + r4)*(3/2*(1/(2*r1 + 2*r3)^2*(18*r1 - 6*r2 + 6*r3 - 18*r4) -
3/(r1 + r3))/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
9/(2*r1 + 2*r3))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 3)^2 - 2*(r1 + r2)*(1/2*(1/(2*r1 + 2*r3)^2*(18*r1 - 6*r2 + 6*r3 - 18*r4) -
3/(r1 + r3))/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
3/(2*r1 + 2*r3))*(1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) - ((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 1)^3 - 2*(r3 + r4)*(1/2*(1/(2*r1 + 2*r3)^2*(18*r1 - 6*r2 + 6*r3 - 18*r4) -
3/(r1 + r3))/((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
3/(2*r1 + 2*r3))*(3/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) - 3*((3*r2 + 3*r4)/(r1 + r3) +
1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2))/(((3*r2 + 3*r4)/(r1 + r3) + 1/(2*r1 + 2*r3)^2*(3*r1 - r2 + r3 - 3*r4)^2)^(1/2) -
1/(2*r1 + 2*r3)*(3*r1 - r2 + r3 - 3*r4) + 3)^3
abl<-ablb*sqrt(summe)+be*1/(2*sqrt(summe))*abls
abl
}
Try the trio package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.