%Professor Frank Fisher
%Stevens Institute of Technology
%Me345: Modeling and Simulation

clear
clc

r1 = 115.7 * 10^-9;
r2 = 50.76 * 10^-9;
r3 = 51.77 * 10^-9;
r4 = 9.50 * 10^-9;
r5 = 7.23 * 10^-9;
r6 = 136.35 * 10^-9;

q_left = [5*10^-4; 0; 0; 0; 0; -5*10^-4];

bigK = [r1, -r1, 0, 0, 0, 0;
       -r1, r1+r2+r3, -r2, -r3, 0, 0;
       0, -r2, r2+r4, 0, -r4, 0;
       0, -r3, 0, r3+r5, -r5, 0;
       0, 0, -r4, -r5, r4+r5+r6, -r6;
       0, 0, 0, 0, -r6, r6];

%now set p6 = 0 (relative pressure)... cut the matrix size as would do by
%hand

q_left2=q_left(1:5);
bigK2(1:5,1:5)=bigK(1:5,1:5);

disp('These are the relative pressures with p6 = 0 and q1 = q1')
rel_press = inv(bigK2) * q_left2   %these are relative pressures

disp(' ')
disp('To check, we should be able to multiply our original 6x6 BigK Global');
disp('Stiffness Matrix by our found pressures plus p6 = 0 and get our known flows');
disp(' ')

pressures_all(1:5)=rel_press;       %these are the five pressures we found
pressures_all(6)=0                  %this was our BC
check = bigK * pressures_all'       %'check' should be our known flows;
                                    %need transpose so multiplication
                                    %works!
                                    
% Can now calculate flows in each pipe; Q = (pi*D^4/128*L*mu)*(pi-pj)

q1 = ((3.14159*0.1^4)/(128*70.71*0.3)) * (pressures_all(1)-pressures_all(2))
q2 = ((3.14159*0.075^4)/(128*50.99*0.3)) * (pressures_all(2)-pressures_all(3))
q3 = ((3.14159*0.075^4)/(128*50*0.3)) * (pressures_all(2)-pressures_all(4))
q4 = ((3.14159*0.05^4)/(128*53.85*0.3)) * (pressures_all(3)-pressures_all(5))
q5 = ((3.14159*0.05^4)/(128*70.71*0.3)) * (pressures_all(4)-pressures_all(5))
q6 = ((3.14159*0.1^4)/(128*60*0.3)) * (pressures_all(5)-pressures_all(6))

