;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
;
; @fie-comments 
; extractoin of subdomains of matrixes; 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
;
; @restrictions res=-1 in case of mistake
;
; 
; @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:
;      
;      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