;+
;
; @file_comments
; cut p parallelogram(s) into p*n^2 parallelograms
;
; @categories basic work
;
; @param x0 {in}{required}
; @param y0 {in}{required}
; @param x1 {in}{required}
; @param y1 {in}{required}
; @param x2 {in}{required}
; @param y2 {in}{required}
; @param x3 {in}{required}
; @param y3 {in}{required}
; 1d arrays of p elements, giving the edge positions.
; The edges must be given as in plot to draw the parallelogram. (see example).
;
; @param n {in}{required}
; each parallelogram will be cut in n^2 pieces
;
; @keyword ENDPOINTS
; see outputs
;
; @keyword ONSPHERE
; to specify that the points are located on a
; sphere. In this case, x and y corresponds to longitude and
; latitude in degrees.
;
; @returns
;  - default: a 3d array(2,n^2,p) giving the center position of each
;  piece of the parallelograms
;  - if /ENDPOINTS : a 3d array(2,(n+1)^2,p) giving the edge positions
;  of each piece of the parallelograms
;
; @uses cutsegment.pro
;
; @examples
;
; IDL> x0 = [2,6,2]
; IDL> y0 = [0,2,6]
; IDL> x1 = [3,8,4]
; IDL> y1 = [4,4,6]
; IDL> x2 = [1,6,4]
; IDL> y2 = [5,6,8]
; IDL> x3 = [0,4,2]
; IDL> y3 = [1,4,8]
; IDL> n = 4
; IDL> splot, [0,10], [0,10], xstyle = 1, ystyle = 1,/nodata
; IDL> 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]]
; IDL> res=cutpar(x0, y0, x1, y1, x2, y2, x3, y3, n)
; IDL> for i=0,2 do oplot, [res[0,*,i]], [res[1,*,i]], color = 20+10*i, psym = 1, thick = 3
;
; @history
;           S. Masson (smasson\@lodyc.jussieu.fr)
;           July 5th, 2002
;
; @version $Id$
;
;-
FUNCTION cutpar, x0, y0, x1, y1, x2, y2, x3, y3, n, ENDPOINTS = endpoints, ONSPHERE = onsphere
;
  compile_opt idl2, strictarrsubs
;
; 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