New cubed sphere 2 grid

[mo-jonasganderton]

Description

Initial implementation of a new cubed sphere grid.

@odlomax and @wdeconinck can both have a look over to recommend any changes - I expect this will have changed once everything else has also been implemented.

Using this implementation generates the following lonlat points - n=12, red=tile0, purple=tile5, darkest point of each colour is tij=0, lightest point is tij=143. This matches the sequence of the points in the LFRic cubed sphere.

Impact

Once everything else is also implemented, this makes the cubed sphere more maintainable.

Checklist

• I have performed a self-review of my own code
• I have made corresponding changes to the documentation
• I have run the unit tests before creating the PR

[FussyDuck]

All committers have signed the CLA.

[odlomax]

Looks good to see a much simpler CS grid! I've added some changes. It's mostly style related. I'm sure @wdeconinck can make some more comments/suggestions.

After this, we should probably plumb the grid into the grid builder.

@wdeconinck, what do you think?

[odlomax]

You should add some tests to check the return values of these functions. It should be relatively straightforward to add a new section to your test that checks the basic API.

[odlomax]

I don't believe we need this if, else structure around the code. Also, for completeness we should add the domain spec to grid_spec.

[odlomax]

I think this tij index logic could be a bit more readable and with less duplication if we did:

using CSIndices = std::array<idx_t, 3>;

// Useful to have this a public method
idx_t CubedSphere2::N() const {return N_;}

CSIndices CubedSphere2::get_cs_indices(gidx_t n) const {
   ATLAS_ASSERT(n >= size());
   const idx_t tile_size = N() * N();
   const idx_t t = n / tile_size;
   const idx_t ij = n % tile_size;
   const idx_t j = ij / N();
   const idx_t i = ij % N();
   return {t, i, j};
}

You can then use auto [t, i, j] = get_cs_indicies(...) to set separate t, i and j values if you need to.

[odlomax]

might as well make nTiles_ constexpr.

[odlomax]

This should be more expressive. For example:

using col = int;
using row = std::array<col, 3>;
using matrix = std::array<row, 3>;

std::array<matrix, nTiles_> rotations_ {matrix{row{0, 0, 1},
                                              row{1, 0, 0},
                                              row{0, -1, 0}},
                                       matrix...}

[odlomax]

I think the intent is clearer if we define:

using PointAlphaBeta = Point2;

PointAlphaBeta CubedSphere2::ij_to_curvilinear_coord(idx_t i, idx_t j) {
  const auto get_coord = [&](idx_t idx){
    return M_PI / 2 * (-0.5 + (0.5 + idx) / N());
  };
  return {get_coord(i), get_coord(j)};
}

PointXY CubedSphere2::curvilinear_to_tangent_coord(const PointAlphaBeta& curvi_coord) {
  return {std::tan(curvi_coord[0]), std::tan(curvi_coord[1])};
}

PointXYZ CubedSphere2::tangent_to_xyz_coord(const PointXY& tan_coord) {
  /// Perform coordinate rotations...
}

[odlomax]

You want to call projection_.xy2lonlat(xy) here.

[odlomax]

eckit has cartesian to spherical transforms already implemented here.

Also the degrees to radians constants are already defined here.