geoserver - How properly calculate BBox for mocked GetMap request?

I am trying to imitate GetMap requests for GeoServer and have an application which generates these requests (requests are for cached layers). The problem is that sometimes I get geowebcache-cache-result:HIT and sometimes geowebcache-cache-result:MISS with geowebcache-miss-reason:request does not align to grid(s) 'XXXX'. I understand that for these "miss" requests I have generated invalid BBox parameter. For example (pay attention to marked numbers):

...REQUEST=GetMap&BBOX=55|6|865.8627022989,6104809.344416159,55|8|155.8627022989,6106099.344416159...

results in HIT

when

...REQUEST=GetMap&BBOX=55|5|865.8627022989,6104809.344416159,55|7|155.8627022989,6106099.344416159...

results in MISS with geowebcache-miss-reason:request does not align to grid(s) 'XXXX'.

What would be the formula for generating correct BBox?

NOTE: I have checked and the generated BBox coordinates fit inside layer's and grid extent.

UPDATE: Question also could be formulated as - how to validate if BBox coordinates are aligned to the grid that is used?

Answer

The easy way to do this is to make WMTS requests instead of WMS requests since then all you need to do is specify the zoom level, row and column number to get a hit.

If you have to calculate the tile bounds you should fetch the WMTS capabilities document. Then under the layer you are asking for will be a TileMatrixSetLink that tells you the projection and limits for each level of tiles.

 
    EPSG:4326
    
      
        EPSG:4326:0
        0
        0

        0
        0
      
      
        EPSG:4326:1
        0
        0
        0
        1
      

      
      [.....]

This gives you the row and column ranges for each zoom level, don't request a tile outside these limits. Then look for the TileMatrixSet towards the end of the file that matches this name.

  
  EPSG:4326
  urn:ogc:def:crs:EPSG::4326
  
    EPSG:4326:0
    2.795411320143589E8

    90.0 -180.0
    256
    256
    2
    1
  
  
    EPSG:4326:1
    1.3977056600717944E8
    90.0 -180.0

    256
    256
    4
    2
  
  [......]

Now for any zoom level, you can look up the origin of the grid and then calculate the bounds for any tile as:

pixelSize = scaleDenom * 2.8e-4 
if CRS is degrees then / by 111319 else convert from metres if needed

tileWidth = 256 * pixelSize
tileHeight = 256 * pixelSize

minX = col * tileWidth + originX
maxY = originY - row * tileHeight 
maxX = minX + tileWidth
minY = maxY - tileHeight

You can see a full implementation with axis order handling and all the other exciting issues from GeoTools.