/*
			CIRCUIT SOLVER
	      		~~~~~~~~~~~~~~

 This program carries out a steady state phasor analysis of R-C-L circuits.
 The program is called through 'circuit_solve' which has as arguments
	1) the angular frequency for the analysis.
	2) the component list.
	3) the list of nodes which are to be 'grounded' - otherwise
		all voltages are relative.
	4) the 'Selection' list - see below.

 The circuit is defined by a list of components.  Each component is
 described by the component type, name, value and the nodes to which
 it is connected.  The component type is used to determine the component
 characteristics. 

 This version provides for selective printout capability through the
 'selective_print' predicate.  This can be used to print out all components
 connected to specified nodes, or to print out specific component information.  
 If the 'Selection' argument is null, then all components are printed out,
 otherwise the list of nodes/components is used to print out the required
 information.

 New component types can be added easily in the declarative CLP language.
 This is illustrated by a definition of the transformer, a simple small
 signal model for the transistor (note the ease of adding new transistor
 types), and the diode.  Note that the diode can not be analysed using
 phasors.  However in the absence of capacitors, inductors and
 transformers, our phasor analysis program reduces to DC analysis.
 Hence we can apply it to circuits with diodes, in the absence of
 these other components.  An interesting aspect of the diode modeling is
 the use of piecewise linear approximation.  Such modeling can give useful
 results for more complex components (eg. the modeling of the DC
 characteristics of transistors).  Furthermore it can be easily
 implemented within the CLP language.
*/

circuit_solve(W, L, G, Selection) :-
	get_node_vars(L, NV),
	solve(W, L, NV, Handles, G), 
	format_print(Handles, Selection).

get_node_vars([[Comp, Num, X, Ns]|Ls], NV) :-
	get_node_vars(Ls, NV1),
	insert_list(Ns, NV1, NV).
get_node_vars([], []).

insert_list([N|Ns], NV1, NV3) :-
	insert_list(Ns, NV1, NV2),
	insert(N, NV2, NV3).
insert_list([], NV, NV).

insert(N, [[N, V, I]|NV1], [[N, V, I]|NV1]) :- !.
insert(N, [[N1, V, I]|NV1], [[N1, V, I]|NV2]) :-
	insert(N, NV1, NV2).
insert(N, [], [[N, V, c(0,0)]]).

solve(W, [X|Xs], NV, [H|Hs], G) :-
	addcomp(W, X, NV, NV1,H),
		solve(W, Xs, NV1, Hs, G).
solve(W, [], NV, [], G) :- 
		!,
		zero_currents(NV),
		ground_nodes(NV, G).

zero_currents([[N, V, c(X,Y)]|Ls]) :-
	zero(X),
	zero(Y),
	zero_currents(Ls).
zero_currents([]).

/* There is some redundancy setting all node current sums to zero - We only
 need to some of them to zero.  Setting all of them to zero can cause some
 problems because of numerical instability; hence the following code. */

zero(0) :- !.
zero(X) :- 
	ground(X),
	%printf("Is it ground: % ?\n", [X]),
	EPS = 0.001,
	X > -1*EPS,
	X < EPS.

ground_nodes(Vs,[N|Ns]) :-
	ground_node(Vs, N),
	ground_nodes(Vs, Ns).
ground_nodes(Vs, []).
ground_node([[N, c(0,0), I]|Vs], N) :- !.
ground_node([[N1, V, I]|Vs], N) :- ground_node(Vs, N).
ground_node([], N) :- printf("Error could be: node %doesn't exist\n", [N]).


% ***** CLAUSES TO DEFINE COMPONENT CHARACTERISTICS *****

