; $Id: polar_surface.pro,v 1.15 2004/01/21 15:54:58 scottm Exp $ ; ; Copyright (c) 1992-2004, Research Systems, Inc. All rights reserved. ; Unauthorized reproduction prohibited. function polar_surface, z, r, theta, GRID = grid, SPACING = sp, $ BOUNDS = bounds, QUINTIC = quintic, MISSING = missing ;+ ; NAME: ; POLAR_SURFACE ; ; PURPOSE: ; This function interpolates a surface from polar coordinates ; (R, Theta, Z) to rectangular coordinates (X, Y, Z). ; ; CATEGORY: ; Gridding. ; ; CALLING SEQUENCE: ; Result = POLAR_SURFACE(Z, R, Theta) ; ; INPUTS: ; Z: An array containing the surface value at each point. ; If the data are regularly gridded (GRID=1) in R and ; Theta, Z is a two dimensional array, where Z[i,j] has a ; radius of R[i] and an azimuth of Theta[j]. If GRID is ; not set, R[i] and Theta[i] contain the radius and azimuth ; of each Z[i]. ; R: The radius. If GRID is set, Z[i,j] has a radius of R[i]. ; If GRID is not set, R must have the same number of elements ; as Z, and contains the radius of each point. ; Theta: The azimuth, in radians. If GRID is set, Z[i,j] has an ; azimuth of Theta[j]. If GRID is not set, Theta must ; have the same number of elements as Z, and contains ; the azimuth of each point. ; ; KEYWORD PARAMETERS: ; GRID: Set GRID to indicate that Z is regularly gridded in R ; and Theta. ; SPACING: A two element vector containing the desired grid spacing ; of the resulting array in X and Y. If omitted, the grid ; will be approximately 51 by 51. ; BOUNDS: A four element vector, [X0, Y0, X1, Y1], containing the ; limits of the XY grid of the resulting array. If omitted, ; the extent of input data sets the limits of the grid. ; QUINTIC: If set, the function uses quintic interpolation, which is ; slower but smoother than the default linear interpolation. ; MISSING: A value to use for areas within the grid but not within ; the convex hull of the data points. The default is 0.0. ; ; OUTPUTS: ; This function returns a two-dimensional array of the same type as Z. ; ; PROCEDURE: ; First, each data point is transformed to (X, Y, Z). Then ; the TRIANGULATE and TRIGRID procedures are used to interpolate ; the surface over the rectangular grid. ; ; EXAMPLE: ; r = findgen(50) / 50. ;Radius ; theta = findgen(50) * (2 * !pi / 50.) ;Theta ; z = r # sin(theta) ;Make a function (tilted circle) ; SURFACE, POLAR_SURFACE(z, r, theta, /GRID) ;Show it ; ; MODIFICATION HISTORY: ; DMS Sept, 1992 Written ;- COMPILE_OPT strictarr zeroes = where(r eq 0, count) ;Avoid multiple points with a 0 radius ;Fudge the center point?? if (keyword_set(grid) and count ne 0) OR (count gt 1) then begin range = max(abs(r)) rtemp = r rtemp[zeroes] = range / 1.0e5 ;Fudge a little endif else rtemp = r IF keyword_set(grid) THEN BEGIN ;Regulary gridded? s = size(z) if s[0] ne 2 then message, "Z must be 2d if GRID is set" if n_elements(r) ne s[1] then $ message, "R has wrong number of elements" if n_elements(theta) ne s[2] then $ message, "Theta has wrong number of elements" x = rtemp # cos(theta) y = rtemp # sin(theta) ENDIF ELSE BEGIN if (n_elements(r) ne n_elements(z)) or $ (n_elements(theta) ne n_elements(z)) then $ message, "R and Theta must have as many elements as Z" x = rtemp * cos(theta) y = rtemp * sin(theta) ENDELSE TRIANGULATE, x, y, tr ;Triangulate it if n_elements(bounds) lt 4 then $ bounds = [ min(x, max=x1), min(y, max=y1), x1, y1] if n_elements(sp) lt 2 then $ sp = [ bounds[2] - bounds[0], bounds[3] - bounds[1]] / 50. if n_elements(missing) le 0 then missing = 0 RETURN, TRIGRID(x, y, z, tr, sp, bounds, QUINTIC = KEYWORD_SET(quintic), $ MISSING = missing) END