ps[mm,r] = sum(ps1[,mm,r])
nY2[mm,r] = sum(Y2[,mm,r]^2)
rho[mm,mm,r] = (ps[mm,r]+sqrt(ps[mm,r]^2+4*nY2[mm,r]*gam2[r])) / (2*nY2[mm,r])
ps[mm,r] = sum(ps1[,mm,r])
nY2[mm,r] = sum(Y2[,mm,r]^2)
rho[mm,mm,r] = (ps[mm,r]+sqrt(ps[mm,r]^2+4*nY2[mm,r]*gam2[r])) / (2*nY2[mm,r])