Bivariate Normal cumulative distribution function

function BiVariateNCDF(a , b , rho )
{
 aarray=new Array;
 barray=new Array;
    rho1=0.0;
 rho2=0.0;
 delta=0.0;
    ap=0.0;
 bp=0.0;
    Sum=0;
    i=0;
 j=0;
 Pi=0;
    denum=0.0;
 aarray[0]=0.24840615;
 aarray[1]=0.39233107;
 aarray[2]=0.21141819;
 aarray[3]=0.03324666;
 aarray[4]=0.00082485334;
 barray[0]=0.10024215;
 barray[1]=0.48281397;
 barray[2]=1.0609498;
 barray[3]=1.7797294;
 barray[4]=2.6697604;
 Pi = 3.14159265358979;
    
    ap = a / Math.sqrt(2 * (1 - rho * rho));
    bp = b / Math.sqrt(2 * (1 - rho * rho ));
 BivN=0.0;
 if (a <= 0 && b <= 0 && rho <= 0)
 {
  Sum = 0;
        for (i = 0; i< 5;i++)
   for (j = 0; j< 5;j++)
                Sum = Sum + aarray[i] * aarray[j] * f(barray[i], barray[j], ap, bp, rho);
        BivN = Sum * Math.sqrt(1 - rho * rho) / Pi;
 }
 else if( a <= 0 && b >= 0 && rho >= 0)
  BivN = NCDF2(a) - BiVariateNCDF(a, -b, -rho);
 else if(a >= 0 && b <= 0 && rho >= 0 )
 {
  BivN = NCDF2(b) - BiVariateNCDF(-a, b, -rho);
 }
 else if(a >= 0 && b >= 0 && rho <= 0 )
  BivN = NCDF2(a) + NCDF2(b) - 1 + BiVariateNCDF(-a, -b, rho);
 else if(a * b * rho >= 0)
 {
  denum = Math.sqrt(Math.abs(a*a-2*rho*a*b+b*b));
        rho1 = (rho * a - b) * Sgn(a) / denum;
        rho2 = (rho * b - a) * Sgn(b) / denum;
        delta = (1 - Sgn(a) * Sgn(b)) / 4;
  
        BivN = BiVariateNCDF(a, 0, rho1) + BiVariateNCDF(b, 0, rho2) - delta;
 }
 
 return BivN;
}


function CumulativeBivariateNormal(a, b, rho)
 {
    aarray=new Array;
 aarray[0]= 0.24840615;
 aarray[1]=0.39233107;
 aarray[2]=0.21141819;
 aarray[3]=0.03324666;
 aarray[4]=0.00082485334;
    barray=new Array;
 barray[0]=0.10024215;
 barray[1]=0.48281397;
 barray[2]=1.0609498;
 barray[3]=1.7797294;
 barray[4]=2.6697604;
  //   static Real aarray[4]= {0.325303, 0.4211071, 0.1334425, 0.006374323};
  //static Real barray[4]={0.1337764, 0.6243247, 1.3425378, 2.2626645};
    rho1=0;  rho2=0;  delta =0;
 ap=0;  bp =0;
    Sum=0; Pi=0;
    sgnA=0;sgnB=0;
    i=0; j=0 ;
    BivN=0; denum=0;

    sgnA=Sgn(a);
    sgnB=Sgn(b);
    
    Pi = 3.14159265358979;
    if(Math.pow(rho,2)==1)
        rho=Sgn(rho)*0.99999;
    ap = a / Math.sqrt(2 * (1 - Math.pow(rho,2)));
    bp = b / Math.sqrt(2 * (1 - Math.pow(rho,2)));
     
    if( (a <= 0) && (b <= 0) && (rho <= 0) ){
        Sum = 0;
        for (i = 0 ;i< 5;i++)
           for (j = 0 ;j< 5;j++)
                Sum +=  + aarray[i] * aarray[j] * SubFunctionForBivariateNormal(barray[i], barray[j], ap, bp, rho);
        BivN = Sum * Math.sqrt(1 - Math.pow(rho,2)) / Pi;
    }  
    else if ((a <= 0 )&& (b >= 0) && (rho >= 0)) 
        BivN = NCDF2(a) - CumulativeBivariateNormal(a, -b, -rho); 
    else if ((a >= 0) && (b <= 0) && (rho >= 0)) 
        BivN = NCDF2(b) - CumulativeBivariateNormal(-a, b, -rho);
    else if ((a >= 0) && (b >= 0) && (rho <= 0))
        BivN = NCDF2(a) + NCDF2(b) - 1 + CumulativeBivariateNormal(-a, -b, rho);
    else if (a * b * rho >= 0) {
        denum = Math.sqrt(Math.pow(a,2) - 2 * rho * a * b + Math.pow(b,2));
        rho1 = (rho * a - b) * sgnA / denum;
        rho2 = (rho * b - a) * sgnB / denum;
        delta = (1 - sgnA*sgnB) / 4;
        BivN = CumulativeBivariateNormal(a, 0, rho1) + CumulativeBivariateNormal(b, 0, rho2) - delta;
    }       
    return BivN;
}
 
function SubFunctionForBivariateNormal( X, y, ap, bp, rho){
    r = ap * (2 * X - ap) + bp * (2 * y - bp) + 2 * rho * (X - ap) * (y - bp);
   return Math.exp(r);
}