function [th,F,dI,dD,head,SQ_F,F_DI,SQ_DI]=didd_calc_fun(file,I0,spikes)
% Вычисление вариаций dI,dD по измерениям пяти модулей F
% для GEM dIdD-магнитометра
%
%   [th,F,dI,dD[,head,SQ_F,F_DI,SQ_DI]]=didd_calc_fun(file,I0[,spikes]);
%
% file  - имя файла с данными
%   I0  - среднее значение наклонения (необходимо при вычислении dD)
% spikes- амплитуда импульсной помехи для отбраковки (если 0 - нет
%           отбраковки, по умолчанию = 0)
%   th  - время измерении (в часах)
%   F   - модуль поля (без тока в кольцах)
% dI,dD - вариации наклонения и склонения (в град.)
%  head - заголовок
%  SQ_F - SQ-параметр для "чистого" модуля F
%  F_DI - модули с токами в кольцах (Dp, Dm, Ip, Im)
% SQ_DI - SQ-параметры для F с токами в кольцах
%
% Используется формулы (2-11) из статьи Alldredge L.R., Saldukas I.
% "The automatic standard magnetic observatory" in Technical Bull. - 1966.
% - No.31. - 35p. с учетом того, что колечная система не поворачивается на
% 180 град. и поправка зашита в софт. (?)

% 30.06.2012 - С.Ю.Хомутов, Новосибирск
% 24.11.2012 -  1) обработка для файла без заголовка в начале
%               2) выводяться также модули Dp,Dm,Ip,Im и SQ
% 26.07.2013 - предусмотрена отбраковка импульсных помех в модулях
% 20.04.2016 - наклонение для вычисления dD берется по I0 и вариациям
% 09.08.2017 - добавил чистку выбросов в F-канале

if ~exist(file),
    error(sprintf('!!! Не могу найти файл %s !!!\n',file))
end

if nargin>1,    % задано значение наклонения для вычисления dD
    I=I0;
    if nargin<2,% не задана величина спайка (импульсной помехи)
        spikes=0;
    end
end

% формат строки в файлах данных dIdD
dIdD_fmt='%2d%2d%f %f %d %f %d %f %d %f %d %f %d %f %f %f %f %f %f %f';

fid=fopen(file,'rt');
    str=char(fread(fid))';  % считываем в строку
fclose(fid);

head=[];th=[];F=[];Dp=[];Dm=[];Ip=[];Im=[];SQ_F=[];SQ_DI=[];
kh=findstr(str,'//');       % ищем строки заголовка
% проверяем, есть ли заголовок
if isempty(kh),  % нет строк-заголовков с "//"
    disp(sprintf('... Файл %s нестандартного формата, возможно, нет заголовка ...\n',file))
    head='No head';
    Dbase=0;
    % считываем весь файл, если получится
    C=textscan(str,dIdD_fmt);  % массив cell
    th=dms_deg(double(cell2mat(C(1))),double(cell2mat(C(2))),cell2mat(C(3)));
    F=cell2mat(C(4)); SQ_F=cell2mat(C(5));
    Dp=cell2mat(C(6)); Dm=cell2mat(C(8));
    Ip=cell2mat(C(10)); Im=cell2mat(C(12));
    SQ_DI=[cell2mat(C(7)) cell2mat(C(9)) cell2mat(C(11)) cell2mat(C(13))];
else % считываем заголовок (3 служебных строки)
    Nsect=length(kh)/3;     % кол-во секций с заголовками
    if rem(Nsect,1)>1E-5,   % кол-во строк заголовков не кратно 3
        error('!!! Ошибка в заголовках - строках с "//" !!!')
    end
    if Nsect>1,
        ks=reshape(kh,3,Nsect)';    % в виде мартицы - кол-во строк равно кол-ву секций
    else
        ks=kh(:)';
    end
    for k=1:Nsect,
        strhead=str(ks(k,1):ks(k,2)-2); % 1-я строка заголовка
        head=[head; cellstr(strhead)];
        strhead=str(ks(k,2):ks(k,3)-2); % 2-я строка заголовка
        head=[head; cellstr(strhead)];
        % находим базовое значение склонения
        kd=findstr(strhead,'Declination :');
        if ~isempty(kd),
            Dbase=str2num(strhead(kd+13:length(strhead)));
        else
            Dbase=0;
        end
        strhead=str(ks(k,3):ks(k,3)+137); % 3-я строка заголовка
        head=[head; cellstr(strhead)];
        if k<Nsect,
            k1=ks(k,3)+139; k2=ks(k+1)-1;
        else    % для последней секции - до конца файла
            k1=ks(k,3)+139; k2=length(str);
        end
        C=textscan(str(k1:k2),dIdD_fmt);  % массив cell
        % разбираем на числовые массивы
        th=[th;dms_deg(double(cell2mat(C(1))),double(cell2mat(C(2))),cell2mat(C(3)))];
        F=[F;cell2mat(C(4))]; SQ_F=[SQ_F;cell2mat(C(5))];
        Dp=[Dp;cell2mat(C(6))]; Dm=[Dm;cell2mat(C(8))];
        Ip=[Ip;cell2mat(C(10))]; Im=[Im;cell2mat(C(12))];
        SQ_DI=[SQ_DI; [cell2mat(C(7)) cell2mat(C(9)) cell2mat(C(11)) cell2mat(C(13))]];
    end
end

if spikes>0,    % чистим единичные выбросы (спайки)
    fprintf('F  - '); F=remove_spikes(F,spikes);
    fprintf('Dp - '); Dp=remove_spikes(Dp,spikes);
    fprintf('Dm - '); Dm=remove_spikes(Dm,spikes);
    fprintf('Ip - '); Ip=remove_spikes(Ip,spikes);
    fprintf('Im - '); Im=remove_spikes(Im,spikes);
    fprintf('\n');
end
% 
F_DI=[Dp Dm Ip Im];

% вычисление вариации dI
A=sqrt((Ip.^2+Im.^2)/2-F.^2);
dI=r2d((Ip.^2-Im.^2)./(4*F.*A));

I=I0+dI; % 20.04.2016

% вычисление вариации dD
A=sqrt((Dp.^2+Dm.^2)/2-F.^2);
dD=r2d(((Dp.^2-Dm.^2)./(4*F.*A)./cos(d2r(I))));
% dD=-r2d( (asin((Dp.^2-Dm.^2)./(4*F.*A)/cos(d2r(I)))) );
