wilcox_R: What the package does (short line)

relplot <-
function(xv,yv,C=36,epsilon=.0001,plotit=TRUE,pch='*'){
#   Compute the measures of location, scale and correlation used in the
#   bivariate boxplot of Goldberg and Iglewicz,
#   Technometrics, 1992, 34, 307-320.
#
#   The code in relplot plots the boxplot.
#
#   This code assumes the data are in xv and yv
#
#   This code uses the function biloc, stored in the file biloc.b7 and
#   bivar stored in bivar.b7
#
plotit<-as.logical(plotit)
#
# Do pairwise elimination of missing values
#
temp<-matrix(c(xv,yv),ncol=2)
temp<-elimna(temp)
xv<-temp[,1]
yv<-temp[,2]
tx<-biloc(xv)
ty<-biloc(yv)
sx<-sqrt(bivar(xv))
sy<-sqrt(bivar(yv))
z1<-(xv-tx)/sx+(yv-ty)/sy
z2<-(xv-tx)/sx-(yv-ty)/sy
ee<-((z1-biloc(z1))/sqrt(bivar(z1)))^2+
((z2-biloc(z2))/sqrt(bivar(z2)))^2
w<-(1-ee/C)^2
if(length(w[w==0])>=length(xv)/2)warning("More than half of the w values equal zero")
sumw<-sum(w[ee<C])
tempx<-w*xv
txb<-sum(tempx[ee<C])/sumw
tempy<-w*yv
tyb<-sum(tempy[ee<C])/sumw
tempxy<-w*(xv-txb)*(yv-tyb)
tempx<-w*(xv-txb)^2
tempy<-w*(yv-tyb)^2
sxb<-sum((tempx[ee<C]))/sumw
syb<-sum((tempy[ee<C]))/sumw
rb<-sum(tempxy[ee<C])/(sqrt(sxb*syb)*sumw)
z1<-((xv-txb)/sqrt(sxb)+(yv-tyb)/sqrt(syb))/sqrt(2*(1+rb))
z2<-((xv-txb)/sqrt(sxb)-(yv-tyb)/sqrt(syb))/sqrt(2*(1-rb))
wo<-w
ee<-z1^2+z2^2
w<-(1-ee/C)^2
sumw<-sum(w[ee<C])
tempx<-w*xv
txb<-sum(tempx[ee<C])/sumw
tempy<-w*yv
tyb<-sum(tempy[ee<C])/sumw
tempxy<-w*(xv-txb)*(yv-tyb)
tempx<-w*(xv-txb)^2
tempy<-w*(yv-tyb)^2
sxb<-sum((tempx[ee<C]))/sumw
syb<-sum((tempy[ee<C]))/sumw
rb<-sum(tempxy[ee<C])/(sqrt(sxb*syb)*sumw)
z1<-((xv-txb)/sqrt(sxb)+(yv-tyb)/sqrt(syb))/sqrt(2*(1+rb))
z2<-((xv-txb)/sqrt(sxb)-(yv-tyb)/sqrt(syb))/sqrt(2*(1-rb))
iter<-0
while(iter<=10){
iter<=iter+1
ee<-z1^2+z2^2
w<-(1-ee/C)^2
sumw<-sum(w[ee<C])
tempx<-w*xv
txb<-sum(tempx[ee<C])/sumw
tempy<-w*yv
tyb<-sum(tempy[ee<C])/sumw
tempxy<-w*(xv-txb)*(yv-tyb)
tempx<-w*(xv-txb)^2
tempy<-w*(yv-tyb)^2
sxb<-sum((tempx[ee<C]))/sumw
syb<-sum((tempy[ee<C]))/sumw
rb<-sum(tempxy[ee<C])/(sqrt(sxb*syb)*sumw)
z1<-((xv-txb)/sqrt(sxb)+(yv-tyb)/sqrt(syb))/sqrt(2*(1+rb))
z2<-((xv-txb)/sqrt(sxb)-(yv-tyb)/sqrt(syb))/sqrt(2*(1-rb))
wo<-w
ee<-z1^2+z2^2
w<-(1-ee/C)^2
dif<-w-wo
crit<-sum(dif^2)/(mean(w))^2
if(crit <=epsilon)break
}
if(plotit){
em<-median(sqrt(ee))
r1<-em*sqrt((1+rb)/2)
r2<-em*sqrt((1-rb)/2)
temp<-c(0:179)
thet<-2*3.141593*temp/180
theta1<-r1*cos(thet)
theta2<-r2*sin(thet)
xplot1<-txb+(theta1+theta2)*sqrt(sxb)
yplot1<-tyb+(theta1-theta2)*sqrt(syb)
emax<-max(sqrt(ee[ee<7*em^2]))
r1<-emax*sqrt((1+rb)/2)
r2<-emax*sqrt((1-rb)/2)
theta1<-r1*cos(thet)
theta2<-r2*sin(thet)
xplot<-txb+(theta1+theta2)*sqrt(sxb)
yplot<-tyb+(theta1-theta2)*sqrt(syb)
totx<-c(xv,xplot,xplot1)
toty<-c(yv,yplot,yplot1)
plot(totx,toty,type="n",xlab="x",ylab="y")
points(xv,yv,pch=pch)
points(xplot,yplot,pch=".")
points(xplot1,yplot1,pch=".")
}
list(mest=c(txb,tyb),mvar=c(sxb,syb),mrho=rb)
}