User: Explain this code to me

/* Copyright 2002-2024 CS GROUP
 * Licensed to CS GROUP (CS) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * CS licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *   http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.orekit.models.earth.tessellation;

import java.util.ArrayList;
import java.util.Collection;
import java.util.IdentityHashMap;
import java.util.Iterator;
import java.util.LinkedList;
import java.util.List;
import java.util.Map;
import java.util.NoSuchElementException;
import java.util.Queue;

import org.hipparchus.exception.LocalizedCoreFormats;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.geometry.partitioning.BSPTree;
import org.hipparchus.geometry.partitioning.Hyperplane;
import org.hipparchus.geometry.partitioning.RegionFactory;
import org.hipparchus.geometry.partitioning.SubHyperplane;
import org.hipparchus.geometry.spherical.oned.ArcsSet;
import org.hipparchus.geometry.spherical.twod.Circle;
import org.hipparchus.geometry.spherical.twod.S2Point;
import org.hipparchus.geometry.spherical.twod.Sphere2D;
import org.hipparchus.geometry.spherical.twod.SphericalPolygonsSet;
import org.hipparchus.geometry.spherical.twod.SubCircle;
import org.hipparchus.util.FastMath;
import org.hipparchus.util.MathUtils;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.bodies.OneAxisEllipsoid;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitInternalError;

/** Class used to tessellate an interest zone on an ellipsoid in either
 * {@link Tile tiles} or grids of {@link GeodeticPoint geodetic points}.
 * <p>
 * This class is typically used for Earth Observation missions, in order to
 * create tiles or grids that may be used as the basis of visibility event
 * detectors. Tiles are used when surface-related elements are needed, the
 * tiles created completely cover the zone of interest. Grids are used when
 * point-related elements are needed, the points created lie entirely within
 * the zone of interest.
 * </p>
 * <p>
 * One should note that as tessellation essentially creates a 2 dimensional
 * almost Cartesian map, it can never perfectly fulfill geometrical dimensions
 * because neither sphere nor ellipsoid are developable surfaces. This implies
 * that the tesselation will always be distorted, and distortion increases as
 * the size of the zone to be tessellated increases.
 * </p>
 * @author Luc Maisonobe
 * @since 7.1
 */
public class EllipsoidTessellator {

    /** Safety limit to avoid infinite loops during tesselation due to numerical noise.
     * @since 10.3.1
     */
    private static final int MAX_ITER = 1000;

    /** Number of segments tiles sides are split into for tiles fine positioning. */
    private final int quantization;

    /** Aiming used for orienting tiles. */
    private final TileAiming aiming;

    /** Underlying ellipsoid. */
    private final OneAxisEllipsoid ellipsoid;

    /** Simple constructor.
     * <p>
     * The {@code quantization} parameter is used internally to adjust points positioning.
     * For example when quantization is set to 4, a complete tile that has 4 corner points
     * separated by the tile lengths will really be computed on a grid containing 25 points
     * (5 rows of 5 points, as each side will be split in 4 segments, hence will have 5
     * points). This quantization allows rough adjustment to balance margins around the
     * zone of interest and improves geometric accuracy as the along and across directions
     * are readjusted at each points.
     * </p>
     * <p>
     * It is recommended to use at least 2 as the quantization parameter for tiling. The
     * rationale is that using only 1 for quantization would imply all points used are tiles
     * vertices, and hence would lead small zones to generate 4 tiles with a shared vertex
     * inside the zone and the 4 tiles covering the four quadrants at North-West, North-East,
     * South-East and South-West. A quantization value of at least 2 allows to shift the
     * tiles so the center point is an inside point rather than a tile vertex, hence allowing
     * a single tile to cover the small zone. A value even greater like 4 or 8 would allow even
     * finer positioning to balance the tiles margins around the zone.
     * </p>
     * @param ellipsoid underlying ellipsoid
     * @param aiming aiming used for orienting tiles
     * @param quantization number of segments tiles sides are split into for tiles fine positioning
     */
    public EllipsoidTessellator(final OneAxisEllipsoid ellipsoid, final TileAiming aiming,
                                final int quantization) {
        this.ellipsoid    = ellipsoid;
        this.aiming       = aiming;
        this.quantization = quantization;
    }

