FUNCTION vu_gettime, hdr, ff, $
        ut      = ut        , $
        count   = count     , $
        ff_ut   = ff_ut     , $
        fixgrid = fixgrid   , $
        bracket = bracket   , $
        nowrap  = nowrap    , $
        silent  = silent
    @compile_opt.pro            ; On error, return to caller

is_time = IsTime(hdr)
is_hdr  = IsType(hdr, /struct) AND NOT is_time

count = n_elements(hdr)

CASE count OF
0: BEGIN

    CASE 1 OF
    is_hdr : hdr_ut = {vu_header}
    is_time: hdr_ut = !TheTime
    ELSE   : hdr_ut = 0*hdr[0]
    ENDCASE

END
ELSE: BEGIN

    InitVar, fixgrid, /key

    utt = Carrington(ut, /get_time)

    CASE 1 OF
    is_hdr : time = vu_get(hdr, /uttime)
    is_time: time = hdr
    ELSE   : time = Carrington(hdr)
    ENDCASE

    uday = TimeUnit(/day)

    IF fixgrid THEN BEGIN                   ; Change time to nearest grid time

        dt = TimeOp(/subtract, time, utt, uday)

        bracket = (where(dt GE 0))[0]

        IF bracket EQ -1 THEN       $
            bracket = count-1       $
        ELSE IF bracket GT 0 THEN   $
            bracket = bracket-(abs(dt[bracket-1]) LT dt[bracket])

        IF dt[bracket] NE 0 THEN BEGIN
            say, threshold=0, silent=silent, 'shifted by'+strcompress(dt[bracket])+' days to nearest gridpoint'
            utt = time[bracket]
        ENDIF

    ENDIF

    dt = TimeOp(/subtract, time, utt, uday)

    bracket = (where(dt EQ 0))[0]

    IF bracket EQ -1 THEN BEGIN
        bracket = (where(dt GT 0))[0]
        IF bracket EQ -1 THEN bracket = count
        bracket = ((bracket+[-1,0]) > 0) < (count-1)
        IF bracket[0] EQ bracket[1] THEN bracket = bracket[0]
    ENDIF

    CASE 1 OF
    is_hdr: BEGIN

        ; The index pair 'bracket' brackets the time 'ut'.
        ; Interpolate between the neigbouring times to get the matrix at time ut

        CASE n_elements(bracket) OF     ; ut is grid time: no interpolation needed
        1: BEGIN
            hdr_ut = hdr[bracket]
            ff_ut  = ff[*,*,*,*,bracket]
        END
        2: BEGIN                        ; Interpolate between times bracketing ut

            ; Interpolate on time array index as a function of Carrington time
            ; to calculate the density matrix at the proper (interpolated) time.

            w1 = bracket[0]
            w2 = bracket[1]

            w  = TimeInterpol(indgen(count), time, utt, quadratic=count gt 2)
            ;print, 'vu_gettime:', w2-w,w-w1,Carrington(vu_get(hdr[w1], /uttime)),  $
            ;   Carrington(utt), Carrington(vu_get(hdr[w2], /uttime)), format='(A,5F10.4)'

            hdr_ut = vu_mean(hdr[w1:w2], ff[*,*,*,*,w1:w2], [w2-w,w-w1], tt_weight=utt, $
                ff_mean=ff_ut, nowrap=nowrap)   ; /setbad)

        END
        ENDCASE

    END
    IsTime(ut): hdr_ut = utt
    ELSE      : hdr_ut = Carrington(utt, /get_variable)
    ENDCASE

    count = 1

END
ENDCASE

RETURN, hdr_ut  &  END