# -*- coding: utf-8 -*-
"""
Created on Wed Sep 23 19:14:52 2020

@author: jbobowsk
"""

# Created using Spyder(Python 3.7)
# Solving systems of equations

# To solve a system of equations, we will use the Python module
# 'sympy' and the function 'linsolve'.  We first have to tell Python which 
# variables we will be using. In this first trivial example, we will use "x" 
# as the variable.

# Import sympy:
from sympy import *

# Define our variable:
x = Symbol('x')

# Next, we define our equation(s).  It is assumed that the expression is 
# equal to zero.  For example, the RHS of the equation is x - 2 and the LHS
# is zero.
eqn1 = x - 2

# Now we can implement linsolve([...], [...]). In the first set of 
# square brackets, list the equations.  In the second set of square brackets,
# list the variable that you want to solve for.
soln = linsolve([eqn1], [x])
print(soln)

# Python returns the solution as a `FiniteSet'.  We can covert the FiniteSet
# to list via:
solnList = list(soln)[0]
print(solnList)

# Finally, we can index the list to get the final solution
solnFinal = solnList[0]
print(solnFinal)

# Here's a less trivial system of two equations and two unknowns.
y = Symbol('y')

eqn1 = 10*x - 3*y - 5
eqn2 = -2*x - 4*y - 7

# We now implement linsolve.
soln = list(linsolve([eqn1, eqn2], [x, y]))[0];
print('x =', soln[0], 'and y =', soln[1])

# It is also possible to algebraically solve a system of equations in terms
# of other variables.  Below is a system of 6 complex equations and 6 unknowns.

# The example below solves for the currents in an all-pass filter.
# (Optional experiment in PHYS 231.)

# We must define the six currents that we wish to solve for as variables as
# well as the other relevant variables (R is resistance, C is capacitance,
# L is inductance, v is voltage amplitude, and w is angular freqeuncy.
v = Symbol('v')
R = Symbol('R')
w = Symbol('w')
L = Symbol('L')
C = Symbol('C')

# Here is a method to define all of six of the i0 to i5 symbols in a single line
i0, i1, i2, i3, i4, i5 = symbols('i0:6')

# Before we write our system of equations, first note that in Python you make
# a number complex by following it with a 'j' (no space).  For example,
# to enter z = x + j*y we would type 'z = x + y*1j'.  In these expressions
# j represents sqrt(-1).  Below we evaluate the square of of j which should
# be equal to -1.  
print(1j**2)

# To continue an input on a new line, use the backslash notation.
eqns = [i0 - i1 - i2, i1 - i3 - i4, i2 + i4 - i5,\
        i0*R + 1j*w*i1*L + i4*R + 1j*w*i5*L - v,\
        i0*R + i2/(1j*w*C) + 1j*w*i5*L - v, i3/(1j*w*C) - i4*R - 1j*w*i5*L]

# For some reason, 'linsolve' doesn't work here.  However, we can use 'solve'.
soln = solve(eqns, [i0, i1, i2, i3, i4, i5])

# As you can see the solutions is very complicated!  
print(soln)

# The output of 'solve' is an object called a 'Python dictionary'.
# The individual solutions can be accessed using soln[i0], soln[i1], ...
# Here is the solution for i0:
print(soln[i0])

# It is possible to use 'pprint()' to display the expression a little more nicely:
pprint(soln[i0])

# We can use 'subs([],[],[],...)' to enter in numerical values for the various
# symbols.  I'll import 'math' too so that we can access pi.  Here is a numerical
# value of i0.
import numpy as np
i0num = soln[i0].subs([(v, 1), (C, 10e-9), (L, 10e-3), (R, 1000), (w, 2*np.pi*20e3)])
print(i0num)

# We can force Python to tidy up the number using 'N()'.
print(N(i0num))

# We can use an option with N() to specify how many sig figs to keep.  In the 
# output, the real part of the numerical value shows only 1 sig. fig. because 
# it is exact.
print(N(i0num, 3))

# Here are numerical values for all six of the currents.
print('i0 =', N(soln[i0].subs([(v, 1), (C, 10e-9), (L, 10e-3), (R, 1000), (w, 2*np.pi*20e3)]), 3))
print('i1 =', N(soln[i1].subs([(v, 1), (C, 10e-9), (L, 10e-3), (R, 1000), (w, 2*np.pi*20e3)]), 3))
print('i2 =', N(soln[i2].subs([(v, 1), (C, 10e-9), (L, 10e-3), (R, 1000), (w, 2*np.pi*20e3)]), 3))
print('i3 =', N(soln[i3].subs([(v, 1), (C, 10e-9), (L, 10e-3), (R, 1000), (w, 2*np.pi*20e3)]), 3))
print('i4 =', N(soln[i4].subs([(v, 1), (C, 10e-9), (L, 10e-3), (R, 1000), (w, 2*np.pi*20e3)]), 3))
print('i5 =', N(soln[i5].subs([(v, 1), (C, 10e-9), (L, 10e-3), (R, 1000), (w, 2*np.pi*20e3)]), 3))