    /** Tessellate a zone of interest into tiles.
     * <p>
     * The created tiles will completely cover the zone of interest.
     * </p>
     * <p>
     * The distance between a vertex at a tile corner and the vertex at the same corner
     * in the next vertex are computed by subtracting the overlap width (resp. overlap length)
     * from the full width (resp. full length). If for example the full width is specified to
     * be 55 km and the overlap in width is specified to be +5 km, successive tiles would span
     * as follows:
     * </p>
     * <ul>
     *   <li>tile 1 covering from   0 km to  55 km</li>
     *   <li>tile 2 covering from  50 km to 105 km</li>
     *   <li>tile 3 covering from 100 km to 155 km</li>
     *   <li>...</li>
     * </ul>
     * <p>
     * In order to achieve the same 50 km step but using a 5 km gap instead of an overlap, one would
     * need to specify the full width to be 45 km and the overlap to be -5 km. With these settings,
     * successive tiles would span as follows:
     * </p>
     * <ul>
     *   <li>tile 1 covering from   0 km to  45 km</li>
     *   <li>tile 2 covering from  50 km to  95 km</li>
     *   <li>tile 3 covering from 100 km to 155 km</li>
     *   <li>...</li>
     * </ul>
     * @param zone zone of interest to tessellate
     * @param fullWidth full tiles width as a distance on surface, including overlap (in meters)
     * @param fullLength full tiles length as a distance on surface, including overlap (in meters)
     * @param widthOverlap overlap between adjacent tiles (in meters), if negative the tiles
     * will have a gap between each other instead of an overlap
     * @param lengthOverlap overlap between adjacent tiles (in meters), if negative the tiles
     * will have a gap between each other instead of an overlap
     * @param truncateLastWidth if true, the first tiles strip will be started as close as
     * possible to the zone of interest, and the last tiles strip will have its width reduced
     * to also remain close to the zone of interest; if false all tiles strip will have the
     * same {@code fullWidth} and they will be balanced around zone of interest
     * @param truncateLastLength if true, the first tile in each strip will be started as close as
     * possible to the zone of interest, and the last tile in each strip will have its length reduced
     * to also remain close to the zone of interest; if false all tiles in each strip will have the
     * same {@code fullLength} and they will be balanced around zone of interest
     * @return a list of lists of tiles covering the zone of interest,
     * each sub-list corresponding to a part not connected to the other
     * parts (for example for islands)
     */
    public List<List<Tile>> tessellate(final SphericalPolygonsSet zone,
                                       final double fullWidth, final double fullLength,
                                       final double widthOverlap, final double lengthOverlap,
                                       final boolean truncateLastWidth, final boolean truncateLastLength) {

        final double                  splitWidth  = (fullWidth  - widthOverlap)  / quantization;
        final double                  splitLength = (fullLength - lengthOverlap) / quantization;
        final Map<Mesh, List<Tile>>   map         = new IdentityHashMap<>();
        final RegionFactory<Sphere2D> factory     = new RegionFactory<>();
        SphericalPolygonsSet          remaining   = (SphericalPolygonsSet) zone.copySelf();
        S2Point                       inside      = getInsidePoint(remaining);

        int count = 0;
        while (inside != null) {

            if (++count > MAX_ITER) {
                throw new OrekitException(LocalizedCoreFormats.MAX_COUNT_EXCEEDED, MAX_ITER);
            }

            // find a mesh covering at least one connected part of the zone
            final List<Mesh.Node> mergingSeeds = new ArrayList<>();
            Mesh mesh = new Mesh(ellipsoid, zone, aiming, splitLength, splitWidth, inside);
            mergingSeeds.add(mesh.getNode(0, 0));
            List<Tile> tiles = null;
            while (!mergingSeeds.isEmpty()) {

                // expand the mesh around the seed
                neighborExpandMesh(mesh, mergingSeeds, zone);

                // extract the tiles from the mesh
                // this further expands the mesh so tiles dimensions are multiples of quantization,
                // hence it must be performed here before checking meshes independence
                tiles = extractTiles(mesh, zone, lengthOverlap, widthOverlap, truncateLastWidth, truncateLastLength);

                // check the mesh is independent from existing meshes
                mergingSeeds.clear();
                for (final Map.Entry<Mesh, List<Tile>> entry : map.entrySet()) {
                    if (!factory.intersection(mesh.getCoverage(), entry.getKey().getCoverage()).isEmpty()) {
                        // the meshes are not independent, they intersect each other!

                        // merge the two meshes together
                        mesh = mergeMeshes(mesh, entry.getKey(), mergingSeeds);
                        map.remove(entry.getKey());
                        break;

                    }
                }

            }

            // remove the part of the zone covered by the mesh
            remaining = (SphericalPolygonsSet) factory.difference(remaining, mesh.getCoverage());
            inside    = getInsidePoint(remaining);

            map.put(mesh, tiles);

        }

        // concatenate the lists from the independent meshes
        final List<List<Tile>> tilesLists = new ArrayList<>(map.size());
        for (final Map.Entry<Mesh, List<Tile>> entry : map.entrySet()) {
            tilesLists.add(entry.getValue());
        }

        return tilesLists;

    }

    /** Sample a zone of interest into a grid sample of {@link GeodeticPoint geodetic points}.
     * <p>
     * The created points will be entirely within the zone of interest.
     * </p>
     * @param zone zone of interest to sample
     * @param width grid sample cells width as a distance on surface (in meters)
     * @param length grid sample cells length as a distance on surface (in meters)
     * @return a list of lists of points sampling the zone of interest,
     * each sub-list corresponding to a part not connected to the other
     * parts (for example for islands)
     */
    public List<List<GeodeticPoint>> sample(final SphericalPolygonsSet zone,
                                            final double width, final double length) {

        final double                         splitWidth  = width  / quantization;
        final double                         splitLength = length / quantization;
        final Map<Mesh, List<GeodeticPoint>> map         = new IdentityHashMap<>();
        final RegionFactory<Sphere2D>        factory     = new RegionFactory<>();
        SphericalPolygonsSet                 remaining   = (SphericalPolygonsSet) zone.copySelf();
        S2Point                              inside      = getInsidePoint(remaining);

        int count = 0;
        while (inside != null) {

            if (++count > MAX_ITER) {
                throw new OrekitException(LocalizedCoreFormats.MAX_COUNT_EXCEEDED, MAX_ITER);
            }

            // find a mesh covering at least one connected part of the zone
            final List<Mesh.Node> mergingSeeds = new ArrayList<>();
            Mesh mesh = new Mesh(ellipsoid, zone, aiming, splitLength, splitWidth, inside);
            mergingSeeds.add(mesh.getNode(0, 0));
            List<GeodeticPoint> sample = null;
            while (!mergingSeeds.isEmpty()) {

                // expand the mesh around the seed
                neighborExpandMesh(mesh, mergingSeeds, zone);

                // extract the sample from the mesh
                // this further expands the mesh so sample cells dimensions are multiples of quantization,
                // hence it must be performed here before checking meshes independence
                sample = extractSample(mesh);

                // check the mesh is independent from existing meshes
                mergingSeeds.clear();
                for (final Map.Entry<Mesh, List<GeodeticPoint>> entry : map.entrySet()) {
                    if (!factory.intersection(mesh.getCoverage(), entry.getKey().getCoverage()).isEmpty()) {
                        // the meshes are not independent, they intersect each other!

                        // merge the two meshes together
                        mesh = mergeMeshes(mesh, entry.getKey(), mergingSeeds);
                        map.remove(entry.getKey());
                        break;

                    }
                }

            }

            // remove the part of the zone covered by the mesh
            remaining = (SphericalPolygonsSet) factory.difference(remaining, mesh.getCoverage());
            inside    = getInsidePoint(remaining);

            map.put(mesh, sample);

        }

        // concatenate the lists from the independent meshes
        final List<List<GeodeticPoint>> sampleLists = new ArrayList<>(map.size());
        for (final Map.Entry<Mesh, List<GeodeticPoint>> entry : map.entrySet()) {
            sampleLists.add(entry.getValue());
        }

        return sampleLists;

    }

    /** Get an inside point from a zone of interest.
     * @param zone zone to mesh
     * @return a point inside the zone or null if zone is empty or too thin
     */
    private S2Point getInsidePoint(final SphericalPolygonsSet zone) {

        final InsidePointFinder finder = new InsidePointFinder(zone);
        zone.getTree(false).visit(finder);
        return finder.getInsidePoint();

    }

    /** Expand a mesh so it surrounds at least one connected part of a zone.
     * <p>
     * This part of mesh expansion is neighbors based. It includes the seed
     * node neighbors, and their neighbors, and the neighbors of their
     * neighbors until the path-connected sub-parts of the zone these nodes
     * belong to are completely surrounded by the mesh taxicab boundary.
     * </p>
     * @param mesh mesh to expand
     * @param seeds seed nodes (already in the mesh) from which to start expansion
     * @param zone zone to mesh
     */
    private void neighborExpandMesh(final Mesh mesh, final Collection<Mesh.Node> seeds,
                                    final SphericalPolygonsSet zone) {

        // mesh expansion loop
        boolean expanding = true;
        final Queue<Mesh.Node> newNodes = new LinkedList<>(seeds);
        int count = 0;
        while (expanding) {

            if (++count > MAX_ITER) {
                throw new OrekitException(LocalizedCoreFormats.MAX_COUNT_EXCEEDED, MAX_ITER);
            }

            // first expansion step: set up the mesh so that all its
            // inside nodes are completely surrounded by at least
            // one layer of outside nodes
            while (!newNodes.isEmpty()) {

                // retrieve an active node
                final Mesh.Node node = newNodes.remove();

                if (node.isInside()) {
                    // the node is inside the zone, the mesh must contain its 8 neighbors
                    addAllNeighborsIfNeeded(node, mesh, newNodes);
                }

            }

            // second expansion step: check if the loop of outside nodes
            // completely surrounds the zone, i.e. there are no peaks
            // pointing out of the loop between two nodes
            expanding = false;
            final List<Mesh.Node> boundary = mesh.getTaxicabBoundary(false);
            if (boundary.size() > 1) {
                Mesh.Node previous = boundary.get(boundary.size() - 1);
                for (final Mesh.Node node : boundary) {
                    if (meetInside(previous.getS2P(), node.getS2P(), zone)) {
                        // part of the mesh boundary is still inside the zone!
                        // the mesh must be expanded again
                        addAllNeighborsIfNeeded(previous, mesh, newNodes);
                        addAllNeighborsIfNeeded(node,     mesh, newNodes);
                        expanding = true;
                    }
                    previous = node;
                }
            }

        }

    }

    /** Extract tiles from a mesh.
     * @param mesh mesh from which tiles should be extracted
     * @param zone zone covered by the mesh
     * @param lengthOverlap overlap between adjacent tiles
     * @param widthOverlap overlap between adjacent tiles
     * @param truncateLastWidth true if we can reduce last tile width
     * @param truncateLastLength true if we can reduce last tile length
     * @return extracted tiles
     */
    private List<Tile> extractTiles(final Mesh mesh, final SphericalPolygonsSet zone,
                                    final double lengthOverlap, final double widthOverlap,
                                    final boolean truncateLastWidth, final boolean truncateLastLength) {

        final List<Tile>      tiles = new ArrayList<>();
        final List<RangePair> rangePairs = new ArrayList<>();

        final int minAcross = mesh.getMinAcrossIndex();
        final int maxAcross = mesh.getMaxAcrossIndex();
        for (Range acrossPair : nodesIndices(minAcross, maxAcross, truncateLastWidth)) {

            int minAlong = mesh.getMaxAlongIndex() + 1;
            int maxAlong = mesh.getMinAlongIndex() - 1;
            for (int c = acrossPair.lower; c <= acrossPair.upper; ++c) {
                minAlong = FastMath.min(minAlong, mesh.getMinAlongIndex(c));
                maxAlong = FastMath.max(maxAlong, mesh.getMaxAlongIndex(c));
            }

            for (Range alongPair : nodesIndices(minAlong, maxAlong, truncateLastLength)) {

                // get the base vertex nodes
                final Mesh.Node node0 = mesh.addNode(alongPair.lower, acrossPair.lower);
                final Mesh.Node node1 = mesh.addNode(alongPair.upper, acrossPair.lower);
                final Mesh.Node node2 = mesh.addNode(alongPair.upper, acrossPair.upper);
                final Mesh.Node node3 = mesh.addNode(alongPair.lower, acrossPair.upper);

                // apply tile overlap
                final S2Point s2p0 = node0.move(new Vector3D(-0.5 * lengthOverlap, node0.getAlong(),
                                                             -0.5 * widthOverlap,  node0.getAcross()));
                final S2Point s2p1 = node1.move(new Vector3D(+0.5 * lengthOverlap, node1.getAlong(),
                                                             -0.5 * widthOverlap,  node1.getAcross()));
                final S2Point s2p2 = node2.move(new Vector3D(+0.5 * lengthOverlap, node2.getAlong(),
                                                             +0.5 * widthOverlap,  node2.getAcross()));
                final S2Point s2p3 = node3.move(new Vector3D(-0.5 * lengthOverlap, node2.getAlong(),
                                                             +0.5 * widthOverlap,  node2.getAcross()));

                // create a quadrilateral region corresponding to the candidate tile
                final SphericalPolygonsSet quadrilateral =
                        new SphericalPolygonsSet(zone.getTolerance(), s2p0, s2p1, s2p2, s2p3);

                if (!new RegionFactory<Sphere2D>().intersection(zone.copySelf(), quadrilateral).isEmpty()) {

                    // the tile does cover part of the zone, it contributes to the tessellation
                    tiles.add(new Tile(toGeodetic(s2p0), toGeodetic(s2p1), toGeodetic(s2p2), toGeodetic(s2p3)));
                    rangePairs.add(new RangePair(acrossPair, alongPair));

                }

            }
        }

        // ensure the taxicab boundary follows the built tile sides
        // this is done outside of the previous loop in order
        // to avoid one tile changing the min/max indices of the
        // neighboring tile as they share some nodes that will be enabled here
        for (final RangePair rangePair : rangePairs) {
            for (int c = rangePair.across.lower; c < rangePair.across.upper; ++c) {
                mesh.addNode(rangePair.along.lower, c + 1).setEnabled();
                mesh.addNode(rangePair.along.upper, c).setEnabled();
            }
            for (int l = rangePair.along.lower; l < rangePair.along.upper; ++l) {
                mesh.addNode(l,     rangePair.across.lower).setEnabled();
                mesh.addNode(l + 1, rangePair.across.upper).setEnabled();
            }
        }

        return tiles;

    }

    /**
     * Extract a sample of points from a mesh.
     *
     * @param mesh mesh from which grid should be extracted
     * @return extracted grid
     */
    private List<GeodeticPoint> extractSample(final Mesh mesh) {

        // find how to select sample points taking quantization into account
        // to have the largest possible number of points while still
        // being inside the zone of interest
        int selectedAcrossModulus = -1;
        int selectedAlongModulus  = -1;
        int selectedCount         = -1;
        for (int acrossModulus = 0; acrossModulus < quantization; ++acrossModulus) {
            for (int alongModulus = 0; alongModulus < quantization; ++alongModulus) {

                // count how many points would be selected for the current modulus
                int count = 0;
                for (int across = mesh.getMinAcrossIndex() + acrossModulus;
                        across <= mesh.getMaxAcrossIndex();
                        across += quantization) {
                    for (int along = mesh.getMinAlongIndex() + alongModulus;
                            along <= mesh.getMaxAlongIndex();
                            along += quantization) {
                        final Mesh.Node  node = mesh.getNode(along, across);
                        if (node != null && node.isInside()) {
                            ++count;
                        }
                    }
                }

                if (count > selectedCount) {
                    // current modulus are better than the selected ones
                    selectedAcrossModulus = acrossModulus;
                    selectedAlongModulus  = alongModulus;
                    selectedCount         = count;
                }
            }
        }

        // extract the sample points
        final List<GeodeticPoint> sample = new ArrayList<>(selectedCount);
        for (int across = mesh.getMinAcrossIndex() + selectedAcrossModulus;
                across <= mesh.getMaxAcrossIndex();
                across += quantization) {
            for (int along = mesh.getMinAlongIndex() + selectedAlongModulus;
                    along <= mesh.getMaxAlongIndex();
                    along += quantization) {
                final Mesh.Node  node = mesh.getNode(along, across);
                if (node != null && node.isInside()) {
                    sample.add(toGeodetic(node.getS2P()));
                }
            }
        }

        return sample;

    }

    /** Merge two meshes together.
     * @param mesh1 first mesh
     * @param mesh2 second mesh
     * @param mergingSeeds collection where to put the nodes created during the merge
     * @return merged mesh (really one of the instances)
     */
    private Mesh mergeMeshes(final Mesh mesh1, final Mesh mesh2,
                             final Collection<Mesh.Node> mergingSeeds) {

        // select the way merge will be performed
        final Mesh larger;
        final Mesh smaller;
        if (mesh1.getNumberOfNodes() >= mesh2.getNumberOfNodes()) {
            // the larger new mesh should absorb the smaller existing mesh
            larger  = mesh1;
            smaller = mesh2;
        } else {
            // the larger existing mesh should absorb the smaller new mesh
            larger  = mesh2;
            smaller = mesh1;
        }

        // prepare seed nodes for next iteration
        for (final Mesh.Node insideNode : smaller.getInsideNodes()) {

            // beware we cannot reuse the node itself as the two meshes are not aligned!
            // we have to create new nodes around the previous location
            Mesh.Node node = larger.getClosestExistingNode(insideNode.getV());

            while (estimateAlongMotion(node, insideNode.getV()) > +mesh1.getAlongGap()) {
                // the node is before desired index in the along direction
                // we need to create intermediates nodes up to the desired index
                node = larger.addNode(node.getAlongIndex() + 1, node.getAcrossIndex());
            }

            while (estimateAlongMotion(node, insideNode.getV()) < -mesh1.getAlongGap()) {
                // the node is after desired index in the along direction
                // we need to create intermediates nodes up to the desired index
                node = larger.addNode(node.getAlongIndex() - 1, node.getAcrossIndex());
            }

            while (estimateAcrossMotion(node, insideNode.getV()) > +mesh1.getAcrossGap()) {
                // the node is before desired index in the across direction
                // we need to create intermediates nodes up to the desired index
                node = larger.addNode(node.getAlongIndex(), node.getAcrossIndex() + 1);
            }

            while (estimateAcrossMotion(node, insideNode.getV()) < -mesh1.getAcrossGap()) {
                // the node is after desired index in the across direction
                // we need to create intermediates nodes up to the desired index
                node = larger.addNode(node.getAlongIndex(), node.getAcrossIndex() - 1);
            }

            // now we are close to the inside node,
            // make sure the four surrounding nodes are available
            final int otherAlong  = (estimateAlongMotion(node, insideNode.getV()) < 0.0) ?
                                    node.getAlongIndex()  - 1 : node.getAlongIndex() + 1;
            final int otherAcross = (estimateAcrossMotion(node, insideNode.getV()) < 0.0) ?
                                    node.getAcrossIndex()  - 1 : node.getAcrossIndex() + 1;
            addNode(node.getAlongIndex(), node.getAcrossIndex(), larger, mergingSeeds);
            addNode(node.getAlongIndex(), otherAcross,           larger, mergingSeeds);
            addNode(otherAlong,           node.getAcrossIndex(), larger, mergingSeeds);
            addNode(otherAlong,           otherAcross,           larger, mergingSeeds);

        }

        return larger;

    }

