How to structure my sage worksheet to make it understandable?

edit

I work towards a worksheet that will allow me to develope a physical model of a real world problem. I expect it to grow and collect additional special cases and factor in some more physical phenomena.

I would like this worksheet to remain readable and make the flow of thoughts transparent and comprehensible. Also it should be easy to see that it is actually correct.

What would be a "unit test" that the process so far has been correct? Is there support for switching such tests on and off?

Below, you can see just the start of the cases and it looks confusing already. How would YOU strucure this? What guidelines do you have to achive something like this?

x,z, d_e, e_m, e_s, v_ges, Q, U, a =var('x z d_e e_m e_s v_ges Q U a')
z_c = var('z_c', latex_name=r'\varDelta z_{dc}')
z_m = var('z_m', latex_name=r'\varDelta z_{ds}')

#========= dielectric constants ==================================
e_0 = var('e_0', latex_name=r'\varepsilon_0')
# sediment
e_s = var('e_s', latex_name=r'\varepsilon_s')
# seawater
e_m = var('e_m', latex_name=r'\varepsilon_m')

#e_Ersatz = var('e_Ersatz', latex_name=r'\varepsilon_{Ersatz}')
e_Ersatz(e_1,e_2,d_1,d_2) = e_1*e_2*(d_1+d_2)/(d_1*e_2+d_2*e_1)
print( SR.symbol('e_Ersatz') == e_Ersatz)
#latex(SR.symbol('e_Ersatz') == e_Ersatz)
#============ electric potential field ===============================
#p = var ('p', latex_name=r'\varphi')
p_ges = var('p_ges', latex_name=r'\varphi_{ges}')
e = var('e', latex_name=r'\varepsilon')

p(e, Q, r) = 1/(4*pi*e)*Q/r
print( SR.symbol ('p') == p)

p_1 = p(  heaviside( z_s -z )                           * e_0*e_m
        + heaviside( z- z_s ) * heaviside(z_c + z_s -z) * e_0*e_Ersatz(e_m, e_s, z_s            , z-z_s)
        + heaviside( z -(z_c + z_s))                    * e_0*e_Ersatz(e_m, e_s, (z-2*(z_c-z_s)), 2*(z_c-z_s)),
          1, sqrt(x^2 + z^2)
)
p_2 = p(  heaviside( z_s -z )                           * e_0*e_m
        + heaviside( z- z_s ) * heaviside(z_c + z_s -z) * e_0*e_Ersatz(e_m, e_s, z_s            , z-z_s)
        + heaviside( z -(z_c + z_s))                    * e_0*e_Ersatz(e_m, e_s, ( z -2*(z_c-z_s)), 2*(z_c-z_s)),
        -1, sqrt((d_e-x)^2 + z^2)
)
p_3 = p(    heaviside( z_s -z)                            * e_0*e_Ersatz(e_m, e_s, (-z +2*(z_c-z_s)), 2*(z_c-z_s))
          + heaviside( z- z_s ) * heaviside(z_c + z_s -z) * e_0*e_Ersatz(e_m, e_s, z_s              , 2*z_c-z_s -z)
          + heaviside( z -(z_c + z_s))                    * e_0*e_m,
          -1, sqrt((2*z_c-z)^2+x^2))
p_4 = p(    heaviside( z_s -z)                            * e_0*e_Ersatz(e_m, e_s, (-z +2*(z_c-z_s)), 2*(z_c-z_s))
          + heaviside( z- z_s ) * heaviside(z_c + z_s -z) * e_0*e_Ersatz(e_m, e_s, z_s              , 2*z_c-z_s -z)
          + heaviside( z -(z_c + z_s))                    * e_0*e_m,    
            1, sqrt((2*z_c-z)^2+(d_e -x)^2))
p_fix_rhs = p_1 + p_2 + p_3 + p_4

d = var ('d')
d = U/p_fix_rhs

print( "Ladung, so dass am Punkt a die Spannung U ist")
show( Q==d )