;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
;+
; NAME:extrait
;
; PURPOSE:extraction de sous domaines de matrices. Meme si le sous
; domaine est "troue". cf. : l''exemple
; Par defaut IDL peut faire des extractions de sous domaines:
;
;      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
; mais 
;      IDL> print, a[[0,2],[3,4]] 
;            15      22
; alors que
;      IDL> print, extrait(a,[0,2],[3,4])  
;            15      17
;            20      22
;
; CATEGORY:bidouille avec les matrices
;
; CALLING SEQUENCE:res=extrait(tab,indicex[,indicey[,indicez,[indicet]]])
;
; INPUTS:
;        tab: un tableau 1,2,3 ou 4 d
;
;        indicex: indicex peut avoir deux formes:
;              1)un vecteur contenant les indices des lignes a garder
;              2)le string '*' dans ce cas touts les lignes sont
;              gardees
;
;        indicey/z/t: la meme chose que indicex mais pour les
;        dimensions 2,3 et 4.
;
;        rq: il faut autant de vecteurs indice que tab a de dimensions
;
; KEYWORD PARAMETERS:
;
; OUTPUTS:
;         res: une matice 1,2,3 ou 4d extraite a partir de tab
;
; COMMON BLOCKS:
;
; SIDE EFFECTS:res=-1 en cas d''erreur
;
; RESTRICTIONS:
;
; EXAMPLE:
; j''ai une matrice A de dim 2. je veux en extraire une petite matrice
; 2d interscetion de la ligne 2,3 et 7 et de la colonne 0 et 1
;      
;      res=extrait(A,[2,3,7],[0,1])
;
;autre ex:
;      IDL> print, a 
;      a b c
;      d e f
;      g h i
;      IDL> print, extrait(a,[0,2],[0,2])        
;      a c
;      g i
;
; MODIFICATION HISTORY: Sebastien Masson (smasson@lodyc.jussieu.fr)
;                       12/1/1999
;                       29/4/1999: correction d''un bug et complement
;                       de l''en-tete
;-
;------------------------------------------------------------
;------------------------------------------------------------
;------------------------------------------------------------
FUNCTION extrait, tab, indicex, indicey, indicez, indicet
;------------------------------------------------------------
   taille = size(tab)
;------------------------------------------------------------
; test du nombre de parametres
; et de la nature de indice (pour LE cas 'x')
;------------------------------------------------------------
   if n_params() NE taille[0]+1 THEN $
    return, report('il faut autant d''indices que de dimension du tableau')
   IF n_params() GE 5 THEN BEGIN 
      if size(indicet,/type) EQ 7 then indicet = lindgen(taille[4]) $
      ELSE indicet = long(indicet)
      nt = n_elements(indicet) 
   ENDIF 
   IF n_params() GE 4 THEN BEGIN 
      if size(indicez,/type) EQ 7 then indicez = lindgen(taille[3]) $
      ELSE indicez = long(indicez)
      nz = n_elements(indicez) 
   ENDIF 
   IF n_params() GE 3 THEN BEGIN 
      if size(indicey,/type) EQ 7 then indicey = lindgen(taille[2]) $
      ELSE indicey = long(indicey)
      ny = n_elements(indicey) 
   ENDIF
   IF n_params() GE 2 THEN BEGIN 
      if size(indicex,/type) EQ 7 then indicex = lindgen(taille[1]) $
      ELSE indicex = long(indicex)
      nx = n_elements(indicex) 
   ENDIF
   
;------------------------------------------------------------
; construction du tableau d''indice et du resultat suivant la taille
; de tab
;------------------------------------------------------------
  case taille[0] of
      1:res = tab[indicex]
      2:BEGIN
         indice = indicex#replicate(1, ny)+taille[1]*replicate(1, nx)#indicey
         res = tab[indice]
      END 
      3:BEGIN
         indice = indicex#replicate(1, ny)+taille[1]*replicate(1, nx)#indicey
         indice = temporary(indice[*])#replicate(1, nz) $
          +taille[1]*taille[2]*replicate(1, nx*ny)#indicez
         res = tab[reform(indice, nx, ny, nz, /over)]
      END 
      4:BEGIN
         indice = indicex#replicate(1, ny)+taille[1]*replicate(1, nx)#indicey
         indice = temporary(indice[*])#replicate(1, nz) $
          +taille[1]*taille[2]*replicate(1, nx*ny)#indicez
         indice = temporary(indice[*])#replicate(1, nt) $
          +taille[1]*taille[2]*taille[3]*replicate(1, nx*ny*nz)#indicet
         res = tab[reform(indice, nx, ny, nz, nz, /over)]
      END 
   endcase


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