Dear yt users, I am trying to use the pf.h.surface routine, for example on density: *surface=pf.h.surface(primary,"Density",external_density)* and then using the sum of the triangles area (via surface.triangles) as an estimate of the surface area. But this method is giving me a wrong value in a test dataset, where I already know the right value for the surface area. I would like hence to ask if there is a better way to estimate the surface area of an iso-quantity surface, or if the process I am following is wrong. Below are the commands I am using to estimate the surface: * #computes the surface of the equidensity surface at the desired density* * #obtains the vertices of the triangles composing the surface tassellation approximation* * surface=pf.h.surface(primary,"Density",density)* * surface_density = surface['Density']* * surface_trg_vertices = surface.triangles* * #rearranges the arrays with the vertex of the triangles (coordinates names: a, b, c) to simplify the vectorial computation* * a_x = surface_trg_vertices[:,0][:,0]* * a_y = surface_trg_vertices[:,0][:,1]* * a_z = surface_trg_vertices[:,0][:,2]* * b_x = surface_trg_vertices[:,1][:,0]* * b_y = surface_trg_vertices[:,1][:,1]* * b_z = surface_trg_vertices[:,1][:,2]* * c_x = surface_trg_vertices[:,2][:,0]* * c_y = surface_trg_vertices[:,2][:,1]* * c_z = surface_trg_vertices[:,2][:,2]* * #finds the sides lengths of the triangles (an array for each set of side)* * side_ab = ((a_x - b_x)**2.0 + (a_y - b_y)**2.0 + (a_z - b_z)**2.0)**0.5* * side_bc = ((b_x - c_x)**2.0 + (b_y - c_y)**2.0 + (b_z - c_z)**2.0)**0.5* * side_ca = ((c_x - a_x)**2.0 + (c_y - a_y)**2.0 + (c_z - a_z)**2.0)**0.5* * #computes the altitude of the triangles wrt the ab side* * s = 0.5 * (side_ab + side_bc + side_ca) #semiperimeter* * altitude_side_ab = (2.0 * (s * (s - side_ab) * (s - side_bc) * (s - side_ca))**0.5) / side_ab* * #computes the area of the triangles using ab sides and altitudes * * trg_surface_area = (side_ab * altitude_side_ab) / 2.0* * #computes the total surface area* * total_surface_area = np.sum(trg_surface_area)* Thanks, Roberto