[yt-users] Computing the surface area of an iso-quantity surface

[trobolo dinni]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.spacepope.org/pipermail/yt-users-spacepope.org/attachments/20140203/4b245bf3/attachment.htm>