code: func (s *Surface) GradientAt(point *Point) *Point {
  	if !s.bound.Contains(point) {
  		return NewPoint(0, 0)
  	}
  
  	xi, yi, deltaX, deltaY := s.gridCoordinate(point)
  
  	xi1 := xi + 1
  	if limit := s.Width - 1; xi1 > limit {
  		xi = limit - 1
  		xi1 = limit
  		deltaX = 1.0
  	}
  
  	yi1 := yi + 1
  	if limit := s.Height - 1; yi1 > limit {
  		yi = limit - 1
  		yi1 = limit
  		deltaY = 1.0
  	}
  
  	u1 := s.Grid[xi][yi]*(1-deltaX) + s.Grid[xi1][yi]*deltaX
  	u2 := s.Grid[xi][yi1]*(1-deltaX) + s.Grid[xi1][yi1]*deltaX
  
  	w1 := (1 - deltaY) * (s.Grid[xi1][yi] - s.Grid[xi][yi])
  	w2 := deltaY * (s.Grid[xi1][yi1] - s.Grid[xi][yi1])
  
  	return NewPoint((w1+w2)/s.gridBoxWidth(), (u2-u1)/s.gridBoxHeight())
  }
docs: // GradientAt returns the surface gradient at the given point.
  // Bilinearlly interpolates the grid cell to find the gradient.