pro hsi_matadd, a, adiag, b, bdiag, result, resultdiag

resultdiag = adiag AND bdiag

if (adiag EQ bdiag) then begin
    result = a+b
    return
endif

nbins = min([n_elements(a),n_elements(b)])  ;Whichever is a vector

if adiag AND NOT bdiag then begin
    result = b
    for i=0,nbins-1 do result[i,i] = result[i,i]+a[i]
    return
endif

if NOT adiag AND bdiag then begin
    result = a
    for i=0,nbins-1 do result[i,i] = result[i,i]+b[i]
    return
endif

end