% Jacobian of boundary value problem
%
% ============================================

function result = BVP_HomT_jac(x,p,YS,YU,J)

global homTds

ups = reshape(x,homTds.nphase,homTds.npoints);
%p = num2cell(p);
% Allocate space for sparse jacobian
ks=homTds.nphase*homTds.npoints+(homTds.nphase-homTds.nu)*homTds.nu+(homTds.nphase-homTds.ns)*homTds.ns;
result = spalloc(ks,ks+1,ks/2);
n=homTds.nphase;
N=homTds.npoints;
x1=ups(:,1); 
% Component 1 (the initial fixed point)
% ===========

A1=homT_jac(x1,p,J);
result(1:n, 1:n) = A1-eye(n);
%Derivatives w.r.t active parameter
%jp=homjacp(x1,p,J);
%result(1:n, ks+1) = jp;
%result(1:2,:),pause


% Component 2 (the iteration conditions)
% ===========
for j=3:N
result((j-2)*n+1:(j-1)*n, (j-2)*n+1:(j-1)*n)=homT_jac(ups(:,j-1),p,J);

result((j-2)*n+1:(j-1)*n, (j-1)*n+1:j*n)=-eye(n);
end

%Derivatives w.r.t active parameter
% for j=3:N
% jp=homjacp(ups(:,j-1),p,J);
% result((j-2)*n+1:(j-1)*n, ks+1)=jp;%jp(:,homTds.ActiveParams);
% end


% Component 3 (UNSTABLE)
% Ricatti blocks from unstable eigenspace
% F(Y_U)=R22Y_U-Y_UR11+E21-Y_UR12Y_U;
% ===========

% D1=R_22Y_U
Q0U = homTds.Q0;
[R11, R12, E21, R22] = homT_RicattiCoeff(Q0U,A1,homTds.nu);

% derivatives of D1 = R22 * YU w.r.t (Y_U)ij
l=n*(N-1); h=N*n;
for j=1:n-homTds.nu
    for i=1:homTds.nu
         idx1=l+i+(j-1)*homTds.nu;
         idx2=h+1+(i-1)*(n-homTds.nu); 
         idx3=h+(n-homTds.nu)+(i-1)*(n-homTds.nu);
         result(idx1,idx2:idx3)=result(idx1,idx2:idx3)+R22(j,1:n-homTds.nu);
    end
end


% derivatives of D2 = YU * R11 w.r.t (Y_U)ij

for j=1:n-homTds.nu
    for s=1:homTds.nu
        idx1=l+1+(j-1)*homTds.nu;
        idx2=l+homTds.nu+(j-1)*homTds.nu;
        idx3=h+j+(s-1)*(n-homTds.nu);
       % result(idx1:idx2,idx3)=result(idx1:idx2,idx3)-(R11(1:homTds.nu,i))';
       result(idx1:idx2,idx3)=result(idx1:idx2,idx3)-(R11(s,1:homTds.nu))';
    end
end


% derivatives of D3 = YU * R12*YU w.r.t (Y_U)ij
D31=-homTds.YU*R12;
for j=1:n-homTds.nu
    for i=1:homTds.nu
        idx1=l+i+(j-1)*homTds.nu;
        idx2=h+1+(i-1)*(n-homTds.nu);
        idx3=h+n-homTds.nu+(i-1)*(n-homTds.nu);
        result(idx1,idx2:idx3)=result(idx1,idx2:idx3)+D31(j,:);
    end
end
D32=-R12*homTds.YU;
for j=1:n-homTds.nu
    for s=1:homTds.nu
        idx1=l+i+(j-1)*homTds.nu;
        idx2=l+homTds.nu+(j-1)*homTds.nu;
        idx3=h+j+(s-1)*(n-homTds.nu);
        result(idx1:idx2,idx3)=result(idx1:idx2,idx3)+D31(j,:);
    end
end

% Component 4
% derivatives of F(YU) w.r.t  x1 

hess=HomT_hess(x1,p,J);
l=n*(N-1); Q0=homTds.Q0;

for i=1:n
    D=Q0'*hess(:,:,i)*Q0;
    D1=D(1:homTds.nu,1:homTds.nu);
    D2=D(1:homTds.nu,homTds.nu+1:n);
    D3=D(homTds.nu+1:n,1:homTds.nu);
    D4=D(homTds.nu+1:n,homTds.nu+1:n);
    for j=1:n-homTds.nu
        for s=1:homTds.nu;
            idx=l+s+(j-1)*homTds.nu;
            result(idx,i)=result(idx,i)+D4(j,:)*homTds.YU(:,s)-homTds.YU(j,:)*D1(:,s)+D3(j,s);
            for k=1:homTds.nu
                result(idx,i)=result(idx,i)-homTds.YU(j,k)*(D2(k,:)*YU(:,s));
            end
        end           
    end
