Circuit Solver

I finally got a basic circuit solver going! What has been done so far:

1. Only voltage source, resistor and inductor have been included. Capacitor has not yet been tested.
2. Ammeters can be introduced in the circuit wherever needed.

So here is the code (click "view raw" below the code box to see the code in a new window):

	#! /usr/bin/env python
	import sys, math, matrix
	from network_reader import *
	from solver import *
	def csv_tuple(csv_elem):
	""" Convert a cell position from spreadsheet form
	to [row, tuple] form. """
	csv_elem.upper()
	# Create a dictionary of alphabets
	csv_col="A B C D E F G H I J K L M N O P Q R S T U V W X Y Z"
	csv_dict={}
	csv_col_list=csv_col.split(" ")
	# Make the alphabets correspond to integers
	for c1 in range(1, 27):
	csv_dict[csv_col_list[c1-1]]=c1
	# The cell position starts with a number
	flag="number"
	c1=0
	while flag=="number":
	# When conversion to int fails
	# it means the element is an alphabet
	try:
	int(csv_elem[c1])
	except ValueError:
	flag="alphabet"
	else:
	c1+=1
	# Split them up into numbers and alphabets
	pol_row=int(csv_elem[0:c1])
	pol_col=csv_elem[c1:]
	elem_tuple=[pol_row-1, 0]
	# Convert the alphabets to number
	# Similar to converting binary to decimal
	for c1 in range(len(pol_col)-1, -1, -1):
	if len(pol_col)-1-c1>0:
	elem_tuple[1]+=26*(len(pol_col)-1-c1)*csv_dict[pol_col[c1]]
	else:
	elem_tuple[1]+=csv_dict[pol_col[c1]]-1
	return elem_tuple
	def reading_params(param_file):
	""" Read a file. Ramove additional quotes and
	carriage returns. Remove leading spaces. """
	from_file=[]
	for line in param_file:
	from_file.append(line.split(","))
	for c1 in range(len(from_file)):
	for c2 in range(len(from_file[c1])-1, -1, -1):
	# Remove additional quotes and carriage returns
	if from_file[c1][c2]:
	scrub_elements(from_file, c1, c2)
	# Remove blank spaces and null elements
	if from_file[c1][c2]==" " or from_file[c1][c2]=="":
	del from_file[c1][c2]
	return from_file
	class Resistor:
	def __init__(self, res_index, res_pos, res_tag):
	self.res_number=res_index
	self.res_pos=res_pos
	self.res_tag=res_tag
	self.resistor=100.0
	self.current=0.0
	self.voltage=0.0
	def display(self):
	print "Resistor is ",
	print self.res_tag,
	print "= %f" %self.resistor,
	print " located at ",
	print self.res_pos
	def ask_values(self, x_list, ckt_mat):
	res_params=["Resistor"]
	res_params.append(self.res_tag)
	res_params.append(self.res_pos)
	res_params.append(self.resistor)
	x_list.append(res_params)
	def get_values(self, x_list, ckt_mat):
	self.resistor=float(x_list[0])
	def transfer_to_sys(self, sys_loops, mat_e, mat_a, mat_b, mat_u, source_list):
	""" The matrix A in E.dx/dt=Ax+Bu will be updated by the
	resistor value."""
	for c1 in range(len(sys_loops)):
	for c2 in range(c1, len(sys_loops)):
	# Updating the diagonal element corresponding
	# to a loop.
	if c1==c2:
	if csv_tuple(self.res_pos) in sys_loops[c1][c2]:
	mat_a.data[c1][c2]+=self.resistor
	else:
	# Updating the offdiagonal element depending
	# on the sense of the loops (aiding or opposing)
	for c3 in range(len(sys_loops[c1][c2])):
	if csv_tuple(self.res_pos) in sys_loops[c1][c2][c3]:
	if sys_loops[c1][c2][c3][-1]=="forward":
	mat_a.data[c1][c2]+=self.resistor
	else:
	mat_a.data[c1][c2]-=self.resistor
	# Because the matrices are symmetric
	mat_a.data[c2][c1]=mat_a.data[c1][c2]
	return
	def generate_val(self, source_lst, sys_loops, mat_u, t, dt):
	pass
	def update_val(self, sys_loops, mat_e, mat_a, mat_b, state_vec, mat_u):
	pass
	class Inductor:
	def __init__(self, ind_index, ind_pos, ind_tag):
	self.ind_number=ind_index
	self.ind_pos=ind_pos
	self.ind_tag=ind_tag
	self.inductor=0.001
	self.current=0.0
	self.i_dbydt=0.0
	self.voltage=0.0
	def display(self):
	print "Inductor is ",
	print self.ind_tag,
	print "=%f" %self.inductor,
	print " located at ",
	print self.ind_pos
	def ask_values(self, x_list, ckt_mat):
	ind_params=["Inductor"]
	ind_params.append(self.ind_tag)
	ind_params.append(self.ind_pos)
	ind_params.append(self.inductor)
	x_list.append(ind_params)
	def get_values(self, x_list, ckt_mat):
	self.inductor=float(x_list[0])
	def transfer_to_sys(self, sys_loops, mat_e, mat_a, mat_b, mat_u, source_list):
	""" The matrix E in E.dx/dt=Ax+Bu will be updated by the
	inductor value."""
	for c1 in range(len(sys_loops)):
	for c2 in range(c1, len(sys_loops)):
	# Updating the diagonal element corresponding
	# to a loop.
	if c1==c2:
	if csv_tuple(self.ind_pos) in sys_loops[c1][c2]:
	mat_e.data[c1][c2]+=self.inductor
	else:
	# Updating the offdiagonal element depending
	# on the sense of the loops (aiding or opposing)
	for c3 in range(len(sys_loops[c1][c2])):
	if csv_tuple(self.ind_pos) in sys_loops[c1][c2][c3]:
	if sys_loops[c1][c2][c3][-1]=="forward":
	mat_e.data[c1][c2]+=self.inductor
	else:
	mat_e.data[c1][c2]-=self.inductor
	# Because the matrices are symmetric
	mat_e.data[c2][c1]=mat_e.data[c1][c2]
	return
	def generate_val(self, source_lst, sys_loops, mat_u, t, dt):
	pass
	def update_val(self, sys_loops, mat_e, mat_a, mat_b, state_vec, mat_u):
	pass
	class Capacitor:
	def __init__(self, cap_index, cap_pos, cap_tag):
	self.cap_number=cap_index
	self.cap_pos=cap_pos
	self.cap_tag=cap_tag
	self.capacitor=10.0e-6
	self.current=0.0
	self.voltage=0.0
	self.v_dbydt=0.0
	def display(self):
	print "Capacitor is ",
	print self.cap_tag,
	print "= %f" %self.capacitor,
	print " located at ",
	print self.cap_pos,
	print " with positive polarity towards %s" %(csv_element(self.cap_polrty))
	def ask_values(self, x_list, ckt_mat):
	cap_params=["Capacitor"]
	cap_params.append(self.cap_tag)
	cap_params.append(self.cap_pos)
	cap_params.append(self.capacitor)
	# Looking for a default value of polarity
	# in the neighbouring cells
	self.cap_elem=csv_tuple(self.cap_pos)
	if self.cap_elem[0]>0:
	if ckt_mat[self.cap_elem[0]-1][self.cap_elem[1]]:
	self.cap_polrty=[self.cap_elem[0]-1, self.cap_elem[1]]
	if self.cap_elem[1]>0:
	if ckt_mat[self.cap_elem[0]][self.cap_elem[1]-1]:
	self.cap_polrty=[self.cap_elem[0], self.cap_elem[1]-1]
	if self.cap_elem[0]<len(ckt_mat)-1:
	if ckt_mat[self.cap_elem[0]+1][self.cap_elem[1]]:
	self.cap_polrty=[self.cap_elem[0]+1, self.cap_elem[1]]
	if self.cap_elem[1]<len(ckt_mat)-1:
	if ckt_mat[self.cap_elem[0]][self.cap_elem[1]+1]:
	self.cap_polrty=[self.cap_elem[0], self.cap_elem[1]+1]
	cap_params.append("Positive polarity towards (cell) = %s" %csv_element(self.cap_polrty))
	x_list.append(cap_params)
	def get_values(self, x_list, ckt_mat):
	self.capacitor=float(x_list[0])
	cap_polrty=x_list[1].split("=")[1]
	# Convert the human readable form of cell
	# to [row, column] form
	while cap_polrty[0]==" ":
	cap_polrty=cap_polrty[1:]
	self.cap_polrty=csv_tuple(cap_polrty)
	if not ckt_mat[self.cap_polrty[0]][self.cap_polrty[1]]:
	print "Polarity incorrect. Branch does not exist at %s" %csv_element(self.cap_polrty)
	def transfer_to_sys(self, sys_loops, mat_e, mat_a, mat_b, mat_u, source_list):
	""" The matrix B in E.dx/dt=Ax+Bu will be updated by the
	polarity of the capacitor."""
	for c1 in range(len(sys_loops)):
	if csv_tuple(self.cap_pos) in sys_loops[c1][c1]:
	# If the positive polarity appears before the capacitor position
	# it means as per KVL, we are moving from +ve to -ve
	# and so the capacitor voltage will be taken negative
	if sys_loops[c1][c1].index(self.cap_polrty)<sys_loops[c1][c1].index(csv_tuple(self.cap_pos)):
	mat_b.data[c1][source_list.index(self.cap_pos)]=-1.0
	else:
	mat_b.data[c1][source_list.index(self.cap_pos)]=1.0
	return
	def generate_val(self, source_lst, sys_loops, mat_u, t, dt):
	""" The capacitor voltage is updated in the matrix u in
	E.dx/dt=Ax+Bu ."""
	self.v_dbydt=self.current/self.capacitor
	self.voltage+=self.v_dbydt*dt
	mat_u.data[source_lst.index(self.cap_pos)][0]=self.voltage
	def update_val(self, sys_loops, mat_e, mat_a, mat_b, state_vec, mat_u):
	""" The capacitor current is calculated as a result of the KVL."""
	c1=0
	while c1<len(sys_loops):
	# The first occurance of the capacitor
	if csv_tuple(self.cap_pos) in sys_loops[c1][c1]:
	# Check for polarity. If positive polarity is before the capacitor
	# it means the current is entering the positive terminal and so
	# current is positive.
	if sys_loops[c1][c1].index(self.cap_polrty)<sys_loops[c1][c1].index(csv_tuple(self.cap_pos)):
	self.current=state_vec.data[c1][0]
	else:
	self.current=-state_vec.data[c1][0]
	# Check for capacitor in other loops with the capacitor
	# as common branch
	for c2 in range(c1+1, len(sys_loops)):
	# These lists are nested lists as two loops
	# may have separate branches in common
	for c3 in range(len(sys_loops[c1][c2])):
	# If capacitor is found in one of the branches.
	# Can only be in one branch at most.
	if csv_tuple(self.cap_pos) in sys_loops[c1][c2][c3]:
	# Check for polarity again
	if sys_loops[c1][c2][c3].index(self.cap_polrty)<sys_loops[c1][c2][c3].index(csv_tuple(self.cap_pos)):
	# Then check is the loop is aiding or opposing
	# the main loop.
	if sys_loops[c1][c2][c3][-1]=="forward":
	self.current+=state_vec.data[c2][0]
	else:
	self.current-=state_vec.data[c2][0]
	else:
	if sys_loops[c1][c2][c3][-1]=="forward":
	self.current-=state_vec.data[c2][0]
	else:
	self.current+=state_vec.data[c2][0]
	# This is to jump out of the while
	# because looking along a row of the
	# time capacitor is found is enough.
	c1=len(sys_loops)
	# Increment the counter if capacitor
	# is not found.
	c1+=1
	class Voltage_Source:
	def __init__(self, volt_index, volt_pos, volt_tag):
	self.volt_number=volt_index
	self.volt_pos=volt_pos
	self.volt_tag=volt_tag
	self.v_peak=120.0
	self.v_freq=60.0
	self.v_phase=0.0
	self.voltage=0.0
	self.current=0.0
	self.op_value=0.0
	self.v_polrty=[-1, -1]
	def display(self):
	print "Voltage Source is ",
	print self.volt_tag,
	print "of %f V(peak), %f Hz(frequency) and %f (degrees phase shift)" %(self.v_peak, self.v_freq, self.v_phase),
	print " located at ",
	print self.volt_pos,
	print " with positive polarity towards %s" %(csv_element(self.v_polrty))
	def ask_values(self, x_list, ckt_mat):
	volt_params=["VoltageSource"]
	volt_params.append(self.volt_tag)
	volt_params.append(self.volt_pos)
	volt_params.append("Peak (Volts) = %f" %self.v_peak)
	volt_params.append("Frequency (Hertz) = %f" %self.v_freq)
	volt_params.append("Phase (degrees) = %f" %self.v_phase)
	if self.v_polrty==[-1, -1]:
	# Looking for a default value of polarity
	# in the neighbouring cells
	self.volt_elem=csv_tuple(self.volt_pos)
	if self.volt_elem[0]>0:
	if ckt_mat[self.volt_elem[0]-1][self.volt_elem[1]]:
	self.v_polrty=[self.volt_elem[0]-1, self.volt_elem[1]]
	if self.volt_elem[1]>0:
	if ckt_mat[self.volt_elem[0]][self.volt_elem[1]-1]:
	self.v_polrty=[self.volt_elem[0], self.volt_elem[1]-1]
	if self.volt_elem[0]<len(ckt_mat)-1:
	if ckt_mat[self.volt_elem[0]+1][self.volt_elem[1]]:
	self.v_polrty=[self.volt_elem[0]+1, self.volt_elem[1]]
	if self.volt_elem[1]<len(ckt_mat)-1:
	if ckt_mat[self.volt_elem[0]][self.volt_elem[1]+1]:
	self.v_polrty=[self.volt_elem[0], self.volt_elem[1]+1]
	volt_params.append("Positive polarity towards (cell) = %s" %csv_element(self.v_polrty))
	x_list.append(volt_params)
	def get_values(self, x_list, ckt_mat):
	self.v_peak=float(x_list[0].split("=")[1])
	self.v_freq=float(x_list[1].split("=")[1])
	self.v_phase=float(x_list[2].split("=")[1])
	volt_polrty=x_list[3].split("=")[1]
	# Convert the human readable form of cell
	# to [row, column] form
	while volt_polrty[0]==" ":
	volt_polrty=volt_polrty[1:]
	self.v_polrty=csv_tuple(volt_polrty)
	if not ckt_mat[self.v_polrty[0]][self.v_polrty[1]]:
	print "Polarity incorrect. Branch does not exist at %s" %csv_element(self.v_polrty)
	def transfer_to_sys(self, sys_loops, mat_e, mat_a, mat_b, mat_u, source_list):
	""" The matrix B in E.dx/dt=Ax+Bu will be updated by the
	polarity of the voltage source."""
	for c1 in range(len(sys_loops)):
	if csv_tuple(self.volt_pos) in sys_loops[c1][c1]:
	# If the positive polarity appears before the voltage position
	# it means as per KVL, we are moving from +ve to -ve
	# and so the voltage will be taken negative
	if sys_loops[c1][c1].index(self.v_polrty)<sys_loops[c1][c1].index(csv_tuple(self.volt_pos)):
	mat_b.data[c1][source_list.index(self.volt_pos)]=-1.0
	else:
	mat_b.data[c1][source_list.index(self.volt_pos)]=1.0
	return
	def generate_val(self, source_lst, sys_loops, mat_u, t, dt):
	""" The source voltage is updated in the matrix u in
	E.dx/dt=Ax+Bu ."""
	self.voltage=self.v_peak*math.sin(2*math.pi*self.v_freq*t+self.v_phase)
	mat_u.data[source_lst.index(self.volt_pos)][0]=self.voltage
	self.op_value=self.voltage
	def update_val(self, sys_loops, mat_e, mat_a, mat_b, state_vec, mat_u):
	pass
	class Ammeter:
	def __init__(self, amm_index, amm_pos, amm_tag):
	self.amm_number=amm_index
	self.amm_pos=amm_pos
	self.amm_tag=amm_tag
	self.current=0.0
	self.op_value=0.0
	self.amm_polrty=[-1, -1]
	def display(self):
	print "Ammeter is ",
	print self.amm_tag,
	print " located at ",
	print self.amm_pos,
	print " with positive polarity towards %s" %(csv_element(self.amm_polrty))
	def ask_values(self, x_list, ckt_mat):
	amm_params=["Ammeter"]
	amm_params.append(self.amm_tag)
	amm_params.append(self.amm_pos)
	if self.amm_polrty==[-1, -1]:
	# Looking for a default value of polarity
	# in the neighbouring cells
	self.amm_elem=csv_tuple(self.amm_pos)
	if self.amm_elem[0]>0:
	if ckt_mat[self.amm_elem[0]-1][self.amm_elem[1]]:
	self.amm_polrty=[self.amm_elem[0]-1, self.amm_elem[1]]
	if self.amm_elem[1]>0:
	if ckt_mat[self.amm_elem[0]][self.amm_elem[1]-1]:
	self.amm_polrty=[self.amm_elem[0], self.amm_elem[1]-1]
	if self.amm_elem[0]<len(ckt_mat)-1:
	if ckt_mat[self.amm_elem[0]+1][self.amm_elem[1]]:
	self.amm_polrty=[self.amm_elem[0]+1, self.amm_elem[1]]
	if self.amm_elem[1]<len(ckt_mat)-1:
	if ckt_mat[self.amm_elem[0]][self.amm_elem[1]+1]:
	self.amm_polrty=[self.amm_elem[0], self.amm_elem[1]+1]
	amm_params.append("Positive polarity towards (cell) = %s" %csv_element(self.amm_polrty))
	x_list.append(amm_params)
	def get_values(self, x_list, ckt_mat):
	amm_polrty=x_list[0].split("=")[1]
	# Convert the human readable form of cell
	# to [row, column] form
	while amm_polrty[0]==" ":
	amm_polrty=amm_polrty[1:]
	self.amm_polrty=csv_tuple(amm_polrty)
	if not ckt_mat[self.amm_polrty[0]][self.amm_polrty[1]]:
	print "Polarity incorrect. Branch does not exist at %s" %csv_element(self.amm_polrty)
	def transfer_to_sys(self, sys_loops, mat_e, mat_a, mat_b, mat_u, source_list):
	pass
	def generate_val(self, source_lst, sys_loops, mat_u, t, dt):
	pass
	def update_val(self, sys_loops, mat_e, mat_a, mat_b, state_vec, mat_u):
	""" The ammeter current is calculated as a result of the KVL."""
	c1=0
	while c1<len(sys_loops):
	# The first occurance of the ammeter
	if csv_tuple(self.amm_pos) in sys_loops[c1][c1]:
	# Check for polarity. If positive polarity is after the ammeter
	# it means the current is in the direction of ammeter current and so
	# current is positive.
	if sys_loops[c1][c1].index(self.amm_polrty)>sys_loops[c1][c1].index(csv_tuple(self.amm_pos)):
	self.current=state_vec.data[c1][0]
	else:
	self.current=-state_vec.data[c1][0]
	# Check for ammeter in other loops with the ammeter
	# as common branch
	for c2 in range(c1+1, len(sys_loops)):
	# These lists are nested lists as two loops
	# may have separate branches in common
	for c3 in range(len(sys_loops[c1][c2])):
	if csv_tuple(self.amm_pos) in sys_loops[c1][c2][c3]:
	# Check for polarity again
	if sys_loops[c1][c2][c3].index(self.amm_polrty)>sys_loops[c1][c2][c3].index(csv_tuple(self.amm_pos)):
	# Then check is the loop is aiding or opposing
	# the main loop.
	if sys_loops[c1][c2][c3][-1]=="forward":
	self.current+=state_vec.data[c2][0]
	else:
	self.current-=state_vec.data[c2][0]
	else:
	if sys_loops[c1][c2][c3][-1]=="forward":
	self.current-=state_vec.data[c2][0]
	else:
	self.current+=state_vec.data[c2][0]
	# This is to jump out of the while
	# because looking along a row of the
	# time ammeter is found is enough.
	c1=len(sys_loops)
	# Increment the counter if ammeter
	# is not found.
	c1+=1
	# Since it is a meter, this is its output value
	self.op_value=self.current
	component_list={"resistor":Resistor, "inductor":Inductor, "capacitor":Capacitor, "voltagesource":Voltage_Source, "ammeter":Ammeter}
	nw_input=raw_input("CSV file containing the network layout --> ")
	nw_layout=nw_input+".csv"
	test_ckt=open(nw_layout,"r")
	# Read the circuit into tst_mat
	# Also performs a scrubbing of tst_mat
	tst_mat=csv_reader(test_ckt)
	components_found={}
	for c1 in range(len(tst_mat)):
	for c2 in range(len(tst_mat[0])):
	elem=tst_mat[c1][c2]
	if elem:
	# wire is a zero resistance connection
	if elem.lower()!="wire":
	if len(elem.split("_"))==1:
	jump_det=elem.split("_")[0]
	if len(jump_det)>3:
	if jump_det.lower()[0:4]=="jump":
	pass
	else:
	print "Error! Component at %s does not have a unique name/tag." %csv_element([c1, c2])
	else:
	print "Error! Component at %s does not have a unique name/tag." %csv_element([c1, c2])
	else:
	[elem_name, elem_tag]=elem.split("_")
	elem_type=elem_name.lower()
	if elem_type[0]==" ":
	elem_type=elem_type[1:]
	if elem_tag[0]==" ":
	elem_tag=elem_tag[1:]
	# Check if component exists
	if elem_type in component_list.keys():
	# If found for the first time
	# Create that dictionary element with key
	# as component type
	if elem_type not in components_found:
	components_found[elem_type]=[[csv_element([c1, c2]), elem_tag]]
	else:
	# If already found, append it to
	# dictionary item with that key.
	components_found[elem_type].append([csv_element([c1, c2]), elem_tag])
	else:
	print "Error! Component at %s doesn't exist." %csv_element([c1, c2])
	# Check if a component of the same type has the same tag.
	for items in components_found.keys():
	for c1 in range(len(components_found[items])):
	for c2 in range(len(components_found[items])):
	if c1!=c2:
	if components_found[items][c1][1]==components_found[items][c2][1]:
	print "Duplicate labels found for components of type %s at %s and %s" %(items, components_found[items][c1][0], components_found[items][c2][0])
	component_objects={}
	for items in components_found.keys():
	# Take every type of component found
	# item -> resistor, inductor etc
	for c1 in range(len(components_found[items])):
	# Each component type will be occurring
	# multiple times. Iterate through every find.
	# The list corresponding to each component is
	# the unique cell position in the spreadsheet
	component_objects[components_found[items][c1][0]] = \
	component_list[items](c1+1, components_found[items][c1][0], components_found[items][c1][1])
	parameters_file=nw_input+"_params.csv"
	# Check if the *_params.csv file exists.
	try:
	csv_check_values=open(parameters_file,"r")
	# If not, it has to be created and filled
	# with default values.
	except:
	#param_flag="no"
	pass
	# Check if any of the components with the same
	# tags are present in nw_params.csv. If so, take
	# those parameters from nw_params.csv and replace
	# the default parameters in the component objects.
	else:
	params_from_file=reading_params(csv_check_values)
	for c1 in range(len(params_from_file)):
	# Remove leading spaces if any
	# The first column is the type of element
	if params_from_file[c1][0][0]==" ":
	params_from_file[c1][0]=params_from_file[c1][0][1:]
	name_from_file=params_from_file[c1][0].lower()
	# Check if the component type
	# exists in the new circuit
	try:
	components_found[name_from_file]
	except:
	# If the component doesn't exist, just don't
	# try to read the parameters.
	pass
	else:
	# If it exists, check for the tag.
	for c2 in range(len(components_found[name_from_file])):
	# Remove leading spaces if any
	if params_from_file[c1][1][0]==" ":
	params_from_file[c1][1]=params_from_file[c1][1][1:]
	# Check if the component tag exists in
	# components found so far
	if params_from_file[c1][1]==components_found[name_from_file][c2][1]:
	# If so take the parameters and move them into the object
	# having of that type and having the new cell position
	component_objects[components_found[name_from_file][c2][0]].get_values(params_from_file[c1][3:], tst_mat)
	csv_check_values.close()
	values_to_file=[]
	for items in component_objects.keys():
	# Each component object has a method
	# ask_values that prints in the csv file
	# default values for parameters.
	component_objects[items].ask_values(values_to_file, tst_mat)
	csv_ask_values=open(parameters_file,"w")
	for c1 in range(len(values_to_file)):
	for c2 in range(len(values_to_file[c1])):
	csv_ask_values.write("%s" %values_to_file[c1][c2])
	csv_ask_values.write(", ")
	csv_ask_values.write("\n")
	csv_ask_values.close()
	# Wait for the user to enter parameters before
	# reading the nw_params.csv file.
	cont_ans="n"
	while cont_ans.lower()!="y":
	cont_ans=raw_input("Enter parameters in file %s. When ready press y and enter to continue -> " %parameters_file)
	print
	csv_get_values=open(parameters_file,"r")
	params_from_file=reading_params(csv_get_values)
	csv_get_values.close()
	for c1 in range(len(params_from_file)):
	# Getting rid of the beginning spaces
	# in the component keys
	if params_from_file[c1][2][0]==" ":
	params_from_file[c1][2]=params_from_file[c1][2][1:]
	component_objects[params_from_file[c1][2]].get_values(params_from_file[c1][3:], tst_mat)
	# Just checking the objects
	for items in component_objects.keys():
	component_objects[items].display()
	node_list, branch_map, loop_list, loop_branches, conn_matrix, \
	[number_of_nodes, number_of_branches, loop_count] = network_solver(nw_layout)
	loop_count=len(loop_branches)
	print "*"*50
	print "Number of nodes",
	print number_of_nodes
	print "Number of branches",
	print number_of_branches
	print "Number of loops",
	print loop_count
	print "*"*50
	def human_loop(loop):
	""" Takes a loop as a list of tupes.
	And prints a series of elements in spreadsheet format. """
	for c1 in range(len(loop)):
	print csv_element(loop[c1]),
	return
	for c1 in range(len(loop_branches)):
	human_loop(loop_branches[c1])
	print
	print "*"*50
	def comm_elem_in_loop(loop1, loop2):
	""" Takes two loops and returns a list
	which has all the elements (tuples) that are
	common between the loops. """
	loop_comm=[]
	# Check every element of loop1
	# w.r.t to every element of loop2
	for c1 in range(len(loop1)):
	for c2 in range(len(loop2)):
	# Check if elements are equal
	if loop1[c1]==loop2[c2]:
	# Check if either of the elements
	# are the last of the loops
	if c1<len(loop1)-1 and c2<len(loop2)-1:
	# Check if they have already been
	# identified as common elements
	if loop1[c1] not in loop_comm:
	loop_comm.append(loop1[c1])
	elif loop2[c2-1]==loop_comm[-1]:
	# This is a special condition.
	# The first and last element of
	# every loop are the same.
	# Therefore, it will fail the condition that
	# the element should not be found before.
	# But, it may be possible that the segment of the
	# loops that is common does not contain the last
	# element. So, the check is:
	# If the latest element to be found is loop2[c2-1]
	# i.e is the second last element of loop2, in that
	# case, the last element i.e loop2[c2] should
	# also be appended to the common segment
	loop_comm.append(loop2[c2])
	return loop_comm
	def comm_branch_in_loop(loop1, loop2):
	""" Takes the common elements (loop1) found between
	two loops (out of which one is loop2) and break these
	elements up into separate branches."""
	# The collection of branches
	loop_comm=[]
	# Each branch
	loop_segment=[]
	# starting element
	prev_elem=loop1[0]
	# Iterate from second to last element
	for c1 in range(1, len(loop1)):
	# Check if the index between this current element
	# and the previous element is less than or equal to 1
	# This means it is a continuation of a branch
	if abs(loop2.index(loop1[c1])-loop2.index(prev_elem))<=1:
	loop_segment.append(prev_elem)
	# If not, it means it is a new branch
	else:
	# Complete the branch with the previous element
	loop_segment.append(prev_elem)
	# If it is the final element of the loop
	# Don't leave it out but add that too
	if c1==len(loop1)-1:
	loop_segment.append(loop1[c1])
	# Add that to the collection of branches
	loop_comm.append(loop_segment)
	loop_segment=[]
	# Refresh the previous element
	prev_elem=loop1[c1]
	# This is a special condition.
	# If there is only one common branch, the main
	# condition will fail because a new branch will not be found
	# In that case, the addition function needs to be repeated.
	if not loop_comm:
	loop_segment.append(prev_elem)
	loop_comm.append(loop_segment)
	return loop_comm
	# A temporary list that stores the nodes
	# common to two loops
	nodes_in_loop=[]
	# A array of all the loops of the system
	# including common branches between loops.
	system_loops=[]
	for c1 in range(loop_count):
	row_vector=[]
	for c2 in range(loop_count):
	row_vector.append([])
	system_loops.append(row_vector)
	for c1 in range(len(loop_branches)):
	# The diagonal elements of system_loops
	# will be the loops themselves.
	for c2 in range(len(loop_branches[c1])):
	system_loops[c1][c1].append(loop_branches[c1][c2])
	# The system_loops array will be symmetric
	for c2 in range(c1+1, len(loop_branches)):
	# Find the nodes in loop1
	for c3 in range(len(node_list)):
	if node_list[c3] in loop_branches[c1]:
	nodes_in_loop.append(node_list[c3])
	# Find out the nodes common to loop1
	# and loop2.
	for c3 in range(len(nodes_in_loop)-1, -1, -1):
	if nodes_in_loop[c3] not in loop_branches[c2]:
	del nodes_in_loop[c3]
	# If there are two or more nodes common
	# between loop1 and loop2, there are
	# common elements.
	if len(nodes_in_loop)>1:
	comm_seg_in_loop=comm_elem_in_loop(loop_branches[c1], loop_branches[c2])
	# The list of common branches between
	# loop c1 and loop c2
	sys_loop_off_diag=comm_branch_in_loop(comm_seg_in_loop, loop_branches[c1])
	for c3 in range(len(sys_loop_off_diag)):
	#for c4 in range(len(sys_loop_off_diag[c3])):
	system_loops[c1][c2].append(sys_loop_off_diag[c3])
	system_loops[c2][c1].append(sys_loop_off_diag[c3])
	# Determining the direction of common
	# branches between loops c1 and c2
	for c3 in range(len(sys_loop_off_diag)):
	st_node=sys_loop_off_diag[c3][0]
	end_node=sys_loop_off_diag[c3][-1]
	# Map the start and end nodes of common
	# branch c3 to loops c1 and c2
	loop1_st_node=loop_branches[c1].index(st_node)
	loop1_end_node=loop_branches[c1].index(end_node)
	loop2_st_node=loop_branches[c2].index(st_node)
	loop2_end_node=loop_branches[c2].index(end_node)
	# This check is because the first node of a
	# loop is the same as the last node.
	# So this is to make sure, the correct start
	# and end nodes are mapped, the difference between
	# the indices is compared to the length of the common
	# branch. If not equal, this means, the last node of
	# a loop needs to be taken instead of the first node.
	if (abs(loop1_end_node-loop1_st_node)!= len(sys_loop_off_diag[c3])):
	if (len(loop_branches[c1])-loop1_st_node==len(sys_loop_off_diag[c3])):
	loop1_end_node=len(loop_branches[c1])-1
	if (len(loop_branches[c1])-loop1_end_node==len(sys_loop_off_diag[c3])):
	loop1_st_node=len(loop_branches[c1])-1
	if (abs(loop2_end_node-loop2_st_node)!= len(sys_loop_off_diag[c3])):
	if (len(loop_branches[c2])-loop2_st_node==len(sys_loop_off_diag[c3])):
	loop2_end_node=len(loop_branches[c2])-1
	if (len(loop_branches[c2])-loop2_end_node==len(sys_loop_off_diag[c3])):
	loop2_st_node=len(loop_branches[c2])-1
	if (loop1_st_node<loop1_end_node):
	loop1_dir="ahead"
	else:
	loop1_dir="back"
	if (loop2_st_node<loop2_end_node):
	loop2_dir="ahead"
	else:
	loop2_dir="back"
	if loop1_dir=="ahead" and loop2_dir=="ahead":
	sys_loop_off_diag[c3].append("forward")
	if loop1_dir=="back" and loop2_dir=="back":
	sys_loop_off_diag[c3].append("forward")
	if loop1_dir=="ahead" and loop2_dir=="back":
	sys_loop_off_diag[c3].append("reverse")
	if loop1_dir=="back" and loop2_dir=="ahead":
	sys_loop_off_diag[c3].append("reverse")
	#print c1, c2
	#human_loop(comm_seg_in_loop)
	#print
	#for item in sys_loop_off_diag:
	#print "[",
	#human_loop(item[:-1]),
	#print item[-1],
	#print "]",
	#print
	nodes_in_loop=[]
	system_size=len(loop_branches)
	sys_mat_a=matrix.Matrix(system_size, system_size)
	sys_mat_e=matrix.Matrix(system_size, system_size)
	curr_state_vec=matrix.Matrix(system_size)
	next_state_vec=matrix.Matrix(system_size)
	source_list=[]
	try:
	components_found["voltagesource"]
	except:
	pass
	else:
	for c1 in range(len(components_found["voltagesource"])):
	source_list.append(components_found["voltagesource"][c1][0])
	try:
	components_found["capacitor"]
	except:
	pass
	else:
	for c1 in range(len(components_found["capacitor"])):
	source_list.append(components_found["capacitor"][c1][0])
	if source_list:
	sys_mat_b=matrix.Matrix(system_size, len(source_list))
	sys_mat_u=matrix.Matrix(len(source_list))
	else:
	sys_mat_b=matrix.Matrix(system_size)
	sys_mat_u=0.0
	meter_list=[]
	try:
	components_found["ammeter"]
	except:
	pass
	else:
	for c1 in range(len(components_found["ammeter"])):
	meter_list.append(components_found["ammeter"][c1][0])
	for comps in components_found.keys():
	for c1 in range(len(components_found[comps])):
	comp_pos=components_found[comps][c1][0]
	component_objects[comp_pos].transfer_to_sys(system_loops, sys_mat_e, sys_mat_a, sys_mat_b, sys_mat_u, source_list)
	##mat_ode_reduce(e, a, b)
	validity_flag="no"
	while validity_flag=="no":
	dt=raw_input("Simulation time step in seconds --> ")
	dt=float(dt)
	if dt>0.0:
	validity_flag="yes"
	else:
	print "Error. Invalid time step."
	validity_flag="no"
	validity_flag="no"
	while validity_flag=="no":
	t_limit=raw_input("Duration of simulation in seconds --> ")
	t_limit=float(t_limit)
	if t_limit>0.0:
	validity_flag="yes"
	else:
	print "Error. Invalid time duration."
	validity_flag="no"
	op_dat_file=raw_input("Output data file name (.dat extension will be added) --> ")
	op_dat_file=op_dat_file+".dat"
	f=open(op_dat_file,"w")
	ode_k1=matrix.Matrix(system_size)
	ode_k2=matrix.Matrix(system_size)
	ode_k3=matrix.Matrix(system_size)
	ode_k4=matrix.Matrix(system_size)
	ode_dbydt=matrix.Matrix(system_size)
	ode_var=[ode_k1, ode_k2, ode_k3, ode_k4, ode_dbydt]
	t=0.0
	while t<t_limit:
	for c1 in range(len(source_list)):
	component_objects[source_list[c1]].generate_val(source_list, system_loops, sys_mat_u, t, dt)
	mat_ode(sys_mat_e, sys_mat_a, sys_mat_b, [curr_state_vec, next_state_vec], sys_mat_u, dt, ode_var)
	for comps in component_objects.keys():
	component_objects[comps].update_val(system_loops, sys_mat_e, sys_mat_a, sys_mat_b, next_state_vec, sys_mat_u)
	curr_state_vec=next_state_vec
	#f.write("%s \t %s \t %s \t %s \n" %(str(t), str(sys_mat_u.data[0][0]), str(next_state_vec.data[0][0]), str(next_state_vec.data[1][0])))
	f.write("%s " %str(t),)
	for c1 in range(len(source_list)):
	f.write("%s " %component_objects[source_list[c1]].op_value,)
	for c1 in range(len(meter_list)):
	f.write("%s " %component_objects[meter_list[c1]].op_value,)
	f.write("\n")
	t=t+dt
	print sys_mat_a
	print sys_mat_e
	print sys_mat_b
	f.close()

view raw circuit_solver.py hosted with ❤ by GitHub

The main comments are:

1. The loop currents are calculated at every simulation time step. Currents in every branch and every element will not be calculated unless necessary. For example, at a later stage, saturation in the inductor or heating in the resistor can be included. So then, current will be calculated and used.

2. If you want to know the currents, you connect an ammeter. Voltmeters need to be coded. Ammeter was simple as I need to manipulate the loop currents. But voltmeters, I need to think how to use them. If I connect a voltmeter across two nodes, it will be seen as a branch. So it will appear in loop equations. This means, to make sure the currents don't change a large resistance will have to be chosen. However, a simpler method will be to look for the voltmeter tag and exclude that branch from the system matrix. I can then back calculate later.

A sample circuit: