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

;+
;
; @file_comments cut p segments into p*n equal parts
;
; @categories basic work
;
; @examples  
; res = cutsegment(x0, y0, x1, y1, n)
;
;         @param x0,y0 and x1,y1  {in}{required}  1d arrays of p elements, the coordinates of
;         the endpoints of the p segmements
;         @param  n {in}{required}  the number of pieces we want to cut each segment
;
;
; @keyword         /endpoints see ouputs
;
; @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
;        defaut: a 3d array (2,n,p) that gives the coordinates of the
;        middle of the cutted 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
;-
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