; $Id: array_indices.pro,v 1.3 2004/01/21 15:54:48 scottm Exp $
; Copyright (c) 2002-2004, Research Systems, Inc.  All rights reserved.
;       Unauthorized reproduction prohibited.
;+
; NAME:
;   ARRAY_INDICES
;
; PURPOSE:
;   Given an input array, this function converts one-dimensional
;   subscripts back into the corresponding multi-dimensional subscripts.
;
;
; CALLING SEQUENCE:
;   Result = ARRAY_INDICES(Array, Index)
;
;
; RETURN VALUE:
;   If Index is a scalar, then the Result will be a vector containing
;   the M multi-dimensional subscripts. If Index is a vector containing
;   N elements, then the Result will be a (M x N) array, with each row
;   containing the multi-dimensional subscripts corresponding to that Index.
;
;
; INPUTS:
;   Array: An array of any type, whose dimensions should be used in
;          converting the subscripts.
;
;   Index: A scalar or vector containing the subscript(s) to be converted.
;
;
; KEYWORD PARAMETERS:
;   None.
;
;
; EXAMPLE:
;
;
; MODIFICATION HISTORY:
;   Written by:  CT, RSI, October 2002
;   Modified:
;
;-

function array_indices, array, indices

    compile_opt idl2

    ON_ERROR, 2

    ; Error checking.
    if (N_PARAMS() ne 2) then $
        MESSAGE, 'Incorrect number of arguments.'


    ; Check for valid types.
    switch SIZE(indices, /TYPE) of

        6:      ; complex
        7:      ; string
        8:      ; struct
        9:      ; dcomplex
        10:     ; pointer
        11: $   ; objref
            MESSAGE, 'Index must be an integer.'
        else: ; okay, do nothing

    endswitch


    ; Check for valid index range.
    mn = MIN(indices, MAX=mx)

    if ((mn lt 0) or (mx ge N_ELEMENTS(array))) then $
        MESSAGE, 'Index out of range.'


    ndim = SIZE(array, /N_DIMENSIONS)

    ; If we have a scalar or vector, we're done. Just return the indices.
    if (ndim le 1) then $
        return, indices


    ni = N_ELEMENTS(indices)
    dim = SIZE(array, /DIMENSIONS)

    ; Result type is either Long or Long64 depending upon platform.
    result = MAKE_ARRAY(ndim, ni, /NOZERO, TYPE=SIZE(dim, /TYPE))

    if (ndim eq 2) then begin   ; Simple 2D case.

        result[0, *] = indices mod dim[0]  ; first dimension
        result[1, *] = indices/dim[0]      ; second dimension

    endif else begin  ; Multidimensional case.

        ; Product of all "previous" dimensions.
        dimProduct = PRODUCT(dim, /CUMULATIVE)

        temp = indices   ; make a working copy

        for d=ndim-1, 1, -1 do begin
            ; Remove indices from even higher dimensions.
            if (d lt ndim-1) then $
                temp mod= dimProduct[d]
            ; Indices for higher dimensions.
            result[d, *] = temp/dimProduct[d-1]
        endfor

        ; Indices for first dimension.
        result[0, *] = temp mod dim[0]

    endelse

    return, result

end