Finding intersections between an implicit fill pattern and an outer layer for formation of a 3D printing slice

I am sure there are lots of solutions on the internet but it is pretty hard for me to understand someone else's code. Therefore, I want to ask a question about my own case. Basically, I want to put an infill pattern (which is defined by a implicit function such as cos(x)+cos(y)+cos(z)=0) into a slice. How can I find the intersections and plot the slice accordingly ?

(I have also heard and searched about marching-cubes algorithm but I'm not sure how to implement it in this case.)

Here is my code:

import numpy as np
from stl import mesh
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

class Triangle:
    
    def Triangle_builder(triangle_number):
        px=stl.x[triangle_number][:]
        px=np.append(px, px[0])
        py=stl.y[triangle_number][:]
        py=np.append(py, py[0])
        pz=stl.z[triangle_number][:]
        pz=np.append(pz, pz[0])
    
        # ax.plot3D(px,py,pz)
        
        p0=stl.points[triangle_number][0:3]
        p1=stl.points[triangle_number][3:6]
        p2=stl.points[triangle_number][6:9]
        p3=stl.points[triangle_number][0:3]
        
        return p0,p1,p2,p3
          
def isect(p0,p1,plane_p,plane_n):
    vector_x=p1[0]-p0[0]
    vector_y=p1[1]-p0[1]
    vector_z=p1[2]-p0[2]
    vector=np.array([vector_x,vector_y,vector_z])
    
    denom=np.dot(vector,plane_n)
    if denom==0:
        return None #parallel
    else:
        isect_p=p0+vector*(np.dot((plane_p-p0),plane_n)/denom)
        if isect_p[2]>max(p0[2],p1[2]) or isect_p[2]<min(p0[2],p1[2]):
            return None
        else:
            isect_px=isect_p[0]
            isect_py=isect_p[1]
            isect_pz=isect_p[2]
            
    return isect_px,isect_py,isect_pz
       
def plotter(hdy):
    for i in range(len(hdy)-1):
        # print(hdy[i])
        # print(hdy[i+1])
        ax.plot3D([hdy[i][0],hdy[i+1][0]],[hdy[i][1],hdy[i+1][1]],[hdy[i][2],hdy[i+1][2]])
    return ax

def filter(list_):
    for i in range(len(list_)):
        if list_[i]==None:
            list_.pop(i)
    return list_
    
def hyp_part1(x,y,z):
    return np.cos(x)+np.cos(y)+np.cos(z)
    
stl = mesh.Mesh.from_file("sphere.stl")
row,col=stl.z.shape
z=np.arange(stl.min_[2],stl.max_[2],1)
for k in z:
    plane_p=np.array([0,0,k])
    plane_n=np.array([0,0,1])
    isect_list=[]
    i=0
    while i<row:
        j=0
        p_list=Triangle.Triangle_builder(i)
        i=i+1
        isect_list=[]
        while j<3:
            p0=p_list[j]
            p1=p_list[j+1]
            isect_list.append(isect(p0, p1, plane_p, plane_n))
            j=j+1
            filter(isect_list)
            # print(isect_list)
            if len(isect_list)==2:
                ax=plotter(isect_list)

x_range = np.arange(0,4,5) 
y_range = np.arange(0,4,5)
X,Y = np.meshgrid(x_range,y_range)
A = np.linspace(-10, 10, 8)

A1,A2 = np.meshgrid(A,A)    

for z in A: 
    X,Y = A1, A2
    Z = hyp_part1(X,Y,z)
    ax.contour(X, Y, Z+z, [z], zdir='z')

for y in A: 
    X,Z= A1, A2
    Y = hyp_part1(X,y,Z)
    ax.contour(X, Y+y, Z, [y], zdir='y')

for x in A:
    Y,Z = A1, A2 
    X = hyp_part1(x,Y,Z)
    ax.contour(X+x, Y, Z, [x], zdir='x')

Slices and the pattern with excess lines Slices alone