function sol2knut(sol,pars,period,ntst,ncol,output)
% This m-function translates a periodic solution as simulated by dde23
% within MATLAB into a starting file suitable for PDDE-CONT. 
% The input parameters are
% 1. The periodic solution
% 2. The parameters from the simulation
% 3. An guess for the period
% 4. The number of mesh points for the solution
% 5. The number of collocation points
% 6. The name of the output mat-file contains several matrices described in the
% PDDE-CONT/KNUT manual:
% Matrix name   Rows Columns    Content
% knut_npoints  1       1       Number of points computed in the branch of solutions.
% knut_par      NPAR+3 STEPS    Values of parameters.
% knut_mul      NMUL   STEPS    The Floquet multipliers.
% knut_ntrivmul 1       1       Number of trivial multipliers.
% knut_ndim     1      STEPS    NDIM: dimesion of the right-hand-side% function f
% knut_elem     NDEG+1 STEPS    Mesh points in a collocation interval.
% knut_mesh     NINT+1 STEPS    Mesh points of the solution profile.
% knut_prof     NDIM*(NINT*NDEG+1) STEPS The solution profile.
% Of course, we might use STEPS=1, but the software requires the
% zero-starting point to be specified in the cpp-file and all others with
% higher labels.

knut_npoints=2;
knut_ntrivmul = [1;0;0];
knut_mul=complex(zeros(3,knut_npoints,'double'));
knut_ndim=size(sol.y,1);
knut_elem=[(0:ncol); (0:ncol)]'/ncol;
knut_mesh=[(0:ntst); (0:ntst)]'/ntst;

% Finding an 'optimal' period
times=find(abs(sol.x-sol.x(end)+period)<.1*period);
dist=zeros(size(times,2),1);
for i=1:size(times,2)
  dist(i)=norm(sol.y(:,times(i))-sol.y(:,end));
end
[Y,i]=min(dist);starttime=sol.x(times(i));
period=sol.x(end)-sol.x(times(i));
knut_par(:,1)=[period,pars,0,0,0];
knut_par(:,2)=[period,pars,0,0,0];

%Use linear approximation to construct a starting profile on an equidistant mesh.
knut_prof=zeros(knut_ndim*(ntst*ncol+1),2);
knut_prof(1:knut_ndim,1)=sol.y(:,end);
knut_prof(end-knut_ndim+1:end,1)=sol.y(:,end);
for i=1:ntst
  for j=1:ncol
    t3=starttime+(knut_mesh(i)+knut_elem(j)/ntst)*period;
    t1=find(sol.x-t3<0,1,'last');
    t2=find(sol.x-t3>=0,1,'first');
    knut_prof(((i-1)*ncol+j)*knut_ndim+[1:knut_ndim],1)=...
      sol.y(:,t1)+(sol.y(:,t2)-sol.y(:,t1))*(t3-sol.x(t1))/(sol.x(t2)-sol.x(t1));
  end
end
%On this first mesh we want to minimize interpolation error, this gives the
%second point and the initial cycle for the continuer.
xx=reshape(knut_prof(:,1),knut_ndim,ntst*ncol+1);xx=xx(:,1:ncol:end);
cumulerror=[0,cumsum(sqrt(sum(diff(xx,1,2).^2,1)))];
idealcumulerror = (0:ntst)*cumulerror(end)/ntst;

% idealcumulerror contains the ideal cumulative error to have at each mesh
% point (over each mesh interval). We go over each of these points.
for j=1:length(idealcumulerror)
    % At each point, find the index of the cumulative error-vector of the
    % old situation, where you just passed the ideal cumulative error...
    m = find(cumulerror > idealcumulerror(j),1,'first');
    if isempty(m)
        m = length(cumulerror);
    end
    % ...and go subtract 1 to get the index right before reaching the ideal
    % cumulative error
    m = m-1;%if(m==0) m=1;end
    
    % Now we will find the new mesh point by interpolating between two old
    % mesh points: the one right before, and the one right after reaching
    % the ideal cumulative error.
    t = (idealcumulerror(j) - cumulerror(m)) / (cumulerror(m+1) - cumulerror(m));
    knut_mesh(j,2) = (1-t)*knut_mesh(m,1) + t*knut_mesh(m+1,1);
end
for i=1:ntst
  for j=1:ncol
    t3=starttime+(knut_mesh(i,2)+knut_elem(j+1)*(knut_mesh(i+1,2)-knut_mesh(i,2)))*period;
    t1=find(sol.x-t3<0,1,'last');
    t2=find(sol.x-t3>=0,1,'first');
    knut_prof(((i-1)*ncol+j)*knut_ndim+[1:knut_ndim],2)=...
      sol.y(:,t1)+(sol.y(:,t2)-sol.y(:,t1))*(t3-sol.x(t1))/(sol.x(t2)-sol.x(t1));
  end
end
knut_prof(1:knut_ndim,2)=sol.y(:,end);
knut_prof(end-knut_ndim+1:end,2)=sol.y(:,end);
knut_ndim=[size(sol.y,1) size(sol.y,1)];

% knut_prof(:,2)=knut_prof(:,1);
% knut_mesh(:,2)=knut_mesh(:,1);
knut_parnames=zeros(size(knut_par,1),13);

save(output, 'knut_npoints', 'knut_ntrivmul', 'knut_elem', 'knut_ndim', ...
      'knut_par', 'knut_mul', 'knut_mesh', 'knut_prof', 'knut_parnames', '-V4')

end
