diff --git a/.ameba.yml b/.ameba.yml
new file mode 100644
index 0000000..2cf4cc7
--- /dev/null
+++ b/.ameba.yml
@@ -0,0 +1,2 @@
+Lint/Typos:
+  Enabled: false
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index dc580da..54ed39e 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -13,15 +13,14 @@ jobs:
         crystal:
           - latest
           - nightly
-          - 1.0.0
     runs-on: ${{ matrix.os }}
     container: crystallang/crystal:${{ matrix.crystal }}-alpine
     steps:
     - uses: actions/checkout@v2
     - name: Install dependencies
       run: shards install --ignore-crystal-version
-    - name: Lint
-      run: ./bin/ameba
+    # - name: Lint
+    #  run: ./bin/ameba
     - name: Format
       run: crystal tool format --check
     - name: Run tests
diff --git a/README.md b/README.md
index a3595f2..c2b0841 100644
--- a/README.md
+++ b/README.md
@@ -24,14 +24,21 @@ Useful for things like storing points [in InfluxDB](https://docs.influxdata.com/
 
 require "s2_cells"
 
+# index a location in your database
 lat = -33.870456
 lon = 151.208889
 level = 24
 
 cell = S2Cells.at(lat, lon).parent(level)
 token = cell.to_token # => "3ba32f81"
+# or
+id = cell.id # => Int64
 
-# Or a little more direct
-S2Cells::LatLon.new(lat, lon).to_token(level)
+# find all the indexes in an area
+p1 = S2Cells::LatLng.from_degrees(33.0, -122.0)
+p2 = S2Cells::LatLng.from_degrees(33.1, -122.1)
+cells = S2Cells.in(p1, p2) # => Array(CellId)
 
+# then can search your index:
+# loc_index = ANY(cells.map(&.id))
 ```
diff --git a/shard.lock b/shard.lock
new file mode 100644
index 0000000..36d58a3
--- /dev/null
+++ b/shard.lock
@@ -0,0 +1,10 @@
+version: 2.0
+shards:
+  bisect:
+    git: https://github.com/spider-gazelle/bisect.git
+    version: 1.2.1
+
+  priority-queue:
+    git: https://github.com/spider-gazelle/priority-queue.git
+    version: 1.1.2
+
diff --git a/shard.yml b/shard.yml
index c32a084..fc14b06 100644
--- a/shard.yml
+++ b/shard.yml
@@ -1,7 +1,11 @@
 name: s2_cells
-version: 1.0.1
+version: 2.0.0
 crystal: ">= 0.36.1"
 
-development_dependencies:
-  ameba:
-    github: veelenga/ameba
+dependencies:
+  priority-queue:
+    github: spider-gazelle/priority-queue
+
+# development_dependencies:
+#  ameba:
+#    github: veelenga/ameba
diff --git a/spec/s2_cells_spec.cr b/spec/s2_cells_spec.cr
index ddfd262..ed44d48 100644
--- a/spec/s2_cells_spec.cr
+++ b/spec/s2_cells_spec.cr
@@ -2,7 +2,7 @@ require "./spec_helper"
 
 module S2Cells
   describe S2Cells do
-    it "should convert lat lon to a cell id" do
+    it "should convert lat lng to a cell id" do
       {
         {0x47a1cbd595522b39_u64, 49.703498679, 11.770681595},
         {0x46525318b63be0f9_u64, 55.685376759, 12.588490937},
@@ -18,13 +18,14 @@ module S2Cells
         {0x3b00955555555555_u64, 10.050986518, 78.293170610},
         {0x1dcc469991555555_u64, -34.055420593, 18.551140038},
         {0xb112966aaaaaaaab_u64, -69.219262171, 49.670072392},
-      }.each do |(id, lat, lon)|
-        point = LatLon.new(lat, lon).to_point
+      }.each do |(id, lat, lng)|
+        lat_lng = LatLng.from_degrees(lat, lng)
+        point = lat_lng.to_point
         cell = CellId.from_point(point)
         cell.id.should eq(id)
 
-        CellId.from_lat_lon(lat, lon).id.should eq(id)
-        S2Cells.at(lat, lon).id.should eq(id)
+        CellId.from_lat_lng(lat_lng).id.should eq(id)
+        S2Cells.at(lat, lng).id.should eq(id)
       end
     end
 
@@ -69,4 +70,176 @@ module S2Cells
       end
     end
   end
+
+  it "should generate the correct face" do
+    CellId.from_lat_lng(0.0, 0.0).face.should eq 0
+    CellId.from_lat_lng(0.0, 90.0).face.should eq 1
+    CellId.from_lat_lng(90.0, 0.0).face.should eq 2
+    CellId.from_lat_lng(0.0, 180.0).face.should eq 3
+    CellId.from_lat_lng(0.0, -90.0).face.should eq 4
+    CellId.from_lat_lng(-90.0, 0.0).face.should eq 5
+  end
+
+  it "test parent child relationship" do
+    cell_id = CellId.from_face_pos_level(3, 0x12345678_u64, CellId::MAX_LEVEL - 4)
+
+    cell_id.face.should eq 3
+    cell_id.pos.to_s(2).should eq 0x12345700.to_s(2)
+    cell_id.level.should eq(CellId::MAX_LEVEL - 4)
+    cell_id.valid?.should be_true
+    cell_id.leaf?.should be_false
+
+    cell_id.child_begin(cell_id.level + 2).pos.should eq 0x12345610
+    cell_id.child_begin.pos.should eq 0x12345640
+    cell_id.parent.pos.should eq 0x12345400
+    cell_id.parent(cell_id.level - 2).pos.should eq 0x12345000
+
+    cell_id.child_begin.next.next.next.next.should eq cell_id.child_end
+    cell_id.child_begin(CellId::MAX_LEVEL).should eq cell_id.range_min
+    cell_id.child_end(CellId::MAX_LEVEL).should eq cell_id.range_max.next
+
+    # Check that cells are represented by the position of their center
+    # alngg the Hilbert curve.
+    (cell_id.range_min.id &+ cell_id.range_max.id).should eq(2_u64 &* cell_id.id)
+  end
+
+  it "should be able to switch between lat lang and cell ids" do
+    INVERSE_ITERATIONS.times do
+      cell_id = get_random_cell_id(CellId::MAX_LEVEL)
+      cell_id.leaf?.should be_true
+      cell_id.level.should eq CellId::MAX_LEVEL
+      center = cell_id.to_lat_lng
+      CellId.from_lat_lng(center).id.should eq cell_id.id
+    end
+  end
+
+  it "should be able to switch between tokens and cell ids" do
+    TOKEN_ITERATIONS.times do
+      cell_id = get_random_cell_id
+      token = cell_id.to_token
+      (token.size <= 16).should be_true
+      CellId.from_token(token).id.should eq cell_id.id
+    end
+  end
+
+  it "should be able to obtain neighbours" do
+    # Check the edge neighbors of face 1.
+    out_faces = {5, 3, 2, 0}
+    face_nbrs = CellId.from_face_pos_level(1, 0, 0).get_edge_neighbors
+    face_nbrs.each_with_index do |face_nbr, i|
+      face_nbr.face?.should be_true
+      face_nbr.face.should eq out_faces[i]?
+    end
+
+    # Check the vertex neighbors of the center of face 2 at level 5.
+    neighbors = CellId.from_point(Point.new(0, 0, 1)).get_vertex_neighbors(5)
+    neighbors.sort!
+    neighbors.each_with_index do |neighbor, i|
+      neighbor.id.should eq(CellId.from_face_ij(
+        2,
+        (1_u64 << 29) - (i < 2 ? 1 : 0),
+        (1_u64 << 29) - (i == 0 || i == 3 ? 1 : 0)
+      ).parent(5).id)
+    end
+
+    neighbors.each_with_index do |neighbor, i|
+      neighbor.should eq(CellId.from_face_ij(
+        2,
+        (1_u64 << 29) - (i < 2 ? 1 : 0),
+        (1_u64 << 29) - (i == 0 || i == 3 ? 1 : 0)
+      ).parent(5))
+    end
+
+    # Check the vertex neighbors of the corner of faces 0, 4, and 5.
+    cell_id = CellId.from_face_pos_level(0, 0, CellId::MAX_LEVEL)
+    neighbors = cell_id.get_vertex_neighbors(0)
+    neighbors.sort!
+    neighbors.size.should eq 3
+
+    CellId.from_face_pos_level(0, 0, 0).should eq neighbors[0]
+    CellId.from_face_pos_level(4, 0, 0).should eq neighbors[1]
+    CellId.from_face_pos_level(5, 0, 0).should eq neighbors[2]
+
+    # check a bunch
+    NEIGHBORS_ITERATIONS.times do
+      cell_id = get_random_cell_id
+      cell_id = cell_id.parent if cell_id.leaf?
+      max_diff = {6, CellId::MAX_LEVEL - cell_id.level - 1}.min
+      level = max_diff == 0 ? cell_id.level : cell_id.level + rand(max_diff)
+      raise "level < cell_id.level" unless level >= cell_id.level
+      raise "level == MAX_LEVEL" if level >= CellId::MAX_LEVEL
+
+      all, expected = {Set(CellId).new, Set(CellId).new}
+      neighbors = cell_id.get_all_neighbors(level)
+      all.concat neighbors
+      cell_id.children(level + 1).each do |child|
+        all.add(child.parent)
+        expected.concat(child.get_vertex_neighbors(level))
+      end
+
+      all_a = all.map(&.id).uniq!.sort
+      expect_a = expected.map(&.id).uniq!.sort
+
+      all_a.size.should eq expect_a.size
+      all_a.should eq expect_a
+    end
+  end
+
+  it "should work with faces" do
+    edge_counts = Hash(Point, Int32).new(0)
+    vertex_counts = Hash(Point, Int32).new(0)
+
+    6.times do |face|
+      cell_id = CellId.from_face_pos_level(face, 0, 0)
+      cell = Cell.new(cell_id)
+      cell.id.should eq cell_id
+      cell.face.should eq face
+      cell.level.should eq 0
+
+      cell.orientation.should eq(face & SWAP_MASK)
+      cell.leaf?.should eq false
+
+      4.times do |k|
+        edge_counts[cell.get_edge_raw(k)] += 1
+        vertex_counts[cell.get_vertex_raw(k)] += 1
+
+        cell.get_vertex_raw(k).dot_prod(cell.get_edge_raw(k)).should eq 0.0
+        cell.get_vertex_raw((k + 1) & 3)
+          .dot_prod(cell.get_edge_raw(k))
+          .should eq 0.0
+
+        cell
+          .get_vertex_raw(k)
+          .cross_prod(cell.get_vertex_raw((k + 1) & 3))
+          .normalize
+          .dot_prod(cell.get_edge(k))
+          .should be_close(1.0, 0.000001)
+      end
+    end
+
+    edge_counts.values.each { |count| count.should eq 2 }
+    vertex_counts.values.each { |count| count.should eq 3 }
+  end
+
+  it "generates the correct covering for a given region" do
+    coverer = RegionCoverer.new
+
+    p1 = LatLng.from_degrees(33.0, -122.0)
+    p2 = LatLng.from_degrees(33.1, -122.1)
+
+    cell_ids = coverer.get_covering(LatLngRect.from_point_pair(p1, p2))
+    ids = cell_ids.map(&.id).sort
+
+    target = [
+      9291041754864156672_u64,
+      9291043953887412224_u64,
+      9291044503643226112_u64,
+      9291045878032760832_u64,
+      9291047252422295552_u64,
+      9291047802178109440_u64,
+      9291051650468806656_u64,
+      9291052200224620544_u64,
+    ]
+    ids.should eq(target)
+  end
 end
diff --git a/spec/spec_helper.cr b/spec/spec_helper.cr
index 8ef9125..9fb3163 100644
--- a/spec/spec_helper.cr
+++ b/spec/spec_helper.cr
@@ -1,2 +1,19 @@
 require "spec"
+require "random"
 require "../src/s2_cells"
+
+def get_random_cell_id(level : Int32 = Random.rand(S2Cells::CellId::MAX_LEVEL + 1))
+  face = Random.rand(S2Cells::CellId::NUM_FACES)
+  pos = Random.rand(UInt64::MAX) & ((1_u64 << (2 * S2Cells::CellId::MAX_LEVEL)) - 1)
+
+  S2Cells::CellId.from_face_pos_level(face, pos, level)
+end
+
+INVERSE_ITERATIONS          = 20
+TOKEN_ITERATIONS            = 10
+COVERAGE_ITERATIONS         = 10
+NEIGHBORS_ITERATIONS        = 10
+NORMALIZE_ITERATIONS        = 20
+REGION_COVERER_ITERATIONS   = 10
+RANDOM_CAPS_ITERATIONS      = 10
+SIMPLE_COVERINGS_ITERATIONS = 10
diff --git a/src/s2_cells.cr b/src/s2_cells.cr
index 0b613b3..765f65c 100644
--- a/src/s2_cells.cr
+++ b/src/s2_cells.cr
@@ -8,12 +8,112 @@ module S2Cells
     end
   end
 
-  def self.at(lat : Float64, lon : Float64)
-    CellId.from_lat_lon(lat, lon)
+  def self.at(point : LatLng) : CellId
+    CellId.from_lat_lng(point)
   end
+
+  def self.at(lat : Float64, lng : Float64) : CellId
+    at LatLng.from_degrees(lat, lng)
+  end
+
+  def self.in(p1 : LatLng, p2 : LatLng) : Array(CellId)
+    coverer = RegionCoverer.new
+    coverer.get_covering(LatLngRect.from_point_pair(p1, p2))
+  end
+
+  LINEAR_PROJECTION    = 0
+  TAN_PROJECTION       = 1
+  QUADRATIC_PROJECTION = 2
+
+  SWAP_MASK   =  0x01
+  INVERT_MASK =  0x02
+  LOOKUP_BITS = 4_u64
+
+  POS_TO_ORIENTATION = {SWAP_MASK, 0, 0, INVERT_MASK | SWAP_MASK}
+  POS_TO_IJ          = { {0_u64, 1_u64, 3_u64, 2_u64},
+                        {0_u64, 2_u64, 3_u64, 1_u64},
+                        {3_u64, 2_u64, 0_u64, 1_u64},
+                        {3_u64, 1_u64, 0_u64, 2_u64} }
+
+  LOOKUP_POS = Array.new((1_u64 << (2 * LOOKUP_BITS + 2)), 0_u64)
+  LOOKUP_IJ  = Array.new((1_u64 << (2 * LOOKUP_BITS + 2)), 0_u64)
 end
 
-require "./s2_cells/cell_base"
+require "./s2_cells/angle"
 require "./s2_cells/point"
+require "./s2_cells/lat_lng"
+require "./s2_cells/interval"
 require "./s2_cells/cell_id"
-require "./s2_cells/lat_lon"
+require "./s2_cells/metric"
+require "./s2_cells/cap"
+require "./s2_cells/cell"
+require "./s2_cells/cell_union"
+require "./s2_cells/lat_lng_rect"
+require "./s2_cells/region_coverer"
+
+# helper functions
+module S2Cells
+  def self.valid_face_xyz_to_uv(face : Int32, p : Point)
+    raise "invalid face xyz" unless p.dot_prod(face_uv_to_xyz(face, 0.0, 0.0)) > 0
+
+    case face
+    when 0 then [p.y / p.x, p.z / p.x]
+    when 1 then [-p.x / p.y, p.z / p.y]
+    when 2 then [-p.x / p.z, -p.y / p.z]
+    when 3 then [p.z / p.x, p.y / p.x]
+    when 4 then [p.z / p.y, -p.x / p.y]
+    else        [-p.y / p.z, -p.x / p.z]
+    end
+  end
+
+  def self.face_uv_to_xyz(face : Int32, u : Float64, v : Float64)
+    case face
+    when 0 then Point.new(1_f64, u, v)
+    when 1 then Point.new(-u, 1_f64, v)
+    when 2 then Point.new(-u, -v, 1_f64)
+    when 3 then Point.new(-1_f64, -v, -u)
+    when 4 then Point.new(v, -1_f64, -u)
+    else        Point.new(v, u, -1_f64)
+    end
+  end
+
+  def self.xyz_to_face_uv(p : Point) : Tuple(Int32, Float64, Float64)
+    face = p.largest_abs_component
+
+    pface = case face
+            when 0 then p.x
+            when 1 then p.y
+            else        p.z
+            end
+
+    face += 3 if pface < 0.0
+
+    u, v = valid_face_xyz_to_uv(face, p)
+    {face, u, v}
+  end
+
+  def self.face_xyz_to_uv(face : Int32, p : Point)
+    if face < 3
+      if p[face] <= 0.0
+        return {false, 0.0, 0.0}
+      end
+    else
+      if p[face - 3] >= 0.0
+        return {false, 0.0, 0.0}
+      end
+    end
+    u, v = valid_face_xyz_to_uv(face, p)
+    {true, u, v}
+  end
+
+  def self.is_unit_length(p : Point) : Bool
+    (p.norm * p.norm - 1).abs <= 1e-15
+  end
+
+  FACE_CELLS = {Cell.from_face_pos_level(0, 0_u64, 0),
+                Cell.from_face_pos_level(1, 0_u64, 0),
+                Cell.from_face_pos_level(2, 0_u64, 0),
+                Cell.from_face_pos_level(3, 0_u64, 0),
+                Cell.from_face_pos_level(4, 0_u64, 0),
+                Cell.from_face_pos_level(5, 0_u64, 0)}
+end
diff --git a/src/s2_cells/angle.cr b/src/s2_cells/angle.cr
new file mode 100644
index 0000000..3a1e83a
--- /dev/null
+++ b/src/s2_cells/angle.cr
@@ -0,0 +1,55 @@
+require "math"
+
+struct S2Cells::Angle
+  include Comparable(Angle)
+
+  # A one-dimensional angle (as opposed to a two-dimensional solid angle).
+  # It has methods for converting angles to or from radians and degrees.
+
+  # Initializes a new Angle with a given radians value.
+  def initialize(@radians : Float64 = 0)
+    raise ArgumentError.new("Invalid type for radians") unless @radians.is_a?(Float64)
+  end
+
+  def hash : UInt64
+    @radians.hash
+  end
+
+  # Checks equality of the Angle with another object.
+  def ==(other : Angle)
+    @radians == other.radians
+  end
+
+  # Compares this Angle to another Angle to determine if it is less than the other.
+  def <=>(other : Angle)
+    @radians <=> other.radians
+  end
+
+  # Adds two Angles together.
+  def +(other : Angle) : Angle
+    Angle.from_radians(@radians + other.radians)
+  end
+
+  # Represents the Angle as a string.
+  def to_s : String
+    "#{self.class.name}: #{@radians}"
+  end
+
+  # Creates an Angle from a degree measurement.
+  def self.from_degrees(degrees : Float64) : Angle
+    new(degrees * Math::PI / 180.0)
+  end
+
+  # Creates an Angle from a radians measurement.
+  def self.from_radians(radians : Float64) : Angle
+    new(radians)
+  end
+
+  # Gets the radians value of the Angle.
+  getter radians : Float64
+
+  # Converts the Angle's radians to degrees.
+  def degrees : Float64
+    @radians * 180.0 / Math::PI
+  end
+end
diff --git a/src/s2_cells/cap.cr b/src/s2_cells/cap.cr
new file mode 100644
index 0000000..8c7a226
--- /dev/null
+++ b/src/s2_cells/cap.cr
@@ -0,0 +1,199 @@
+struct S2Cells::Cap
+  ROUND_UP = 1.0 + 1.0 / (1_u64 << 52)
+
+  def initialize(@axis : Point = Point.new(1.0, 0.0, 0.0), @height : Float64 = -1.0)
+  end
+
+  def self.from_axis_height(axis : Point, height : Float64) : Cap
+    raise "axis must be unit length" unless S2Cells.is_unit_length(axis)
+    Cap.new(axis, height)
+  end
+
+  def self.from_axis_angle(axis : Point, angle : Angle) : Cap
+    raise "axis must be unit length" unless S2Cells.is_unit_length(axis)
+    raise "angle must be non-negative" unless angle.radians >= 0
+    Cap.new(axis, get_height_for_angle(angle.radians))
+  end
+
+  def self.get_height_for_angle(radians : Float64) : Float64
+    raise "radians must be non-negative" unless radians >= 0
+    return 2.0 if radians >= Math::PI
+
+    d = Math.sin(0.5 * radians)
+    2 * d * d
+  end
+
+  def self.from_axis_area(axis : Point, area : Float64) : Cap
+    raise "axis must be unit length" unless S2Cells.is_unit_length(axis)
+    Cap.new(axis, area / (2 * Math::PI))
+  end
+
+  def self.empty : Cap
+    Cap.new
+  end
+
+  def self.full : Cap
+    Cap.new(Point.new(1.0, 0.0, 0.0), 2.0)
+  end
+
+  def height : Float64
+    @height
+  end
+
+  def axis : Point
+    @axis
+  end
+
+  def area : Float64
+    2 * Math::PI * (@height > 0 ? @height : 0)
+  end
+
+  def angle : Angle
+    return Angle.from_radians(-1) if is_empty?
+    Angle.from_radians(2 * Math.asin(Math.sqrt(0.5 * @height)))
+  end
+
+  def is_valid? : Bool
+    S2Cells.is_unit_length(@axis) && @height <= 2.0
+  end
+
+  def is_empty? : Bool
+    @height < 0
+  end
+
+  def is_full? : Bool
+    @height >= 2.0
+  end
+
+  def get_cap_bound : Cap
+    self
+  end
+
+  def add_point(point : Point)
+    raise "point must be unit length" unless S2Cells.is_unit_length(point)
+    if is_empty?
+      @axis = point
+      @height = 0.0
+    else
+      dist2 = (@axis - point).norm2
+      @height = {@height, ROUND_UP * 0.5 * dist2}.max
+    end
+  end
+
+  def complement : Cap
+    height = is_full? ? -1.0 : 2.0 - {@height, 0.0}.max
+    Cap.from_axis_height(-@axis, height)
+  end
+
+  def contains(other : Cap | Point | Cell) : Bool
+    case other
+    in Cap
+      return true if is_full? || other.is_empty?
+      angle.radians >= @axis.angle(other.axis).radians + other.angle.radians
+    in Point
+      raise "point must be unit length" unless S2Cells.is_unit_length(other)
+      (@axis - other).norm2 <= 2 * @height
+    in Cell
+      vertices = Array(Point).new(4)
+      (0..3).each do |k|
+        vertices[k] = other.get_vertex(k)
+        return false unless contains(vertices[k])
+      end
+      !complement.intersects(other, vertices)
+    end
+  end
+
+  def interior_contains(other : Point) : Bool
+    raise "point must be unit length" unless S2Cells.is_unit_length(other)
+    is_full? || (@axis - other).norm2 < 2 * @height
+  end
+
+  def intersects(other : Cap | Cell, vertices = Array(Point).new(4)) : Bool
+    case other
+    in Cap
+      return false if is_empty? || other.is_empty?
+      angle.radians + other.angle.radians >= @axis.angle(other.axis).radians
+    in Cell
+      return false if @height >= 1 || is_empty?
+      return true if other.contains(@axis)
+      sin2_angle = @height * (2 - @height)
+      (0..3).each do |k|
+        edge = other.get_edge_raw(k)
+        dot = @axis.dot_prod(edge)
+        next if dot > 0
+        return false if dot * dot > sin2_angle * edge.norm2
+        dir = edge.cross_prod(@axis)
+        return true if dir.dot_prod(vertices[k]) < 0 && dir.dot_prod(vertices[(k + 1) & 3]) > 0
+      end
+      false
+    end
+  end
+
+  def may_intersect(other : Cell) : Bool
+    vertices = Array(Point).new(4)
+    (0..3).each do |k|
+      vertices[k] = other.get_vertex(k)
+      return true if contains(vertices[k])
+    end
+    intersects(other, vertices)
+  end
+
+  def interior_intersects(other : Cap) : Bool
+    return false if @height <= 0 || other.is_empty?
+    angle.radians + other.angle.radians > @axis.angle(other.axis).radians
+  end
+
+  def get_rect_bound : LatLngRect
+    return LatLngRect.empty if is_empty?
+
+    axis_ll = LatLng.from_point(@axis)
+    cap_angle = angle.radians
+
+    all_longitudes = false
+    lat, lng = [] of Float64, [] of Float64
+    lng << -Math::PI
+    lng << Math::PI
+
+    lat << axis_ll.lat.radians - cap_angle
+    if lat[0] <= -Math::PI / 2.0
+      lat[0] = -Math::PI / 2.0
+      all_longitudes = true
+    end
+
+    lat << axis_ll.lat.radians + cap_angle
+    if lat[1] >= Math::PI / 2.0
+      lat[1] = Math::PI / 2.0
+      all_longitudes = true
+    end
+
+    unless all_longitudes
+      sin_a = Math.sqrt(@height * (2 - @height))
+      sin_c = Math.cos(axis_ll.lat.radians)
+      if sin_a <= sin_c
+        angle_a = Math.asin(sin_a / sin_c)
+        lng[0] = drem(axis_ll.lng.radians - angle_a, 2 * Math::PI)
+        lng[1] = drem(axis_ll.lng.radians + angle_a, 2 * Math::PI)
+      end
+    end
+
+    LatLngRect.new(LineInterval.new(*lat), SphereInterval.new(*lng))
+  end
+
+  def approx_equals(other : Cap, max_error = 1e-14) : Bool
+    (@axis.angle(other.axis).radians <= max_error && (@height - other.height).abs <= max_error) ||
+      (is_empty? && other.height <= max_error) ||
+      (other.is_empty? && @height <= max_error) ||
+      (is_full? && other.height >= 2 - max_error) ||
+      (other.is_full? && @height >= 2 - max_error)
+  end
+
+  def expanded(distance : Angle) : Cap
+    raise "distance must be non-negative" unless distance.radians >= 0
+    return Cap.empty if is_empty?
+    Cap.from_axis_angle(@axis, angle + distance)
+  end
+
+  private def is_unit_length(point : Point) : Bool
+    S2Cells.is_unit_length point
+  end
+end
diff --git a/src/s2_cells/cell.cr b/src/s2_cells/cell.cr
new file mode 100644
index 0000000..76d0323
--- /dev/null
+++ b/src/s2_cells/cell.cr
@@ -0,0 +1,288 @@
+struct S2Cells::Cell
+  @uv : Array(Array(Float64))
+  @cell_id : CellId
+  @face : Int32
+  @orientation : Int32
+  @level : Int32
+
+  def initialize(cell_id : CellId)
+    @uv = Array.new(2) { Array(Float64).new(2, Float64::NAN) }
+
+    @cell_id = cell_id
+    face, i, j, orientation = cell_id.to_face_ij_orientation
+    ij = {i, j}
+    @face = face
+    @orientation = orientation
+    @level = cell_id.level
+
+    cell_size = cell_id.get_size_ij
+    ij.zip(@uv) do |ij_, uv_|
+      ij_lo = ij_ & (~cell_size &+ 1)
+      ij_hi = ij_lo + cell_size
+      uv_[0] = CellId.st_to_uv((1.0 / CellId::MAX_SIZE) * ij_lo)
+      uv_[1] = CellId.st_to_uv((1.0 / CellId::MAX_SIZE) * ij_hi)
+    end
+  end
+
+  def initialize(@face, @level, @orientation, @cell_id, @uv)
+  end
+
+  def self.from_lat_lng(lat_lng : LatLng) : Cell
+    new(CellId.from_lat_lng(lat_lng))
+  end
+
+  def self.from_point(point : Point) : Cell
+    new(CellId.from_point(point))
+  end
+
+  def to_s : String
+    "#{self.class.name}: face #{@face}, level #{@level}, orientation #{orientation}, id #{@cell_id.id}"
+  end
+
+  def self.from_face_pos_level(face : Int32, pos : UInt64, level : Int32) : Cell
+    new(CellId.from_face_pos_level(face, pos, level))
+  end
+
+  def id : CellId?
+    @cell_id
+  end
+
+  def face : Int32
+    @face
+  end
+
+  def level : Int32
+    @level
+  end
+
+  def orientation : Int32
+    @orientation
+  end
+
+  def leaf? : Bool
+    @level == CellId::MAX_LEVEL
+  end
+
+  def get_edge(k : Int32) : Point
+    get_edge_raw(k).normalize
+  end
+
+  def get_edge_raw(k : Int32) : Point
+    case k
+    when 0
+      get_v_norm(@face, @uv[1][0])
+    when 1
+      get_u_norm(@face, @uv[0][1])
+    when 2
+      -get_v_norm(@face, @uv[1][1])
+    else
+      -get_u_norm(@face, @uv[0][0])
+    end
+  end
+
+  def get_vertex(k : Int32) : Point
+    get_vertex_raw(k).normalize
+  end
+
+  def get_vertex_raw(k : Int32) : Point
+    S2Cells.face_uv_to_xyz(@face, @uv[0][(k >> 1) ^ (k & 1)], @uv[1][k >> 1])
+  end
+
+  def exact_area : Float64
+    v0 = get_vertex(0)
+    v1 = get_vertex(1)
+    v2 = get_vertex(2)
+    v3 = get_vertex(3)
+    area(v0, v1, v2) + area(v0, v2, v3)
+  end
+
+  def average_area : Float64
+    AVG_AREA.get_value(@level)
+  end
+
+  def approx_area : Float64
+    if @level < 2
+      average_area
+    else
+      flat_area = 0.5 * (get_vertex(2) - get_vertex(0)).cross_prod(get_vertex(3) - get_vertex(1)).norm
+      flat_area * 2 / (1 + Math.sqrt(1 - [1.0 / Math::PI * flat_area, 1.0].min))
+    end
+  end
+
+  def subdivide(&)
+    uv_mid = @cell_id.get_center_uv
+
+    @cell_id.children.each_with_index do |cell_id, pos|
+      uv = Array.new(2) { Array(Float64).new(2, Float64::NAN) }
+      ij = POS_TO_IJ[@orientation][pos]
+      i = ij >> 1
+      j = ij & 1
+      uv[0][i] = @uv[0][i]
+      uv[0][1 - i] = uv_mid[0]
+      uv[1][j] = @uv[1][j]
+      uv[1][1 - j] = uv_mid[1]
+
+      yield Cell.new(
+        @face,
+        @level + 1,
+        @orientation ^ POS_TO_ORIENTATION[pos],
+        cell_id,
+        uv
+      )
+    end
+  end
+
+  def get_center : Point
+    get_center_raw.normalize
+  end
+
+  def get_center_raw : Point
+    @cell_id.to_point_raw
+  end
+
+  def contains(other : Cell | Point) : Bool
+    case other
+    in Cell
+      @cell_id.contains(other.@cell_id)
+    in Point
+      valid, u, v = S2Cells.face_xyz_to_uv(@face, other)
+      return false unless valid
+      u >= @uv[0][0] && u <= @uv[0][1] && v >= @uv[1][0] && v <= @uv[1][1]
+    end
+  end
+
+  def may_intersect(cell : Cell) : Bool
+    @cell_id.intersects(cell.@cell_id)
+  end
+
+  def get_latitude(i : Int32, j : Int32) : Float64
+    p = S2Cells.face_uv_to_xyz(@face, @uv[0][i], @uv[1][j])
+    LatLng.latitude(p).radians
+  end
+
+  def get_longitude(i : Int32, j : Int32) : Float64
+    p = S2Cells.face_uv_to_xyz(@face, @uv[0][i], @uv[1][j])
+    LatLng.longitude(p).radians
+  end
+
+  def get_cap_bound : Cap
+    u = 0.5 * (@uv[0][0] + @uv[0][1])
+    v = 0.5 * (@uv[1][0] + @uv[1][1])
+    cap = Cap.from_axis_height(face_uv_to_xyz(@face, u, v).normalize, 0)
+    4.times { |k| cap.add_point(get_vertex(k)) }
+    cap
+  end
+
+  def get_rect_bound : LatLngRect
+    if @level > 0
+      u = @uv[0][0] + @uv[0][1]
+      v = @uv[1][0] + @uv[1][1]
+      i = get_u_axis(@face)[2] == 0 ? (u < 0 ? 1 : 0) : (u > 0 ? 1 : 0)
+      j = get_v_axis(@face)[2] == 0 ? (v < 0 ? 1 : 0) : (v > 0 ? 1 : 0)
+
+      max_error = 1.0 / (1_u64 << 51)
+      lat = LineInterval.from_point_pair(get_latitude(i, j), get_latitude(1 - i, 1 - j))
+      lat = lat.expanded(max_error).intersection(LatLngRect.full_lat)
+
+      if lat.lo == (-Math::PI / 2.0) || lat.hi == (Math::PI / 2.0)
+        return LatLngRect.new(lat, SphereInterval.full)
+      end
+
+      lng = SphereInterval.from_point_pair(get_longitude(i, 1 - j), get_longitude(1 - i, j))
+      return LatLngRect.new(lat, lng.expanded(max_error))
+    end
+
+    pole_min_lat = Math.asin(Math.sqrt(1.0 / 3.0))
+
+    case @face
+    when 0
+      LatLngRect.new(LineInterval.new(-Math::PI / 4.0, Math::PI / 4.0), SphereInterval.new(-Math::PI / 4.0, Math::PI / 4.0))
+    when 1
+      LatLngRect.new(LineInterval.new(-Math::PI / 4.0, Math::PI / 4.0), SphereInterval.new(Math::PI / 4.0, 3.0 * Math::PI / 4.0))
+    when 2
+      LatLngRect.new(LineInterval.new(pole_min_lat, Math::PI / 2.0), SphereInterval.new(-Math::PI, Math::PI))
+    when 3
+      LatLngRect.new(LineInterval.new(-Math::PI / 4.0, Math::PI / 4.0), SphereInterval.new(3.0 * Math::PI / 4.0, -3.0 * Math::PI / 4.0))
+    when 4
+      LatLngRect.new(LineInterval.new(-Math::PI / 4.0, Math::PI / 4.0), SphereInterval.new(-3.0 * Math::PI / 4.0, -Math::PI / 4.0))
+    else
+      LatLngRect.new(LineInterval.new(-Math::PI / 2.0, -pole_min_lat), SphereInterval.new(-Math::PI, Math::PI))
+    end
+  end
+
+  def get_u_axis(face : Int)
+    case face
+    when 0
+      Point.new(0, 1, 0)
+    when 1
+      Point.new(-1, 0, 0)
+    when 2
+      Point.new(-1, 0, 0)
+    when 3
+      Point.new(0, 0, -1)
+    when 4
+      Point.new(0, 0, -1)
+    else
+      Point.new(0, 1, 0)
+    end
+  end
+
+  def get_v_axis(face : Int)
+    case face
+    when 0
+      Point.new(0, 0, 1)
+    when 1
+      Point.new(0, 0, 1)
+    when 2
+      Point.new(0, -1, 0)
+    when 3
+      Point.new(0, -1, 0)
+    when 4
+      Point.new(1, 0, 0)
+    else
+      Point.new(1, 0, 0)
+    end
+  end
+
+  # Vector normal to the positive v-axis and the plane through the origin.
+  #
+  # The vector is normal to the positive v-axis and a plane that contains the
+  # origin and the v-axis.
+  def get_u_norm(face : Int, u : Float64)
+    case face
+    when 0
+      Point.new(u, -1, 0)
+    when 1
+      Point.new(1, u, 0)
+    when 2
+      Point.new(1, 0, u)
+    when 3
+      Point.new(-u, 0, 1)
+    when 4
+      Point.new(0, -u, 1)
+    else
+      Point.new(0, -1, -u)
+    end
+  end
+
+  # Vector normal to the positive u-axis and the plane through the origin.
+  #
+  # The vector is normal to the positive u-axis and a plane that contains the
+  # origin and the u-axis.
+  def get_v_norm(face, v)
+    case face
+    when 0
+      Point.new(-v, 0, 1)
+    when 1
+      Point.new(0, -v, 1)
+    when 2
+      Point.new(0, -1, -v)
+    when 3
+      Point.new(v, -1, 0)
+    when 4
+      Point.new(1, v, 0)
+    else
+      Point.new(1, 0, v)
+    end
+  end
+end
diff --git a/src/s2_cells/cell_base.cr b/src/s2_cells/cell_base.cr
deleted file mode 100644
index 9d1e555..0000000
--- a/src/s2_cells/cell_base.cr
+++ /dev/null
@@ -1,44 +0,0 @@
-module S2Cells::CellBase
-  LINEAR_PROJECTION    = 0
-  TAN_PROJECTION       = 1
-  QUADRATIC_PROJECTION = 2
-
-  MAX_LEVEL   = 30
-  NUM_FACES   =  6
-  POS_BITS    = 2 * MAX_LEVEL + 1
-  MAX_SIZE    = 1_u64 << MAX_LEVEL
-  SWAP_MASK   = 0x01_u64
-  INVERT_MASK = 0x02_u64
-  LOOKUP_BITS =    4_u64
-  POS_TO_OR   = {SWAP_MASK, 0_u64, 0_u64, INVERT_MASK | SWAP_MASK}
-  POS_TO_IJ   = { {0_u64, 1_u64, 3_u64, 2_u64},
-                 {0_u64, 2_u64, 3_u64, 1_u64},
-                 {3_u64, 2_u64, 0_u64, 1_u64},
-                 {3_u64, 1_u64, 0_u64, 2_u64} }
-
-  LOOKUP_POS = Array.new((1 << (2 * LOOKUP_BITS + 2)), 0_u64)
-  LOOKUP_IJ  = Array.new((1 << (2 * LOOKUP_BITS + 2)), 0_u64)
-
-  def self.lookup_cells(level, i, j, orig_orientation, pos, orientation)
-    return lookup_bits(i, j, orig_orientation, pos, orientation) if level == LOOKUP_BITS
-
-    r = POS_TO_IJ[orientation]
-    4.times do |index|
-      lookup_cells(
-        level + 1, (i << 1) + (r[index] >> 1), (j << 1) + (r[index] & 1),
-        orig_orientation, (pos << 2) + index, orientation ^ POS_TO_OR[index]
-      )
-    end
-  end
-
-  def self.lookup_bits(i, j, orig_orientation, pos, orientation)
-    ij = (i << LOOKUP_BITS) + j
-    LOOKUP_POS[(ij << 2) + orig_orientation] = (pos << 2) + orientation
-    LOOKUP_IJ[(pos << 2) + orig_orientation] = (ij << 2) + orientation
-  end
-
-  lookup_cells(0_u64, 0_u64, 0_u64, 0_u64, 0_u64, 0_u64)
-  lookup_cells(0_u64, 0_u64, 0_u64, SWAP_MASK, 0_u64, SWAP_MASK)
-  lookup_cells(0_u64, 0_u64, 0_u64, INVERT_MASK, 0_u64, INVERT_MASK)
-  lookup_cells(0_u64, 0_u64, 0_u64, SWAP_MASK | INVERT_MASK, 0_u64, SWAP_MASK | INVERT_MASK)
-end
diff --git a/src/s2_cells/cell_id.cr b/src/s2_cells/cell_id.cr
index 386214d..57b95ec 100644
--- a/src/s2_cells/cell_id.cr
+++ b/src/s2_cells/cell_id.cr
@@ -1,42 +1,216 @@
-class S2Cells::CellId
-  include CellBase
+struct S2Cells::CellId
+  include Comparable(CellId)
+
   getter id : UInt64
 
   def initialize(@id)
   end
 
-  # ToToken returns a hex-encoded string of the uint64 cell id, with leading
-  # zeros included but trailing zeros stripped.
-  def to_token : String
-    token = @id.to_s(16).rjust(16, '0').rstrip('0')
-    return "X" if token.size == 0
-    token
+  def hash : UInt64
+    id
   end
 
-  # returns a cell given a hex-encoded string of its uint64 ID
-  def self.from_token(token : String)
-    raise ArgumentError.new("token size was #{token.bytesize}, max size is 16 bytes") if token.bytesize > 16
-    # pad to 16 characters
-    self.new(token.ljust(16, '0').to_u64(16))
+  def ==(other : CellId)
+    @id == other.id
   end
 
-  def parent(level)
-    new_lsb = lsb_for_level(level)
-    s2 = CellId.new((@id & ((~new_lsb) &+ 1)) | new_lsb)
+  def <=>(other : CellId)
+    @id <=> other.id
+  end
 
-    raise InvalidLevel.new(level) unless valid?(s2.id)
-    s2
+  MAX_LEVEL   = 30
+  NUM_FACES   =  6
+  FACE_BITS   =  3
+  POS_BITS    = 2 * MAX_LEVEL + 1
+  MAX_SIZE    = 1_u64 << MAX_LEVEL
+  MAX_SIZE_I  = MAX_SIZE.to_i128
+  WRAP_OFFSET = NUM_FACES << POS_BITS
+
+  def self.lookup_cells(level, i, j, orig_orientation, pos, orientation)
+    return lookup_bits(i, j, orig_orientation, pos, orientation) if level == LOOKUP_BITS
+
+    r = POS_TO_IJ[orientation]
+    4.times do |index|
+      lookup_cells(
+        level + 1, (i << 1) + (r[index] >> 1), (j << 1) + (r[index] & 1),
+        orig_orientation, (pos << 2) + index, orientation ^ POS_TO_ORIENTATION[index]
+      )
+    end
   end
 
-  def prev
-    CellId.new(@id - (lsb << 1))
+  def self.lookup_bits(i, j, orig_orientation, pos, orientation)
+    ij = (i << LOOKUP_BITS) + j
+    LOOKUP_POS[(ij << 2) + orig_orientation] = (pos << 2) + orientation
+    LOOKUP_IJ[(pos << 2) + orig_orientation] = (ij << 2) + orientation
   end
 
-  def next
-    CellId.new(@id + (lsb << 1))
+  lookup_cells(0_u64, 0_u64, 0_u64, 0_u64, 0_u64, 0_u64)
+  lookup_cells(0_u64, 0_u64, 0_u64, SWAP_MASK, 0_u64, SWAP_MASK)
+  lookup_cells(0_u64, 0_u64, 0_u64, INVERT_MASK, 0_u64, INVERT_MASK)
+  lookup_cells(0_u64, 0_u64, 0_u64, SWAP_MASK | INVERT_MASK, 0_u64, SWAP_MASK | INVERT_MASK)
+
+  def self.from_point(p : Point)
+    face, u, v = S2Cells.xyz_to_face_uv(p)
+    i = st_to_ij(uv_to_st(u))
+    j = st_to_ij(uv_to_st(v))
+
+    from_face_ij(face, i, j)
+  end
+
+  def self.from_lat_lng(lat : Float64, lng : Float64)
+    from_point LatLng.from_degrees(lat, lng).to_point
   end
 
-  def level
+  def self.from_lat_lng(lat : Angle, lng : Angle)
+    from_point LatLng.from_angles(lat, lng).to_point
+  end
+
+  def self.from_lat_lng(lat_lng : LatLng)
+    from_point lat_lng.to_point
+  end
+
+  def self.from_face_pos_level(face : Int, pos : UInt64, level : Int) : CellId
+    face = face.to_u64
+    raise ArgumentError.new("Invalid face: #{face} (must be in range 0..5)") unless (0_u64..5_u64).includes?(face)
+
+    # Ensure that the position fits within the bits allocated by the level
+    pos_length = level * 2
+    max_pos = (1_u64 << pos_length) - 1_u64
+    pos &= max_pos # Mask the position to fit the level's precision
+
+    # Shift the face to the top three bits of the UInt64
+    id = face << POS_BITS
+    id = id + (pos | 1_u64)
+
+    lsb_on_level = lsb_for_level(level)
+    CellId.new((id & (~lsb_on_level &+ 1)) | lsb_on_level)
+  end
+
+  def self.from_face_ij(face : Int32, i : Int128 | UInt64, j : Int128 | UInt64) : CellId
+    face = face.to_u64
+    n = face << (POS_BITS - 1)
+    bits = face & SWAP_MASK
+
+    7.downto(0) do |k|
+      mask = (1_u64 << LOOKUP_BITS) - 1
+      bits += (((i >> (k * LOOKUP_BITS)) & mask) << (LOOKUP_BITS + 2))
+      bits += (((j >> (k * LOOKUP_BITS)) & mask) << 2)
+      bits = LOOKUP_POS[bits]
+      n |= (bits >> 2) << (k * 2 * LOOKUP_BITS)
+      bits &= SWAP_MASK | INVERT_MASK
+    end
+
+    new(n * 2 + 1)
+  end
+
+  # Conversion of from_face_ij_wrap method
+  def self.from_face_ij_wrap(face : Int32, i : Int128, j : Int128) : CellId
+    # Convert i and j to the coordinates of a leaf cell just beyond the
+    # boundary of this face.  This prevents 32-bit overflow in the case
+    # of finding the neighbors of a face cell
+    i = {-1_i128, {MAX_SIZE_I, i}.min}.max
+    j = {-1_i128, {MAX_SIZE_I, j}.min}.max
+
+    # Convert (i, j) to (u, v)
+    scale = 1.0 / MAX_SIZE.to_f
+    u = scale * ((i << 1) + 1 - MAX_SIZE_I)
+    v = scale * ((j << 1) + 1 - MAX_SIZE_I)
+
+    # Convert from (u, v) back to face, (s, t) and then to (i, j)
+    face, u, v = S2Cells.xyz_to_face_uv(S2Cells.face_uv_to_xyz(face, u, v))
+    from_face_ij(
+      face,
+      st_to_ij(0.5 * (u + 1)),
+      st_to_ij(0.5 * (v + 1)),
+    )
+  end
+
+  def self.from_face_ij_same(face : Int32, i : Int128, j : Int128, same_face : Bool)
+    same_face ? from_face_ij(face, i, j) : from_face_ij_wrap(face, i, j)
+  end
+
+  def to_s(io : IO)
+    io << "CellId: "
+    io << to_token
+  end
+
+  def self.st_to_ij(s : Float64) : UInt64
+    {0_u64, {MAX_SIZE - 1_u64, (MAX_SIZE * s).floor.to_u64}.min}.max
+  end
+
+  def self.lsb_for_level(level : Int)
+    1_u64 << (2 * (MAX_LEVEL - level))
+  end
+
+  def parent : CellId
+    raise "face cells don't have a parent" if face?
+    new_lsb = lsb << 2
+    self.class.new((@id & (~new_lsb &+ 1)) | new_lsb)
+  end
+
+  def parent(level : Int)
+    current_level = self.level
+    raise "invalid level: #{level}" unless (0...current_level).includes?(level)
+    new_lsb = self.class.lsb_for_level(level)
+    self.class.new((@id & (~new_lsb &+ 1)) | new_lsb)
+  end
+
+  # def child(pos : Int)
+  #  raise "Invalid cell id" unless valid?
+  #  raise "Child position out of range" if leaf?
+  #  new_lsb = lsb >> 2
+  #  self.class.new(@id &+ (2 * pos + 1 - 4) &* new_lsb)
+  # end
+  #
+  def contains(other : CellId) : Bool
+    raise "Invalid cell id" unless valid?
+    raise "Invalid cell id" unless other.valid?
+    other >= range_min && other <= range_max
+  end
+
+  def intersects(other : CellId) : Bool
+    raise "Invalid cell id" unless valid?
+    raise "Invalid cell id" unless other.valid?
+    other.range_min <= range_max && other.range_max >= range_min
+  end
+
+  def face?
+    (@id & (self.class.lsb_for_level(0) - 1)) == 0
+  end
+
+  def self.valid?(cell_id : CellId)
+    # 6 is an invalid face
+    return false unless cell_id.face < NUM_FACES
+    (cell_id.lsb & 0x1555555555555555_u64) != 0
+  end
+
+  private def valid?(cell_id : UInt64)
+    self.class.valid? CellId.new(cell_id)
+  end
+
+  def valid?
+    self.class.valid? self
+  end
+
+  def lsb
+    @id & (~@id &+ 1) # This is equivalent to (@id & -@id) for signed integers
+  end
+
+  def face : Int32
+    # Shift right by 61 bits to move the top 3 bits to the least significant bit position.
+    # Mask with binary 111 to isolate these three bits.
+    (@id >> POS_BITS).to_i & 0b111
+  end
+
+  def pos : UInt64
+    @id & (UInt64::MAX >> FACE_BITS)
+  end
+
+  def leaf?
+    (@id & 1) != 0
+  end
+
+  def level : Int32
     return MAX_LEVEL if leaf?
 
     level = -1
@@ -49,7 +223,7 @@ class S2Cells::CellId
     end
 
     # 2s compliment
-    x &= ((~x) &+ 1)
+    x &= (~x &+ 1)
 
     level += 8 unless (x & 0x00005555_u64).zero?
     level += 4 unless (x & 0x00550055_u64).zero?
@@ -58,98 +232,335 @@ class S2Cells::CellId
     level
   end
 
-  def self.from_lat_lon(lat : Float64, lon : Float64)
-    from_point LatLon.new(lat, lon).to_point
+  # ToToken returns a hex-encoded string of the uint64 cell id, with leading
+  # zeros included but trailing zeros stripped.
+  def to_token : String
+    token = @id.to_s(16).rjust(16, '0').rstrip('0')
+    return "X" if token.size == 0
+    token
   end
 
-  def self.from_point(p : Point)
-    face, u, v = xyz_to_face_uv(p)
-    i = st_to_ij(uv_to_st(u))
-    j = st_to_ij(uv_to_st(v))
+  # returns a cell given a hex-encoded string of its uint64 ID
+  def self.from_token(token : String)
+    raise ArgumentError.new("token size was #{token.bytesize}, max size is 16 bytes") if token.bytesize > 16
+    # pad to 16 characters
+    self.new(token.ljust(16, '0').to_u64(16))
+  end
 
-    CellId.new(from_face_ij(face, i, j))
+  def child_begin
+    old_lsb = lsb
+    self.class.new(@id &- old_lsb &+ (old_lsb >> 2))
   end
 
-  def self.from_face_ij(face : UInt64, i : UInt64, j : UInt64)
-    n = face << (POS_BITS - 1)
-    bits = face & SWAP_MASK
+  def child_begin(level : Int)
+    raise "invalid level: #{level}" unless (0..30).includes?(level)
+    self.class.new(@id &- lsb &+ self.class.lsb_for_level(level))
+  end
 
-    7.downto(0).each do |k|
-      mask = (1_u64 << LOOKUP_BITS) - 1
-      bits += (((i >> (k * LOOKUP_BITS)) & mask) << (LOOKUP_BITS + 2))
-      bits += (((j >> (k * LOOKUP_BITS)) & mask) << 2)
-      bits = LOOKUP_POS[bits]
-      n |= (bits >> 2) << (k * 2 * LOOKUP_BITS)
-      bits &= (SWAP_MASK | INVERT_MASK)
+  def child_end
+    old_lsb = lsb
+    self.class.new(@id &+ old_lsb &+ (old_lsb >> 2))
+  end
+
+  def child_end(level : Int)
+    raise "invalid level: #{level}" unless (0..30).includes?(level)
+    self.class.new(@id &+ lsb &+ self.class.lsb_for_level(level))
+  end
+
+  def prev
+    CellId.new(@id &- (lsb << 1))
+  end
+
+  def next
+    CellId.new(@id &+ (lsb << 1))
+  end
+
+  def children(level : Int32? = nil, &)
+    if level
+      cell_id = child_begin(level)
+      ending = child_end(level)
+    else
+      cell_id = child_begin
+      ending = child_end
     end
 
-    n * 2 + 1
+    while cell_id != ending
+      yield cell_id
+      cell_id = cell_id.next
+    end
   end
 
-  def self.xyz_to_face_uv(p : Point) : Tuple(UInt64, Float64, Float64)
-    face = p.largest_abs_component
+  def children(level : Int32? = nil) : Array(CellId)
+    cells = Array(CellId).new(4)
+    children(level) { |child| cells << child }
+    cells
+  end
 
-    pface = case face
-            when 0_u64 then p.x
-            when 1_u64 then p.y
-            else            p.z
-            end
+  def range_min
+    self.class.new(@id &- (lsb &- 1))
+  end
 
-    face += 3_u64 if pface < 0.0
+  def range_max
+    self.class.new(@id &+ (lsb &- 1))
+  end
+
+  def self.range_begin(level : Int)
+    from_face_pos_level(0_u64, 0_u64, 0).child_begin(level)
+  end
+
+  def self.range_end(level : Int)
+    from_face_pos_level(5_u64, 0_u64, 0).child_end(level)
+  end
+
+  def self.walk(level : Int, &)
+    begin_cell = range_begin(level)
+    cellid_int = begin_cell.id
+    endid_int = range_end(level).id
 
-    u, v = valid_face_xyz_to_uv(face, p)
-    {face, u, v}
+    # Doubling the lsb yields the increment between positions at a certain
+    # level as 64-bit IDs. See CellId docstring for bit encoding.
+    increment = begin_cell.lsb << 1
+
+    while cellid_int != endid_int
+      yield new(cellid_int)
+      cellid_int += increment
+    end
+  end
+
+  def self.none
+    new
   end
 
+  # TODO:: prev_wrap, next_wrap, advance_wrap, advance
+
   def self.uv_to_st(u : Float64)
     return 0.5 * Math.sqrt(1.0 + 3.0 * u) if u >= 0.0
     1.0 - 0.5 * Math.sqrt(1.0 - 3.0 * u)
   end
 
-  def self.st_to_ij(s : Float64) : UInt64
-    {0_u64, {MAX_SIZE - 1_u64, (MAX_SIZE * s).floor.to_u64}.min}.max
+  def prev
+    self.class.new(@id &- (lsb << 1))
+  end
+
+  def next
+    self.class.new(@id &+ (lsb << 1))
+  end
+
+  def to_lat_lng : LatLng
+    LatLng.from_point(self.to_point_raw)
+  end
+
+  def to_point_raw : Point
+    face, si, ti = self.get_center_si_ti
+    S2Cells.face_uv_to_xyz(
+      face,
+      self.class.st_to_uv((0.5 / MAX_SIZE) * si),
+      self.class.st_to_uv((0.5 / MAX_SIZE) * ti),
+    )
   end
 
-  def self.valid_face_xyz_to_uv(face : UInt64, p : Point)
-    raise "invalid face xyz" unless p.dot_prod(face_uv_to_xyz(face, 0.0, 0.0)) > 0
+  def to_point : Point
+    to_point_raw.normalize
+  end
+
+  def get_center_si_ti
+    face, i, j, orientation = self.to_face_ij_orientation
 
-    case face
-    when 0 then [p.y / p.x, p.z / p.x]
-    when 1 then [-p.x / p.y, p.z / p.y]
-    when 2 then [-p.x / p.z, -p.y / p.z]
-    when 3 then [p.z / p.x, p.y / p.x]
-    when 4 then [p.z / p.y, -p.x / p.y]
-    else        [-p.y / p.z, -p.x / p.z]
+    if self.leaf?
+      delta = 1
+    elsif ((i ^ (self.id >> 2)) & 1) != 0_u64
+      delta = 2
+    else
+      delta = 0
     end
+
+    {face, 2_u64 &* i &+ delta, 2_u64 &* j &+ delta}
+  end
+
+  def get_center_uv
+    face, si, ti = get_center_si_ti
+    {
+      self.class.st_to_uv((0.5 / MAX_SIZE) * si),
+      self.class.st_to_uv((0.5 / MAX_SIZE) * ti),
+    }
+  end
+
+  def to_face_ij_orientation : Tuple(Int32, UInt64, UInt64, Int32)
+    i, j = 0_u64, 0_u64
+    face = self.face
+    bits = face & SWAP_MASK
+
+    7.downto(0) do |k|
+      if k == 7
+        nbits = MAX_LEVEL - 7 * LOOKUP_BITS
+      else
+        nbits = LOOKUP_BITS
+      end
+
+      bits += (
+        self.id >> (k * 2 * LOOKUP_BITS + 1) &
+        ((1 << (2 * nbits)) - 1)
+      ) << 2
+      bits = LOOKUP_IJ[bits]
+      i += (bits >> (LOOKUP_BITS + 2)) << (k * LOOKUP_BITS)
+      j += ((bits >> 2) & ((1 << LOOKUP_BITS) - 1)) << (k * LOOKUP_BITS)
+      bits &= SWAP_MASK | INVERT_MASK
+    end
+
+    raise "Assertion failed" unless POS_TO_ORIENTATION[2] == 0
+    raise "Assertion failed" unless SWAP_MASK == POS_TO_ORIENTATION[0]
+    if (self.lsb & 0x1111111111111110_u64) != 0
+      bits ^= SWAP_MASK
+    end
+    orientation = bits
+
+    {face, i, j, orientation.to_i}
+  end
+
+  def get_edge_neighbors : Array(CellId)
+    level = self.level
+    size = get_size_ij(level)
+    face, i, j, orientation = to_face_ij_orientation
+
+    i = i.to_i128
+    j = j.to_i128
+
+    [
+      (self.class.from_face_ij_same(face, i, j - size, j - size >= 0).parent(level)),
+      (self.class.from_face_ij_same(face, i + size, j, i + size < MAX_SIZE_I).parent(level)),
+      (self.class.from_face_ij_same(face, i, j + size, j + size < MAX_SIZE_I).parent(level)),
+      (self.class.from_face_ij_same(face, i - size, j, i - size >= 0).parent(level)),
+    ]
+  end
+
+  def get_vertex_neighbors(level : Int32) : Array(CellId)
+    # "level" must be strictly less than this cell's level so that we can
+    # determine which vertex this cell is closest to.
+    raise "Invalid level" unless level < self.level
+
+    face, i, j, orientation = self.to_face_ij_orientation
+    i = i.to_i128
+    j = j.to_i128
+
+    # Determine the i- and j-offsets to the closest neighboring cell in
+    # each direction. This involves looking at the next bit of "i" and
+    # "j" to determine which quadrant of this->parent(level) this cell
+    # lies in.
+    halfsize = self.get_size_ij(level + 1).to_i128
+    size = halfsize << 1
+    if (i & halfsize) != 0
+      ioffset = size
+      isame = (i + size) < MAX_SIZE_I
+    else
+      ioffset = -size
+      isame = (i - size) >= 0
+    end
+
+    if (j & halfsize) != 0
+      joffset = size
+      jsame = (j + size) < MAX_SIZE_I
+    else
+      joffset = -size
+      jsame = (j - size) >= 0
+    end
+
+    neighbors = [] of CellId
+    neighbors << self.parent(level)
+    neighbors << self.class.from_face_ij_same(face, i + ioffset, j, isame).parent(level)
+    neighbors << self.class.from_face_ij_same(face, i, j + joffset, jsame).parent(level)
+    if isame || jsame
+      neighbors << self.class.from_face_ij_same(face, i + ioffset, j + joffset, isame && jsame).parent(level)
+    end
+
+    neighbors
   end
 
-  def self.face_uv_to_xyz(face : UInt64, u : Float64, v : Float64)
-    case face
-    when 0 then Point.new(1_f64, u, v)
-    when 1 then Point.new(-u, 1_f64, v)
-    when 2 then Point.new(-u, -v, 1_f64)
-    when 3 then Point.new(-1_f64, -v, -u)
-    when 4 then Point.new(v, -1_f64, -u)
-    else        Point.new(v, u, -1_f64)
+  def get_all_neighbors(nbr_level : Int32 = self.level) : Array(CellId)
+    face, i, j, orientation = self.to_face_ij_orientation
+    i = i.to_i128
+    j = j.to_i128
+
+    # Find the coordinates of the lower left-hand leaf cell. Normalize (i, j).
+    size = self.get_size_ij.to_i128
+    i &= -size
+    j &= -size
+
+    nbr_size = self.get_size_ij(nbr_level).to_i128
+    raise "Invalid neighborhood size" unless nbr_size <= size
+
+    neighbors = [] of CellId
+
+    # Compute the N-S, E-W, and diagonal neighbors in one pass.
+    k = -nbr_size
+    loop do
+      if k < 0
+        same_face = (j + k >= 0)
+      elsif k >= size
+        same_face = (j + k < MAX_SIZE_I)
+      else
+        same_face = false
+        # North and South neighbors
+        neighbors << self.class.from_face_ij_same(face, i + k, j - nbr_size, j - size >= 0).parent(nbr_level)
+        neighbors << self.class.from_face_ij_same(face, i + k, j + size, j + size < MAX_SIZE).parent(nbr_level)
+      end
+
+      # East and West neighbors
+      neighbors << self.class.from_face_ij_same(face, i - nbr_size, j + k, same_face && i - size >= 0).parent(nbr_level)
+      neighbors << self.class.from_face_ij_same(face, i + size, j + k, same_face && i + size < MAX_SIZE).parent(nbr_level)
+
+      break if k >= size
+      k += nbr_size
     end
+
+    neighbors
+  end
+
+  def get_size_ij(level = self.level)
+    1_u64 << (MAX_LEVEL - level)
   end
 
-  private def leaf?
-    @id & 1 != 0
+  def self.max_edge
+    LengthMetric.new(max_angle_span.deriv)
   end
 
-  private def valid?(s2_id : UInt64)
-    face = s2_id >> POS_BITS
-    # We're performing 2s compliment
-    lsb = s2_id & ((~s2_id) &+ 1)
-    (face < NUM_FACES) && ((lsb & 0x1555555555555555_u64) != 0)
+  def self.max_angle_span
+    # LINEAR_PROJECTION
+    # LengthMetric.new(2)
+
+    # TAN_PROJECTION
+    # LengthMetric.new(Math::PI / 2)
+
+    # QUADRATIC_PROJECTION
+    LengthMetric.new(1.704897179199218452)
   end
 
-  private def lsb
-    @id & -@id
+  def self.max_diag
+    # LINEAR_PROJECTION
+    # LengthMetric.new(2 * Math.sqrt(2))
+
+    # TAN_PROJECTION
+    # LengthMetric.new(Math::PI * Math.sqrt(2.0 / 3.0))
+
+    # QUADRATIC_PROJECTION
+    LengthMetric.new(2.438654594434021032)
   end
 
-  private def lsb_for_level(level)
-    1_u64 << (2 * (MAX_LEVEL - level))
+  def self.min_width
+    # LINEAR_PROJECTION
+    # LengthMetric.new(Math.sqrt(2))
+
+    # TAN_PROJECTION
+    # LengthMetric.new(Math::PI / 2 * Math.sqrt(2))
+
+    # QUADRATIC_PROJECTION
+    LengthMetric.new(2 * Math.sqrt(2) / 3)
+  end
+
+  def self.st_to_uv(s : Float64) : Float64
+    if s >= 0.5
+      (1.0 / 3.0) * (4 * s * s - 1)
+    else
+      (1.0 / 3.0) * (1 - 4 * (1 - s) * (1 - s))
+    end
   end
 end
diff --git a/src/s2_cells/cell_union.cr b/src/s2_cells/cell_union.cr
new file mode 100644
index 0000000..bb7abb6
--- /dev/null
+++ b/src/s2_cells/cell_union.cr
@@ -0,0 +1,251 @@
+struct S2Cells::CellUnion
+  getter cell_ids : Array(CellId)
+
+  def initialize(cell_ids : Array(CellId)? = nil, raw : Bool = true)
+    if cell_ids.nil?
+      @cell_ids = [] of CellId
+    else
+      @cell_ids = cell_ids.map do |cell_id|
+        cell_id.is_a?(CellId) ? cell_id : CellId.new(cell_id)
+      end
+      normalize if raw
+    end
+  end
+
+  def ==(other : CellUnion) : Bool
+    other.is_a?(CellUnion) && @cell_ids == other.cell_ids
+  end
+
+  def hash : UInt64
+    @cell_ids.hash
+  end
+
+  def to_s : String
+    "#{self.class.name}: #{@cell_ids}"
+  end
+
+  def self.get_union(x : CellUnion, y : CellUnion) : CellUnion
+    new(x.cell_ids + y.cell_ids)
+  end
+
+  def self.get_intersection(cell_union : CellUnion, cell_id : CellId) : CellUnion
+    if cell_union.contains(cell_id)
+      new([cell_id])
+    else
+      index = cell_union.cell_ids.bsearch_index { |id| id >= cell_id.range_min } || 0
+      idmax = cell_id.range_max
+      intersected_cell_ids = [] of CellId
+      while index < cell_union.cell_ids.size && cell_union.cell_ids[index] <= idmax
+        intersected_cell_ids << cell_union.cell_ids[index]
+        index += 1
+      end
+      new(intersected_cell_ids)
+    end
+  end
+
+  def self.get_intersection(x : CellUnion, y : CellUnion) : CellUnion
+    i, j = 0, 0
+    cell_ids = [] of CellId
+    while i < x.num_cells && j < y.num_cells
+      imin = x.cell_ids[i].range_min
+      jmin = y.cell_ids[j].range_min
+      if imin > jmin
+        if x.cell_ids[i] <= y.cell_ids[j].range_max
+          cell_ids << x.cell_ids[i]
+          i += 1
+        else
+          j = y.cell_ids.bsearch_index { |id| id >= imin } || (j + 1)
+          j -= 1 if j > 0 && x.cell_ids[i] <= y.cell_ids[j].range_max
+        end
+      elsif jmin > imin
+        if y.cell_ids[j] <= x.cell_ids[i].range_max
+          cell_ids << y.cell_ids[j]
+          j += 1
+        else
+          i = x.cell_ids.bsearch_index { |id| id >= jmin } || (i + 1)
+          i -= 1 if i > 0 && y.cell_ids[j] <= x.cell_ids[i].range_max
+        end
+      else
+        if x.cell_ids[i] < y.cell_ids[j]
+          cell_ids << x.cell_ids[i]
+          i += 1
+        else
+          cell_ids << y.cell_ids[j]
+          j += 1
+        end
+      end
+    end
+
+    cell_union = new(cell_ids)
+    cell_union.normalize
+    cell_union
+  end
+
+  def expand(level : Int32)
+    output = [] of CellId
+    level_lsb = CellId.lsb_for_level(level)
+    i = num_cells - 1
+    while i >= 0
+      cell_id = @cell_ids[i]
+      if cell_id.lsb < level_lsb
+        cell_id = cell_id.parent(level)
+        while i > 0 && cell_id.contains(@cell_ids[i - 1])
+          i -= 1
+        end
+      end
+      output << cell_id
+      cell_id.append_all_neighbors(level, output)
+      i -= 1
+    end
+    @cell_ids = output
+  end
+
+  def expand(min_radius : Angle, max_level_diff : Int32)
+    min_level = CellId::MAX_LEVEL
+    @cell_ids.each do |cell_id|
+      min_level = {min_level, cell_id.level}.min
+    end
+
+    radius_level = CellId.min_width.get_max_level(min_radius.radians)
+    if radius_level == 0 && min_radius.radians > CellId.min_width.get_value(0)
+      expand(0)
+    end
+    expand({min_level + max_level_diff, radius_level}.min)
+  end
+
+  def self.get_difference(x : CellUnion, y : CellUnion) : CellUnion
+    cell_ids = [] of CellId
+    x.cell_ids.each do |cell_id|
+      __get_difference(cell_id, y, cell_ids)
+    end
+
+    cell_union = new(cell_ids)
+    cell_union.normalize
+    cell_union
+  end
+
+  private def self.__get_difference(cell_id : CellId, y : CellUnion, cell_ids : Array(CellId))
+    if !y.intersects(cell_id)
+      cell_ids << cell_id
+    elsif !y.contains(cell_id)
+      cell_id.children.each do |child|
+        __get_difference(child, y, cell_ids)
+      end
+    end
+  end
+
+  def num_cells : Int32
+    @cell_ids.size
+  end
+
+  def cell_id(i : Int32) : CellId
+    @cell_ids[i]
+  end
+
+  def normalize
+    @cell_ids.sort!
+    output = [] of CellId
+    @cell_ids.each do |cell_id|
+      next if output.any? && output.last.contains(cell_id)
+
+      while output.any? && cell_id.contains(output.last)
+        output.pop
+      end
+
+      while output.size >= 3
+        if (output[-3].id ^ output[-2].id ^ output[-1].id) != cell_id.id
+          break
+        end
+
+        mask = cell_id.lsb << 1
+        mask = ~(mask + (mask << 1))
+        id_masked = (cell_id.id & mask)
+        if (output[-3].id & mask) != id_masked || (output[-2].id & mask) != id_masked || (output[-1].id & mask) != id_masked || cell_id.face?
+          break
+        end
+
+        output.pop
+        output.pop
+        output.pop
+        cell_id = cell_id.parent
+      end
+
+      output << cell_id
+    end
+
+    if output.size < num_cells
+      @cell_ids = output
+      return true
+    end
+    false
+  end
+
+  def denormalize(min_level : Int32, level_mod : Int32) : Array(CellId)
+    raise "Invalid min_level" unless min_level >= 0 && min_level <= CellId::MAX_LEVEL
+    raise "Invalid level_mod" unless level_mod >= 1 && level_mod <= 3
+
+    cell_ids = [] of CellId
+    @cell_ids.each do |cell_id|
+      level = cell_id.level
+      new_level = [min_level, level].max
+      if level_mod > 1
+        new_level += ((CellId::MAX_LEVEL - (new_level - min_level)) % level_mod)
+        new_level = [CellId::MAX_LEVEL, new_level].min
+      end
+      if new_level == level
+        cell_ids << cell_id
+      else
+        cell_id.children(new_level).each do |child|
+          cell_ids << child
+        end
+      end
+    end
+    cell_ids
+  end
+
+  def contains(other : CellUnion | CellId | Cell | Point) : Bool
+    case other
+    in Cell
+      contains(other.id)
+    in CellId
+      cell_id = other
+      index = @cell_ids.bsearch_index { |id| id >= cell_id } || 0
+      if index < @cell_ids.size && @cell_ids[index].range_min <= cell_id
+        return true
+      end
+      index != 0 && @cell_ids[index - 1].range_max >= cell_id
+    in Point
+      contains(CellId.from_point(other))
+    in CellUnion
+      other.cell_ids.each do |cell_id|
+        return false unless contains(cell_id)
+      end
+      true
+    end
+  end
+
+  def intersects(other : CellUnion | CellId) : Bool
+    case other
+    in CellId
+      cell_id = other
+      index = @cell_ids.bsearch_index { |id| id >= cell_id } || 0
+      if index != @cell_ids.size && @cell_ids[index].range_min <= cell_id.range_max
+        return true
+      end
+      index != 0 && @cell_ids[index - 1].range_max >= cell_id.range_min
+    in CellUnion
+      other.cell_ids.each do |cell_id|
+        return true if intersects(cell_id)
+      end
+      false
+    end
+  end
+
+  def get_rect_bound : LatLngRect
+    bound = LatLngRect.empty
+    @cell_ids.each do |cell_id|
+      bound = bound.union(Cell.new(cell_id).get_rect_bound)
+    end
+    bound
+  end
+end
diff --git a/src/s2_cells/interval.cr b/src/s2_cells/interval.cr
new file mode 100644
index 0000000..8db83e3
--- /dev/null
+++ b/src/s2_cells/interval.cr
@@ -0,0 +1,404 @@
+require "math"
+
+module S2Cells
+  abstract struct Interval
+    # Initialize with low and high bounds
+    def initialize(lo : Float | Int, hi : Float | Int)
+      @lo = lo.to_f
+      @hi = hi.to_f
+    end
+
+    # String representation of the Interval object
+    def to_s : String
+      "#{self.class.name}: (#{@lo}, #{@hi})"
+    end
+
+    # Accessor methods for lo and hi
+    getter lo : Float64
+    getter hi : Float64
+
+    # Method to get bound by index
+    def bound(i : Int) : Float64
+      case i
+      when 0 then @lo
+      when 1 then @hi
+      else        raise "Index out of bounds"
+      end
+    end
+
+    # Method to get both bounds as a tuple
+    def bounds : Tuple(Float64, Float64)
+      {@lo, @hi}
+    end
+
+    # Class method to return an empty Interval
+    def self.empty : Interval
+      Interval.new(0.0, 0.0)
+    end
+  end
+
+  struct LineInterval < Interval
+    def initialize(lo : Float | Int = 1.0, hi : Float | Int = 0.0)
+      super(lo, hi)
+    end
+
+    def ==(other : LineInterval) : Bool
+      (self.lo == other.lo && self.hi == other.hi) || (self.empty? && other.empty?)
+    end
+
+    def hash : UInt64
+      {@lo, @hi}.hash
+    end
+
+    def self.from_point_pair(a : Float | Int, b : Float | Int) : LineInterval
+      if a <= b
+        new(a, b)
+      else
+        new(b, a)
+      end
+    end
+
+    def contains(other : LineInterval | Float | Int) : Bool
+      case other
+      when LineInterval
+        return true if other.empty?
+        other.lo >= self.lo && other.hi <= self.hi
+      else
+        other = other.to_f
+        other >= self.lo && other <= self.hi
+      end
+    end
+
+    def interior_contains(other : LineInterval | Float | Int) : Bool
+      case other
+      when LineInterval
+        return true if other.empty?
+        other.lo > self.lo && other.hi < self.hi
+      else
+        other = other.to_f
+        other > self.lo && other < self.hi
+      end
+    end
+
+    def intersects(other : LineInterval) : Bool
+      if self.lo <= other.lo
+        other.lo <= self.hi && other.lo <= other.hi
+      else
+        self.lo <= other.hi && self.lo <= self.hi
+      end
+    end
+
+    def interior_intersects(other : LineInterval) : Bool
+      other.lo < self.hi && self.lo < other.hi && self.lo < self.hi && other.lo <= other.hi
+    end
+
+    def union(other : LineInterval) : LineInterval
+      return other if self.empty?
+      return self if other.empty?
+      LineInterval.new({self.lo, other.lo}.min, {self.hi, other.hi}.max)
+    end
+
+    def intersection(other : LineInterval) : LineInterval
+      LineInterval.new({self.lo, other.lo}.max, {self.hi, other.hi}.min)
+    end
+
+    def expanded(radius : Float | Int) : LineInterval
+      raise "Radius must be non-negative" if radius.negative?
+      return self if self.empty?
+      radius = radius.to_f
+      LineInterval.new(self.lo - radius, self.hi + radius)
+    end
+
+    def get_center : Float64
+      0.5 * (self.lo + self.hi)
+    end
+
+    def get_length : Float64
+      self.hi - self.lo
+    end
+
+    def empty? : Bool
+      self.lo > self.hi
+    end
+
+    def approx_equals?(other : LineInterval, max_error : Float64 = 1e-15) : Bool
+      return other.get_length.to_f <= max_error if self.empty?
+      return self.get_length.to_f <= max_error if other.empty?
+      (other.lo - self.lo).abs + (other.hi - self.hi).abs <= max_error
+    end
+  end
+
+  struct SphereInterval < Interval
+    def initialize(lo : Float | Int = Math::PI, hi : Float | Int = -Math::PI, args_checked : Bool = false)
+      if args_checked
+        super(lo, hi)
+      else
+        clamped_lo, clamped_hi = lo, hi
+        if lo == -Math::PI && hi != Math::PI
+          clamped_lo = Math::PI
+        end
+        if hi == -Math::PI && lo != Math::PI
+          clamped_hi = Math::PI
+        end
+        super(clamped_lo, clamped_hi)
+      end
+      raise "Invalid interval" unless valid?
+    end
+
+    def ==(other : SphereInterval) : Bool
+      self.lo == other.lo && self.hi == other.hi
+    end
+
+    def self.from_point_pair(a : Float64, b : Float64) : SphereInterval
+      raise "Value out of bounds" unless a.abs <= Math::PI && b.abs <= Math::PI
+      a = Math::PI if a == -Math::PI
+      b = Math::PI if b == -Math::PI
+      if positive_distance(a, b) <= Math::PI
+        new(a, b, args_checked: true)
+      else
+        new(b, a, args_checked: true)
+      end
+    end
+
+    def self.positive_distance(a : Float64, b : Float64) : Float64
+      d = b - a
+      d >= 0 ? d : (b + Math::PI) - (a - Math::PI)
+    end
+
+    def self.full : SphereInterval
+      new(-Math::PI, Math::PI, args_checked: true)
+    end
+
+    def full? : Bool
+      (self.hi - self.lo) == 2 * Math::PI
+    end
+
+    def valid? : Bool
+      self.lo.abs <= Math::PI &&
+        self.hi.abs <= Math::PI &&
+        !(self.lo == -Math::PI && self.hi != Math::PI) &&
+        !(self.hi == -Math::PI && self.lo != Math::PI)
+    end
+
+    def inverted? : Bool
+      self.lo > self.hi
+    end
+
+    def empty? : Bool
+      self.lo - self.hi == 2 * Math::PI
+    end
+
+    def get_center : Float64
+      center = 0.5 * (self.lo + self.hi)
+      if !inverted?
+        center
+      elsif center <= 0
+        center + Math::PI
+      else
+        center - Math::PI
+      end
+    end
+
+    def get_length : Float64
+      length = self.hi - self.lo
+      if length >= 0
+        length
+      else
+        length += 2 * Math::PI
+        length > 0.0 ? length : -1.0
+      end
+    end
+
+    def complement : SphereInterval
+      self.lo == self.hi ? self.class.full : self.class.new(self.hi, self.lo)
+    end
+
+    def approx_equals(other : SphereInterval, max_error = 1e-15) : Bool
+      if self.empty?
+        other.get_length <= max_error
+      elsif other.empty?
+        self.get_length <= max_error
+      else
+        (other.lo - self.lo).modulo(2 * Math::PI).abs + (other.hi - self.hi).modulo(2 * Math::PI).abs <= max_error
+      end
+    end
+
+    def fast_contains(other : Float64) : Bool
+      if self.inverted?
+        (other >= self.lo || other <= self.hi) && !self.empty?
+      else
+        other >= self.lo && other <= self.hi
+      end
+    end
+
+    def contains(other : SphereInterval | Float64) : Bool
+      case other
+      when SphereInterval
+        if self.inverted?
+          if other.inverted?
+            other.lo >= self.lo && other.hi <= self.hi
+          else
+            (other.lo >= self.lo || other.hi <= self.hi) && !self.empty?
+          end
+        else
+          if other.inverted?
+            self.full? || other.empty?
+          else
+            other.lo >= self.lo && other.hi <= self.hi
+          end
+        end
+      else
+        raise "Value out of bounds" unless other.abs <= Math::PI
+        other = Math::PI if other == -Math::PI
+        fast_contains(other)
+      end
+    end
+
+    def interior_contains(other : SphereInterval | Float64) : Bool
+      case other
+      when SphereInterval
+        if self.inverted?
+          if !other.inverted?
+            other.lo > self.lo || other.hi < self.hi
+          else
+            (other.lo > self.lo && other.hi < self.hi) || other.empty?
+          end
+        else
+          if other.inverted?
+            self.full? || other.empty?
+          else
+            (other.lo > self.lo && other.hi < self.hi) || self.full?
+          end
+        end
+      else
+        raise "Value out of bounds" unless other.abs <= Math::PI
+        other = Math::PI if other == -Math::PI
+        if self.inverted?
+          other > self.lo || other < self.hi
+        else
+          (other > self.lo && other < self.hi) || self.full?
+        end
+      end
+    end
+
+    def intersects(other : SphereInterval) : Bool
+      return false if self.empty? || other.empty?
+      if self.inverted?
+        other.inverted? || other.lo <= self.hi || other.hi >= self.lo
+      else
+        if other.inverted?
+          other.lo <= self.hi || other.hi >= self.lo
+        else
+          other.lo <= self.hi && other.hi >= self.lo
+        end
+      end
+    end
+
+    def interior_intersects(other : SphereInterval) : Bool
+      return false if self.empty? || other.empty? || self.lo == self.hi
+      if self.inverted?
+        other.inverted? || other.lo < self.hi || other.hi > self.lo
+      else
+        if other.inverted?
+          other.lo < self.hi || other.hi > self.lo
+        else
+          (other.lo < self.hi && other.hi > self.lo) || self.full?
+        end
+      end
+    end
+
+    def union(other : SphereInterval) : SphereInterval
+      return self if other.empty?
+
+      if fast_contains(other.lo)
+        if fast_contains(other.hi)
+          return self.contains(other) ? self : self.class.full
+        end
+        return self.class.new(self.lo, other.hi, args_checked: true)
+      end
+
+      if fast_contains(other.hi)
+        return self.class.new(other.lo, self.hi, args_checked: true)
+      end
+
+      if self.empty? || other.fast_contains(self.lo)
+        return other
+      end
+
+      dlo = self.class.positive_distance(other.hi, self.lo)
+      dhi = self.class.positive_distance(self.hi, other.lo)
+      if dlo < dhi
+        self.class.new(other.lo, self.hi, args_checked: true)
+      else
+        self.class.new(self.lo, other.hi, args_checked: true)
+      end
+    end
+
+    def intersection(other : SphereInterval) : SphereInterval
+      return self.class.empty if other.empty?
+      if fast_contains(other.lo)
+        if fast_contains(other.hi)
+          return other.get_length < self.get_length ? other : self
+        end
+        return self.class.new(other.lo, self.hi, args_checked: true)
+      end
+
+      if fast_contains(other.hi)
+        return self.class.new(self.lo, other.hi, args_checked: true)
+      end
+
+      if other.fast_contains(self.lo)
+        return self
+      end
+      raise "Intervals do not intersect" unless intersects(other)
+      self.class.empty
+    end
+
+    def expanded(radius : Float64) : SphereInterval
+      raise "Radius must be non-negative" if radius < 0
+      return self if self.empty?
+
+      two_pi = 2 * Math::PI
+      return self.class.full if (self.get_length + 2 * radius) >= (two_pi - 1e-15)
+
+      lo = (self.lo - radius).remainder two_pi
+      hi = (self.hi + radius).remainder two_pi
+      lo = Math::PI if lo <= -Math::PI
+
+      self.class.new(lo, hi)
+    end
+
+    def get_complement_center : Float64
+      if self.lo != self.hi
+        return complement.get_center
+      else
+        self.hi <= 0 ? self.hi + Math::PI : self.hi - Math::PI
+      end
+    end
+
+    def get_directed_hausdorff_distance(other : SphereInterval) : Float64
+      return 0.0 if other.contains(self)
+      return Math::PI if other.empty?
+
+      other_complement_center = other.get_complement_center
+      if self.contains(other_complement_center)
+        return self.class.positive_distance(other.hi, other_complement_center)
+      else
+        hi_hi = if self.class.new(other.hi, other_complement_center).contains(self.hi)
+                  self.class.positive_distance(other.hi, self.hi)
+                else
+                  0
+                end
+
+        lo_lo = if self.class.new(other_complement_center, other.lo).contains(self.lo)
+                  self.class.positive_distance(self.lo, other.lo)
+                else
+                  0
+                end
+
+        raise "Invalid distance calculation" if hi_hi <= 0 && lo_lo <= 0
+        {hi_hi, lo_lo}.max
+      end
+    end
+  end
+end
diff --git a/src/s2_cells/lat_lng.cr b/src/s2_cells/lat_lng.cr
new file mode 100644
index 0000000..ac43d3c
--- /dev/null
+++ b/src/s2_cells/lat_lng.cr
@@ -0,0 +1,131 @@
+require "math"
+
+struct S2Cells::LatLng
+  # A point on a sphere in latitude-longitude coordinates.
+
+  # Creates a LatLng from degrees.
+  def self.from_degrees(lat : Float64, lng : Float64) : LatLng
+    from_angles(Angle.from_degrees(lat), Angle.from_degrees(lng))
+  end
+
+  # Creates a LatLng from radians.
+  def self.from_radians(lat : Float64, lng : Float64) : LatLng
+    new(lat, lng)
+  end
+
+  # Creates a LatLng from a Point object.
+  def self.from_point(point : Point) : LatLng
+    new(LatLng.latitude(point).radians, LatLng.longitude(point).radians)
+  end
+
+  # Creates a LatLng from two angles.
+  def self.from_angles(lat : Angle, lng : Angle) : LatLng
+    new(lat.radians, lng.radians)
+  end
+
+  # Creates a default LatLng (0, 0).
+  def self.default : LatLng
+    new(0.0, 0.0)
+  end
+
+  # Creates an invalid LatLng.
+  def self.invalid : LatLng
+    new(Math::PI, 2.0 * Math::PI)
+  end
+
+  # Initializes a LatLng with latitude and longitude in radians.
+  def initialize(@lat : Float64, @lng : Float64)
+  end
+
+  # Checks equality of LatLng with another LatLng.
+  def ==(other : LatLng) : Bool
+    @lat == other.lat && @lng == other.lng
+  end
+
+  # Checks inequality of LatLng with another LatLng.
+  def !=(other : LatLng) : Bool
+    !(self == other)
+  end
+
+  # Returns the hash code of LatLng.
+  def hash : Int32
+    {@lat, @lng}.hash
+  end
+
+  # Returns the string representation of LatLng.
+  def to_s : String
+    "LatLng: #{Math.degrees(@lat)},#{Math.degrees(@lng)}"
+  end
+
+  # Adds two LatLng objects.
+  def +(other : LatLng) : LatLng
+    self.class.new(@lat + other.lat, @lng + other.lng)
+  end
+
+  # Subtracts one LatLng from another.
+  def -(other : LatLng) : LatLng
+    self.class.new(@lat - other.lat, @lng - other.lng)
+  end
+
+  # Multiplies LatLng by a scalar.
+  def *(scalar : Float64) : LatLng
+    self.class.new(scalar * @lat, scalar * @lng)
+  end
+
+  # Static method to calculate latitude from a Point.
+  def self.latitude(point : Point) : Angle
+    Angle.from_radians(Math.atan2(point.z, Math.sqrt(point.x**2 + point.y**2)))
+  end
+
+  # Static method to calculate longitude from a Point.
+  def self.longitude(point : Point) : Angle
+    Angle.from_radians(Math.atan2(point.y, point.x))
+  end
+
+  # Returns the latitude as an Angle.
+  def lat : Angle
+    Angle.from_radians(@lat)
+  end
+
+  # Returns the longitude as an Angle.
+  def lng : Angle
+    Angle.from_radians(@lng)
+  end
+
+  # Checks if the LatLng is valid.
+  def valid? : Bool
+    @lat.abs <= Math::PI / 2 && @lng.abs <= Math::PI
+  end
+
+  # Converts LatLng to a Point.
+  def to_point : Point
+    phi = @lat
+    theta = @lng
+    cosphi = Math.cos(phi)
+    Point.new(Math.cos(theta) * cosphi, Math.sin(theta) * cosphi, Math.sin(phi))
+  end
+
+  # Normalizes the LatLng.
+  def normalized : LatLng
+    self.class.new(@lat.clamp(-Math::PI / 2, Math::PI / 2), Math.drem(@lng, 2 * Math::PI))
+  end
+
+  # Checks if two LatLng objects are approximately equal.
+  def approx_equals(other : LatLng, max_error : Float64 = 1e-15) : Bool
+    Math.abs(@lat - other.lat) < max_error && Math.abs(@lng - other.lng) < max_error
+  end
+
+  # Calculates the great-circle distance between two LatLng objects.
+  def get_distance(other : LatLng) : Angle
+    raise "Invalid LatLng" unless valid? && other.valid?
+
+    from_lat = @lat
+    to_lat = other.lat
+    from_lng = @lng
+    to_lng = other.lng
+    dlat = Math.sin(0.5 * (to_lat - from_lat))
+    dlng = Math.sin(0.5 * (to_lng - from_lng))
+    x = dlat**2 + dlng**2 * Math.cos(from_lat) * Math.cos(to_lat)
+    Angle.from_radians(2 * Math.atan2(Math.sqrt(x), Math.sqrt([0.0, 1.0 - x].max)))
+  end
+end
diff --git a/src/s2_cells/lat_lng_rect.cr b/src/s2_cells/lat_lng_rect.cr
new file mode 100644
index 0000000..c73007e
--- /dev/null
+++ b/src/s2_cells/lat_lng_rect.cr
@@ -0,0 +1,316 @@
+struct S2Cells::LatLngRect
+  getter lat : LineInterval
+  getter lng : SphereInterval
+
+  def initialize
+    @lat = LineInterval.empty
+    @lng = SphereInterval.empty
+  end
+
+  def initialize(lo : LatLng, hi : LatLng)
+    @lat = LineInterval.new(lo.lat.radians, hi.lat.radians)
+    @lng = SphereInterval.new(lo.lng.radians, hi.lng.radians)
+  end
+
+  def initialize(@lat : LineInterval, @lng : SphereInterval)
+  end
+
+  def ==(other : LatLngRect) : Bool
+    self.lat == other.lat && self.lng == other.lng
+  end
+
+  def to_s : String
+    "#{self.class.name}: #{self.lat}, #{self.lng}"
+  end
+
+  def lat_lo : Angle
+    Angle.from_radians(self.lat.lo)
+  end
+
+  def lat_hi : Angle
+    Angle.from_radians(self.lat.hi)
+  end
+
+  def lng_lo : Angle
+    Angle.from_radians(self.lng.lo)
+  end
+
+  def lng_hi : Angle
+    Angle.from_radians(self.lng.hi)
+  end
+
+  def lo : LatLng
+    LatLng.from_angles(self.lat_lo, self.lng_lo)
+  end
+
+  def hi : LatLng
+    LatLng.from_angles(self.lat_hi, self.lng_hi)
+  end
+
+  def self.from_center_size(center : LatLng, size : LatLng) : LatLngRect
+    self.from_point(center).expanded(0.5 * size)
+  end
+
+  def self.from_point(p : LatLng) : LatLngRect
+    raise "Invalid point" unless p.valid?
+    new(p, p)
+  end
+
+  def self.from_point_pair(a : LatLng, b : LatLng) : LatLngRect
+    raise "Invalid point" unless a.valid? && b.valid?
+    LatLngRect.new(
+      LineInterval.from_point_pair(a.lat.radians, b.lat.radians),
+      SphereInterval.from_point_pair(a.lng.radians, b.lng.radians)
+    )
+  end
+
+  def self.full_lat : LineInterval
+    LineInterval.new(-Math::PI / 2.0, Math::PI / 2.0)
+  end
+
+  def self.full_lng : SphereInterval
+    SphereInterval.full
+  end
+
+  def self.full : LatLngRect
+    new(self.full_lat, self.full_lng)
+  end
+
+  def full? : Bool
+    self.lat == self.class.full_lat && self.lng.full?
+  end
+
+  def valid? : Bool
+    (self.lat.lo.abs <= Math::PI / 2.0 &&
+      self.lat.hi.abs <= Math::PI / 2.0 &&
+      self.lng.valid? &&
+      self.lat.empty? == self.lng.empty?)
+  end
+
+  def self.empty : LatLngRect
+    new
+  end
+
+  def get_center : LatLng
+    LatLng.from_radians(self.lat.get_center, self.lng.get_center)
+  end
+
+  def get_size : LatLng
+    LatLng.from_radians(self.lat.get_length, self.lng.get_length)
+  end
+
+  def get_vertex(k : Int32) : LatLng
+    LatLng.from_radians(self.lat.bound(k >> 1), self.lng.bound((k >> 1) ^ (k & 1)))
+  end
+
+  def area : Float64
+    return 0.0 if self.empty?
+    self.lng.get_length * (Math.sin(self.lat_hi.radians) - Math.sin(self.lat_lo.radians)).abs
+  end
+
+  def empty? : Bool
+    self.lat.empty?
+  end
+
+  def point? : Bool
+    self.lat.lo == self.lat.hi && self.lng.lo == self.lng.hi
+  end
+
+  def convolve_with_cap(angle : Angle) : LatLngRect
+    cap = Cap.from_axis_angle(Point.new(1, 0, 0), angle)
+    r = self
+    4.times do |k|
+      vertex_cap = Cap.from_axis_height(self.get_vertex(k).to_point, cap.height)
+      r = r.union(vertex_cap.get_rect_bound)
+    end
+    r
+  end
+
+  def contains(other : LatLngRect | LatLng | Point | Cell) : Bool
+    case other
+    in Point
+      self.contains(LatLng.from_point(other))
+    in LatLng
+      raise "Invalid LatLng" unless other.valid?
+      self.lat.contains(other.lat.radians) && self.lng.contains(other.lng.radians)
+    in LatLngRect
+      self.lat.contains(other.lat) && self.lng.contains(other.lng)
+    in Cell
+      self.contains(other.get_rect_bound)
+    end
+  end
+
+  def interior_contains(other : LatLngRect | LatLng | Point) : Bool
+    case other
+    in Point
+      self.interior_contains(LatLng.new(other))
+    in LatLng
+      raise "Invalid LatLng" unless other.valid?
+      self.lat.interior_contains(other.lat.radians) && self.lng.interior_contains(other.lng.radians)
+    in LatLngRect
+      self.lat.interior_contains(other.lat) && self.lng.interior_contains(other.lng)
+    end
+  end
+
+  def may_intersect(cell : Cell) : Bool
+    self.intersects(cell.get_rect_bound)
+  end
+
+  def intersects(other : LatLngRect | Cell) : Bool
+    case other
+    in LatLngRect
+      self.lat.intersects(other.lat) && self.lng.intersects(other.lng)
+    in Cell
+      return false if self.empty?
+      return true if self.contains(other.get_center_raw)
+      return true if other.contains(self.get_center.to_point)
+      return false unless self.intersects(other.get_rect_bound)
+
+      cell_v = [] of Point
+      cell_ll = [] of LatLng
+      4.times do |i|
+        cell_v << other.get_vertex(i)
+        cell_ll << LatLng.from_point(cell_v[i])
+        return true if self.contains(cell_ll[i])
+        return true if other.contains(self.get_vertex(i).to_point)
+      end
+
+      4.times do |i|
+        edge_lng = SphereInterval.from_point_pair(cell_ll[i].lng.radians, cell_ll[(i + 1) & 3].lng.radians)
+        next unless self.lng.intersects(edge_lng)
+
+        a = cell_v[i]
+        b = cell_v[(i + 1) & 3]
+        if edge_lng.contains(self.lng.lo)
+          return true if self.class.intersects_lng_edge(a, b, self.lat, self.lng.lo)
+        end
+        if edge_lng.contains(self.lng.hi)
+          return true if self.class.intersects_lng_edge(a, b, self.lat, self.lng.hi)
+        end
+        return true if self.class.intersects_lat_edge(a, b, self.lat.lo, self.lng)
+        return true if self.class.intersects_lat_edge(a, b, self.lat.hi, self.lng)
+      end
+      false
+    end
+  end
+
+  def self.intersects_lng_edge(a : Point, b : Point, lat : LineInterval, lng : Float64) : Bool
+    simple_crossing(a, b, LatLng.from_radians(lat.lo, lng).to_point, LatLng.from_radians(lat.hi, lng).to_point)
+  end
+
+  def self.simple_crossing(a : Point, b : Point, c : Point, d : Point) : Bool
+    ab = a.cross_prod(b)
+    acb = -(ab.dot_prod(c))
+    bda = ab.dot_prod(d)
+    return false if acb * bda <= 0
+
+    cd = c.cross_prod(d)
+    cbd = -(cd.dot_prod(b))
+    dac = cd.dot_prod(a)
+    (acb * cbd > 0) && (acb * dac > 0)
+  end
+
+  def self.intersects_lat_edge(a : Point, b : Point, lat : Float64, lng : SphereInterval) : Bool
+    raise "Point not unit length" unless S2Cells.is_unit_length(a) && S2Cells.is_unit_length(b)
+
+    z = robust_cross_prod(a, b).normalize
+    z = -z if z.z < 0
+
+    y = robust_cross_prod(z, Point.new(0, 0, 1)).normalize
+    x = y.cross_prod(z)
+    raise "Point not unit length" unless S2Cells.is_unit_length(x)
+    raise "Invalid X value" unless x.z >= 0
+
+    sin_lat = Math.sin(lat)
+    return false if sin_lat.abs >= x.z
+
+    cos_theta = sin_lat / x.z
+    sin_theta = Math.sqrt(1 - cos_theta * cos_theta)
+    theta = Math.atan2(sin_theta, cos_theta)
+
+    ab_theta = SphereInterval.from_point_pair(
+      Math.atan2(a.dot_prod(y), a.dot_prod(x)),
+      Math.atan2(b.dot_prod(y), b.dot_prod(x))
+    )
+
+    if ab_theta.contains(theta)
+      isect = x * cos_theta + y * sin_theta
+      return true if lng.contains(Math.atan2(isect.y, isect.x))
+    end
+    if ab_theta.contains(-theta)
+      isect = x * cos_theta - y * sin_theta
+      return true if lng.contains(Math.atan2(isect.y, isect.x))
+    end
+    false
+  end
+
+  def interior_intersects(other : LatLngRect) : Bool
+    self.lat.interior_intersects(other.lat) && self.lng.interior_intersects(other.lng)
+  end
+
+  def union(other : LatLngRect) : LatLngRect
+    LatLngRect.new(self.lat.union(other.lat), self.lng.union(other.lng))
+  end
+
+  def intersection(other : LatLngRect) : LatLngRect
+    lat = self.lat.intersection(other.lat)
+    lng = self.lng.intersection(other.lng)
+    return LatLngRect.empty if lat.empty? || lng.empty?
+    LatLngRect.new(lat, lng)
+  end
+
+  def expanded(margin : LatLng) : LatLngRect
+    raise "Invalid margin" unless margin.lat.radians > 0 && margin.lng.radians > 0
+    LatLngRect.new(
+      self.lat.expanded(margin.lat.radians).intersection(self.full_lat),
+      self.lng.expanded(margin.lng.radians)
+    )
+  end
+
+  def approx_equals(other : LatLngRect, max_error = 1e-15) : Bool
+    self.lat.approx_equals(other.lat, max_error) && self.lng.approx_equals(other.lng, max_error)
+  end
+
+  def get_cap_bound : Cap
+    return Cap.empty if self.empty?
+
+    pole_z, pole_angle = if (self.lat.lo + self.lat.hi) < 0
+                           {-1.0, Math::PI / 2.0 + self.lat.hi}
+                         else
+                           {1.0, Math::PI / 2.0 - self.lat.lo}
+                         end
+
+    pole_cap = Cap.from_axis_angle(Point.new(0.0, 0.0, pole_z), Angle.from_radians(pole_angle))
+    lng_span = self.lng.hi - self.lng.lo
+
+    if (lng_span % (2 * Math::PI)) >= 0 && lng_span < (2 * Math::PI)
+      mid_cap = Cap.from_axis_angle(self.get_center.to_point, Angle.from_radians(0))
+      4.times { |k| mid_cap.add_point(self.get_vertex(k).to_point) }
+      return mid_cap if mid_cap.height < pole_cap.height
+    end
+    pole_cap
+  end
+
+  def self.ortho(a : Point) : Point
+    k = a.largest_abs_component - 1
+    k = 2 if k < 0
+    temp = case k
+           when 0
+             Point.new(1, 0.0053, 0.00457)
+           when 1
+             Point.new(0.012, 1, 0.00457)
+           else
+             Point.new(0.012, 0.0053, 1)
+           end
+    a.cross_prod(temp).normalize
+  end
+
+  def self.robust_cross_prod(a : Point, b : Point) : Point
+    raise "Input must be unit length" unless S2Cells.is_unit_length(a) && S2Cells.is_unit_length(b)
+
+    x = (b + a).cross_prod(b - a)
+    return x unless x == Point.new(0, 0, 0)
+
+    ortho(a)
+  end
+end
diff --git a/src/s2_cells/lat_lon.cr b/src/s2_cells/lat_lon.cr
deleted file mode 100644
index 3027d22..0000000
--- a/src/s2_cells/lat_lon.cr
+++ /dev/null
@@ -1,22 +0,0 @@
-class S2Cells::LatLon
-  def initialize(lat_degrees : Float64, lon_degrees : Float64)
-    @lat = lat_degrees * Math::PI / 180.0
-    @lon = lon_degrees * Math::PI / 180.0
-  end
-
-  @lat : Float64
-  @lon : Float64
-
-  def to_point
-    phi = @lat
-    theta = @lon
-    cosphi = Math.cos(phi)
-    Point.new(Math.cos(theta) * cosphi, Math.sin(theta) * cosphi, Math.sin(phi))
-  end
-
-  def to_token(level = 30)
-    CellId.from_point(to_point)
-      .parent(level)
-      .to_token
-  end
-end
diff --git a/src/s2_cells/metric.cr b/src/s2_cells/metric.cr
new file mode 100644
index 0000000..056ff06
--- /dev/null
+++ b/src/s2_cells/metric.cr
@@ -0,0 +1,85 @@
+require "math"
+
+module S2Cells
+  abstract struct Metric
+    def initialize(@deriv : Float64, @dim : Int32)
+    end
+
+    def deriv : Float64
+      @deriv
+    end
+
+    def get_value(level : Int) : Float64
+      Math.ldexp(@deriv, -@dim * level)
+    end
+
+    def get_closest_level(value : Float64) : Int32
+      factor = if @dim == 1
+                 Math.sqrt(2.0)
+               else
+                 2.0
+               end
+      get_min_level(factor * value)
+    end
+
+    def get_min_level(value : Float64) : Int32
+      return CellId::MAX_LEVEL if value <= 0
+
+      m, x = Math.frexp(value / @deriv)
+      level = {0, {CellId::MAX_LEVEL, -((x - 1) >> (@dim - 1))}.min}.max
+      raise "Invalid level" unless level == CellId::MAX_LEVEL || get_value(level) <= value
+      raise "Invalid level" unless level == 0 || get_value(level - 1) > value
+      level
+    end
+
+    def get_max_level(value : Float64) : Int32
+      return CellId::MAX_LEVEL if value <= 0
+
+      m, x = Math.frexp(@deriv / value)
+      level = {0, {CellId::MAX_LEVEL, (x - 1) >> (@dim - 1)}.min}.max
+      raise "Invalid level" unless level == 0 || get_value(level) >= value
+      raise "Invalid level" unless level == CellId::MAX_LEVEL || get_value(level + 1) < value
+      level
+    end
+
+    private def max(a : Int32, b : Int32) : Int32
+      a > b ? a : b
+    end
+
+    private def min(a : Int32, b : Int32) : Int32
+      a < b ? a : b
+    end
+  end
+
+  struct LengthMetric < Metric
+    def initialize(deriv : Float64)
+      super(deriv, 1)
+    end
+  end
+
+  struct AreaMetric < Metric
+    def initialize(deriv : Float64)
+      super(deriv, 2)
+    end
+  end
+
+  AVG_ANGLE_SPAN = LengthMetric.new(Math::PI / 2)         # true for all projections
+  MIN_ANGLE_SPAN = LengthMetric.new(4.0 / 3.0)            # quadratic projection
+  MAX_ANGLE_SPAN = LengthMetric.new(1.704897179199218452) # quadratic projection
+
+  AVG_EDGE = LengthMetric.new(1.459213746386106062)     # quadratic projection
+  MIN_EDGE = LengthMetric.new(2 * Math.sqrt(2.0) / 3.0) # quadratic projection
+  MAX_EDGE = LengthMetric.new(MAX_ANGLE_SPAN.deriv)     # true for all projections
+
+  AVG_DIAG = LengthMetric.new(2.060422738998471683)     # quadratic projection
+  MIN_DIAG = LengthMetric.new(8 * Math.sqrt(2.0) / 9.0) # quadratic projection
+  MAX_DIAG = LengthMetric.new(2.438654594434021032)     # quadratic projection
+
+  AVG_WIDTH = LengthMetric.new(1.434523672886099389)     # quadratic projection
+  MIN_WIDTH = LengthMetric.new(2 * Math.sqrt(2.0) / 3.0) # quadratic projection
+  MAX_WIDTH = LengthMetric.new(MAX_ANGLE_SPAN.deriv)     # true for all projections
+
+  AVG_AREA = AreaMetric.new(4 * Math::PI / 6.0)       # Average cell area for all projections
+  MIN_AREA = AreaMetric.new(8 * Math.sqrt(2.0) / 9.0) # Minimum cell area for quadratic projections
+  MAX_AREA = AreaMetric.new(2.635799256963161491)     # Maximum cell area for quadratic projections
+end
diff --git a/src/s2_cells/point.cr b/src/s2_cells/point.cr
index e4390e0..5c09b44 100644
--- a/src/s2_cells/point.cr
+++ b/src/s2_cells/point.cr
@@ -1,26 +1,112 @@
-class S2Cells::Point
-  getter x : Float64
-  getter y : Float64
-  getter z : Float64
+struct S2Cells::Point
+  @point : Tuple(Float64, Float64, Float64)
 
-  def initialize(@x, @y, @z)
+  def initialize(x : Float64, y : Float64, z : Float64)
+    @point = {x, y, z}
+  end
+
+  def [](index : Int)
+    @point[index]
+  end
+
+  def -
+    self.class.new(-@point[0], -@point[1], -@point[2])
+  end
+
+  def ==(other : Point)
+    @point == other.@point
+  end
+
+  def hash
+    @point.hash
+  end
+
+  def to_s
+    "Point: #{@point}"
+  end
+
+  def +(other : Point)
+    self.class.new(@point[0] + other[0],
+      @point[1] + other[1],
+      @point[2] + other[2])
+  end
+
+  def -(other : Point)
+    self.class.new(@point[0] - other[0],
+      @point[1] - other[1],
+      @point[2] - other[2])
+  end
+
+  def *(other : Number)
+    self.class.new(@point[0] * other,
+      @point[1] * other,
+      @point[2] * other)
+  end
+
+  def x : Float64
+    @point[0]
+  end
+
+  def y : Float64
+    @point[1]
+  end
+
+  def z : Float64
+    @point[2]
   end
 
   def abs
-    {@x.abs, @y.abs, @z.abs}
+    self.class.new(@point[0].abs,
+      @point[1].abs,
+      @point[2].abs)
   end
 
-  def largest_abs_component : UInt64
+  def largest_abs_component
     temp = abs
-
     if temp[0] > temp[1]
-      temp[0] > temp[2] ? 0_u64 : 2_u64
+      temp[0] > temp[2] ? 0 : 2
     else
-      temp[1] > temp[2] ? 1_u64 : 2_u64
+      temp[1] > temp[2] ? 1 : 2
     end
   end
 
-  def dot_prod(o : Point)
-    @x * o.x + @y * o.y + @z * o.z
+  def angle(other : Point)
+    Math.atan2(cross_prod(other).norm, dot_prod(other))
+  end
+
+  def cross_prod(other : Point)
+    x, y, z = @point
+    ox, oy, oz = other.@point
+    self.class.new(y * oz - z * oy,
+      z * ox - x * oz,
+      x * oy - y * ox)
+  end
+
+  def dot_prod(other : Point)
+    x, y, z = @point
+    ox, oy, oz = other.@point
+    x * ox + y * oy + z * oz
+  end
+
+  def norm2
+    x, y, z = @point
+    x * x + y * y + z * z
+  end
+
+  def norm
+    Math.sqrt(norm2)
+  end
+
+  def normalize
+    n = norm
+    n = 1.0 / n if n != 0
+    self.class.new(@point[0] * n, @point[1] * n, @point[2] * n)
+  end
+end
+
+# Define * for Number * Point
+struct Number
+  def *(point : Point)
+    Point.new(self * point[0], self * point[1], self * point[2])
   end
 end
diff --git a/src/s2_cells/region_coverer.cr b/src/s2_cells/region_coverer.cr
new file mode 100644
index 0000000..4d50501
--- /dev/null
+++ b/src/s2_cells/region_coverer.cr
@@ -0,0 +1,230 @@
+require "priority-queue"
+
+struct S2Cells::RegionCoverer
+  struct Candidate
+    getter children : Array(Candidate)
+    property cell : Cell
+    property? terminal : Bool
+
+    def initialize(@cell, @terminal)
+      @children = [] of Candidate
+    end
+
+    def num_children : Int32
+      @children.size
+    end
+
+    def <(other : Candidate) : Bool
+      my_cell = @cell
+      other_cell = other.cell
+      raise "NotImplementedError" unless my_cell && other_cell
+      my_cell.id < other_cell.id
+    end
+  end
+
+  getter min_level : Int32
+  getter max_level : Int32
+  getter level_mod : Int32
+  getter max_cells : Int32
+
+  getter! region : LatLngRect
+
+  @result : Array(CellId)
+  @pq : Priority::Queue(Candidate)
+
+  def initialize(
+    @min_level = 0,
+    @max_level = CellId::MAX_LEVEL,
+    @level_mod = 1,
+    @max_cells = 8,
+    @region = nil,
+    @result = [] of CellId,
+    @pq = Priority::Queue(Candidate).new
+  )
+  end
+
+  def min_level=(value : Int32)
+    raise "Invalid min_level" unless value >= 0 && value <= CellId::MAX_LEVEL
+    @min_level = value
+  end
+
+  def max_level=(value : Int32)
+    raise "Invalid max_level" unless value >= 0 && value <= CellId::MAX_LEVEL
+    @max_level = value
+  end
+
+  def level_mod=(value : Int32)
+    raise "Invalid level_mod" unless value >= 1 && value <= 3
+    @level_mod = value
+  end
+
+  def max_cells=(value : Int32)
+    @max_cells = value
+  end
+
+  def get_covering(region : LatLngRect)
+    @result.clear
+    tmp_union = __get_cell_union(region)
+    tmp_union.denormalize(@min_level, @level_mod)
+  end
+
+  def get_interior_covering(region : LatLngRect)
+    @result.clear
+    tmp_union = __get_interior_cell_union(region)
+    tmp_union.denormalize(@min_level, @level_mod)
+  end
+
+  private def __new_candidate(cell : Cell) : Candidate?
+    return nil unless region.may_intersect(cell)
+    is_terminal = false
+    if cell.level >= @min_level
+      if @interior_covering
+        if region.contains(cell)
+          is_terminal = true
+        elsif cell.level + @level_mod > @max_level
+          return nil
+        end
+      else
+        if cell.level + @level_mod > @max_level || region.contains(cell)
+          is_terminal = true
+        end
+      end
+    end
+
+    Candidate.new cell, is_terminal
+  end
+
+  private def __max_children_shift : Int32
+    2 * @level_mod
+  end
+
+  private def __expand_children(candidate : Candidate, cell : Cell, num_levels : Int32) : Int32
+    num_levels -= 1
+    num_terminals = 0
+    cell.subdivide do |child_cell|
+      if num_levels > 0
+        if region.may_intersect(child_cell)
+          num_terminals += __expand_children(candidate, child_cell, num_levels)
+        end
+        next
+      end
+      child = __new_candidate(child_cell)
+      if child
+        candidate.children << child
+        num_terminals += 1 if child.terminal?
+      end
+    end
+    num_terminals
+  end
+
+  private def __add_candidate(candidate : Candidate?)
+    return unless candidate
+
+    if candidate.terminal?
+      @result << candidate.cell.id
+      return
+    end
+
+    raise "Candidate should have no children" if candidate.num_children > 0
+
+    num_levels = @level_mod
+    num_levels = 1 if candidate.cell.level < @min_level
+    num_terminals = __expand_children(candidate, candidate.cell, num_levels)
+
+    if candidate.num_children == 0
+      # Not needed due to GC
+    elsif !@interior_covering && num_terminals == (1 << __max_children_shift) && candidate.cell.level >= @min_level
+      candidate.terminal = true
+      __add_candidate(candidate)
+    else
+      priority = ((candidate.cell.level << __max_children_shift) + candidate.num_children) << __max_children_shift + num_terminals
+      @pq.push priority, candidate
+    end
+  end
+
+  private def __get_initial_candidates
+    if @max_cells >= 4
+      cap = region.get_cap_bound
+      level = {CellId.min_width.get_max_level(2 * cap.angle.radians), {@max_level, CellId::MAX_LEVEL - 1}.min}.min
+
+      if @level_mod > 1 && level > @min_level
+        level -= (level - @min_level) % @level_mod
+      end
+
+      if level > 0
+        cell_id = CellId.from_point(cap.axis)
+        vertex_neighbors = cell_id.get_vertex_neighbors(level)
+        vertex_neighbors.each do |neighbor|
+          __add_candidate(__new_candidate(Cell.new(neighbor)))
+        end
+        return
+      end
+    end
+
+    6.times do |face|
+      __add_candidate(__new_candidate(FACE_CELLS[face]))
+    end
+  end
+
+  private def __get_covering(region : LatLngRect)
+    raise "Priority queue not empty" unless @pq.empty?
+    raise "Result not empty" unless @result.empty?
+    @region = region
+
+    __get_initial_candidates
+    while !@pq.empty? && (!@interior_covering || @result.size < @max_cells)
+      candidate = @pq.shift.value
+
+      result_size = @interior_covering ? 0 : @pq.size
+      if candidate.cell.level < @min_level || candidate.num_children == 1 || @result.size + result_size + candidate.num_children <= @max_cells
+        candidate.children.each do |child|
+          __add_candidate(child)
+        end
+      elsif @interior_covering
+        # Do nothing here
+      else
+        candidate.terminal = true
+        __add_candidate(candidate)
+      end
+    end
+
+    @pq.clear
+    @region = nil
+  end
+
+  private def __get_cell_union(region) : CellUnion
+    @interior_covering = false
+    __get_covering(region)
+    CellUnion.new(@result)
+  end
+
+  private def __get_interior_cell_union(region) : CellUnion
+    @interior_covering = true
+    __get_covering(region)
+    CellUnion.new(@result)
+  end
+
+  def self.flood_fill(region, start)
+    all_nbrs = Set(CellId).new
+    frontier = [] of CellId
+    all_nbrs.add(start)
+    frontier << start
+    while !frontier.empty?
+      cell_id = frontier.pop
+      next unless region.may_intersect(Cell.new(cell_id))
+      yield cell_id
+
+      neighbors = cell_id.get_edge_neighbors
+      neighbors.each do |nbr|
+        if !all_nbrs.includes?(nbr)
+          all_nbrs.add(nbr)
+          frontier << nbr
+        end
+      end
+    end
+  end
+
+  def self.get_simple_covering(region, start, level)
+    flood_fill(region, CellId.from_point(start).parent(level))
+  end
+end