; $Id: int_tabulated_2d.pro,v 1.7 2004/01/21 15:54:55 scottm Exp $ ; ; Copyright (c) 1998-2004, Research Systems, Inc. All rights reserved. ; Unauthorized reproduction prohibited. ;+ ; NAME: ; INT_TABULATED_2D ; ; PURPOSE: ; Integrate the volume formed by a surface defined by tabulated ; data, (x[i], y[i], f[i]), and the XY plane, over the area ; defined by the convex hull of the set of points (x,y). ; ; CATEGORY: ; Numerical Analysis. ; ; CALLING SEQUENCE: ; Result = INT_TABULATED_2D(X, Y, F) ; ; INPUTS: ; X: The array of tabulated X-value data. This data may be ; irregularly gridded and in random order. ; Y: The tabulated Y-value data, Y[i] corresponds to X[i]. ; F: The tabulated F-value data. F[i] is the function's value ; at the point (X[i],Y[i]).; ; X Y, and F must be of floating point or double precision type. ; ; OUTPUTS: ; This fuction returns the integral of F computed from the tabulated ; data over the convex hull defined by the points in (X,Y). ; Linear interpolation is used. ; ; PROCEDURE: ; The (X,Y) points are triangulated using TRIANGULATE. Then for ; each triangle, in the convex hull, formed by points X[[i,j,k]] ; and Y[[i,j,k]], the volume of the triangular cylinder formed ; by the 6 points: ; (X[i], Y[i], Z[i]), (X[j], Y[j], Z[j]), (X[i], Y[k], Z[k]), ; (X[i], Y[i], 0), (X[j], Y[j], 0), and (X[k], Y[k], 0), is ; computed and summed. ; ; EXAMPLE: ; Compute the volume between the surface f=x^2+y^2 and the XY ; plane, over the interval of x=-1 to +1, y=-1 to 1, with a grid ; increment of 0.1: ; n = 21 ; x = (findgen(n) * (2./(n-1)) - 1.0) # replicate(1.0, n) ; y = transpose(x) ; f = x^2 + y^2 ; print, int_tabulated_2d(x,y,f) ; 2.68000 ; (The correct answer computed symbolically is 8/3 = 2.66667) ; ; MODIFICATION HISTORY: ; DMS, August, 1998. ;- Function TetraVolume, x, y, z ; Return the signed volume of a tetrahedron with 4 vertices contained in the ; x, y, and z arrays. If the vertices are ordered counter-clockwise, ; looking from the side opposite vertex [0], then the volume is ; positive, otherwise its negative. xx = x[1:3]-x[0] ;Translate one point to the origin yy = y[1:3]-y[0] zz = z[1:3]-z[0] s = (xx[0] * (yy[1]*zz[2] - yy[2]*zz[1]) + $ ;Compute the 3x3 determinant yy[0] * (zz[1]*xx[2] - zz[2]*xx[1]) + $ zz[0] * (xx[1]*yy[2] - xx[2]*yy[1])) return, s/6. end ; Old method, now coded in-line. ; Function TriCyl_volume, x, y, z ; ; Determine the volume of a triangular cylinder... Base is the ; ; triangle with vertices index [0,1,2], all with Z value = 0. Top is ; ; given by the vertices [3,4,5]. X and Y values for vertices 0 and 3, ; ; 1 and 4, 2 and 5 are identical. ; s = 0.0 ; ; Composed of 3 tetrahedra with these vertices. ; v = [[0,1,2,4],[0,2,5,4],[0,4,5,3]] ; for i=0,2 do begin ; t = v[*,i] ; s = s + TetraVolume(x[t],y[t],z[t]) ; endfor ; return, s ; end Function INT_TABULATED_2D, x, y, z triangulate, x, y, tr ; Subscript indices for the 3 tetrahedra formed with vertex 0 as one ; point. The faces must be order CCW when looking from the side away ; from vertex 0. s0 = [0,1,2,1] s1 = [0,2,2,1] s2 = [0,1,2,0] vol = 0.0 for i=0, n_elements(tr)/3-1 do begin ;Sum each triangle t = tr[*,i] x1 = x[t] ;Save vertices of this triangle y1 = y[t] z1 = z[t] ; Sum the three tetrahedra that make this triangular cylinder. p = TetraVolume(x1[s0], y1[s0], [0,0,0,z1[1]]) + $ TetraVolume(x1[s1], y1[s1], [0,0,z1[2],z1[1]]) + $ TetraVolume(x1[s2], y1[s2], [0,z1[1],z1[2],z1[0]]) vol = vol + p endfor return, vol end