    /** Ensure all 8 neighbors of a node are in the mesh.
     * @param base base node
     * @param mesh complete mesh containing nodes
     * @param newNodes queue where new node must be put
     */
    private void addAllNeighborsIfNeeded(final Mesh.Node base, final Mesh mesh,
                                         final Collection<Mesh.Node> newNodes) {
        addNode(base.getAlongIndex() - 1, base.getAcrossIndex() - 1, mesh, newNodes);
        addNode(base.getAlongIndex() - 1, base.getAcrossIndex(),     mesh, newNodes);
        addNode(base.getAlongIndex() - 1, base.getAcrossIndex() + 1, mesh, newNodes);
        addNode(base.getAlongIndex(),     base.getAcrossIndex() - 1, mesh, newNodes);
        addNode(base.getAlongIndex(),     base.getAcrossIndex() + 1, mesh, newNodes);
        addNode(base.getAlongIndex() + 1, base.getAcrossIndex() - 1, mesh, newNodes);
        addNode(base.getAlongIndex() + 1, base.getAcrossIndex(),     mesh, newNodes);
        addNode(base.getAlongIndex() + 1, base.getAcrossIndex() + 1, mesh, newNodes);
    }

    /** Add a node to a mesh if not already present.
     * @param alongIndex index in the along direction
     * @param acrossIndex index in the across direction
     * @param mesh complete mesh containing nodes
     * @param newNodes queue where new node must be put
     */
    private void addNode(final int alongIndex, final int acrossIndex,
                         final Mesh mesh, final Collection<Mesh.Node> newNodes) {

        final Mesh.Node node = mesh.addNode(alongIndex, acrossIndex);

        if (!node.isEnabled()) {
            // enable the node
            node.setEnabled();
            newNodes.add(node);
        }

    }

    /** Convert a point on the unit 2-sphere to geodetic coordinates.
     * @param point point on the unit 2-sphere
     * @return geodetic point (arbitrarily set at altitude 0)
     */
    protected GeodeticPoint toGeodetic(final S2Point point) {
        return new GeodeticPoint(0.5 * FastMath.PI - point.getPhi(), point.getTheta(), 0.0);
    }

    /** Build a simple zone (connected zone without holes).
     * <p>
     * In order to build more complex zones (not connected or with
     * holes), the user should directly call Hipparchus
     * {@link SphericalPolygonsSet} constructors and
     * {@link RegionFactory region factory} if set operations
     * are needed (union, intersection, difference ...).
     * </p>
     * <p>
     * Take care that the vertices boundary points must be given <em>counterclockwise</em>.
     * Using the wrong order defines the complementary of the real zone,
     * and will often result in tessellation failure as the zone is too
     * wide.
     * </p>
     * @param tolerance angular separation below which points are considered
     * equal (typically 1.0e-10)
     * @param points vertices of the boundary, in <em>counterclockwise</em>
     * order, each point being a two-elements arrays with latitude at index 0
     * and longitude at index 1
     * @return a zone defined on the unit 2-sphere
     */
    public static SphericalPolygonsSet buildSimpleZone(final double tolerance,
                                                       final double[]... points) {
        final S2Point[] vertices = new S2Point[points.length];
        for (int i = 0; i < points.length; ++i) {
            vertices[i] = new S2Point(points[i][1], 0.5 * FastMath.PI - points[i][0]);
        }
        return new SphericalPolygonsSet(tolerance, vertices);
    }

