% Calculation the mean LOD of 365 solar days from aphelions NOON 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-02-3000 BC'  ); et2=cspice_str2et ( '12-12-2000' );  
      st_win = cspice_wninsd (et1 , et2);
      target  = 'EARTH'; abcorr  = 'NONE'; observer  = 'SUN'; relate  = 'LOCMAX';
      refval  = 0; adjust  = 0.; step    = 0.2 * 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(:,2)= et_aphelions (:,2) + cspice_spd; % выбор одной колонки, там где колонки равны
      et_aphelions=reshape(et_aphelions',1,[])';
                  
      % the moments of the azimuths of the sun = 0
      % ищем моменты азимута солнца = 0 т афелия (полночь)
      refval = 180;
      [observer, topo_frame ]=set_observer('CAIRO');  % широта, долгота, высота - LAN, LON, ALT
      [lat, lon, alt, sslat, sslon] = get_observer (observer);
      target = 'Sun'; crdsys = 'LATITUDINAL';coord = 'LONGITUDE'; 
      relate='='; adjust=0; step = 0.22 * cspice_spd;
      et_noon = cspice_gfposc( target, topo_frame, abcorr, observer, crdsys, coord, relate, refval, ...
                                                              adjust, step,  nintvls, et_aphelions );
      et_noon=(reshape(et_noon,2,[])'); % преобразование окон времени в 2 колонки 1) начало события 2) конец события
      et_noon= et_noon (:,1); % выбор одной колонки, там где колонки равны
                                                         
      n=numel(et_noon); 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_noon(ii), 'DD.MM.YYYY ERA JCAL HR:MN ::RND::JCAL::UTC UTC')) ;
%       [~, ~, ~, time, ~] = cspice_et2lst( et_noon(ii), 399, lon * cspice_rpd, 'PLANETOGRAPHIC');
%       fprintf ( 'LST   : %s \n', time );
%       
      % the ~367 days period from an aphelion (2 days will get off)
      st_win(1,:) = et_noon(ii) - 0.48 * cspice_spd + cspice_tyear;
      st_win(2,:) = et_noon(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='='; adjust=0; step = 0.22 * cspice_spd;
      res = cspice_gfposc( target, topo_frame, abcorr, observer, crdsys, coord, relate, refval, ...
                                                              adjust, step,  nintvls, st_win );
                                                          
      a365th_day(ii)= res(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_noon)/365; 
a=1;
% 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