How to color a cube in perspective projection and make lines around it

I made a turtle python program for projecting a cube in perspective but when want to color it there will be some issues because there is no transparency in turtle, so i need to find a way to avoid filling square on top of each other

In the first try i made the turtle know the furthest point and don't fill any of squares that have this point, there will be two ore four points that have the same z value but that applies on them too, this works well in orthographic projection but don't work on perspective, there will be a squares that is behind another square and get filled as shown in the picture an image that shows the difference

For solving that i made the turtle to fill the largest square over the rest by making it the last of the list, but that will cause some problems if i applied different Angel on the rotation x and rotation y in the same time , it is hard to see but it is there.

For making the lines i tried the same as squares and that worked for orthographic but didn't work for the perspective for the same problem.

That is the code

#imports
import turtle
import time
from math import *

#turtle and some variables
t = turtle.Turtle()
turtle.tracer(0, 0)
turtle.bgcolor("lightblue")
t.width(20)
t.hideturtle()
angle1 = 0
angle2 = 0
scale = 200

#matrices multiplication
def mult(matrix1, matrix2):
    result = []
    if len(matrix1[0]) == len(matrix2):
        for _ in range(len(matrix1)):
            result += [[0]*len(matrix2[0])]
        for i in range(len(matrix1)):
            for j in range(len(matrix2[0])):
                for k in range(len(matrix2)):
                    result[i][j] += matrix1[i][k] * matrix2[k][j]
        return result           

#forms of points
def point_matrix(point_as_3Dpvector):
    p = point_as_3Dpvector
    return [[p[0]], [p[1]], [p[2]]]
    
def point_3Dpvector(point_as_matrix):
    p = point_as_matrix
    return (p[0][0],p[1][0], p[2][0])
    
def point_2Dpvector(point_as_3Dpvector):
    p = point_as_3Dpvector
    return (p[0], p[1])
        
#make a line
def line(ln):
    t.up()
    t.goto(point_2Dpvector(ln[0]))
    t.down()
    t.goto(point_2Dpvector(ln[1]))
    t.up()
    
#make multiple lines from a list
def line_a_list(list):
    for i in list:
        line(i)
    
#fill a square
def fill(s):
    t.color(s[4])
    t.goto(point_2Dpvector(s[0]))
    t.begin_fill()
    t.goto(point_2Dpvector(s[1]))
    t.goto(point_2Dpvector(s[2]))
    t.goto(point_2Dpvector(s[3]))
    t.end_fill()
    t.up()
    
#fill muliple squares from a list
def fill_a_list(list):
    for i in list:
        fill(i)

#knows if a point is in a list or not
def isthere(point, list):
    r = 0
    for i in list:
        if point == i:
            r += 1
    if r > 0:
        return True
    else:
        return False
        
#calculates the area of a square
def area(se):
    a1 = se[0][0]
    a2 = se[1][0]
    a3 = se[2][0]
    b1 = se[0][1]
    b2 = se[1][1]
    b3 = se[2][1]
    a = a2 * b3 - b2 * a3
    b = -a1 * b3 + b1 * a3
    c = a1 * b2 - b1 * a2
    return abs(a + b + c)

#tell the turtle what squares to fill and the ones not
#by knowing the farest point by the z value 
#and do not fill any square that have this point 
def do_squares(list_squares, list_points):
    
    #variables
    furthest_point = list_points[0]
    do_not_points = []
    do_not_squares = []
    large_square = ()
    do_squares = []
    
    #find the farest point by its z value
    for i in list_points:
        if i[2] < furthest_point[2]:
            furthest_point = i
            
    #find the other points that have the same z value
    for i in list_points:
        if i[2] == furthest_point[2]:
            do_not_points += [i]    
                    
    #find the squares that have this points
    for i in list_squares:
        for j in do_not_points:
            if isthere(j, i) == True:
                do_not_squares += [i]
                
    #remove the squares that have the points from
    #the squares list and leave the rest there
    for i in list_squares:
        if isthere(i, do_not_squares) == False:
            do_squares += [i]
    
    #find the largest square and make it the last
    large = do_squares[0]
    for i in do_squares:
        if area(i) > area(large):
            large = i
            
    large_square = large
    do_squares.remove(large)
    do_squares += [large_square]
    
    return do_squares

