function [P,T0]=lafkin(file,ncol,pmin,pmax,ncycl,Ydir);
%
%
%
% Function call: [P,T0]=lafkin(file,ncol,pmin,pmax,ncycl,Ydir);
%                                    OR
%                   lafkin(file,ncol,pmin,pmax,ncycl,Ydir);
%
%  Datafile structure:
%  1st column            - JD 
%  2 ,3 ,4, ..., columns - data to process (light, color, rad. vel. curves)
%
% Input:  file      - filename (write as 'V001.DAT')
%         ncol      - column number (with data)
%         pmin      - minimal search period
%         pmax      - maximal search period
%         ncycl     - the number of probes (constant step by frequency)
%         Ydir      - graph mode ('normal' or 'revers' Y-axis) for a phase curve 
%
% Output: P         - period
%         T0        - moment of Max ('normal' Ydir)/ Min ('revers' Ydir)
%
%
% Examples:           [P,T0]=lafkin(''V001.dat,2,0.3,3.0,100000,'revers');
%                            lafkin(''V001.dat,2,0.3,3.0,100000,'revers');
%  
%
%
%% Author: A.S.Rastorguev, MSU faculty of physics, Sternberg astronomical Institute, 2013
%%


%% Program body
clf
a=load(file,'-ascii');
%
[nr,nc]=size(a);
%
% Step size (in cycles):
dfr=(1/pmin-1/pmax)/(ncycl-1);
% 
fr=1/pmax-dfr;
%
for i=1:ncycl
% Probing frequency:    
    fr=fr+dfr;
% Corresponding period:
    Per(i,1)=1/fr;
% Calculating phases:    
    for j=1:nr
        ph(j,1)=a(j,1)*fr;
        ph(j,1)=ph(j,1)-floor(ph(j,1));
    end
% Ranging by phases in ascending mode:
    [y,I]=sort(ph(:,1),'ascend');
% Now y(:,1:2) - data from source file, sorted by phase     
    y(:,2)=a(I(:,1),ncol);
% Calculating Lafler-Kinman statistics:
    LK(i,1)=0;
    for j=1:(nr-1)
        LK(i,1)=LK(i,1)+(y(j+1,2)-y(j,2))^2;
    end
end
%% LK statistics maximun:
[Y,I] = min(LK(:,1));
P=Per(I);
%
%% Searcing T0 as Max (or Min) value moment
if Ydir=='revers',
    [yy,ii]=min(a(:,ncol));
else
    [yy,ii]=max(a(:,ncol));
end
T0=a(ii,1);
%
%% Building phase curve for "best" period  
for i=1:nr
    ph(i,1)=(a(i,1)-T0)/P;
    ph(i,1)=ph(i,1)-floor(ph(i,1));
end
%% Plotting...
% ... LK spectrum (1 / LK as power spectrum indice):
subplot(211);
plot(Per,1./LK);grid
%plot(Per,-LK);grid
xlabel('Period');ylabel('1 / LK');
title('Power spectrum');
%
% ... phase curve:
subplot(212);
if Ydir=='revers',
    plot(ph,a(:,ncol),'+');grid
    set(gca,'XDir','normal','YDir','reverse');
else
    plot(ph,a(:,ncol),'+');grid
end
xlabel('Phase');ylabel('Value');
title('Phase curve');
%
%% Screen output:
disp('                          ')
disp('-------------------------')
disp('     Period       T0     ')
disp('-------------------------')
disp(sprintf('  %10.6f   %10.4f',P,T0));