- //ps[mm,r] = sum(ps1[,mm,r])
- Real sumPs1 = 0.;
- for (int u=0; u<n; u++)
- sumPs1 += ps1[ai(u,mm,r,n,m,k)];
- ps[mi(mm,r,m,k)] = sumPs1;
- //nY2[mm,r] = sum(Y2[,mm,r]^2)
- Real sumY2 = 0.;
- for (int u=0; u<n; u++)
- sumY2 += Y2[ai(u,mm,r,n,m,k)] * Y2[ai(u,mm,r,n,m,k)];
- nY2[mi(mm,r,m,k)] = sumY2;
- //rho[mm,mm,r] = (ps[mm,r]+sqrt(ps[mm,r]^2+4*nY2[mm,r]*(gam2[r]))) / (2*nY2[mm,r])
- rho[ai(mm,mm,r,m,m,k)] = ( ps[mi(mm,r,m,k)] + sqrt( ps[mi(mm,r,m,k)]*ps[mi(mm,r,m,k)]
- + 4*nY2[mi(mm,r,m,k)] * gam2[r] ) ) / (2*nY2[mi(mm,r,m,k)]);
+ //rho[mm,mm,r] = (ps+sqrt(ps^2+4*nY2*gam2[r])) / (2*nY2)
+ rho[ai(mm,mm,r,m,m,k)] = (ps + sqrt(ps*ps + 4*nY2 * gam2[r])) / (2*nY2);