#' @title Distribution function of a non-central squared copula
#'
#' @description This function computes the distribution function a non-central squared copula
#'
#' @param u         (nx2) data matrix of pseudo-observations.
#' @param family    'Gaussian' , 't' , 'Clayton' , 'Frank' , 'Gumbel'.
#' @param param   c(a1,a2,tau)  where  a1,a2  are the non-negative non-centrality
#                     parameters, tau is Kendall's tau of the copula family.
#' @param dof     degrees of freedom of the Student copula (if needed).
#' @return \item{cdf}{Non-central squared copula evaluated at points u.}
#' @author Bouchra R. Nasri, August 14, 2019

#'
#'
#'
#' @examples param = c(0.8,2.5,0.7);
#' u = matrix(c(0.2,0.6,0.3,0.5,0.7,0.9),ncol=2,byrow=TRUE);
#' cdf=NCSCopCdf(u,'Clayton',param);
#'
#' @export

NCSCopCdf <- function(u,family, param,dof=NULL)
{

a1  = param[1];
a2  = param[2];
KendallTau= param[3];

switch(family,
"Gaussian" = {
theta =  copula::iTau(copula::normalCopula(),  tau = KendallTau)
},

"t" = {
theta =  copula::iTau(copula::normalCopula(),  tau = KendallTau)
},

"Clayton" = {
theta =  copula::iTau(copula::claytonCopula(),  tau = KendallTau)
},

"Frank" = {
theta =  copula::iTau(copula::frankCopula(),  tau = KendallTau)
},

"Gumbel" = {
theta =  copula::iTau(copula::gumbelCopula(),  tau = KendallTau)
}
)

d1 = a1^2;
d2 = a2^2;

U1 = u[,1];
U2 = u[,2];

G1 = sqrt(qchisq(U1,1,d1));  #Ga inverse
G2 = sqrt(qchisq(U2,1,d2));

ha = matrix(c((G1-a1),(-G1-a1) , (G2-a2), (-G2-a2)), ncol=4, byrow=FALSE);

htilde = pnorm(ha)+ 1E-10;
w0     = dnorm(ha);

U11 = matrix(c(htilde[,1], htilde[,3]), ncol=2,byrow=FALSE);
U12 = matrix(c(htilde[,1], htilde[,4]), ncol=2,byrow=FALSE);
U21 = matrix(c(htilde[,2], htilde[,3]), ncol=2,byrow=FALSE);
U22 = matrix(c(htilde[,2], htilde[,4]), ncol=2,byrow=FALSE);

switch(family,
"Gaussian" =
{
F11 = copula::pCopula(U11,copula::normalCopula(theta, dim = 2))+ 1E-20;
F12 = copula::pCopula(U12,copula::normalCopula(theta, dim = 2))+ 1E-20;
F21 = copula::pCopula(U21,copula::normalCopula(theta, dim = 2))+ 1E-20;
F22 = copula::pCopula(U22,copula::normalCopula(theta, dim = 2))+ 1E-20;

},

"t" = {
F11 = copula::pCopula(U11, copula::tCopula(theta, dim = 2,df = dof))+ 1E-20;
F12 = copula::pCopula(U12, copula::tCopula(theta, dim = 2,df = dof))+ 1E-20;
F21 = copula::pCopula(U21, copula::tCopula(theta, dim = 2,df = dof))+ 1E-20;
F22 = copula::pCopula(U22, copula::tCopula(theta, dim = 2,df = dof))+ 1E-20;
},

"Clayton" = {
F11 = copula::pCopula(U11,copula::claytonCopula(theta, dim = 2))+ 1E-20;
F12 = copula::pCopula(U12,copula::claytonCopula(theta, dim = 2))+ 1E-20;
F21 = copula::pCopula(U21,copula::claytonCopula(theta, dim = 2))+ 1E-20;
F22 = copula::pCopula(U22,copula::claytonCopula(theta, dim = 2))+ 1E-20;
},

"Frank" = {
F11 = copula::pCopula(U11,copula::frankCopula(theta, dim = 2))+ 1E-20;
F12 = copula::pCopula(U12,copula::frankCopula(theta, dim = 2))+ 1E-20;
F21 = copula::pCopula(U21,copula::frankCopula(theta, dim = 2))+ 1E-20;
F22 = copula::pCopula(U22,copula::frankCopula(theta, dim = 2))+ 1E-20;
},

"Gumbel" = {
F11 = copula::pCopula(U11,copula::gumbelCopula(theta, dim = 2))+ 1E-20;
F12 = copula::pCopula(U12,copula::gumbelCopula(theta, dim = 2))+ 1E-20;
F21 = copula::pCopula(U21,copula::gumbelCopula(theta, dim = 2))+ 1E-20;
F22 = copula::pCopula(U22,copula::gumbelCopula(theta, dim = 2))+ 1E-20;
}
)

cdf = F11-F12-F21+F22;

return(cdf)
}