addcomp(W, [Comp2, Num, X, [N1, N2]], NV, NV2, [Comp2, Num, X, [N1, V1, I1], [N2, V2, I2]]):-  
/* For two node components of value X

               o  N1
               |              
         +    ---              
              | |   |          
         V    | |   | I      The relationship between V, I, X is depends
              | |   v        by the component type and is defined by the
         -    ---            'iv_reln' predicate. 
               |              
               o  N2
                      
*/
	c_neg(I1, I2),
	iv_reln(Comp2, I1, V, X, W),
	c_add(V, V2, V1),
	subst([N1, V1, Iold1], [N1, V1, Inew1], NV, NV1),
	subst([N2, V2, Iold2], [N2, V2, Inew2], NV1, NV2),
	c_add(I1,Iold1,Inew1),
	c_add(I2,Iold2,Inew2).

% Specific current/voltage relations for the two terminal components
iv_reln(resistor, I, V, R, W) :- 
	c_mult(I,c(R,0),V).
iv_reln(voltage_source, I, V, V, W). 
iv_reln(isource, I, V, I, W). 
iv_reln(capacitor, I, V, C, W) :- 
	c_mult( c(0,W * C),V,I).
iv_reln(inductor, I, V, L, W) :- 
	c_mult(c(0,W * L),I,V).
iv_reln(connection, I, c(0,0), L, W).
iv_reln(open, c(0,0), V, L, W).
iv_reln(diode, I, V, D, W) :- diode(D, I, V).

diode(in914, c(I,0), c(V, 0)) :-
	V < -100,
	DV = V + 100,
	I = 10*DV.
diode(in914, c(I,0), c(V, 0)) :-
	V >= -100,
	V < 0.6,
	I = 0.001*V.
diode(in914, c(I,0), c(V, 0)) :-
	V >= 0.6,
	DV = V - 0.6,
	I = 100*DV.

addcomp(W, [transistor, Num, X, [N1, N2, N3]], NV, NV3, [transistor, Num, X, [N1, V1, I1], [N2, V2, I2], [N3, V3, I3]]):-  
/* A small signal transistor equivalent for transistor of type X 

    N1  o--------------        ----------------o N3
        +             |        |      <--
                      >        ^       I
       Vin        R   >       <I>                I = (Vin/R) * Gain
                      >        v                   
        -             |        |
    N2  ---------------------------------------o
*/

	transistor(X, R, Gain),
	c_add(I1, I3, IT),
	c_neg(I2, IT),
	c_add(Vin, V2, V1),
	c_mult(I1, c(R, 0), Vin),
	c_mult(I1, c(Gain, 0), I3),
	subst([N1, V1, Iold1], [N1, V1, Inew1], NV, NV1),
	subst([N2, V2, Iold2], [N2, V2, Inew2], NV1, NV2),
	subst([N3, V3, Iold3], [N3, V3, Inew3], NV2, NV3),
	subst([N4, V4, Iold4], [N4, V4, Inew4], NV3, NV4),
	c_add(I1,Iold1,Inew1),
	c_add(I2,Iold2,Inew2),
	c_add(I3,Iold3,Inew3),
	c_add(I4,Iold4,Inew4).

transistor(bc108, 1000, 100).

addcomp(W, [transformer, Num, X, [N1, N2, N3, N4]], NV, NV4, [transformer, Num, X, [N1, V1, I1], [N2, V2, I2], [N3, V3, I3], [N4, V4, I4]]):-  
/* For transformers of X units
        _________  
    N1           >   ____
        +        >  <      N3
                 >  <   +	   Vin:Vout = X:1
       Vin       >  <  Vout
                 >  <   -
        -        >  <____  N4
    N2  _________>
		    X
*/

	c_neg(I1, I2),
	c_neg(I3, I4),
	c_add(Vin, V2, V1),
	c_add(Vout, V4, V3),
	c_mult(Vout, c(X, 0), Vin),
	c_mult(I1, c(X, 0), I4),
	subst([N1, V1, Iold1], [N1, V1, Inew1], NV, NV1),
	subst([N2, V2, Iold2], [N2, V2, Inew2], NV1, NV2),
	subst([N3, V3, Iold3], [N3, V3, Inew3], NV2, NV3),
	subst([N4, V4, Iold4], [N4, V4, Inew4], NV3, NV4),
	c_add(I1,Iold1,Inew1),
	c_add(I2,Iold2,Inew2),
	c_add(I3,Iold3,Inew3),
	c_add(I4,Iold4,Inew4).


