source: trunk/SRC/Interpolation/cutsegment.pro @ 330

;+
;
; @file_comments
; cut p segments into p*n equal parts
;
; @categories
; basic work
;
; @param x0 {in}{required}
; @param y0 {in}{required}
; @param x1 {in}{required}
; @param y1 {in}{required}
; 1d arrays of p elements, the coordinates of the endpoints of the p segments
;
; @param n {in}{required}
; the number of pieces we want to cut each segment
;
; @keyword ENDPOINTS
; see returns
;
; @keyword ONSPHERE
; to specify that the points are located on a sphere.
; In this case, x and y correspond to longitude and latitude in degrees.
;
; @returns
;  - default: a 3d array (2,n,p) that gives the coordinates of the
;  middle of the cut segments.
;  - if /ENDPOINTS, a 3d array (2,n+1,p) that gives the
;  coordinates of the endpoints of the cutted segments.
;
; @examples
; IDL> x0=[2,5]
; IDL> y0=[5,1]
; IDL> x1=[9,3]
; IDL> y1=[1,8]
; IDL> res=cutsegment(x0, y0, x1, y1, 10)
; IDL> splot, [0,10], [0,10], xstyle = 1, ystyle = 1,/nodata
; IDL> oplot, [x0[0], x1[0]], [y0[0], y1[0]]
; IDL> oplot, [res[0,*,0]], [res[1,*,0]], color = 20, psym = 1, thick = 3
; IDL> oplot, [x0[1], x1[1]], [y0[1], y1[1]]
; IDL> oplot, [res[0,*,1]], [res[1,*,1]], color = 40, psym = 1, thick = 3
;
; @history
;           S. Masson (smasson\@lodyc.jussieu.fr)
;           July 5th, 2002
;
; @version
; $Id$
;
;-
FUNCTION cutsegment, x0, y0, x1, y1, n $
                   , ENDPOINTS=endpoints, ONSPHERE=onsphere
;
  compile_opt idl2, strictarrsubs
;
; number of segment
  nseg = n_elements(x0)
; number of point to find on each segment
  n2find = n+keyword_set(endpoints)
;
  IF keyword_set(onsphere) THEN BEGIN
; save the inputs arrays
    x0in = temporary(x0)
    y0in = temporary(y0)
    x1in = temporary(x1)
    y1in = temporary(y1)
    sp_cood = [transpose(x0in[*]),transpose(y0in[*]),replicate(1, 1, nseg)]
    rect_coord = CV_COORD(FROM_SPHERE = temporary(sp_cood), /TO_RECT, /DEGREES)
    x0 = rect_coord[0, *]
    y0 = rect_coord[1, *]
    z0 = rect_coord[2, *]
    rect_coord = -1 ;free memory
    sp_cood = [transpose(x1in[*]),transpose(y1in[*]),replicate(1, 1, nseg)]
    rect_coord = CV_COORD(FROM_SPHERE = temporary(sp_cood), /TO_RECT, /DEGREES)
    x1 = rect_coord[0, *]
    y1 = rect_coord[1, *]
    z1 = rect_coord[2, *]
    rect_coord = -1 ;free memory
  ENDIF
;
  resx = replicate(1, n2find)#x0[*]
  resx = temporary(resx)+(1./n*(findgen(n2find) $
                                +.5*(1-keyword_set(endpoints))))#(x1-x0)[*]
  resx = (temporary(resx))[*]
;
  resy = replicate(1, n2find)#y0[*]
  resy = temporary(resy)+(1./n*(findgen(n2find) $
                                +.5*(1-keyword_set(endpoints))))#(y1-y0)[*]
  resy = (temporary(resy))[*]

  IF keyword_set(onsphere) THEN BEGIN
    resz = replicate(1, n2find)#z0[*]
    resz = temporary(resz)+(1./n*(findgen(n2find) $
                                  +.5*(1-keyword_set(endpoints))))#(z1-z0)[*]
    resz = (temporary(resz))[*]

    rec_cood = [transpose(temporary(resx)), transpose(temporary(resy)) $
                , transpose(temporary(resz))]
    res = CV_COORD(FROM_RECT = temporary(rec_cood), /TO_SPHERE, /DEGREES)
; restore the input arrays
    x0 = temporary(x0in)
    y0 = temporary(y0in)
    x1 = temporary(x1in)
    y1 = temporary(y1in)
  ENDIF ELSE res = [transpose(temporary(resx)), transpose(temporary(resy))]

  res = reform(res[0:1, *], 2, n2find, nseg, /overwrite)

  RETURN, res
END