//  Copyright (c) CNES  2008
//
//  This software is part of CelestLab, a CNES toolbox for Scilab
//
//  This software is governed by the CeCILL  license under French law and
//  abiding by the rules of distribution of free software.  You can  use,
//  modify and/ or redistribute the software under the terms of the CeCILL
//  license as circulated by CEA, CNRS and INRIA at the following URL
//  'http://www.cecill.info'.

// ----------------------------------------------------------- 
// Equation of time
//> Equation of time
//> Represents the time difference: (E) = UT minus 'Greenwich solar time' in minutes. 
//> So when it is noon (Greenwich solar time), it is noon+(E) UT. 
// Auteur : A. Lamy
// -----------------------------------------------------------


// eqt = mean local time minus true local time (hours)
function [eqt] = equatemps(an, trel)
  t0 = CL_dat_cal2cjd(an,01,01,0,0,0);
  tab_t = t0 + trel; 
  hlocmoy = 12*ones(tab_t); // mean loc time = 12h
  hlocvraie = CL_op_locTime(tab_t ,'mlh',hlocmoy,'tlh'); // true 
  eqt = hlocmoy - hlocvraie; 
endfunction 

tab_eqt = []; 
trel = 0:365; 
an_ref = 2000; 

for an = 1950:2:2050
  eqt = equatemps(an, trel);
  tab_eqt = [tab_eqt ; eqt ]; 
end

eqt_ref = equatemps(an_ref, trel); 

eqt_min = min(tab_eqt , 'r'); 
eqt_max = max(tab_eqt , 'r'); 
ecmax = max(eqt_max - eqt_min); // max difference: (in case)

// ------------------------------------------
// plot
// ------------------------------------------
f = scf(); 
f.visible = "on"; 
f.immediate_drawing = "off"; 
a = gca();

plot(trel*12/365, eqt_min *60, 'b'); 
plot(trel*12/365, eqt_max *60, 'r'); 
plot(trel*12/365, eqt_ref *60, 'k'); 

CL_g_tag(a,0); 
h=CL_g_select(a, "Polyline");
h.thickness = 2; 

legend(['min between 1950 and 2050', 'max between 1950 and 2050', '2000']); 

// grid
for i=0:12
  plot2d([i,i], [-20,15], style = 1); 
end
CL_g_tag(a,1);
h = CL_g_select(a, "Polyline", 1); 
h.line_style = 3; 


a.title.text = 'Equation of time'; 
a.y_label.text = 'Mean local time minus true local time (mn)'; 

month_names = ['J', 'F', 'M', 'A', 'M', 'J', 'J', 'A', 'S', 'O', 'N', 'D']; 

a = gca(); 
a.x_ticks = tlist(['ticks', 'locations', 'labels'], 0.5:11.5, month_names); 
CL_g_stdaxes(a);
a.sub_ticks(1) = 0; // no substicks for a-axis 
a.grid = [-1,1]; 
a.grid_position = "background"; 


f.immediate_drawing = "on"; 
f.visible = "on";