subst(X, Y, [X|L1], [Y|L1]) :- !.
subst(X, Y, [Z|L1], [Z|L2]) :- subst(X,Y,L1,L2).
subst(X, Y, [], L2) :- printf("Node list incomplete\n", []).


% ***** COMPLEX NUMBER ARITHMETIC *****

c_mult(c(Re1,Im1),c(Re2,Im2),c(Re3,Im3)) :- 
	Re3 = Re1*Re2 + -1 * Im1*Im2,
	Im3 = Re1*Im2 + Re2*Im1.

c_add(c(Re1,Im1),c(Re2,Im2),c(Re3,Im3)) :- 
	Re3 = Re1 + Re2,
	Im3 = Im1 + Im2.

c_neg(c(Re,Im),c(Re1,Im1)) :-  
	Re1 = -1 * Re, Im1 = -1 *Im .
c_eq(c(Re1,Im1),c(Re2,Im2)) :- 
	Re1 = Re2, Im1 = Im2 .

c_real(c(Re,Im),Re).
c_imag(c(Re,Im),Im).


% ****** PRINTOUT ROUTINES ******

format_print(H, []) :- !, all_print(H).  % If no selection is given then
					 % print out all nodes.
format_print(H, Selection) :-
	selective_print(H, Selection).
					 % Otherwise print the selection.

selective_print(Ls, [N|Ns]) :-
	print_nodes(Ls, N, 0),
	print_comps(Ls, N),
	selective_print(Ls, Ns).
selective_print(Ls, []).

print_nodes([[Comp, Num, X|Nodes]|L], N1, Heading_flag_1) :-
	member(N1, Nodes), !,
	Heading_flag_2 = Heading_flag_1 + 1,
	heading(N1, Heading_flag_2),
	all_print([[Comp, Num, X|Nodes]]),
	print_nodes(L, N1, Heading_flag_2).
print_nodes([[Comp, Num, X|Nodes]|L], N1, Heading_flag) :-
	print_nodes(L, N1, Heading_flag).
print_nodes([], N1, Heading_flag).

print_comps([[Comp, Num, X|Nodes]|L], Num) :- 
	all_print([[Comp, Num, X|Nodes]]),
	print_comps(L, Num).

print_comps([[Comp, Num, X|Nodes]|L], Num1) :-
	print_comps(L, Num1).
print_comps([], Num).

heading(N, 1) :-
	!, printf("\nCOMPONENT CONNECTIONS TO NODE %\n", [N]).
heading(N, X).

member(N1, [[N1, V, I]|Ls]).
member(N1, [[N2, V, I]|Ls]) :-
	member(N1, Ls).

write_units(resistor, X) :-
	printf("% Ohms",[X]).
write_units(capacitor, X) :-
	printf("% Farads",[X]).
write_units(inductor, X) :-
	printf("% Henrys", [X]).
write_units(current_source, X) :-
	printf("% Ampere", [X]).
write_units(voltage_source, X) :-
	printf("% Volts", [X]).
write_units(diode, X) :-
	printf("type %", [X]).
write_units(transistor, X) :-
	printf("type %    (base, emitter, collector)", [X]).
write_units(transformer, X) :-
	printf("ratio of %", [X]).

all_print([[Comp, Num, X|Nodes]|L]) :-
	printf("% %: ", [Comp, Num]),
	write_units(Comp, X),
	printf("\n", []),
	pr_nodes(Nodes),
	all_print(L).

all_print([]).

pr_nodes([[N1, V1, I1]|X]) :-
	printf("	Node %\n", [N1]),
	printf("		Voltage %\n", [V1]),
	printf("		Current %\n", [I1]),
	pr_nodes(X).

pr_nodes([]) :-
	printf("\n", []).