    /** Build a simple zone (connected zone without holes).
     * <p>
     * In order to build more complex zones (not connected or with
     * holes), the user should directly call Hipparchus
     * {@link SphericalPolygonsSet} constructors and
     * {@link RegionFactory region factory} if set operations
     * are needed (union, intersection, difference ...).
     * </p>
     * <p>
     * Take care that the vertices boundary points must be given <em>counterclockwise</em>.
     * Using the wrong order defines the complementary of the real zone,
     * and will often result in tessellation failure as the zone is too
     * wide.
     * </p>
     * @param tolerance angular separation below which points are considered
     * equal (typically 1.0e-10)
     * @param points vertices of the boundary, in <em>counterclockwise</em>
     * order
     * @return a zone defined on the unit 2-sphere
     */
    public static SphericalPolygonsSet buildSimpleZone(final double tolerance,
                                                       final GeodeticPoint... points) {
        final S2Point[] vertices = new S2Point[points.length];
        for (int i = 0; i < points.length; ++i) {
            vertices[i] = new S2Point(points[i].getLongitude(),
                                      0.5 * FastMath.PI - points[i].getLatitude());
        }
        return new SphericalPolygonsSet(tolerance, vertices);
    }

    /** Estimate an approximate motion in the along direction.
     * @param start node at start of motion
     * @param end desired point at end of motion
     * @return approximate motion in the along direction
     */
    private double estimateAlongMotion(final Mesh.Node start, final Vector3D end) {
        return Vector3D.dotProduct(start.getAlong(), end.subtract(start.getV()));
    }

    /** Estimate an approximate motion in the across direction.
     * @param start node at start of motion
     * @param end desired point at end of motion
     * @return approximate motion in the across direction
     */
    private double estimateAcrossMotion(final Mesh.Node start, final Vector3D end) {
        return Vector3D.dotProduct(start.getAcross(), end.subtract(start.getV()));
    }

    /** Check if an arc meets the inside of a zone.
     * @param s1 first point
     * @param s2 second point
     * @param zone zone to check arc against
     * @return true if the arc meets the inside of the zone
     */
    private boolean meetInside(final S2Point s1, final S2Point s2,
                               final SphericalPolygonsSet zone) {
        final Circle  circle = new Circle(s1, s2, zone.getTolerance());
        final double alpha1  = circle.toSubSpace(s1).getAlpha();
        final double alpha2  = MathUtils.normalizeAngle(circle.toSubSpace(s2).getAlpha(),
                                                        alpha1 + FastMath.PI);
        final SubCircle sub  = new SubCircle(circle,
                                             new ArcsSet(alpha1, alpha2, zone.getTolerance()));
        return recurseMeetInside(zone.getTree(false), sub);

    }

    /** Check if an arc meets the inside of a zone.
     * <p>
     * This method is heavily based on the Characterization class from
     * Hipparchus library, also distributed under the terms
     * of the Apache Software License V2.
     * </p>
     * @param node spherical zone node
     * @param sub arc to characterize
     * @return true if the arc meets the inside of the zone
     */
    private boolean recurseMeetInside(final BSPTree<Sphere2D> node, final SubHyperplane<Sphere2D> sub) {

        if (node.getCut() == null) {
            // we have reached a leaf node
            if (sub.isEmpty()) {
                return false;
            } else {
                return (Boolean) node.getAttribute();
            }
        } else {
            final Hyperplane<Sphere2D> hyperplane = node.getCut().getHyperplane();
            final SubHyperplane.SplitSubHyperplane<Sphere2D> split = sub.split(hyperplane);
            switch (split.getSide()) {
                case PLUS:
                    return recurseMeetInside(node.getPlus(),  sub);
                case MINUS:
                    return recurseMeetInside(node.getMinus(), sub);
                case BOTH:
                    if (recurseMeetInside(node.getPlus(), split.getPlus())) {
                        return true;
                    } else {
                        return recurseMeetInside(node.getMinus(), split.getMinus());
                    }
                default:
                    // this should not happen
                    throw new OrekitInternalError(null);
            }
        }
    }

