%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%      Compute all reachable / observable spaces
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [R, U] = reachUnobserv(A,B,C,J)

p=size(A,2); 
%%       reachable subspaces
R{1} = sym(myinvariant(A{1}, B{1}));
for k=2:p   
    local_R{k} = sym(myinvariant(A{k}, B{k})); 
    R_invariance{k} = sym(myinvariant(A{k}, J{k-1}*R{k-1}));  
    R{k} = mycolspace(sym([local_R{k},R_invariance{k}])); 
end

 %%    unobservable subspaces
O{p}= sym(myinvariant(A{p}', C{p}'));
U{p} = sym(mynull(O{p}')); 

 for k=p-1:-1:1           
    O{k}= sym(myinvariant(A{k}', C{k}')); 
    local_U{k} = sym(mynull(O{k}'));  
   
    part2{k}=getPreImage(J{k}, U{k+1});
    N{k} = myintersect(local_U{k}, part2{k});  % unobsv
    M{k} = mynull(N{k}');      %  obsrv 
    U2_invariance{k} =sym(myinvariant(A{k},M{k})); 
    U{k} = sym(mynull(U2_invariance{k}'));            % unobsv
  end
%  print("done")
end




