"d2norm" <-
function (x,y,m1=0,m2=0,sd1=1,sd2=1,rho=0) 
{
# Dichte der zweidimensionalen Normalverteilung
f<-(1/(2*pi*sd1*sd2*sqrt(1-rho^2)))*exp((-1/(2*(1-rho^2)))*((x-m1)^2/sd1^2
   -2*rho*(x-m1)*(y-m2)/(sd1*sd2)+(y-m2)^2/sd2^2))
f
}