function [T Y PRC]=prc_fm(sys_dinam,jakobianas,method,per_points,prad_sal,sprend_tiksl)
% v.004.0_en_1
% Function prc_fm() returns 3 vector columns. In the first column T are 
% saved time values where PRC is calculated, in the second column Y are saved 
% limit cycle values where PRC is calculated and third column has PRC values.
% method can be 'simple' or 'advanced'.
% per_points is vector where [precise_period interval_number].
% prad_sal is initial conditions on limit cycle.
% sprend_tiksl is solution accuracy given in the same way as in function
% ode45.


N=size(prad_sal,2);
period=per_points(1);
if strcmpi(method,'advanced') dep_lev=floor(log(per_points(2))/log(2)); n=2^dep_lev; end;
if strcmpi(method,'simple') n=per_points(2); end;
sk=period/n;
mt=eye(N);
prad_sal1=mt(:)';
aux=zeros(N,N,n);
if strcmpi(method,'advanced') aux_level=zeros(N,N,n,dep_lev+1); end;
greit=zeros(N,n);
evol=zeros(N,N,n);
T=zeros(n+1,1);
Y=zeros(n+1,N);
PRC=zeros(n+1,N);


% Here we construct function which returns dinamical system and 
% linearised dinamical system 
    function dydt=sys_dinam_linear(t,y)      
        dydt(1:N,1)=sys_dinam(t,y(1:N));  
            for m=1:N             
        dydt(m*N+1:(m+1)*N,1) = jakobianas(y(1:N))*y(m*N+1:(m+1)*N);
            end
    end
% From this loop we get velocity,limit cycle and time    

for i=1:n
greit(:,i)=sys_dinam(0,prad_sal);
Y(i,:)=prad_sal;
sol = ode45(@sys_dinam_linear,[0 sk],[prad_sal prad_sal1],sprend_tiksl);
prad_sal=(sol.y(1:N,end))';
aux(:,:,i)=reshape((sol.y(N+1:end,end))',N,N);
T(i,:)=sk*(i-1);
end;
Y(n+1,:)=Y(1,:);
T(n+1,:)=period;



% Here we multiply auxiliary matrices in special order

if strcmpi(method,'advanced')
aux_level(:,:,:,1)=aux;
for lev=2:dep_lev+1
    for k=1:2^(dep_lev-lev+1)
        aux_level(:,:,k,lev)=aux_level(:,:,2*k,lev-1)*aux_level(:,:,2*k-1,lev-1);
    end;
end;

for i=1:n
j=i-1;
PART1=eye(N);
while j>0
    bi=fliplr(bitget(uint32(j),dep_lev:-1:1));
    for pp=1:dep_lev
        if bi(pp)==1 PART1=PART1*aux_level(:,:,j/2^(pp-1),pp); j=j-2^(pp-1); break; end;
    end;
end;

j=i-1;
if j==0 PART2=aux_level(:,:,1,dep_lev+1); else PART2=eye(N); end;
while j>0
    bi=fliplr(bitget(uint32(j),dep_lev:-1:1));
    for pp=1:dep_lev
        if bi(pp)==1 PART2=aux_level(:,:,j/2^(pp-1)+1,pp)*PART2; j=mod(j+2^(pp-1),n); break; end;
    end;
end;
evol(:,:,i)=PART1*PART2;
end;

end;

if strcmpi(method,'simple')
for i=1:n
    
    az=eye(N);
    for j=1:n
    k=i-j;
    if k<1 
        k=k+n;
    end;
    az=az*aux(:,:,k);
    end;
    evol(:,:,i)=az;  
end;
clear az;

end;



D=NaN(1,N);
% Here we calculate left eigenvector and PRC
for i=1:n
  for r=1:N
      eilute=evol(:,:,i);
      eilute(r,:)=[];
      eilute(:,r)=[];
      D(r)=det(eilute-eye(N-1));
  end;
[~,nr]=max(abs(D));
eilute=evol(:,:,i);
eilute(:,nr)=[];
up=eilute(nr,:);
eilute(nr,:)=[];
l2=-up/(eilute-eye(N-1));   
L1=[l2(1:nr-1) 1 l2(nr:end)];
PRC(i,:)=L1/(L1*greit(:,i));    
end; 
PRC(n+1,:)=PRC(1,:);
end

