function [M,K] = makeMK2(h,r,branches,Length,rl);

rb = find(branches==0);                 % roots of dendrites
D = length(rb);                         % number of dendrites

% process first dendrite

disp(' ')
disp('*** PROCESSING DENDRITE 1 ***')

if 1 < D,
   ind = 1:rb(2)-1;
else
   ind = 1:length(branches);
end

sind = length(ind);

disp(['It has ' num2str(sind) ' branches']);

%disp([' ind = ' num2str(ind)])

disp(['branches = ' num2str(branches(ind))]);
[M,K] = makeMK1(h,r(ind),branches(ind),Length(ind),rl(ind));

disp(['M(1,1) = ' num2str(M(1,1))]);
disp(['K(1,1) = ' num2str(K(1,1))]);

% process remaining dendrites

for d = 2:D,

    disp(' ')
    disp(['*** PROCESSING DENDRITE ' num2str(d) ' ***'])

    if d < D
       ind = rb(d):rb(d+1)-1;
    else
       ind = rb(d):length(branches);
    end

    disp(['It has ' num2str(length(ind)) ' branch(es)']);
%    disp([' ind = ' num2str(ind)])

    disp(['branches = ' num2str(branches(ind))]);
    [Mtemp,Ktemp] = makeMK1(h,r(ind),branches(ind),Length(ind),rl(ind));

    low = size(M,1)+1;
    Nd = size(Mtemp,1);

    disp(['Nd = ' num2str(Nd) ' low = ' num2str(low)]);

    M(1,1) = Mtemp(1,1) + M(1,1);
    M(1,low) = Mtemp(1,2);
    M(low,1) = M(1,low);

    K(1,1) = Ktemp(1,1) + K(1,1);
    K(1,low) = Ktemp(1,2);
    K(low,1) = K(1,low);

    disp(['M(1,1) = ' num2str(M(1,1))]);
    disp(['K(1,1) = ' num2str(K(1,1))]);

    M(low:low+Nd-2,low:low+Nd-2) = Mtemp(2:Nd,2:Nd);
    K(low:low+Nd-2,low:low+Nd-2) = Ktemp(2:Nd,2:Nd);
   
end

figure(2)
spy(M)
grid minor
disp(' ')

%

function answer = branchID(branch_number,BranchIndex)

%   A parent can be a child to a parent or it is connected to the soma
%   A child can be a parent or an ending
%   Therefore any branch can be:
%    answer =  1 if a child and an ending 
%              2 if a parent to children and a child to a parent 
%              3 if a parent connected to a soma
    
    answer=0;
    button = 1;

    for j = 1:length(BranchIndex),
    
        if floor(BranchIndex(j)/3) == branch_number,
            button = 0;
        end
    end

    if button,
        answer = 1;
    elseif branch_number == 0,
        answer = 3;
    elseif ~button,
        answer = 2;
    end

    return

%
%   [M,K] = makeMK1(h,r,branches,length,r1)
%
%   h is the space step along the dendrite
%   r is the cell array of radii of points on all dendrites
%   As is the surface areas of the cell
%   branches is the vector of branch numbers in base 3
%   length is the length of each dendrite
%   rl are the lengths between points of radii on all dendrites

function [M,K] = makeMK1(h,r,branches,Length,rl)

%   To Do list
%   1. Convert r from cell to integer before processing (done)
%   2. Handle the entire cell:
%       Break up all passed variables by branches
%   3. Determine new MaxValue Index (done)
%   4. Connect the Soma to all branches and all branches to each other

counter = 1;

numberOfBranches = length(branches);

% check for lone branch (only possible as root)

if numberOfBranches == 1,

   [MaxValue,MaxValueIndex] = max(branches);

   parent = find(branches==floor(MaxValue/3));
  
   r1 = str2num(char(r{parent}));
  
   lr1 = Length(parent)/h;
  
   r1 = interp1(str2num(char(rl{parent})),r1,0:h:lr1*h);
