% Calculation the mean LOD of 365 solar days from aphelions for a period
      % расчет и вывод в книгу Excel средних длин 365 солнечных суток с момента афелия
      % и девиаций между соседними периодами

      clear all
      cspice_furnsh( 'd:\Documents and Settings\user\Мои документы\MATLAB\astro\mice\data\standard.tm' );      
      TFMT='DD.MM.YYYY HR:MN:SC ERA::RND::MCAL::TDT'; TFMTH='HR:MN:SC.## ::RND::TDT' ; 
      % the path to output - путь вывода
      save_path = 'd:\Documents and Settings\user\Мои документы\MATLAB\astro\mice\test\NC\';
      
           
      % the epochs of aphelions for period
      % моменты афелиев Земли за период
      et1=cspice_str2et ( '01-07-0000 AD'  ); et2=cspice_str2et ( '2013-07-01' ); 
      st_win = cspice_wninsd (et1 , et2);
      target  = 'EARTH'; abcorr  = 'NONE'; observer  = 'SUN'; relate  = 'LOCMAX';
      refval  = 0; adjust  = 0.; step    = 0.45 * cspice_tyear; nintvls = 1000000;
 
      et_aphelions = cspice_gfdist( target, abcorr, observer, relate, refval, adjust, step,  nintvls, st_win);
      et_aphelions=(reshape(et_aphelions,2,[])'); % преобразование окон времени в 2 колонки 1) начало события 2) конец события
      et_aphelions= et_aphelions (:,1); % выбор одной колонки, там где колонки равны
                  
                          
      % 365 equal azimuths of Sun from the epoch of an given aphelion
      % 365 одинковых азимутов Солнца c момента афелия 
      target = '399';
      [observer, topo_frame ]=set_observer('CAIRO');  % широта, долгота, высота - LAN, LON, ALT
      n=numel(et_aphelions); az=zeros(365,1); mean_365_days = zeros(n,1); a365th_day=zeros(n,1);
      
      for ii=1:n
%       fprintf ( '\nan aphelion, %s\n', cspice_timout( et_aphelions(ii), 'DD.MM.YYYY ERA JCAL HR:MN ::RND::JCAL::UTC UTC')) ;
      % determine an azimuth of the sun on aphelions epoch    
      % определяем азимут солнца на момент афелия
      [pos, ~] = cspice_spkpos ( target, et_aphelions(ii), topo_frame, abcorr, observer); 
      [~, az(ii), ~] = cspice_recrad(pos);
      
      % the ~367 days period from an aphelion (2 days will get off)
      st_win(1,:) = et_aphelions(ii) - 0.48 * cspice_spd + cspice_tyear;
      st_win(2,:) = et_aphelions(ii) + 0.48 * cspice_spd + cspice_tyear; 
            
      % the moments of the repetitive azimuths of the sun for the next 365 days from an aphelion
      % ищем моменты повторений того же азимута солнца начиная с момента данного афелия
      
      target = 'Sun'; crdsys = 'LATITUDINAL';coord = 'LONGITUDE'; 
      relate='='; refval = az(ii); adjust=0; step = 0.2 * cspice_spd;
      res = cspice_gfposc( target, topo_frame, abcorr, observer, crdsys, coord, relate, refval, ...
                                                              adjust, step,  nintvls, st_win );
                                                          
      a365th_day(ii)= res(1); 
      a=et_aphelions(ii)-res(1);
      b=1;
        
%        n=numel(res);                                                    
%        for kk=1:2:n-2 
%            res(kk+1)=res(kk+2)-res(kk);
%        end
%        res=(reshape(res,2,[])'); % преобразование окон времени в 2 колонки 1) начало события 2) конец события
%        res= res (1:365,2); % выбор одной колонки, там где колонки равны
        
      end
mean_365_days= (a365th_day-et_aphelions)/365; 
mean_365_days= mean_365_days(2:numel(mean_365_days)-1,1);
% plot(mean_365_days);     
a=(1:numel(mean_365_days));
cftool(a,mean_365_days);
title('the mean LOD of 365 solar days from aphelions');  
      
% xlswrite(horzcat(save_path, 'earth_lst_deviations.xls'), data);
% fprintf(horzcat('the result saved in file ', 'earth_lst_deviations.xls' ));

      cspice_kclear  
      clear all