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

function result = BVP_Hom_jac(mapfile,x,p,YS,YU,J)

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

A1=homjac(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)=homjac(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(:,homds.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 = homds.Q0;
[R11, R12, E21, R22] = hom_RicattiCoeff(Q0U,A1,homds.nu);

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


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

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


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

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

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

for i=1:n
    D=Q0'*hess(:,:,i)*Q0;
    D1=D(1:homds.nu,1:homds.nu);
    D2=D(1:homds.nu,homds.nu+1:n);
    D3=D(homds.nu+1:n,1:homds.nu);
    D4=D(homds.nu+1:n,homds.nu+1:n);
    for j=1:n-homds.nu
        for s=1:homds.nu;
            idx=l+s+(j-1)*homds.nu;
            result(idx,i)=result(idx,i)+D4(j,:)*homds.YU(:,s)-homds.YU(j,:)*D1(:,s)+D3(j,s);
            for k=1:homds.nu
                result(idx,i)=result(idx,i)-homds.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 = homds.Q1;
[R11, R12, E21, R22] = hom_RicattiCoeff(Q1S,A1,homds.ns);
% derivatives of D1 = R22 * Y_S w.r.t (YS)ij
l=n*(N-1)+(n-homds.nu)*homds.nu; h=N*n+(n-homds.nu)*homds.nu;
for j=1:n-homds.ns
    for i=1:homds.ns
        idx1=l+i+(j-1)*homds.ns;
        idx2=h+1+(i-1)*(n-homds.ns);
        idx3=h+(n-homds.ns)+(i-1)*(n-homds.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-homds.ns
    for s=1:homds.ns
        idx1=l+i+(j-1)*homds.ns;
        idx2=h+homds.nu+(j-1)*homds.ns;
        idx3=h+j+(s-1)*(n-homds.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=-homds.YS*R12;
for j=1:n-homds.ns
    for i=1:homds.ns
        idx1=l+i+(j-1)*homds.ns;
        idx2=h+1+(i-1)*(n-homds.ns);
        idx3=h+n-homds.ns+(i-1)*(n-homds.ns);
        result(idx1,idx2:idx3)=result(idx1,idx2:idx3)+D31(j,:); 
    end
    
end

D32=-R12*homds.YS;
for j=1:n-homds.ns
    for s=1:homds.ns
        idx1=l+1+(j-1)*homds.ns;
        idx2=l+homds.ns+(j-1)*homds.ns;
        idx3=h+j+(s-1)*(n-homds.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=homds.Q0;
D=Q0'*hessp*Q0;
D1=D(1:homds.nu,1:homds.nu);
D2=D(1:homds.nu,1:homds.nu+1:n);
D3=D(homds.nu+1:n,1:homds.nu);
D4=D(homds.nu+1:n,homds.nu+1:n);
for j=1:n-homds.nu
    for s=1:homds.nu;
        idx=l+s+(j-1)*homds.nu;
        result(idx,ks+1)=result(idx,ks+1)+D4(j,:)*homds.YU(:,s)-homds.YU(j,:)*D1(:,s)+D3(j,s);
        for k=1:homds.nu
            result(idx,ks+1)=result(idx,ks+1)-homds.YU(j,k)*(D2(k,:)*YU(:,s));          
        end           
    end
end

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

for i=1:n
    D=Q1'*hess(:,:,i)*Q1;
    D1=D(1:homds.nu,1:homds.nu);
    D2=D(1:homds.nu,homds.nu+1:n);
    D3=D(homds.nu+1:n,1:homds.nu);
    D4=D(homds.nu+1:n,homds.nu+1:n);
    for j=1:n-homds.nu
        for s=1:homds.nu;
            idx=l+s+(j-1)*homds.nu;
            result(idx,i)=result(idx,i)+D4(j,:)*homds.YU(:,s)-homds.YU(j,:)*D1(:,s)+D3(j,s);
            for k=1:homds.nu
                result(idx,i)=result(idx,i)-homds.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)+homds.nu*(n-homds.nu);
Q1=homds.Q1;
D=Q1'*hessp*Q1;
D1=D(1:homds.ns,1:homds.ns);
D2=D(1:homds.ns,1:homds.ns+1:n);
D3=D(homds.ns+1:n,1:homds.ns);
D4=D(homds.ns+1:n,homds.ns+1:n);
for j=1:n-homds.ns
    for s=1:homds.ns;
        idx=l+s+(j-1)*homds.ns;
        result(idx,ks+1)=D4(j,:)*homds.YS(:,s)-homds.YS(j,:)*D1(:,s)+D3(j,s);
        for k=1:homds.nu
            result(idx,ks+1)=result(idx,ks+1)-homds.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 = homds.Q0;
    vect = ups(:,2) - x1;
    QU =Q0U*[-YU'; eye(size(YU,1))];
    l=n*(N-1)+(n-homds.nu)*homds.nu+(n-homds.ns)*homds.ns;
    for i=1:n-homds.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-homds.nu)*homds.nu+(n-homds.ns)*homds.ns;
    for i=1:n-homds.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(homds.nu); 
 
 for i=1:1:n-homds.nu
     idx1=h+1+(n-homds.nu)*(i-1);
     idx2=h+homds.nu+(n-homds.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-1

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