    /** Get an iterator over mesh nodes indices.
     * @param minIndex minimum node index
     * @param maxIndex maximum node index
     * @param truncateLast true if we can reduce last tile
     * @return iterator over mesh nodes indices
     */
    private Iterable<Range> nodesIndices(final int minIndex, final int maxIndex, final boolean truncateLast) {

        final int first;
        if (truncateLast) {

            // truncate last tile rather than balance tiles around the zone of interest
            first = minIndex;

        } else {

            // balance tiles around the zone of interest rather than truncate last tile

            // number of tiles needed to cover the full indices range
            final int range      = maxIndex - minIndex;
            final int nbTiles    = (range + quantization - 1) / quantization;

            // extra nodes that must be added to complete the tiles
            final int extraNodes = nbTiles * quantization  - range;

            // balance the extra nodes before min index and after maxIndex
            final int extraBefore = (extraNodes + 1) / 2;

            first = minIndex - extraBefore;

        }

        return new Iterable<Range>() {

            /** {@inheritDoc} */
            @Override
            public Iterator<Range> iterator() {
                return new Iterator<Range>() {

                    private int nextLower = first;

                    /** {@inheritDoc} */
                    @Override
                    public boolean hasNext() {
                        return nextLower < maxIndex;
                    }

                    /** {@inheritDoc} */
                    @Override
                    public Range next() {

                        if (nextLower >= maxIndex) {
                            throw new NoSuchElementException();
                        }
                        final int lower = nextLower;

                        nextLower += quantization;
                        if (truncateLast && nextLower > maxIndex) {
                            // truncate last tile
                            nextLower = maxIndex;
                        }

                        return new Range(lower, nextLower);

                    }

                    /** {@inheritDoc} */
                    @Override
                    public void remove() {
                        throw new UnsupportedOperationException();
                    }

                };
            }
        };

    }

    /** Local class for a range of indices to be used for building a tile. */
    private static class Range {

        /** Lower index. */
        private final int lower;

        /** Upper index. */
        private final int upper;

        /** Simple constructor.
         * @param lower lower index
         * @param upper upper index
         */
        Range(final int lower, final int upper) {
            this.lower = lower;
            this.upper = upper;
        }

    }

    /** Local class for a pair of ranges of indices to be used for building a tile. */
    private static class RangePair {

        /** Across range. */
        private final Range across;

        /** Along range. */
        private final Range along;

        /** Simple constructor.
         * @param across across range
         * @param along along range
         */
        RangePair(final Range across, final Range along) {
            this.across = across;
            this.along  = along;
        }

    }

}

Assistant: This Java code defines a class `EllipsoidTessellator` that is used to tessellate or sample a geographic area on an ellipsoidal body, such as Earth, into smaller, manageable parts. These smaller parts can be either tiles or grids of geodetic points. The code is a part of a library focused on Earth Observation and satellite missions. Let’s break down the primary components through its methods and supporting classes:  
  
### Main Features  
  
1. **Tessellation of Interest Zone**:  
- The `tessellate` method divides a specified zone on the ellipsoid into tiles. These tiles are adjusted to fit the geographic region of interest, using strategies to handle overlaps and to decide whether to truncate the final row or column of tiles to avoid going beyond the zone.  
  
2. **Sampling of Interest Zone**:  
- The `sample` method generates a grid of geodetic points entirely contained within the specified zone. It allows generating a point-based representation of the zone, which is useful for certain Earth observation applications.  
  
3. **Zone Building**:  
- The `buildSimpleZone` methods are static utility functions for constructing simple geographic zones from latitude and longitude points. This enables users to easily define the area they wish to tessellate or sample.  
  
4. **Mesh Management**:  
- The `Mesh` class (not fully shown in this code but referred to) is responsible for managing nodes in the tessellator that correspond to points on the ellipsoid’s surface. Methods like `neighborExpandMesh`, `addAllNeighborsIfNeeded`, and `addNode` help ensure proper mesh creation and management, ensuring all necessary connections and expansions meet the desired zone’s boundaries.  
  
5. **Handling of Indices and Nodes**:  
- The logic around indices (handled in private methods like `nodesIndices`) makes sure that the divisions (tiles or grids) are properly aligned and overlap or truncate as necessary. This ensures efficient geographic representation while minimizing errors due to ellipsoidal distortion.  
  
6. **Geodetic Translations**:  
- The methods `toGeodetic` and related calculations handle translating between the mathematical representation of a point on the ellipsoid and real-world geographic coordinates.  
  
### Error Handling and Limits  
  
- **MAX_ITER**: The constant `MAX_ITER` is a limit to avoid infinite loops in computations due to numerical issues. This acts as a sanity check during mesh expansions and merges.  
  
### Libraries Used  
  
- **Hipparchus and Orekit**: The code imports classes from the Hipparchus and Orekit libraries, known for advanced mathematical and astrodynamics functionalities, respectively. These libraries assist in geometric and trigonometric computations as well as error management.  
  
### Code Documentation and Licensing  
  
- **Apache License**: The code is licensed under the Apache License 2.0, allowing for broad use and distribution under certain conditions.  
- **Documentation**: The code is well-commented, explaining the purpose of major classes and methods, parameters, and describing complex concepts which might not be immediately apparent.  
  
Overall, the `EllipsoidTessellator` is a sophisticated part of geospatial software infrastructure, crucial for managing Earth surface representations in applications where precise geographic coverage is required, such as satellite-based Earth observation and geographic data analysis.
User: Tell me more about how the mesh class works
Assistant: