COORDINATES IDL program by James Preiss

; +
;
; COORDINATES by James Alan Preiss
;
; PURPOSE:
;     generate multidemensional arrays of coordinates
;     to use in vectorized computation of functions.
;
; PARAMETERS:
;
;   MINS: scalar or one-dimensional array containing the minimum
;     desired coordinate value for each dimension.
;
;   MAXES: scalar or one-dimensional array containing the maximum
;     desired coordinate value for each dimension.
;
;   MINS and MAXES must have the same number of elements.
;
;   TYPE: IDL numerical typecode. If TYPE is not given, Long64
;     integers will be used. Currently, COORDINATES does not
;     check if the chosen data type is large enough to hold the
;     range of values that will be generated.
;
; RETURNS:
;   A pointer array, where element I references the multidimensional
;   array containing coordinate values for the Ith dimension.
;
; EXAMPLES:
;
;   In this example, we generate the vector field for a
;   Lorenz Attractor, using typical values for the parameters.
;
;       coords = coordinates([-20,-20,0],[20,20,40],4)
;       x = *coords[0]
;       y = *coords[1]
;       z = *coords[2]
;
;       dxdt = 10.0*(y - x)
;       dydt = x*(28.0 - z) - y
;       dzdt = x*y - (8.0/3)*z
;
; HISTORY:
;   Written April 7, 2010
;+

function coordinates, mins, maxes, type

    ;a type argument exists - make sure it is usable.
    ;notice that this DOES NOT check for overflow!
    if (n_params() gt 2) then begin
        ;types labeled 1 are unusable, 2 are unsigned
        typecds = bytarr(15)
        typecds[[0,1,6,7,8,9,10,11]] = 1
        typecds[[1,12,13,15]] = 2
        if (typecds[type] eq 1) then message, $
            'TYPE must be a non-complex numerical typecode.'
        if (typecds[type] eq 2) then begin
            if (total(mins lt 0) ne 0) then message, $
                'Unsigned data type cannot contain negative values.'
        endif
    endif

    ;check for proper dimensions in mins and maxes
    ;each should be an n-element vector where n is the number
    ;of dimensions for which we are generating coordinates
    mindims = size(mins, /dimensions)
    maxdims = size(maxes, /dimensions)

    if (~array_equal(mindims, maxdims) || $
        n_elements(mindims) gt 1 || n_elements(maxdims) gt 1) then $
            message, 'MINS and MAXES must be scalars ' + $
                     'or one-dimensional arrays of the same length.'

    ;check for nonsensical ranges -
    ;the min must be less than or equal to the max
    sizes = maxes - mins + 1
    if (total(sizes lt 1) ne 0) then message, $
        'Maxes must be less than or equal to Mins.'

    dims = mindims
    if (dims eq 0) then dims = 1


    out = ptrarr(dims)

    for i=0, dims[0]-1 do begin
        ;make a 1D vector containing the indices for this dimension
        out[i] = ptr_new(make_array(sizes[i], /index, type=type))
        ;shift index vector it contains values from min to max, inclusive
        *out[i] +=+ mins[i]
        ;make a vector to feed into reform for dimensional juggling
        rvec = replicate(1,i+1)
        ;place the index vector into the proper dimension
        rvec[i] = sizes[i]
        ;replicate the index vector across all dimensions
        *out[i] = rebin(reform(*out[i], rvec), sizes)
    endfor

    return, out
end