% Pressure and Flow Calculator for flow networks
% Written by Dr. Shayn Peirce, U of Virginia
% Modified by A.Marathe on 04/18/06

%User prompted for the data set they want to use
filename = input('Enter the name of the data set you want to use\n and include .txt extension\n', 's');
%Loads the data file (a text file)
data = load(filename)
% User prompted for number of nodes
nnodes = input('Enter the number of nodes:\n');
Numberofnodes = nnodes;
% User prompted for number of elements
nelements = input('Enter the number of elements:\n');
Numberofelements = nelements;
% User prompted for height of channels
h = input('Enter the height of microchannels(um):\n');
h = h/1000; %UNITS:mm
% User prompted for number of pressures known
nknowp = input('Enter the number of known pressures:\n');
Numberofpressuresknown = nknowp;
%Initialize the P vector
P = zeros(nnodes, 1);
%Fill in the P vector with the known pressures
for j=1:nknowp
    NodeP(j)=input('Enter node of known pressure\n');
    Pressure(j)=input('What is the pressure at that node (Pa)\n');
    P(NodeP(j))=Pressure(j); %UNITS: Pa
end
% User prompted for viscosity
Visc = input('Enter the viscosity (cP)\n');
u=Visc/1000000; %UNITS: kg/mm-s

%Create the global nodal matrix, C
%Initialize the matrix:
C=zeros(nnodes,nnodes);
%Tell the size of the data matrix:
[datarow,datacol]=size(data);
%Create the off-diagonal entries (all negative numbers):
for i=1:datarow;
    d=data(i,3)/1000; %UNITS:mm
    l=data(i,4)/1000; %UNITS:mm
    row=data(i,1);
    column=data(i,2);
    %if h<=d
        C(row,column)=-(h^3*d)/(12*u*l); %UNITS: mm^4-s/kg
        %else
        %C(row,column)=-(h*d^3)/(12*u*l); %UNITS: mm^4-s/kg
        %end
    C(column,row)=C(row,column);
end
%Create the entries on the diagonal of matrix (-sum of off-diagonal elements):
j=1;
while j<=nnodes
    C(j,j)=-sum(C(j,1:nnodes));
    j=j+1;
end

%Create the Q' vector, QQ
%Initialize the Q' (QQ) vector
QQ=zeros(nnodes,1);
%Matrix manipulation to zero out the conductance matrix as defined by the boundary values that were entered
for i=1:nnodes
    if P(i)~=0
        for k=1:nnodes
            if(i==k)
                C(i,k)=C(i,k);
                QQ(i)=C(i,i)*P(i); %UNITS: mm^4-s-Pa/kg
            else   
               C(i,k)=0;
           end
            if (i~=k)
                QQ(k)=QQ(k)-C(k,i)*P(i); %UNITS: mm^4-s-Pa/kg
                C(k,i)=0;
            end
        end
    end           
end

%Solve for the pressures in the network
PP=(C^-1)*QQ ; %UNITS: Pa
for x=1:nnodes
    pressures(x,:)=[x,PP(x)];
end
%User prompted to enter the name of the file to which the data will be saved
pressfile = input('Enter filename to save pressure data to\n', 's');

fid = fopen(pressfile, 'w');
fprintf(fid,'   Pressure in the Network\n\n');
fprintf(fid,'Pressure at node (%g) = %0.5g Pa\n',pressures');
fclose(fid);

%Calculates the flows from the nodal pressures for each segment
for g=1:datarow;
    a(g)=data(g,1);
    b(g)=data(g,2);
    dd=data(g,3)/1000; %UNITS: mm
    l=data(g,4)/1000; %UNITS: mm

%Calculate pressure drop across each element
deltaP=PP(a(g))-PP(b(g)); %UNITS: mmHg
deltaP = deltaP/1000; %UNITS: kPa, or kg/mm-s^2

%Calculate conductance/L for each segment
%if h <=dd
    Conductance=(h^3*dd)/(12*u*l); %UNITS: mm^4-s/kg
    %else
    %Conductance=(h*dd^3)/(12*u*l); %UNITS: mm^4-s/kg
    %end
Conductance;

%Calculate flow from pressure drop and conductance
Q(g)=(deltaP*Conductance)*60; %UNITS: uL/min
Q(g);

flow(g,:)=[a(g),b(g),Q(g)];
end

%User prompted to enter the name of the file to which the data will be saved
storefile = input('Enter filename to save flow data to\n', 's');

fid = fopen(storefile, 'w');
fprintf(fid,'     Flow in the Network\n\n', storefile);
fprintf(fid,'Flow in (%g,%g) = %0.5g uL/min\n',flow');
fclose(fid);

%Display data in the command window
fprintf('\nPressures in the network (Pa):\n\n');
fprintf('Pressure at node (%g) = %0.5f Pa\n',pressures');
fprintf('\nFlow rates in the network (uL/min):\n\n');
fprintf('Flow in (%g,%g) = %0.5g uL/min\n',flow');