source: trunk/Interpolation/cutpar.pro @ 59

;+
; NAME: cutpar
;
; PURPOSE: cut p parallelogram(s) into p*n^2 parallelograms
;
; CATEGORY: basic work
;
; CALLING SEQUENCE:res = cutpar(x0, y0, x1, y1, x2, y2, x3, y3, n)
;
; INPUTS:
;       x0,y0 1d arrays of p elements, giving the edge positions. The
;       edges must be given as in plot to traw the parallelogram. (see
;       example).
;       n: each parallelogram will be cutted in n^2 pieces
;
; KEYWORD PARAMETERS: 
;
;         /endpoints: see outputs
;
;         /onsphere: to specify that the points are located on a
;         sphere. In this case, x and y corresponds to longitude and
;         latitude in degrees.
;
; OUTPUTS:
;        -defaut: 3d array(2,n^2,p) giving the center position of each
;        piece of the parallelograms
;        -/endpoints: 3d array(2,(n+1)^2,p) giving the edge positions
;        of each piece of the parallelograms
;
; COMMON BLOCKS: no
;
; SIDE EFFECTS: need cutsegment.pro
;
; RESTRICTIONS: ?
;
; EXAMPLE:
;
; x0 = [2,6,2]
; y0 = [0,2,6]
; x1 = [3,8,4]
; y1 = [4,4,6]
; x2 = [1,6,4]
; y2 = [5,6,8]
; x3 = [0,4,2]
; y3 = [1,4,8]
; n = 4
; splot, [0,10], [0,10], xstyle = 1, ystyle = 1,/nodata
; for i=0,2 do oplot, [x0[i],x1[i],x2[i],x3[i],x0[i]],[y0[i],y1[i],y2[i],y3[i],y0[i]]
; res=cutpar(x0, y0, x1, y1, x2, y2, x3, y3, n)
; for i=0,2 do oplot, [res[0,*,i]], [res[1,*,i]], color = 20+10*i, psym = 1, thick = 3
;
; MODIFICATION HISTORY:
;           S. Masson (smasson@lodyc.jussieu.fr)
;           July 5th, 2002
;-
FUNCTION cutpar, x0, y0, x1, y1, x2, y2, x3, y3, n, endpoints = endpoints, onsphere = onsphere
; is it a parallelogram?
; eps = 1e-4
; IF total(abs((x0+x2)/2-(x1+x3)/2) GE eps) GT 0 $
;   OR total(abs((y0+y2)/2-(y1+y3)/2) GE eps) GT 0 $
;   THEN stop; print, 'NOT a parallelogram'
; x0(npar)
  npar = n_elements(x0) 
; firstborder(2,n+keyword_set(endpoints),npar)
  firstborder = cutsegment(x0, y0, x1, y1, n $
                           , endpoints = endpoints, onsphere = onsphere)
  thirdborder = cutsegment(x3, y3, x2, y2, n $
                           , endpoints = endpoints, onsphere = onsphere)
; res(2,n+keyword_set(endpoints),(n+keyword_set(endpoints))*npar)
  res = cutsegment(firstborder[0, *, *], firstborder[1, *, *] $
                   , thirdborder[0, *, *], thirdborder[1, *, *] $
                   , n, endpoints = endpoints, onsphere = onsphere)
; free memory
  firstborder = -1
  thirdborder = -1
; reform the result
  res = reform(res, 2, (n+keyword_set(endpoints))^2, npar, /overwrite)

  RETURN, res
END