end
%result(21,1),result(21,2),pause

%%%%%%%%%%%%%%%%%%%%%%%%%
% Component  5% STABLE
% Ricatti blocks from stable eigenspace
% F(Y_S)=R22Y_S-Y_SR11+E21-YR12Y_S;
% =========
% D1=R_22Y_S
Q1S = homTds.Q1;
[R11, R12, E21, R22] = homT_RicattiCoeff(Q1S,A1,homTds.ns);
% derivatives of D1 = R22 * Y_S w.r.t (YS)ij
l=n*(N-1)+(n-homTds.nu)*homTds.nu; h=N*n+(n-homTds.nu)*homTds.nu;
for j=1:n-homTds.ns
    for i=1:homTds.ns
        idx1=l+i+(j-1)*homTds.ns;
        idx2=h+1+(i-1)*(n-homTds.ns);
        idx3=h+(n-homTds.ns)+(i-1)*(n-homTds.ns);
        result(idx1,idx2:idx3)=result(idx1,idx2:idx3)+R22(j,:);
    end
end
% derivatives of D2= Y_S * R11 w.r.t (YS)ij

for j=1:n-homTds.ns
    for s=1:homTds.ns
        idx1=l+i+(j-1)*homTds.ns;
        idx2=h+homTds.nu+(j-1)*homTds.ns;
        idx3=h+j+(s-1)*(n-homTds.ns);
        result(idx1,idx2:idx3)=result(idx1,idx2:idx3)-(R11(s,:))';
    end
end
% derivatives of D3 = YS * R12*YS w.r.t (YS)ij
D31=-homTds.YS*R12;
for j=1:n-homTds.ns
    for i=1:homTds.ns
        idx1=l+i+(j-1)*homTds.ns;
        idx2=h+1+(i-1)*(n-homTds.ns);
        idx3=h+n-homTds.ns+(i-1)*(n-homTds.ns);
        result(idx1,idx2:idx3)=result(idx1,idx2:idx3)+D31(j,:); 
    end
    
end

D32=-R12*homTds.YS;
for j=1:n-homTds.ns
    for s=1:homTds.ns
        idx1=l+1+(j-1)*homTds.ns;
        idx2=l+homTds.ns+(j-1)*homTds.ns;
        idx3=h+j+(s-1)*(n-homTds.ns);
        result(idx1:idx2,idx3)=result(idx1:idx2,idx3)+(D32(s,:))';
        
    end
end

%result(22,24),pause
%%%%%%%%%%%%%%%%%%%%%%%%


% Component 7
% derivatives of F(Y_U)=R22Y_U-Y_UR11+E21-Y_UR12Y_U; w.r.t  the active parameter 

% hessp=Hom_hessp(x1,p,J);
% l=n*(N-1); Q0=homTds.Q0;
% D=Q0'*hessp*Q0;
% D1=D(1:homTds.nu,1:homTds.nu);
% D2=D(1:homTds.nu,1:homTds.nu+1:n);
% D3=D(homTds.nu+1:n,1:homTds.nu);
% D4=D(homTds.nu+1:n,homTds.nu+1:n);
% for j=1:n-homTds.nu
%     for s=1:homTds.nu;
%         idx=l+s+(j-1)*homTds.nu;
%         result(idx,ks+1)=result(idx,ks+1)+D4(j,:)*homTds.YU(:,s)-homTds.YU(j,:)*D1(:,s)+D3(j,s);
%         for k=1:homTds.nu
%             result(idx,ks+1)=result(idx,ks+1)-homTds.YU(j,k)*(D2(k,:)*YU(:,s));          
%         end           
%     end
% end

% Component 7
% derivatives of F(Y_S) w.r.t  x1
hess=HomT_hess(x1,p,J);
l=n*(N-1)+(n-homTds.nu)*homTds.nu;
Q1=homTds.Q1;

for i=1:n
    D=Q1'*hess(:,:,i)*Q1;
    D1=D(1:homTds.nu,1:homTds.nu);
    D2=D(1:homTds.nu,homTds.nu+1:n);
    D3=D(homTds.nu+1:n,1:homTds.nu);
    D4=D(homTds.nu+1:n,homTds.nu+1:n);
    for j=1:n-homTds.nu
        for s=1:homTds.nu;
            idx=l+s+(j-1)*homTds.nu;
            result(idx,i)=result(idx,i)+D4(j,:)*homTds.YU(:,s)-homTds.YU(j,:)*D1(:,s)+D3(j,s);
            for k=1:homTds.nu
                result(idx,i)=result(idx,i)-homTds.YU(j,k)*(D2(k,:)*YU(:,s));
            end
        end           
    end
end

%result(22,1),result(22,2),pause

% Derivatives of F(Y_S)=R22Y_S-Y_SR11+E21-Y_SR12Y_S w.r.t the active parameter 

% hessp=Hom_hessp(x1,p,J); 
% l=n*(N-1)+homTds.nu*(n-homTds.nu);
% Q1=homTds.Q1;
% D=Q1'*hessp*Q1;
% D1=D(1:homTds.ns,1:homTds.ns);
% D2=D(1:homTds.ns,1:homTds.ns+1:n);
% D3=D(homTds.ns+1:n,1:homTds.ns);
% D4=D(homTds.ns+1:n,homTds.ns+1:n);
% for j=1:n-homTds.ns
%     for s=1:homTds.ns;
%         idx=l+s+(j-1)*homTds.ns;
%         result(idx,ks+1)=D4(j,:)*homTds.YS(:,s)-homTds.YS(j,:)*D1(:,s)+D3(j,s);
%         for k=1:homTds.nu
%             result(idx,ks+1)=result(idx,ks+1)-homTds.YS(j,k)*(D2(k,:)*YS(:,s));    
%         end           
%     end
% end
%result,pause


% Component  7 (First vector along unstable eigenspaces)
% ===========
% derivatives w.r.t x1 
    Q0U = homTds.Q0;
    vect = ups(:,2) - x1;
    QU =Q0U*[-YU'; eye(size(YU,1))];
    l=n*(N-1)+(n-homTds.nu)*homTds.nu+(n-homTds.ns)*homTds.ns;
    for i=1:n-homTds.nu
       result(l+i, 1:n )=-QU(:,i)';
    end
   %result(23,1:2),pause
    
    
   % derivatives w.r.t x2
    vect = ups(:,2) - x1;
    l=n*(N-1)+(n-homTds.nu)*homTds.nu+(n-homTds.ns)*homTds.ns;
    for i=1:n-homTds.nu
       result(l+i, n+1:2*n )=QU(:,i)';
    end
   
   % result(23,3:4),pause
    
 % derivatives w.r.t  components of YU_{(nu+i)}   
 H=vect'*Q0U;h=n*N;h=n*N;H=H(homTds.nu); 
 
 for i=1:1:n-homTds.nu
     idx1=h+1+(n-homTds.nu)*(i-1);
     idx2=h+homTds.nu+(n-homTds.nu)*(i-1);
     result(l+i,idx1:idx2)=result(l+i,idx1:idx2)-H;
 end
%result(23,23),pause

 % Component  8 (Last vectors along stable eigenspaces)
% ===========
% derivatives w.r.t xN

    Q1S = homTds.Q1;
    %vect = ups(:,N-1) - x1;
    vect = ups(:,N) - x1;
    QS =Q1S*[-YS'; eye(size(YS,1))];
    l=n*(N-1)+(n-homTds.nu)*homTds.nu+(n-homTds.ns)*homTds.ns+homTds.ns;
    for i=1:n-homTds.ns
       result(l+i, n*(N-1)+1:n*N)=result(l+i, n*(N-1)+1:n*N)+QS(:,i)'; 
    end
   %result(24,21:22),pause
    
   % derivatives w.r.t x1   
    l=n*(N-1)+(n-homTds.nu)*homTds.nu+(n-homTds.ns)*homTds.ns+homTds.ns;
    for i=1:n-homTds.nu
        
       result(l+i, 1:n )=-QS(:,i)'; 
    end
    % result(24,1:2),pause
    
    % derivatives w.r.t  components of YU_{(nu+i)}
     vect = ups(:,N) - x1;
     H=vect'*Q1S;h=n*N+(n-homTds.nu)*homTds.nu;
     for i=1:1:n-homTds.ns
      for j=1:homTds.ns
        result(l+i,h+j+(n-homTds.ns)*(i-1))=-H(j);
      end
     end
   

88,size(result),pause