Документ взят из кэша поисковой машины. Адрес оригинального документа : http://lnfm1.sai.msu.ru/~rastor/Software/lafkin.m Дата изменения: Thu Mar 7 13:58:07 2013 Дата индексирования: Fri Feb 28 00:30:09 2014 Кодировка:

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));