#as the def before but for lines
def do_lines(list_lines, list_points):
    
    furthest_point = list_points[0]
    do_not_points = []
    do_not_lines = []
    do_lines = []
    
    for i in list_points:
        if i[2] < furthest_point[2]:
            furthest_point = i
            
    for i in list_points:
        if i[2] == furthest_point[2]:
            do_not_points += [i]    
                
    for i in list_lines:
        for j in do_not_points:
            if isthere(j, i) == True:
                do_not_lines += [i]
                
    for i in list_lines:
        if isthere(i, do_not_lines) == False:
            do_lines += [i]
            
    return do_lines
            
#the cube's points and points' list
p0 = (-0.5, -0.5, -0.5)
p1 = (0.5, -0.5, -0.5)
p2 = (0.5, 0.5, -0.5)
p3 = (-0.5, 0.5, -0.5)
p4 = (-0.5, -0.5, 0.5)
p5 = (0.5, -0.5, 0.5)
p6 = (0.5, 0.5, 0.5)
p7 = (-0.5, 0.5, 0.5)

points_list = [p0, p1, p2, p3, p4, p5, p6, p7]

#the loop of operations that do the work
while True:
    
    #variables
    t.clear()
    l = []   # list of rotated squares
    angle1 += .012 / 3  #radian
    angle2 += .005 / 3
    
    #rotation matrices
    rotateZ = [
    [cos(angle1), -sin(angle1), 0],
    [sin(angle1), cos(angle1), 0],
    [0, 0, 1]
    ]                      #Z axis
    
    rotateX = [
    [1, 0, 0],
    [0, cos(angle1), -sin(angle1)],
    [0, sin(angle1), cos(angle1)]
    ]                      #X axis
    
    rotateY = [
    [cos(angle2), 0, -sin(angle2)],
    [0, 1, 0],
    [sin(angle2), 0, cos(angle2)]
    ]                      #Y axis
    
    
    for point in points_list:
        
        #rotate points
        rotated = mult(rotateY, point_matrix(point))
        rotated = mult(rotateX, rotated)
        rotated = mult(rotateZ, rotated)
        
        #project points
        
        z = 1 / (2 - rotated[2][0]) * 2   #prespective
        #z = 1   #orthographic
        
        projection =[
        [z*scale, 0, 0],
        [0, z*scale, 0],
        [0, 0, 1]
        ]
    
        projected_points = mult(projection, rotated)
        l += [point_3Dpvector(projected_points)]
        
    # list of squares
    squar = [
    (l[0], l[1], l[2], l[3], "pink"),
    (l[4], l[5], l[6], l[7], "#5676AA"),
    (l[1], l[5], l[6], l[2], "#FF7777"),
    (l[0], l[4], l[7], l[3], "cyan"),
    (l[0], l[1], l[5], l[4], "#FFAA55"),
    (l[3], l[2], l[6], l[7], "#DA90F5")
    ]
    fill_a_list(do_squares(squar,l))
    
    #list of lines
    lines = [
    (l[0], l[1]),
    (l[0], l[4]),
    (l[1], l[5]),
    (l[4], l[5]),
    (l[1], l[2]),
    (l[5], l[6]),
    (l[6], l[2]),
    (l[0], l[3]),
    (l[2], l[3]),
    (l[3], l[7]),
    (l[7], l[4]),
    (l[7], l[6])
    ]
    t.color("black")
    line_a_list(do_lines(lines, l)) #you can turn it of to see the coloring problem
    turtle.update()
        
turtle.done()

Visibility for 3D rendering is not an easy problem. But in your case, with a cube, we can use a method based on the convexity of the mesh: this property guarantees that every face will be either fully visible or not visible at all.

So, to check the visibility status of each face, you can consider any point M on this face (not on an edge, the center of the face would be a good candidate) and calculate the 3D line passing through this point M and the center of the camera C. Then you can compute the intersection P of this line with every other face and check if this intersection exists or not (ie faces overlap or not) and if the intersection P is between M and C (ie the face is hidden) or behind M (ie the face is not hidden by this other face)

Here is a 2D example: if we consider the top face, containing M1, we can see that the corresponding line intersects the left face, but not the other two faces. This intersection P1 is between C and M1, so the top face is hidden. Considering the left face, with point M2, we find an intersection P2 of the corresponding line with the right face, and no intersection with the other faces. This intersection is not between C and M2. As no other intersection exists, it means that the left face is visible.

About the visibility of the edges, once you have a working algorithm for the faces, you can simply check for each edge if it belongs to a visible face or not.