%   disp([' r1 = ' num2str(r1)])
  
   N1 = length(r1);     

   sa = sqrt(h^2+diff(r1).^2);

   Kd(1) = r1(1)^2 + r1(2)^2 + r1(1)*r1(2);
   Ko(1) = -Kd(1);

   Md(1) = sa(1) * (r1(1)/4 + r1(2)/12);
   Mo(1) = sa(1) * (r1(1) + r1(2))/12;

   for j=2:N1-1,

       ra = r1(j-1)^2 + r1(j)^2 + r1(j)*r1(j-1);
       rb = r1(j+1)^2 + r1(j)^2 + r1(j)*r1(j+1);

       Kd(j) = ra + rb;
       Ko(j) = -rb;

       Md(j) = sa(j-1) * (r1(j)/4 + r1(j-1)/12) + sa(j)*(r1(j)/4+r1(j+1)/12);
       Mo(j) = sa(j)*(r1(j)+r1(j+1))/12;

   end
        
   Kd(N1) = (r1(N1-1)^2 + r1(N1)^2 + r1(N1)*r1(N1-1));
   Md(N1) = sa(N1-1)*(r1(N1)/4 + r1(N1-1)/12);

   K = (1/(3*h)) * spdiags([[Ko,0]' Kd' [0,Ko]'], -1:1, N1,N1);

   M = spdiags([[Mo,0]' Md' [0,Mo]'], -1:1, N1,N1);

   return

end

numberOfNodes = (numberOfBranches -1)/2;

BranchIndex = branches;

while counter <= numberOfNodes,
   
   [MaxValue,MaxValueIndex] = max(branches);
      
   parent = find(branches==floor(MaxValue/3));
   
   disp(' ')

   disp(['Parent = ' num2str(floor(MaxValue/3))  ...
         '  Children = ' num2str(MaxValue-1) ' and ' num2str(MaxValue)]);
   
%   bo = [bo floor(MaxValue/3) MaxValue-1 MaxValue];  branch order
   %converts the cell array into vectors    

   r1 = str2num(char(r{parent}));
   r2 = str2num(char(r{MaxValueIndex-1}));
   r3 = str2num(char(r{MaxValueIndex}));
   
   lr1 = Length(parent)/h;
   lr2 = Length(MaxValueIndex-1)/h;
   lr3 = Length(MaxValueIndex)/h;
   
   r1 = interp1(str2num(char(rl{parent})),r1,0:h:lr1*h);
%   disp([' r1 = ' num2str(r1)])
   r2 = interp1(str2num(char(rl{MaxValueIndex-1})),r2,0:h:lr2*h);
%   disp([' r2 = ' num2str(r2)])
   r3 = interp1(str2num(char(rl{MaxValueIndex})),r3,0:h:lr3*h);
%   disp([' r3 = ' num2str(r3)])
   
   N1 = length(r1);	% if r1 is not the root then N1 is too big, use nr
   N2 = length(r2)-1;
   N3 = length(r3)-1;
   
   N = N1 + N2 + N3;
   
   disp(['size of r1, r2, r3 : ' num2str(N1) ' ' num2str(N2) ' ' num2str(N3)]);
       
   [Mtemp,Ktemp] = makeMK(h,r1,r2,r3);
   disp('calling makeMK')
   
   disp(['size of Mtemp = ' num2str(size(Mtemp,2))]);
   
   counter = counter +1;
   branches(MaxValueIndex-1) = -1;
   branches(MaxValueIndex) = -1;
        
   % number of space steps before start of parent and children
        
   lrest1 = sum(Length(1:parent-1))/h;
   nr = 0;
   if lrest1 > 0,		% if parent is not root, fudge
      lrest1 = lrest1 + 1;
      nr = 1;			% use this to correct N1
   end
   lrest2 = 1+sum(Length(1:MaxValueIndex-2))/h;
   lrest3 = 1+sum(Length(1:MaxValueIndex-1))/h;
       
   %  copy branch end of 1 and beginning of 3 into Matrix M and K. 
   %  Since processing is made continous at any node,
   %  it doesn't matter whether or not the children are ends.
                
    % beginning of 3 and end of 1 

    M(lrest1+N1-nr,lrest3+1) = Mtemp(N1,N1+N2+1);
    K(lrest1+N1-nr,lrest3+1) = Ktemp(N1,N1+N2+1);

    figure(2)
    spy(M)
    %pause

     M(lrest3+1,lrest1+N1-nr) = M(lrest1+N1-nr,lrest3+1);
     K(lrest3+1,lrest1+N1-nr) = K(lrest1+N1-nr,lrest3+1);

figure(2)
    spy(M)
    %pause

     
     % beginning of 2 and end of 1
     
     M(lrest1+N1-nr,lrest2+1) = Mtemp(N1,N1+1);
     K(lrest1+N1-nr,lrest2+1) = Ktemp(N1,N1+1);

figure(2)
    spy(M)
    %pause

     
     M(lrest2+1,lrest1+N1-nr) = M(lrest1+N1-nr,lrest2+1);
     K(lrest2+1,lrest1+N1-nr) = K(lrest1+N1-nr,lrest2+1);

figure(2)
    spy(M)
    %pause

     
     % create temporary array to store diagonal info

     Mtemp1 = Mtemp(1+nr:N1,1+nr:N1);
     Mtemp2 = Mtemp(N1+1:N1+N2,N1+1:N1+N2);
     Mtemp3 = Mtemp(N1+N2+1:N1+N2+N3,N1+N2+1:N1+N2+N3);
     
     Ktemp1 = Ktemp(1+nr:N1,1+nr:N1);
     Ktemp2 = Ktemp(N1+1:N1+N2,N1+1:N1+N2);
     Ktemp3 = Ktemp(N1+N2+1:N1+N2+N3,N1+N2+1:N1+N2+N3);
     
     FirstChildID = branchID(MaxValue,BranchIndex);
     SecondChildID = branchID(MaxValue-1,BranchIndex);
     ParentID = branchID(floor(MaxValue/3),BranchIndex);
     
     disp(['lrest1 = ' num2str(lrest1) ...
           '  lrest2 = ' num2str(lrest2) '  lrest3 = ' num2str(lrest3)]);
    
     if FirstChildID == 1 & SecondChildID == 1,
         
         disp('Neither Children have Children');
         
         %  copy the branch info into Matrix M and K

         M(lrest1+1:lrest1+N1-nr,lrest1+1:lrest1+N1-nr) = Mtemp1;
         M(lrest2+1:lrest2+N2,lrest2+1:lrest2+N2) = Mtemp2;
         M(lrest3+1:lrest3+N3,lrest3+1:lrest3+N3) = Mtemp3;
         
         K(lrest1+1:lrest1+N1-nr,lrest1+1:lrest1+N1-nr) = Ktemp1;
         K(lrest2+1:lrest2+N2,lrest2+1:lrest2+N2) = Ktemp2;
         K(lrest3+1:lrest3+N3,lrest3+1:lrest3+N3) = Ktemp3;
        
     elseif FirstChildID == 1,
         
         disp(['Branch number ' num2str(MaxValue-1) ' has children']);
         
         M(lrest3+1:lrest3+N3,lrest3+1:lrest3+N3) = Mtemp3;
         
         K(lrest3+1:lrest3+N3,lrest3+1:lrest3+N3) = Ktemp3;
         
         % place elements (upper half) from children into corresponding
         % block in M&K
         
         for j=1:N2,
             for i = 1:N2-j+1,
                 M(lrest2+i,lrest2+j) = Mtemp2(i,j);
                 K(lrest2+i,lrest2+j) = Ktemp2(i,j);
             end
         end
       
         % place lower half of parent
         for j=1:N1-nr,
             for i = N1-nr-j+1:N1-nr,
                 M(lrest1+i,lrest1+j)=Mtemp1(i,j);
                 K(lrest1+i,lrest1+j)=Ktemp1(i,j);
             end
         end
         
     elseif SecondChildID ==1,
         
         disp(['Branch number ' num2str(MaxValue) ' has children']);
         
         M(lrest2+1:lrest2+N2,lrest2+1:lrest2+N2) = Mtemp2;
         K(lrest2+1:lrest2+N2,lrest2+1:lrest2+N2) = Ktemp2;
         
         % place lower half of parent

         for j=1:N1-nr,
             for i = N1-nr-j+1:N1-nr,
                 M(lrest1+i,lrest1+j)=Mtemp1(i,j);
                 K(lrest1+i,lrest1+j)=Ktemp1(i,j);
             end
         end
         
         for j=1:N3,
             for i = 1:N3-j+1,
                 M(lrest3+i,lrest3+j) = Mtemp3(i,j);
                 K(lrest3+i,lrest3+j) = Ktemp3(i,j);
             end
         end
         
     else,
         
         disp('Both Children have children');
         
         % place lower half of parent
         for j=1:N1-nr,
             for i = N1-nr-j+1:N1-nr,
                 M(lrest1+i,lrest1+j)=Mtemp1(i,j);
                 K(lrest1+i,lrest1+j)=Ktemp1(i,j);
             end
         end
         
         for j=1:N2,
             for i = 1:N2-j+1,
                 M(lrest2+i,lrest2+j) = Mtemp2(i,j);
                 K(lrest2+i,lrest2+j) = Ktemp2(i,j);
             end
         end
         
         for j=1:N3,
             for i = 1:N3-j+1,
                 M(lrest3+i,lrest3+j) = Mtemp3(i,j);
                 K(lrest3+i,lrest3+j) = Ktemp3(i,j);
             end
         end
         
     end
     
     
     if ParentID == 3,
         
         disp([num2str(floor(MaxValue/3)) ' is a root']);
         full(Mtemp1);
         
         for j=1:N1-nr,
             for i = 1:N1-nr-j+1,
                 M(lrest1+i,lrest1+j) = Mtemp1(i,j);
                 K(lrest1+i,lrest1+j) = Ktemp1(i,j);
             end
         end

     end

     %keyboard

 end
 
 return
 
%
%	make M & K for the wishbone, convention
%
%      r1(n1) = r2(1) = r3(1)	where n1 = length(r1)
%

function [M,K] = makeMK(h,r1,r2,r3)

n1 = length(r1);
n2 = length(r2);
n3 = length(r3);

N = n1+n2+n3-2;
Kd = zeros(1,N);
Md = Kd;
Ko = zeros(1,N-1);
Mo = Ko;

r = r1;
n = n1;

Kd(1) = r(1)^2 + r(2)^2 + r(1)*r(2);
Ko(1) = -Kd(1);

sa = sqrt(h^2+diff(r).^2);
sa1 = sa;
sa2 = sqrt(h^2+diff(r2).^2);
sa3 = sqrt(h^2+diff(r3).^2);

Md(1) = sa(1) * (r(1)/4 + r(2)/12);
Mo(1) = sa(1) * (r(1) + r(2))/12;

for j=2:n1-1,

    ra = r(j-1)^2 + r(j)^2 + r(j)*r(j-1);
    rb = r(j+1)^2 + r(j)^2 + r(j)*r(j+1);

    Kd(j) = ra + rb;
    Ko(j) = -rb;

    Md(j) = sa(j-1) * (r(j)/4 + r(j-1)/12) + sa(j)*(r(j)/4+r(j+1)/12);
    Mo(j) = sa(j)*(r(j)+r(j+1))/12;

end

Kd(n) = r(n-1)^2 + r(n)^2 + r(n-1)*r(n) + ...
        r2(1)^2 + r2(2)^2 + r2(1)*r2(2) + ...
        r3(1)^2 + r3(2)^2 + r3(1)*r3(2) ;

Ko(n) = -(r2(1)^2 + r2(2)^2 + r2(1)*r2(2));

Md(n) = sa(n-1)*(r(n)/4 + r(n-1)/12) + ...
        sa2(1)*(r2(1)/4 + r2(2)/12) + ...
        sa3(1)*(r3(1)/4 + r3(2)/12) ;

r = r2;
sa = sa2;
Mo(n) = sa(1)*(r2(1)+r2(2))/12;

for j=2:n2-1,

    ra = r(j-1)^2 + r(j)^2 + r(j)*r(j-1);
    rb = r(j+1)^2 + r(j)^2 + r(j)*r(j+1);

    Kd(n1+j-1) = ra + rb;
    Ko(n1+j-1) = -rb;

    Md(n1+j-1) = sa(j-1) * (r(j)/4 + r(j-1)/12) + sa(j)*(r(j)/4+r(j+1)/12);
    Mo(n1+j-1) = sa(j)*(r(j)+r(j+1))/12;

end
        
Kd(n1+n2-1) = (r(n2-1)^2 + r(n2)^2 + r(n2)*r(n2-1));
Md(n1+n2-1) = sa(n2-1)*(r(n2)/4 + r(n2-1)/12);

r = r3;
sa = sqrt(h^2+diff(r).^2);

for j=2:n3-1,

    ra = r(j-1)^2 + r(j)^2 + r(j)*r(j-1);
    rb = r(j+1)^2 + r(j)^2 + r(j)*r(j+1);

    Kd(n1+n2+j-2) = ra + rb;
    Ko(n1+n2+j-2) = -rb;

    Md(n1+n2+j-2) = sa(j-1)*(r(j)/4+r(j-1)/12)+sa(j)*(r(j)/4+r(j+1)/12);
    Mo(n1+n2+j-2) = sa(j)*(r(j)+r(j+1))/12;

end
        
Kd(n1+n2+n3-2) = r(n3-1)^2 + r(n3)^2 + r(n3)*r(n3-1);
Md(n1+n2+n3-2) = sa(n3-1)*(r(n3)/4 + r(n3-1)/12);

n = n1 + n2 + n3 - 2;

K = (1/(3*h)) * spdiags([[Ko,0]' Kd' [0,Ko]'], -1:1, n,n);

K(n1,n1+n2) = -(r3(1)^2 + r3(2)^2 + r3(1)*r3(2))/3/h;
K(n1+n2,n1) = K(n1,n1+n2);

M = spdiags([[Mo,0]' Md' [0,Mo]'], -1:1, n,n);

M(n1,n1+n2) = sa3(1)*(r3(1)+r3(2))/12;
M(n1+n2,n1) = M(n1,n1+n2);

return