//  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'.

// -----------------------------------------------------------
//> Sun elevation from any location on Earth 
//> - Elevation is the angle from the location between the Sun direction and its 
//>   projection onto the horizontal plane (tangent to Earth surface). 
//> - The location is on the ground and is defined by its geodetic longitude 
//>   and latitude.  
//> - The plotted parameter can be the elevation of the Sun or the cosine of the solar 
//>   zenith angle (=  sine of elevation)
//>
//> NB:
//> - UT1 is here considered the same as UTC 
//
// Author: A. Lamy
// -----------------------------------------------------------


// Sun elevation 
// lon, lat: geodetic longitude and latitude (rad) (1x1 or 1xN)
// t: UTC time (1x1 or 1xN)
function [elev] = sun_elev(lon, lat, t)

  N = max(size(lon,2), size(lat,2), size(t,2)); 
  lon = lon .* ones(1:N); 
  lat = lat .* ones(1:N); 
  t = t .* ones(1:N); 

  ra = CL_op_locTime(t, 'lon', lon, 'ra'); 

  pos = CL_co_ell2car([ra; lat; %CL_eqRad * ones(lat)]); 
  vert = CL_co_sph2car([ra; lat; ones(lat)]); 

  pos_sun = CL_eph_sun(t); 

  elev = %pi/2 - CL_vectAngle(vert, pos_sun-pos);  

endfunction


// -----------------------------------------------------------
// Parameters input
// -----------------------------------------------------------

lon = 0 * %CL_deg2rad; // geodetic longitude
lat = 45 * %CL_deg2rad; // geodetic latitude
opt = 1; // option (see below)

desc_param = list(..
   CL_defParam("Geodetic longitude of location", lon, units=['rad', 'deg'], valid='$x >= 0 & $x <= 360'),..
   CL_defParam("Geodetic latitude of location", lat, units=['rad', 'deg'], valid='$x >= -90 & $x <= 90'),..
   CL_defParam("Option: 1=elevation, 2=cos(zenith angle)", opt, accv = [1,2])..
);

[lon, lat, opt] = CL_inputParam(desc_param);


// -----------------------------------------------------------
// computations
// -----------------------------------------------------------

date_now = CL_dat_now("cal");
year = date_now(1);  // year number (current year)

t0 = CL_dat_cal2cjd(year,1,1); 
tf = CL_dat_cal2cjd(year+1,1,1); 
yearlen = tf-t0; 
t = t0 : tf; 

s = (t - t0) / yearlen; 

utmin = 0;
utmax = 24; 
ut = linspace(utmin, utmax, 200); 

[T, Ut] = ndgrid(t, ut); 
elev = matrix(sun_elev(lon, lat, T(:)'+Ut(:)'/24), size(T));  


// -----------------------------------------------------------
// plot
// -----------------------------------------------------------

f = scf(); 
f.visible="on";
f.immediate_drawing="off";

nbcol = 50;
f.color_map = 0.5 + 0.5*jetcolormap(nbcol);

col_bk = addcolor([0,0,0]); 
col_g1 = addcolor([1,1,1] * 0); 
col_g2 = addcolor([1,1,1] * 0.3); 

a=gca();
CL_g_tag(a, 0);

x = s; // abscissa for plot
y = ut; // ordinate for plot

if (opt == 1)
  z = elev * %CL_rad2deg; 
else 
  z = sin(elev); 
end

zmin = min(z); 
zmax = max(z); 

// Zmin, Zmax: forced range for colorbar
if (opt == 1) // case "elevation"
  Zmin = -90; 
  Zmax = 90; 
else // case "cos(sza)"
  Zmin = -1; 
  Zmax = 1; 
end

colorbar(Zmin,Zmax,colminmax=[1,nbcol],fmt="%.1f");

Sgrayplot(x, y, z, colminmax=[1,nbcol], ...
          zminmax=[Zmin, Zmax], colout=[0,0]);

[levels, sublevels] = CL_autoLevels(zmin, zmax);

contour2d(x, y, z, levels, style=col_g1*ones(levels));
CL_g_tag(a,1);

contour2d(x, y, z, sublevels, style=col_g2*ones(sublevels));
CL_g_tag(a,2);


// plot adjustments

h = CL_g_select(a, "Text", 2);
CL_g_delete(h);

h = CL_g_select(a, "Text");
CL_g_set(h, "text", string(strtod(h.text)));
h.font_foreground=col_g1;
h.font_size=2;
h.font_style=8;

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


// titles

a.x_label.text = "Day / Month (" + string(year) + ")"; 
a.y_label.text = "Time (UT)"; 

if (opt == 1) 
  str = "Sun elevation (deg)"; 
else
  str = "Cos(zenith angle)"; 
end

a.title.text = sprintf("%s - lon = %.1f deg, lat = %.1f deg", str, lon*%CL_rad2deg, lat*%CL_rad2deg); 


// ticks
m = [1:12, 1]; 
d = ones(m);
y = [year * ones(1:12), year+1]; 
t_ticks = CL_dat_cal2cjd(y, m, d); 
labels = "1/" + string(m); 
x_ticks = (t_ticks - t0) / yearlen; 

a.x_ticks = tlist(['ticks', 'locations', 'labels'], x_ticks, labels);

y_ticks = 0:2:24; 
a.y_ticks = tlist(['ticks', 'locations', 'labels'], y_ticks, string(y_ticks));

a.data_bounds = [0, 0; 1, 24]; 
a.tight_limits = "on"; 

CL_g_stdaxes(a, colg=col_g2);
a.sub_ticks = [0,1]; // no substicks 


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