I am trying to solve seven coupled differential equations using odeint

Question

I tried to slove this using the following code.I am just a beginner in programming trying to learn. So, i took the help of chatgpt. The code is syntactically correct but there seems to be a logical error The problem i am facing is with the term 'phi'. phi can take values 0 and 1 depending on the time of the day. So, it should change every 12 hours. EX: first 12 hrs should be 1 next 12 hrs should be 0, similarly, it keep oscillating like this. i dont know how define a condition to make this happen.

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# Define the system of coupled differential equations
def model(y, t, params):
    CLm, CLc, CLn, CTm, CTc, CTn, CPn = y
    n1, a, g1, m1, p1, r1, r2, m2, k1, k2, m3, k3, n2, g2, b, m4, k4, p2, r3, r4, m5, k5, m6, k6, p3, m7, k7, q1, q2 = params
    
    # Determine the value of phi based on time
    phi = 1 if (t // 12) % 2 == 0 else 0
    
    
    # Differential equations for other variables
    dCLmdt = (q1 * CPn * phi) + (n1 * CTn**a) / (g1**a + CTn**a) - ((m1 * CLm)/(k1+CLm))
    dCLcdt = (p1 * CLm) - (r1 * CLc) + (r2 * CLn) - ((m2 * CLc) / (k2 + CLc))
    dCLndt = (r1 * CLc) -(r2 * CLn) - ((m3 * CLn) / (k3 + CLn))
    dCTmdt = ((n2 * g2**b) / (g2**b + CLn**b)) - ((m4 * CTm) / (k4 + CTm))
    dCTcdt = (p2 * CTm) - (r3 * CTc) + (r4 * CTn) -((m5 * CTc) / (k5 + CTc))
    dCTndt = (r3 * CTc) - (r4 * CTn) - ((m6 * CTc) / (k6 + CTn))
    dCPndt = ((1 - phi) * p3) - ((m7 * CPn) / (k7 + CPn)) - (q2 * phi * CPn)
    
    return [dCLmdt, dCLcdt, dCLndt, dCTmdt, dCTcdt, dCTndt, dCPndt]

# Initial conditions
y0 = [1, 1, 1, 1, 1, 1, 1]

# Time points
t = np.linspace(0, 24*10, 100)  # 15 days simulation with 1000 evenly spaced data points

# Parameters
params = (n1, a, g1, m1, p1, r1, r2, m2, k1, k2, m3, k3, n2, g2, b, m4, k4, p2, r3, r4, m5, k5, m6, k6, p3, m7, k7, q1, q2) = (
    16.9711, 1, 3.0351, 9.3383, 4.9753, 1.4563, 0.8421, 16.9058, 1.3294, 0.8085, 1.0214, 0.1445, 1.3043, 0.386, 2, 1.6859, 0.2089, 1.2947, 0.0451, 0.0018, 0.4242, 0.3187, 0.0484, 0.3505, 0.5, 1.2, 1.2, 2.9741, 1.0
)

# Solve the differential equations
sol = odeint(model, y0, t, args=(params,), rtol=1e-6, atol=1e-6)

# Extract results
CLm, CLc, CLn, CTm, CTc, CTn, CPn = sol[:, 0], sol[:, 1], sol[:, 2], sol[:, 3], sol[:, 4], sol[:, 5], sol[:, 6]

# Plotting
plt.figure(figsize=(10, 6))
plt.plot(t, CLm, label='CLm')
plt.plot(t, CLc, label='CLc')
plt.plot(t, CLn, label='CLn')
plt.plot(t, CTm, label='CTm')
plt.plot(t, CTc, label='CTc')
plt.plot(t, CTn, label='CTn')
plt.plot(t, CPn, label='CPn')
plt.xlabel('Time')
plt.ylabel('Variable')
plt.legend()
plt.show()
print(t)

I tried a logic where i defined time by taking the remainder of the quotient by dividing it by 2. if it is 0 then the value of phi should be 1 else 0. but this condition is not working properly for all the timepoints. when i printed phi for 240 hrs with 100 timepoints it is printing 12000 phi values. Can someone help me with this?