FUNCTION index2ij,index,sizes

;+
; NAME:
;       index2ij
; CALLING SEQUENCE:
;       coor=index2ij(ind,sizes)
; PURPOSE:
;       returns the pixel postion of index ind in an array
; INPUTS:
;       index = the index for which coordinates are needed (value or 1D)
;	   sizes = Array dimensions
; OPTIONAL INPUT:
;	   None
; OUTPUTS:
;	   coordinates corresponding to the index
; PROCEDURE AND SUBROUTINE USED
;       
; SIDE EFFECTS:
; MODIFICATION HISTORY:
;       written by Jean-Philippe Bernard 09-92
;	modified JPB 11-Sep-92 for 2D arrays
;	modified JPB to allow complete indexes to be handled
;       modified JPB 21-June-94 for large sets of indexes
;       modified JPB 30-July-01 replaced fix by long allow large values of output
; modified 19 Jul 2016 JPB+AL fixed a bug in division remaining affecting first index of large arrays
; EXAMPLE:
;       a=fltarr(12,15,17)
;       a(10,10,5)=1 & a(11,5,16)=1
;       ind=where(a NE 0,count)
;       ij=index2ij(ind,[12,15,17])
;       help,ij
;       print,ij
;       ind2=ij2index(ij,[12,15,17])
;       print,ind,ind2
;-

;on_error,2

IF n_params(0) NE 2 THEN BEGIN
  print,'Calling sequence: coor=index2ij(index,sizes)'
  print,'                  index can be 1D or single value'
  goto,sortie
ENDIF

taille=n_elements(sizes)
si=size(index)

;print,taille
res=lonarr(si[1],taille)
FOR i=0L,si[1]-1 DO BEGIN
  ind_val=index[i]
  FOR tt=taille,2,-1 DO BEGIN
    size_prod=1.
    FOR j=0L,tt-2 DO size_prod=size_prod*sizes[j]
    res[i,tt-1]=long(ind_val/size_prod)
    ind_val=ind_val-long(size_prod*res[i,tt-1]) ;remains in division
  ENDFOR
;  res(i,1)=fix(ind_val/sizes(0))
  res[i,0]=long(ind_val)
  ;stop
ENDFOR  

return,res

sortie:

end