intfun<-function(a,b,rho,sqnp){ ## # Area under the bivariate normal calculated # # analogously to the midpoint method. # # Finds P[X1