I keep getting a divide by zero error in my code and not sure how to fix it.
# coding: utf-8
# In[74]:
import numpy
import math
Theta_ep = [253.15, 283.15, 313.15]
Theta_e = [[0.0, 0.0, 0.0],[0.0, 0.0, 0.0],[0.0, 0.0, 0.0],[0.0, 0.0, 0.0]]
p = [1000, 1050, 850, 500, 100]
Lambda = 3.504
K1 = [0, 0, 0, 0]
K2 = [0, 0, 0, 0]
D = [0, 0, 0, 0,]
R_v = 461 #specific gas constant
Rho = 0.597728
C = 273.15
T_E = [[0, 0, 0],[0, 0, 0],[0, 0, 0],[0, 0, 0]]
T_w = [[0, 0, 0],[0, 0, 0],[0, 0, 0],[0, 0, 0]]
A = 2675.0
e_s = [[0, 0, 0],[0, 0, 0],[0, 0, 0],[0, 0, 0]]
for i in range (0, 4):
for j in range(0, 3):
Pi = (p[i+1] / p[0]) ** (1.0 / Lambda)
K1[i] = -38.5 * (Pi ** (2.0)) + 137.81 * Pi - 53.737
K2[i] = -4.392 * (Pi ** (2.0)) + 56.831 * Pi - 0.384
D[i] = (0.1859 * (p[i+1] / p[0]) + 0.6512) ** (-1.0)
T_E[i][j] = Theta_ep[j] * Pi
e_s[i][j] = 6.112 * math.exp((17.67*T_E[i][j]) / (T_E[i][j] + 243.5))
r = 0.622 * (e_s[i][j] / (p[i+1] * 1000.0 - e_s[i][j]))
Derivative = (4302.65) / ((T_E[i][j] + 243.5) ** 2.0)
if (C / T_E[i][j]) ** (Lambda) > D[i]:
T_w[i][j] = -(T_E[i][j] - C - ((A * r) / (1.0 + A * r * Derivative)))
elif 1.0 <= (C / T_E[i][j]) ** (Lambda) <= D[i]:
T_w[i][j] = -(K1[i] - K2[i] * ((C / T_E[i][j]) ** (Lambda)))
elif 0.4 <= (C / T_E[i][j]) ** (Lambda) < 1.0:
T_w[i][j] = ((K1[i] - 1.21) - (K2[i] - 121) * (C / T_E[i][j]) ** (Lambda))
elif (C / T_E[i][j]) ** (Lambda) < 0.4:
T_w[i][j] = ((K1[i] - 2.66) - (K2[i] - 1.21) * (C / T_E[i][j]) ** (Lambda))
Theta_e[i][j] = T_E[i][j]*((1000.0 / p[i+1]) ** (0.2854 * (1.0-0.28 * 10.0 ** (-3.0 * r)))) * math.exp(((3.376 / T_w[i][j] - 0.00254) * (r * (1.0 + 0.81 * 10.0 ** (-3.0 * r)))))
print ('T_w values:')
print (T_w)
print ('Theta_e values:')
print (Theta_e)
I get a float division by zero error at the first if statement and was wondering if there is a way to fix it?