import numpy as np

from tabulate import tabulate

# Input values for inflow and initial discharge

inflow_variable = 552
primary_discharge = 100

# Initial flow rates based on the continuity equation

outflow_rate = 50
adjusted_outflow = inflow_variable - 8 * outflow_rate
secondary_discharge = primary_discharge + inflow_variable
average_discharge1 = 0.5 * (primary_discharge + adjusted_outflow)
average_discharge2 = 0.5 * (primary_discharge + adjusted_outflow)
final_discharge1 = average_discharge2 + outflow_rate
final_discharge2 = outflow_rate + final_discharge1
discharge3 = 0.5 * (secondary_discharge - outflow_rate)
discharge5 = 0.5 * (secondary_discharge - outflow_rate)
discharge4 = discharge3 - outflow_rate
discharge6 = 0.5 * (discharge5 - outflow_rate)
discharge7 = discharge6 + discharge4 - outflow_rate
discharge8 = 0.5 * (discharge5 - outflow_rate)
discharge9 = discharge7 - final_discharge2 - outflow_rate

# Define loops and their discharges

circuit_1 = [1, 2, 5, 8, 10, primary_discharge, secondary_discharge, discharge5, discharge8,
average_discharge1]
circuit_2 = [5, 3, 4, 6, -discharge5, discharge3, discharge4, -discharge6]
circuit_3 = [8, 6, 7, 9, -discharge8, discharge6, discharge7, discharge9]
circuit_4 = [11, 10, 9, 13, 12, average_discharge2, -average_discharge1, -discharge9,
final_discharge2, final_discharge1]

# Resistance values for branches

branch_resistances = np.array([0.4, 0.4, 0.3, 0.3, 0.3, 0.2, 0.2, 0.2, 0.3, 0.3, 0.3, 0.2, 0.2])

# Connection matrix indicating loops and branches

connectivity_matrix = np.array([
[ 1, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0], # Circuit 1
[ 0, 0, 1, 1, -1, 0, 0, 0, 0, 0, 0, 1, 0], # Circuit 2
[ 0, 0, 0, 0, 0, 1, 1, -1, 0, 0, 0, 0, 0], # Circuit 3
[ 0, 0, 0, 0, 0, 0, 0, 0, -1, -1, 1, 1, 1], # Circuit 4
])

# Initialize flow direction for branches

direction_flow = np.ones(len(branch_resistances))
# Iteratively refine flow rates
for iter_count in range(1, 101):
def calculate_correction(circuit):
numerator, denominator = 0, 0
loop_half = int(0.5 * len(circuit))
for i in range(loop_half):
branch_idx = circuit[i] - 1
flow_rate = circuit[i + loop_half]
abs_flow = abs(flow_rate)
numerator += branch_resistances[branch_idx] * flow_rate * abs_flow
denominator += branch_resistances[branch_idx] * 2 * abs_flow
return -numerator / denominator

# Compute adjustments for each circuit

adjustments = [calculate_correction(circ) for circ in [circuit_1, circuit_2, circuit_3, circuit_4]]

# Update flow rates based on corrections

def update_flow(flow, adjustment, flow_dir, branch_idx):
updated_flow = flow + adjustment
if flow * updated_flow < 0:
flow_dir[branch_idx] *= -1
return updated_flow

# Apply corrections to all flows

primary_discharge = update_flow(primary_discharge, adjustments[0], direction_flow, 0)
secondary_discharge = update_flow(secondary_discharge, adjustments[0], direction_flow, 1)
discharge3 = update_flow(discharge3, adjustments[1], direction_flow, 2)
discharge4 = update_flow(discharge4, adjustments[1], direction_flow, 3)
discharge5 = update_flow(discharge5, adjustments[0] - adjustments[1], direction_flow, 4)
discharge6 = update_flow(discharge6, -adjustments[1] + adjustments[2], direction_flow, 5)
discharge7 = update_flow(discharge7, adjustments[2], direction_flow, 6)
discharge8 = update_flow(discharge8, adjustments[0] - adjustments[2], direction_flow, 7)
discharge9 = update_flow(discharge9, adjustments[2] - adjustments[3], direction_flow, 8)
average_discharge1 = update_flow(average_discharge1, adjustments[0] - adjustments[3],
direction_flow, 9)
average_discharge2 = update_flow(average_discharge2, adjustments[3], direction_flow, 10)
final_discharge1 = update_flow(final_discharge1, adjustments[3], direction_flow, 11)
final_discharge2 = update_flow(final_discharge2, adjustments[3], direction_flow, 12)

# Update circuits with new values

circuit_1 = [1, 2, 5, 8, 10, primary_discharge, secondary_discharge, discharge5, discharge8,
average_discharge1]
circuit_2 = [5, 3, 4, 6, -discharge5, discharge3, discharge4, -discharge6]
circuit_3 = [8, 6, 7, 9, -discharge8, discharge6, discharge7, discharge9]
circuit_4 = [11, 10, 9, 13, 12, average_discharge2, -average_discharge1, -discharge9,
final_discharge2, final_discharge1]

# Exit if corrections are below threshold

if all(abs(adj) < 0.001 for adj in adjustments):
break

# Verify nodal conditions

print("Nodal Condition Verification:")
conditions = [
primary_discharge + inflow_variable - secondary_discharge,
primary_discharge + adjusted_outflow - average_discharge1 - average_discharge2,
average_discharge2 + outflow_rate - final_discharge1,
outflow_rate + discharge3 + discharge5 - secondary_discharge,
outflow_rate + discharge6 + discharge8 - discharge5,
discharge8 + discharge9 - outflow_rate - average_discharge1,
outflow_rate + discharge4 - discharge3,
discharge4 + discharge6 - outflow_rate - discharge7,
discharge9 + outflow_rate - discharge7 + final_discharge2,
outflow_rate + final_discharge1 - final_discharge2,
]
for condition in conditions:
print(f"{condition:.6f}")

# Output final flow rates in tabular format

flow_labels = [
"Primary Discharge (Q1)", "Secondary Discharge (Q2)", "Discharge 3 (Q3)",
"Discharge 4 (Q4)", "Discharge 5 (Q5)", "Discharge 6 (Q6)",
"Discharge 7 (Q7)", "Discharge 8 (Q8)", "Discharge 9 (Q9)",
"Discharge 10 (Q10)", "Discharge 11 (Q11)",
"Discharge 12 (Q12)", "Discharge 13 (Q13)"
]
final_flows = [
primary_discharge, secondary_discharge, discharge3, discharge4,
discharge5, discharge6, discharge7, discharge8, discharge9,
average_discharge1, average_discharge2, final_discharge1, final_discharge2
]

table = tabulate(
[[label, f"{flow:.2f}"] for label, flow in zip(flow_labels, final_flows)],
headers=["Flow Component", "Final Value"],
tablefmt="grid"
)
print("\nFinal Flow Rates in Table Format:")
print(f"Iterations performed: {iter_count}\n")
print(table)