function [EpochMin,Period,Scatter]=FourierPeriodogram(fname,order,N1st,Pmin,Pmax,delta,epoch);
%
% -------------------------------------------------------------------------
% Program calculates Fourier power spectrum of the time series
%
% Program call:
% [EpochMin,Period,Scatter]=FourierPeriodogram(fname,order,N1st,Pmin,Pmax,delta,epoch)
% or
% FourierPeriodogram(fname,order,N1st,Pmin,Pmax,delta,epoch)
%
% Input:
%       fname           file name with the time series
%                       1st column - time
%                       2nd column - data
%       order           Fourier expansion order
%       N1st            the number of garmonics (from the 1st one)
%                       used to estimate the power spectrum by summing
%                       squared Fourier coefficients
%       Pmin,Pmax       minimal and maximal trial period values
%       delta           phase shift between consecutive trial phases
%       epoch           initial epoch (any)
%
% Output (screen and variables):
%       EpochMin        epoch of the minimum
%       Period          period
%       Scatter         scatter of the phase curve
% 
% Programmer: A.S.Rastorguev (Sternberg Astronomical Institute, Lomonosov
%                             Moscow State University, 7(495)9391616)
% -------------------------------------------------------------------------
close all
data=load(fname);
JD(:,1)=data(:,1);
mag(:,1)=data(:,2);
NN=order;
[nr,nc]=size(data);

% Начальная частота:
omega=1/Pmax;
% Начало счёта циклов:
M=0;
while omega <= 1/Pmin,
   M=M+1;
   om(M,1)=omega;
   for i=1:nr
       % Вычисление фаз для текущего пробного периода
       t(i,1)=omega*(JD(i)-epoch);
       t(i,1)=t(i,1)-floor(t(i,1));
   end
   %% Вычисление коэффициентов Фурье-разложения для текущего периода:
   %
   %  y =  x(1)+x(2)*sin (2pi*t*1) + ... + x(NN+1)*sin (2pi*t*N) +
   %       +  x(NN+2)*cos (2pi*t*1) + ...+x(2*NN+1)*cos (2pi*t*N)
   %
   clear A
   A = zeros(nr,2*NN+1);
   for i = 1:nr
       A(i,1) = 1;
   end
   for j = 2:(NN+1)
       A(:,j) = sin (2*pi*t.*(j-1));
   end
   for j = (NN+2):(2*NN+1)
       A(:,j) = cos (2*pi*t.*(j-NN-1));
   end
   % Определение массива коэффициентов:
   x=A\mag;
   % Сумма квадратов всех коэффициентов Фурье-разложения:
   sp_all(M,1)=0;
   for j=2:(2*NN+1)
       sp_all(M,1)=sp_all(M,1)+x(j)^2;
   end
   % Сумма квадратов амплитуд первых N1st гармоник
   sp1st(M,1)=0;
   for j=1:N1st
       sp1st(M,1)=sp1st(M,1)+x(j+1)^2+x(j+NN+1)^2;
   end
   % Восстановленное приближение:
   z=A*x;
   % Рассеяние данных относительно Фурье-разложения:
   ff(M,1)=std(mag-z);
   % Увеличиваем частоту на delta для следующего цикла:
   omega=omega+delta;
end % while
%
% Вывод периодограммы (спектра мощности)
% Размер дисплея:
screen_size=get(0,'MonitorPosition');
% Размер рисунка:
swidth1=screen_size(3)*5/6;
swidth2=screen_size(3)*2/3;
sheight=screen_size(4)*4/5;
h1=figure('Position',[0 screen_size(4)-sheight-100 swidth1 sheight],'NumberTitle','off','Name','Power spectrum, Scatter and 1/Scatter Graphs');
subplot(311)
plot(om,sp1st,'k.');grid;
ylabel('First garmonics');
subplot(312)
plot(om,ff,'k.');grid;
ylabel('Scatter');
subplot(313)
plot(om,1./ff,'k.');grid;
xlabel('Cycles per day');ylabel('1/Scatter');
%
% Определяем пик периодограммы:
[phase_max,I] = min(ff);
P=1/om(I);
% Минимальное рассеяние:
MinScatter=ff(I);
% Рисование фазовой кривой для наилучшего периода:
omega=om(I);
for i=1:nr
   % Вычисление фаз для найденного периода
   t(i,1)=omega*(JD(i)-epoch);
   t(i,1)=t(i,1)-floor(t(i,1));
end
%% Вычисление коэффициентов Фурье-разложения:
%
%  y =  x(1)+x(2)*sin (2pi*t*1) + ... + x(NN+1)*sin (2pi*t*N) +
%       +  x(NN+2)*cos (2pi*t*1) + ...+x(2*NN+1)*cos (2pi*t*N)
%
clear A
A = zeros(nr,2*NN+1);
for i = 1:nr
    A(i,1) = 1;
end
for j = 2:(NN+1)
    A(:,j) = sin (2*pi*t.*(j-1));
end
for j = (NN+2):(2*NN+1)
    A(:,j) = cos (2*pi*t.*(j-NN-1));
end
% Определение массива коэффициентов:
x=A\mag;
% Восстановленное приближение:
z=A*x;
%
% Определение момента минимума по сглаженной кривой (разложению):
[phase_max,J] = min(z);
epoch=JD(J);
% Пересчёт разложения с новой начальной эпохой (минимума):
for i=1:nr
   % Вычисление фаз для найденного периода
   t(i,1)=omega*(JD(i)-epoch);
   t(i,1)=t(i,1)-floor(t(i,1));
end
%% Вычисление коэффициентов Фурье-разложения:
%
%  y =  x(1)+x(2)*sin (2pi*t*1) + ... + x(NN+1)*sin (2pi*t*N) +
%       +  x(NN+2)*cos (2pi*t*1) + ...+x(2*NN+1)*cos (2pi*t*N)
%
clear A
A = zeros(nr,2*NN+1);
for i = 1:nr
    A(i,1) = 1;
end
for j = 2:(NN+1)
    A(:,j) = sin (2*pi*t.*(j-1));
end
for j = (NN+2):(2*NN+1)
    A(:,j) = cos (2*pi*t.*(j-NN-1));
end
% Определение массива коэффициентов:
x=A\mag;
% Восстановленное приближение:
z=A*x;
% Сортировка фаз в порядке возрастания:
[Y,I] = sort(t,'ascend');
%
h2=figure('Position',[0 screen_size(4)-sheight-100 swidth2 sheight],'NumberTitle','off','Name','Light curve');
plot(Y,z(I),'r',Y,mag(I),'k.');grid
xlabel('Phase');ylabel('Magnitude');
%%
disp('       To          Period       Scatter     ')
disp('--------------------------------------------')
disp(sprintf('  %10.6f   %10.6f   %10.6f',epoch,P,MinScatter));
%
% Вывод результатов - присвоение переменных:
EpochMin=epoch;
Period=P;
Scatter=MinScatter;
%
end % function