source: trunk/SRC/Matrix/extrac2.pro @ 230

;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
;
; @file_comments
; extraction of subdomains of matrices;
; Even if the subdomain is "pierced" (see the example)
; By default, IDL can make extractions of subdomain:
;
;      IDL> a=indgen(5,5)
;      IDL> print, a
;             0       1       2       3       4
;             5       6       7       8       9
;            10      11      12      13      14
;            15      16      17      18      19
;            20      21      22      23      24
;      IDL> print, a[[0,2],3]
;            15      17
;      IDL> print, a[[0,2],*]
;             0       2
;             5       7
;            10      12
;            15      17
;            20      22
; but
;      IDL> print, a[[0,2],[3,4]]
;            15      22
; while
;      IDL> print, extrac2(a,[0,2],[3,4])
;            15      17
;            20      22
;
; @categories 
; Utilities
;
; @param array {in}{required}
; a 1,2,3 or 4 dim input array
;
; @param index1 {in}{required}
; can have 2 forms:
; 1)a vector containing indexes of lines we want to keep
; 2)the string '*' if we want to keep all lines.
;
; @param index2 {in}{required}
; the same thing that index1 but for dim 2.
;
; @param index3 {in}{required}
; the same thing that index1 but for dim 3.
;
; @param index4 {in}{required}
; the same thing that index1 but for dim 4.
;
; @returns
; a matrix 1,2,3 or 4d extract from input array
; -1 in case of mistake
;
; @restrictions
;
; @examples
; I have a dim 2 matrix named A. I want extract a small intersection
; matrix 2d of the line 2,3 and 7 and of the column 0 and 1:
; 
; IDL> res=extrac2(A,[2,3,7],[0,1])
;
; other ex:
; IDL> print, a
; a b c
; d e f
; g h i
; IDL> print, extrac2(a,[0,2],[0,2])
; a c
; g i
;
; @history Sebastien Masson (smasson\@lodyc.jussieu.fr)
;                       12/1/1999
;                       29/4/1999: correction of a bug and complement of the heading
;
; @version $Id$
;
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
FUNCTION extrac2, array, index1, index2, index3, index4
;------------------------------------------------------------
;
  compile_opt idl2, strictarrsubs
;
   taille = size(array)
;------------------------------------------------------------
; test of the numbre of parameters
; and of the nature of the index (for THE case 'x')
;------------------------------------------------------------
   if n_params() NE taille[0]+1 THEN $
    return, report('we need as many indexes as the number of dimensions of the input array')
   IF n_params() GE 5 THEN BEGIN
      if size(index4,/type) EQ 7 then index4 = lindgen(taille[4]) $
      ELSE index4 = long(index4)
      nt = n_elements(index4)
   ENDIF
   IF n_params() GE 4 THEN BEGIN
      if size(index3,/type) EQ 7 then index3 = lindgen(taille[3]) $
      ELSE index3 = long(index3)
      nz = n_elements(index3)
   ENDIF
   IF n_params() GE 3 THEN BEGIN
      if size(index2,/type) EQ 7 then index2 = lindgen(taille[2]) $
      ELSE index2 = long(index2)
      ny = n_elements(index2)
   ENDIF
   IF n_params() GE 2 THEN BEGIN
      if size(index1,/type) EQ 7 then index1 = lindgen(taille[1]) $
      ELSE index1 = long(index1)
      nx = n_elements(index1)
   ENDIF

;------------------------------------------------------------
; construction of an array of indexes and of results following the size of array
;------------------------------------------------------------
  case taille[0] of
      1:res = array[index1]
      2:BEGIN
         index = index1#replicate(1, ny)+taille[1]*replicate(1, nx)#index2
         res = array[index]
      END
      3:BEGIN
         index = index1#replicate(1, ny)+taille[1]*replicate(1, nx)#index2
         index = temporary(index[*])#replicate(1, nz) $
          +taille[1]*taille[2]*replicate(1, nx*ny)#index3
         res = array[reform(index, nx, ny, nz, /over)]
      END
      4:BEGIN
         index = index1#replicate(1, ny)+taille[1]*replicate(1, nx)#index2
         index = temporary(index[*])#replicate(1, nz) $
          +taille[1]*taille[2]*replicate(1, nx*ny)#index3
         index = temporary(index[*])#replicate(1, nt) $
          +taille[1]*taille[2]*taille[3]*replicate(1, nx*ny*nz)#index4
         res = array[reform(index, nx, ny, nz, nz, /over)]
      END
   endcase


;------------------------------------------------------------
;------------------------------------------------------------
   return, res
end