diff --git a/experimental/FTheoryTools/src/FTheoryTools.jl b/experimental/FTheoryTools/src/FTheoryTools.jl
index f345d829f8c..a14867a89f7 100644
--- a/experimental/FTheoryTools/src/FTheoryTools.jl
+++ b/experimental/FTheoryTools/src/FTheoryTools.jl
@@ -30,6 +30,7 @@ include("G4Fluxes/constructors.jl")
 include("G4Fluxes/attributes.jl")
 include("G4Fluxes/properties.jl")
 include("G4Fluxes/special_attributes.jl")
+include("G4Fluxes/special-intersection-theory.jl")
 
 include("Serialization/tate_models.jl")
 include("Serialization/weierstrass_models.jl")
diff --git a/experimental/FTheoryTools/src/G4Fluxes/special-intersection-theory.jl b/experimental/FTheoryTools/src/G4Fluxes/special-intersection-theory.jl
new file mode 100644
index 00000000000..daa97e4f6ce
--- /dev/null
+++ b/experimental/FTheoryTools/src/G4Fluxes/special-intersection-theory.jl
@@ -0,0 +1,299 @@
+# ---------------------------------------------------------------------------------------------------------
+# (1) Compute the intersection product of an algebraic cycle with a hypersurface.
+# ---------------------------------------------------------------------------------------------------------
+
+function sophisticated_intersection_product(v::NormalToricVariety, indices::NTuple{4, Int64}, hypersurface_equation::MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}, inter_dict::Dict{NTuple{4, Int64}, ZZRingElem}, s_inter_dict::Dict{String, ZZRingElem}, data::NamedTuple)
+
+  # (A) Have we computed this intersection number in the past? If so, just use that result...
+  if haskey(inter_dict, indices)
+    return inter_dict[indices]
+  end
+  
+  # (B) Get the indices of the variables that we try to intersect and, by virtue of the SR-ideal, intersect trivially
+  for sr_set in data.sr_ideal_pos
+    if is_subset(sr_set, indices)
+      inter_dict[indices] = ZZ(0)
+      return ZZ(0)
+    end
+  end
+
+  # (C) Deal with self-intersection and should-never-happen case.
+  variable_pos = Set(indices)
+  if length(variable_pos) < 4 && length(variable_pos) >= 1
+    return intersection_from_equivalent_cycle(v, indices, hypersurface_equation, inter_dict, s_inter_dict, data)
+  end
+  if length(variable_pos) == 0
+    println("WEIRD! THIS SHOULD NEVER HAPPEN! INFORM THE AUTHORS!")
+    println("")
+  end
+
+
+  # (D) Deal with transverse intersection...
+
+  # D.1 Work out the intersection locus in detail.
+  pt_reduced, gs_reduced, remaining_vars, reduced_scaling_relations = Oscar._reduce_hypersurface_equation(v, hypersurface_equation, indices, data)
+
+  # D.2 If pt == 0, then we are not looking at a transverse intersection. So take an equivalent cycle and try again...
+  if is_zero(pt_reduced)
+    return intersection_from_equivalent_cycle(v, indices, hypersurface_equation, inter_dict, s_inter_dict, data)
+  end
+
+  # D.3 If pt is constant and non-zero, then the intersection is trivial.
+  if is_constant(pt_reduced) && is_zero(pt_reduced) == false
+    inter_dict[indices] = ZZ(0)
+    return ZZ(0)
+  end
+
+  # D.4 Helper function for the cases below
+  function has_one_and_rest_zero(vec)
+    return count(==(1), vec) == 1 && all(x -> x == 0 || x == 1, vec)
+  end
+
+  # C.5 Cover a case that seems to appear frequently for our investigation:
+  # pt_reduced of the form a * x + b * y for non-zero number a,b and remaining variables x, y subject to a reduced SR generator x * y and scaling relation [1,1].
+  # This will thus always give exactly one solution (x = 1, y = -a/b), and so the intersection number is one.
+  if length(gs_reduced) == 1 && length(remaining_vars) == 2
+    mons_list = collect(monomials(pt_reduced))
+    if length(mons_list) == 2
+      if all(x -> x != 0, collect(coefficients(pt_reduced)))
+        exps_list = [collect(exponents(k))[1] for k in mons_list]
+        if has_one_and_rest_zero(exps_list[1]) && has_one_and_rest_zero(exps_list[2])
+          if gs_reduced[1] == remaining_vars[1] * remaining_vars[2]
+            if reduced_scaling_relations == matrix(ZZ, [[1,1]])
+              inter_dict[indices] = ZZ(1)
+              return ZZ(1)
+            end
+          end
+        end
+      end
+    end
+  end
+
+  # C.6 Cover a case that seems to appear frequently for our investigation:
+  # pt_reduced of the form a * x for non-zero number a and remaining variables x, y subject to a reduced SR generator x * y and scaling relation [*, != 0].
+  # This only gives the solution [0:1], so one intersection point.
+  if length(gs_reduced) == 1 && length(remaining_vars) == 2
+    mons_list = collect(monomials(pt_reduced))
+    if length(mons_list) == 1 && collect(coefficients(pt_reduced))[1] != 0
+      list_of_exps = collect(exponents(mons_list[1]))[1]
+      number_of_zeros = count(==(0), list_of_exps)
+      if number_of_zeros == length(list_of_exps) - 1
+        if gs_reduced[1] == remaining_vars[1] * remaining_vars[2]
+          if reduced_scaling_relations[1,1] != 0 && reduced_scaling_relations[1,2] != 0
+            highest_power = list_of_exps[findfirst(x -> x > 0, list_of_exps)]
+            if highest_power == 1
+              inter_dict[indices] = highest_power
+              return highest_power
+            end
+          end
+        end
+      end
+    end
+  end
+    
+  # C.7 Cover a case that seems to appear frequently for our investigation. It looks as follows:
+  # pt_reduced = -5700*w8*w10
+  # remaining_vars = MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}[w8, w10]
+  # gs_reduced = MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}[w8*w10]
+  # reduced_scaling_relations = [1 1]
+  # This gives exactly two solutions, namely [0:1] and [1:0].
+  if length(gs_reduced) == 1 && length(remaining_vars) == 2
+    mons_list = collect(monomials(pt_reduced))
+    if length(mons_list) == 1 && mons_list[1] == remaining_vars[1] * remaining_vars[2]
+      if gs_reduced[1] == remaining_vars[1] * remaining_vars[2]
+        if reduced_scaling_relations == matrix(ZZ, [[1,1]])
+          inter_dict[indices] = 2
+          return 2
+        end
+      end
+    end
+  end
+
+  # C.8 Check if this was covered in our special cases
+  if haskey(s_inter_dict, string([pt_reduced, gs_reduced, remaining_vars, reduced_scaling_relations]))
+    numb = s_inter_dict[string([pt_reduced, gs_reduced, remaining_vars, reduced_scaling_relations])]
+    inter_dict[indices] = numb
+    return numb
+  end
+
+  # C.9 In all other cases, proceed via a rationally equivalent cycle
+  println("")
+  println("FOUND CASE THAT CANNOT YET BE DECIDED!")
+  println("$pt_reduced")
+  println("$remaining_vars")
+  println("$gs_reduced")
+  println("$indices")
+  println("$reduced_scaling_relations")
+  println("TRYING WITH EQUIVALENT CYCLE")
+  println("")
+  numb = intersection_from_equivalent_cycle(v, indices, hypersurface_equation, inter_dict, s_inter_dict, data)
+  s_inter_dict[string([pt_reduced, gs_reduced, remaining_vars, reduced_scaling_relations])] = numb
+  return numb
+
+end
+
+
+
+# ---------------------------------------------------------------------------------------------------------
+# (2) Compute the intersection product from a rationally equivalent cycle.
+# ---------------------------------------------------------------------------------------------------------
+
+function intersection_from_equivalent_cycle(v::NormalToricVariety, indices::NTuple{4, Int64}, hypersurface_equation::MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}, inter_dict::Dict{NTuple{4, Int64}, ZZRingElem}, s_inter_dict::Dict{String, ZZRingElem}, data::NamedTuple)
+  coeffs_list, tuple_list = Oscar._rationally_equivalent_cycle(v, indices, data)
+  intersect_numb = 0
+  for k in 1:length(tuple_list)
+    if !is_zero(coeffs_list[k])
+      intersect_numb += coeffs_list[k] * sophisticated_intersection_product(v, tuple_list[k], hypersurface_equation, inter_dict, s_inter_dict, data)
+    end
+  end
+  @req is_integer(intersect_numb) "Should have expected to find only integer intersection numbers, but found $intersect_numb for $indices"
+  inter_dict[indices] = ZZ(intersect_numb)
+  return ZZ(intersect_numb)
+end
+
+
+
+# ---------------------------------------------------------------------------------------------------------
+# (3) A function to reduce the hypersurface polynomial to {xi = 0} with i in indices
+# ---------------------------------------------------------------------------------------------------------
+
+function _reduce_hypersurface_equation(v::NormalToricVariety, p_hyper::MPolyRingElem, indices::NTuple{4, Int64}, data::NamedTuple)
+
+  # Set variables to zero in the hypersurface equation
+  vanishing_vars_pos = unique(indices)  
+  new_p_hyper = divrem(p_hyper, data.gS[vanishing_vars_pos[1]])[2]
+  for m in 2:length(vanishing_vars_pos)
+    new_p_hyper = divrem(new_p_hyper, data.gS[vanishing_vars_pos[m]])[2]
+  end
+
+  # Is the resulting polynomial constant?
+  if is_constant(new_p_hyper)
+    return [new_p_hyper, [], [], zero_matrix(ZZ, 0, 0)]
+  end
+
+  # Identify the remaining variables
+  remaining_vars_pos = Set(1:length(data.gS))
+  for my_exps in data.sr_ideal_pos
+    len_my_exps = length(my_exps)
+    inter_len = count(idx -> idx in vanishing_vars_pos, my_exps)
+    if len_my_exps == inter_len + 1
+      delete!(remaining_vars_pos, my_exps[findfirst(idx -> !(idx in vanishing_vars_pos), my_exps)])
+    end
+  end
+  set_to_one_list = sort([k for k in 1:length(data.gS) if k ∉ remaining_vars_pos])
+  remaining_vars_pos = setdiff(collect(remaining_vars_pos), vanishing_vars_pos)
+  remaining_vars = [data.gS[k] for k in remaining_vars_pos]
+
+  # Extract remaining Stanley-Reisner ideal relations
+  sr_reduced = Vector{MPolyDecRingElem{QQFieldElem, QQMPolyRingElem}}()
+  for k in 1:length(data.sr_ideal_pos)
+    if isdisjoint(set_to_one_list, data.sr_ideal_pos[k])
+      #for m in data.sr_ideal_pos[k] but not in vanishing_vars_pos
+      push!(sr_reduced, prod(data.gS[m] for m in data.sr_ideal_pos[k] if m ∉ vanishing_vars_pos))
+      #push!(sr_reduced, prod(data.gS[m] for m in data.sr_ideal_pos[k]))
+    end
+  end
+
+  # Identify the remaining scaling relations
+  prepared_scaling_relations = zero_matrix(ZZ, length(data.scalings[1]), length(set_to_one_list) + length(remaining_vars_pos))
+  for k in 1:(length(set_to_one_list) + length(remaining_vars_pos))
+    col = k <= length(set_to_one_list) ? data.scalings[set_to_one_list[k]] : data.scalings[remaining_vars_pos[k - length(set_to_one_list)]]
+    for l in 1:length(col)
+      prepared_scaling_relations[l, k] = col[l]
+    end
+  end
+  prepared_scaling_relations = hnf(prepared_scaling_relations)
+  reduced_scaling_relations = prepared_scaling_relations[length(set_to_one_list) + 1: nrows(prepared_scaling_relations), length(set_to_one_list) + 1 : ncols(prepared_scaling_relations)]
+
+  # Identify the final form of the reduced hypersurface equation, by setting all variables to one that we can
+  images = [k in remaining_vars_pos ? data.gS[k] : one(data.S) for k in 1:length(data.gS)]
+  pt_reduced = evaluate(new_p_hyper, images)
+
+  # Return the result
+  return [pt_reduced, sr_reduced, remaining_vars, reduced_scaling_relations]
+end
+
+
+
+# ---------------------------------------------------------------------------------------------------------
+# (4) A function to find a rationally equivalent algebraic cycle.
+# ---------------------------------------------------------------------------------------------------------
+
+function _rationally_equivalent_cycle(v::NormalToricVariety, indices::NTuple{4, Int64}, data::NamedTuple)
+
+  # Identify positions of the single and triple variable
+  power_variable = indices[rand(1:length(indices))]
+  for k in Set(indices)
+    if count(==(k), indices) > 1
+      power_variable = k
+      break
+    end
+  end
+  other_variables = filter(!=(power_variable), indices)
+  pos_power_variable = findfirst(==(power_variable), indices)
+
+  # Let us simplify the problem by extracting the entries in the columns of single_variables and double_variables of the linear relation matrix
+  simpler_matrix = data.linear_relations[[other_variables..., power_variable], :]
+  b = zero_matrix(QQ, length(other_variables) + 1, 1)
+  b[nrows(b), 1] = 1
+  A = solve(simpler_matrix, b; side =:right)
+  
+  # Now form the relation in case...
+  employed_relation = -sum((data.linear_relations[:, k] .* A[k]) for k in 1:5)
+  employed_relation[power_variable] = 0
+
+  # Populate `coeffs` and `tuples` that form rationally equivalent algebraic cycle
+  coeffs = Vector{QQFieldElem}()
+  tuples = Vector{NTuple{4, Int64}}()
+  for k in 1:length(employed_relation)
+    if employed_relation[k] != 0
+
+      # Form the new tuple
+      new_tuple = (indices[1:pos_power_variable-1]..., k, indices[pos_power_variable+1:end]...)
+      new_tuple = Tuple(sort(collect(new_tuple)))
+
+      pos = findfirst(==(new_tuple), tuples)
+      if pos !== nothing
+        coeffs[pos] += employed_relation[k]
+      else
+        push!(coeffs, employed_relation[k])
+        push!(tuples, Tuple(new_tuple))
+      end
+
+      #=
+      # If there are no repeated entries, just add this tuple...
+      if length(new_tuple) == length(Set(new_tuple))
+        
+        pos = findfirst(==(new_tuple), tuples)
+        if pos !== nothing
+          coeffs[pos] += employed_relation[k]
+        else
+          push!(coeffs, employed_relation[k])
+          push!(tuples, Tuple(new_tuple))
+        end
+      
+      # Otherwise, replace with a rationally equivalent cycle.
+      else
+
+        new_coeffs, new_tuples = Oscar._rationally_equivalent_cycle(v, new_tuple, data)
+        for l in 1:length(new_coeffs)
+
+          pos = findfirst(==(new_tuples[l]), tuples)
+          if pos !== nothing
+            coeffs[pos] += employed_relation[k] * new_coeffs[l]
+          else
+            push!(coeffs, employed_relation[k] * new_coeffs[l])
+            push!(tuples, new_tuples[l])
+          end
+
+        end
+
+      end
+      =#
+
+    end
+  end
+
+  return [coeffs, tuples]
+
+end
diff --git a/experimental/FTheoryTools/src/G4Fluxes/special_attributes.jl b/experimental/FTheoryTools/src/G4Fluxes/special_attributes.jl
index 49e20bc07bc..7a6656bfb13 100644
--- a/experimental/FTheoryTools/src/G4Fluxes/special_attributes.jl
+++ b/experimental/FTheoryTools/src/G4Fluxes/special_attributes.jl
@@ -262,3 +262,225 @@ function ambient_space_models_of_g4_fluxes(m::AbstractFTheoryModel; check::Bool
   return filtered_h22_basis
 
 end
+
+
+# ---------------------------------------------------------------------------------------------------------
+# (2) Internal helper function.
+# (2) Internal helper function.
+# ---------------------------------------------------------------------------------------------------------
+
+# The following is an internal function, that is being used to identify all well-quantized G4-fluxes.
+# Let G4 in H^(2,2)(toric_ambient_space) a G4-flux ambient space candidate, i.e. the physically
+# truly relevant quantity is the restriction of G4 to a hypersurface V(pt) in the toric_ambient_space.
+# To tell if this candidate is well-quantized, we execute a necessary check, namely we verify that
+# the integral of G4 * [pt] * [d1] * [d2] over X_Sigma is an integer for any two toric divisors d1, d2.
+# While we wish to execute this test for any two toric divisors d1, d2 some pairs can be ignored.
+# Say, because their intersection locus is trivial because of the SR-ideal, or because their intersection
+# has empty intersection with the hypersurface. The following method identifies the remaining pairs of
+# toric divisors d1, d2 that we must consider.
+
+function _ambient_space_divisors_to_be_considered(m::AbstractFTheoryModel)::Vector{Tuple{Int64, Int64}}
+
+  if has_attribute(m, :_ambient_space_divisors_to_be_considered)
+    return get_attribute(m, :_ambient_space_divisors_to_be_considered)
+  end
+
+  gS = gens(cox_ring(ambient_space(m)))
+  mnf = Oscar._minimal_nonfaces(ambient_space(m))
+  ignored_sets = Set([Tuple(sort(Vector{Int}(Polymake.row(mnf, i)))) for i in 1:Polymake.nrows(mnf)])
+
+  list_of_elements = Vector{Tuple{Int64, Int64}}()
+  for k in 1:n_rays(ambient_space(m))
+    for l in k:n_rays(ambient_space(m))
+
+      # V(x_k, x_l) = emptyset?
+      (k,l) in ignored_sets && continue
+
+      # Simplify the hypersurface polynomial by setting relevant variables to zero.
+      # If all coefficients of this new polynomial add to sum, then we keep this generator.
+      new_p_hyper = divrem(hypersurface_equation(m), gS[k])[2]
+      if k != l
+        new_p_hyper = divrem(new_p_hyper, gS[l])[2]
+      end
+      if sum(coefficients(new_p_hyper)) == 0
+        push!(list_of_elements, (k,l))
+        continue
+      end
+      
+      # Determine remaining variables, after scaling "away" others.
+      remaining_vars_list = Set(1:length(gS))
+      for my_exps in ignored_sets
+        len_my_exps = length(my_exps)
+        inter_len = count(idx -> idx in [k,l], my_exps)
+        if (len_my_exps == 2 && inter_len == 1) || (len_my_exps == 3 && inter_len == 2)
+          delete!(remaining_vars_list, my_exps[findfirst(idx -> !(idx in [k,l]), my_exps)])
+        end
+      end
+      remaining_vars_list = collect(remaining_vars_list)
+  
+      # If one monomial of `new_p_hyper` has unset positions (a.k.a. new_p_hyper is not constant upon
+      # scaling the remaining variables), then keep this generator.
+      for exps in exponents(new_p_hyper)
+        if any(x -> x != 0, exps[remaining_vars_list])
+          push!(list_of_elements, (k,l))
+          break
+        end
+      end
+
+    end
+  end
+  
+  # Remember this result as attribute and return the findings.
+  set_attribute!(m, :_ambient_space_divisors_to_be_considered, list_of_elements)
+  return list_of_elements
+  
+end
+
+
+@doc raw"""
+    well_quantized_ambient_space_models_of_g4_fluxes(m::AbstractFTheoryModel; check::Bool = true)::Vector{CohomologyClass}
+
+Given an F-theory model $m$ defined as hypersurface in a simplicial and
+complete toric base, this method computes a basis of all well-quantized
+ambient space $G_4$-fluxes.
+
+Note that it can be computationally very demanding to check if a toric variety
+$X$ is complete (and simplicial). The optional argument `check` can be set
+to `false` to skip these tests.
+
+# Examples
+```jldoctest; setup = :(Oscar.LazyArtifacts.ensure_artifact_installed("QSMDB", Oscar.LazyArtifacts.find_artifacts_toml(Oscar.oscardir)))
+julia> B3 = projective_space(NormalToricVariety, 3)
+Normal toric variety
+
+julia> Kbar = anticanonical_divisor_class(B3)
+Divisor class on a normal toric variety
+
+julia> t = literature_model(arxiv_id = "1109.3454", equation = "3.1", base_space = B3, defining_classes = Dict("w"=>Kbar))
+Construction over concrete base may lead to singularity enhancement. Consider computing singular_loci. However, this may take time!
+
+Global Tate model over a concrete base -- SU(5)xU(1) restricted Tate model based on arXiv paper 1109.3454 Eq. (3.1)
+
+julia> g4_amb_list = well_quantized_ambient_space_models_of_g4_fluxes(t)
+2-element Vector{CohomologyClass}:
+ Cohomology class on a normal toric variety given by z^2
+ Cohomology class on a normal toric variety given by y^2
+
+julia> t_res = load("/datadisk/ToDo/1511.03209/FinalData/1511-03209-resolved.mrdi");
+
+# Oscar.__ambient_space_divisors_to_be_considered(t_res)
+# save("/home/i/list_of_divisor_pairs_to_be_considered.mrdi", list_of_divisor_pairs_to_be_considered)
+
+julia> list_of_divisor_pairs_to_be_considered = load("/datadisk/ToDo/1511.03209/FinalData/list_of_divisor_pairs_to_be_considered.mrdi");
+julia> set_attribute!(t_res, :_ambient_space_divisors_to_be_considered, list_of_divisor_pairs_to_be_considered);
+
+# ambient_space_flux_candidates_basis = ambient_space_models_of_g4_fluxes(t_res, check = false)
+# save("/home/i/ambient_space_flux_candidates_basis.mrdi", ambient_space_flux_candidates_basis)
+
+julia> ambient_space_flux_candidates_basis = load("/datadisk/ToDo/1511.03209/FinalData/ambient_space_flux_candidates_basis.mrdi");
+julia> set_attribute!(t_res, :ambient_space_models_of_g4_fluxes, ambient_space_flux_candidates_basis);
+
+julia> m = t_res;
+
+julia> Oscar._reduce_hypersurface_equation(ambient_space(m), hypersurface_equation(m), my_tuple)
+
+julia> Oscar._rationally_equivalent_cycle(ambient_space(m), my_tuple)
+
+julia> intersection_dict = Dict{NTuple{4, Int64}, ZZRingElem}()
+
+julia> special_intersection_dict = Dict{String, ZZRingElem}()
+
+julia> Oscar.sophisticated_intersection_product(ambient_space(m), my_tuple, hypersurface_equation(m), intersection_dict, special_intersection_dict)
+
+julia> res = well_quantized_ambient_space_models_of_g4_fluxes(t_res, check = false)
+
+julia> intersection_dict = Dict{NTuple{4, Int64}, ZZRingElem}()
+
+julia> special_intersection_dict = Dict{String, ZZRingElem}()
+
+julia> sophisticated_intersection_product(ambient_space(m), (173, 175, 101, 102), hypersurface_equation(m), intersection_dict, special_intersection_dict)
+
+julia> my_tuple = (173, 175, 101, 102)
+
+julia> my_tuple = (175, 175, 101, 102)
+
+julia> my_tuple = (175, 175, 175, 102)
+
+julia> my_tuple = (175, 175, 175, 175)
+
+julia> Oscar._rationally_equivalent_cycle(ambient_space(m), my_tuple, data)
+
+julia> return Oscar._reduce_hypersurface_equation(ambient_space(m), hypersurface_equation(m), my_tuple, data)
+
+julia> sophisticated_intersection_product(ambient_space(m), my_tuple, hypersurface_equation(m), inter_dict, s_inter_dict, data)
+```
+"""
+function well_quantized_ambient_space_models_of_g4_fluxes(m::AbstractFTheoryModel; check::Bool = true)
+#function well_quantized_ambient_space_models_of_g4_fluxes(m::AbstractFTheoryModel; check::Bool = true)::Vector{CohomologyClass}
+
+  # (1) Entry checks
+  @req base_space(m) isa NormalToricVariety "Computation of well-quantized G4-fluxes only supported for toric base and ambient spaces"
+  @req dim(ambient_space(m)) == 5 "Computation of well-quantized G4-fluxes only supported for 5-dimensional toric ambient spaces"
+  if check
+    @req is_complete(ambient_space(m)) "Computation of well-quantized G4-fluxes only supported for complete toric ambient spaces"
+    @req is_simplicial(ambient_space(m)) "Computation of well-quantized G4-fluxes only supported for simplicial toric ambient space"
+  end
+  #if has_attribute(m, :well_quantized_ambient_space_models_of_g4_fluxes)
+  #  return get_attribute(m, :well_quantized_ambient_space_models_of_g4_fluxes)
+  #end
+
+  # (2) Compute data, that is frequently used by the sophisticated intersection product below
+  S = cox_ring(ambient_space(m))
+  gS = gens(cox_ring(ambient_space(m)))
+  linear_relations = matrix(QQ, rays(ambient_space(m)))
+  scalings = [c.coeff for c in S.d]
+  mnf = Oscar._minimal_nonfaces(ambient_space(m))
+  sr_ideal_pos = [Vector{Int}(Polymake.row(mnf, i)) for i in 1:Polymake.nrows(mnf)]
+  data = (
+    S = S,
+    gS = gS,
+    linear_relations = linear_relations,
+    scalings = scalings,
+    sr_ideal_pos = sr_ideal_pos
+  )
+  inter_dict = Dict{NTuple{4, Int64}, ZZRingElem}()
+  s_inter_dict = Dict{String, ZZRingElem}()
+
+  my_tuple = (81, 101, 129, 130)
+  
+  return sophisticated_intersection_product(ambient_space(m), my_tuple, hypersurface_equation(m), inter_dict, s_inter_dict, data)
+
+  # (2) Obtain critical information - this may take some time!
+  ambient_space_flux_candidates_basis = ambient_space_models_of_g4_fluxes(m, check = check) # 629
+  ambient_space_flux_candidates_basis_indices = Vector{Tuple{Int64, Int64}}(undef, length(ambient_space_flux_candidates_basis))
+  if arxiv_doi(m) == "10.48550/arXiv.1511.03209"
+    for k in 1:length(ambient_space_flux_candidates_basis)
+      idxs = findall(x -> x != 0, collect(exponents(polynomial(ambient_space_flux_candidates_basis[k]).f))[1])
+      if length(idxs) == 1
+        idxs = [idxs[1], idxs[1]]
+      end
+      ambient_space_flux_candidates_basis_indices[k] = Tuple(idxs)
+    end
+  end
+  list_of_divisor_pairs_to_be_considered = Oscar._ambient_space_divisors_to_be_considered(m) # 1794
+
+  # (3) Work out the relevant intersection products
+  constraint_matrix = Vector{Vector{Int64}}()
+  for i in 1:length(ambient_space_flux_candidates_basis)
+    condition = Vector{ZZRingElem}()
+    for j in 1:length(list_of_divisor_pairs_to_be_considered)
+      println("$i, $j")
+      if arxiv_doi(m) == "10.48550/arXiv.1511.03209"
+        my_tuple = (ambient_space_flux_candidates_basis_indices[i]..., list_of_divisor_pairs_to_be_considered[j]...)
+        push!(condition, sophisticated_intersection_product(ambient_space(m), my_tuple, hypersurface_equation(m), inter_dict, s_inter_dict, data))
+      else
+        # CORRECT CODE TO BE ADDED!
+        push!(condition, 0)
+      end
+    end
+    push!(constraint_matrix, condition)
+  end
+
+  return [matrix(ZZ, hcat(constraint_matrix...)), inter_dict, s_inter_dict]
+
+end
diff --git a/experimental/FTheoryTools/src/exports.jl b/experimental/FTheoryTools/src/exports.jl
index 6a3f2b2acd4..53ce0b2e962 100644
--- a/experimental/FTheoryTools/src/exports.jl
+++ b/experimental/FTheoryTools/src/exports.jl
@@ -217,5 +217,6 @@ export weighted_resolution_generating_sections
 export weighted_resolution_zero_sections
 export weighted_resolutions
 export weights
+export well_quantized_ambient_space_models_of_g4_fluxes
 export zero_section
 export zero_section_class