2nd order derivative

2nd order derivatives implemented through the [Hessian matrix](http://en.wikipedia.org/wiki/Hessian_matrix)

Addition to high-level interface:

    # example from test code
    julia> y = rand(7,8);
    julia> yi = InterpGrid(y, BCnil, InterpQuadratic);
    julia> x = [2.2, 3.1];
    julia> v,g,h = valgradhess(yi, x);  # new!
    julia> h
      2x2 Array{Float64,2}:
       -0.0846558   0.148568 
        0.148568   -0.0781553

Addition to low-level interface:

    julia> yc = copy(y)
    julia> interp_invert!(yc, BCnan, InterpQuadratic)
    julia> ic = InterpGridCoefs(yc, InterpQuadratic)
    julia> set_hessian_coordinate(ic, 1, 2) # new! cross-derivative, d2/dxdy
    julia> set_hessian_coordinate(ic, 2, 2) # new! second-order derivative d2/dy2

*Note:*

In order to get `valgradhess` working correctly, I had to introduce another parameter to `set_position`, as part of the low-level interface. Currently, the change I made is breaking (one needs to provide a second boolean argument, indicating wether to calculate hessian coefficients) but one might want to make that change non-breaking somehow. I don't know if that function is used a lot outside of the library (it is exported), but I couldn't make up a non-breaking way to change that without keeping its api consistent with itself.

A better (more future-proof), but still breaking, change, could be to make `calc_grad` and `calc_hess` keyword arguments, allowing to add further ones if desired.

Todo:

* [x] Documentation.