# A.Azzalini (18.10.2008)
# costruiamo regione di ammissibilita` per CP nel caso d=2
# Nome originario: "check-invertCP-2.R"
#-----------------
gamma1 <- function(delta){
  mu.z <- sqrt(2/pi)*delta
  sigma.z <- sqrt(1-mu.z^2)
  0.5*(4-pi) * (mu.z/sigma.z)^3
  }
# 
cat("For a bivariate SN distribution, this program displays the maximum 
and minimum value of the correlation (rho) as a function of delta[1:2] 
and of gamma1[1:2]. Press <enter> to move from first to the second plot.
The plots can be rotated using the mouse. The library('rgl') is required.\n")
#
library("rgl")
npts <- 51
delta1 <- delta2 <- seq(-0.999, 0.999, length=npts)
max.rho <- min.rho <- matrix(NA, npts, npts)
# uses (2.10) of Azzalini & Dalla Valle (1996) & top formula of p.720
for(i in 1:npts){
  for(j in 1:npts){
    d1 <- delta1[i]
    d2 <- delta2[j]
    max.omega <- d1*d2 + sqrt((1-d1^2)*(1-d2^2))
    min.omega <- d1*d2 - sqrt((1-d1^2)*(1-d2^2))
    den <- sqrt((1-2*d1^2/pi)*(1-2*d2^2/pi))
    max.rho[i,j] <- (max.omega-2*d1*d2/pi)/den
    min.rho[i,j] <- (min.omega-2*d1*d2/pi)/den
    }
}

rgl.open()
# rgl.bbox(color=c("#333377","white"), emission="#333377", 
#           specular="#3333FF", shininess=5, alpha=0.8 )
axes3d(edges = "bbox", labels = TRUE, tick = TRUE, nticks = 5)
rgl.surface(delta1, delta2, max.rho, color="cyan")
rgl.surface(delta1, delta2, min.rho, color="orange")
title3d('min/max rho ','','delta1', 'rho', 'delta2', color="white")
#
readline("<enter>")
rgl.clear()
gamma1.1 <- gamma1(delta1)
gamma1.2 <- gamma1(delta2)
axes3d(edges = "bbox", labels = TRUE, tick = TRUE, nticks = 5)
rgl.surface(gamma1.1, gamma1.2, max.rho, color="cyan")
rgl.surface(gamma1.1, gamma1.2, min.rho, color="orange")
title3d('min/max rho ','','gamma1[1]', 'rho', 'gamma1[2]', color="white")
readline("<enter>")
rgl.close()
#
# persp(gamma1.1, gamma1.2, max.rho, phi=30, theta=40, ticktype="detailed")
# persp(gamma1.1, gamma1.2, min.rho, phi=30, theta=40, ticktype="detailed")
#