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

// -----------------------------------------------------------
//> Mean solar array efficiency: 
//> It is defined as the mean value of max(cos(ang),0) over one orbit,
//> where ang is the angle between the normal of the solar array (SA) 
//> and the direction of the Sun.   
//> 
//> The SA may be rotating or not. 
//> If it is rotating, the SA normal points to the Sun as much as 
//> possible.  
//> 
//> The 2 required directions (defined by spherical coordinates) are: 
//> - The SA normal "at rest"  
//> - The SA rotation axis (if the array is rotating) 
//> 
//> NB: 
//> - The orbit is circular. 
//> - Eclipses are taken into account (spherical Earth). 
//> - The Sun is supposed fixed / orbit plane over one orbit. 
//
// Author: A. Lamy
// -----------------------------------------------------------

alt = 700.e3; 
sa_type = 1; 
// SA normal - alpha, delta in qsw
sa_n_ad = [0, 0] * %CL_deg2rad;  
// SA rot axis - alpha, delta in qsw
sa_axis_ad = [0, 90] * %CL_deg2rad;  
beta_min = 0 * %CL_deg2rad; 
beta_max = 90 * %CL_deg2rad; 
opt = 1; 

desc_param = list(..
  CL_defParam("Altitude (sma minus equ. radius)", alt, units=['m','km'], valid='$x>=0'),..
  CL_defParam("Solar array type: 1=fixed 2=rotating", sa_type, accv=[1,2]),..
  CL_defParam("SA normal / qsw frame (rot angle = 0) - sph. coord. ", sa_n_ad, units=['rad', 'deg'], dim=2),..
  CL_defParam("SA rotation axis / qsw frame - sph. coord. - if rotating", sa_axis_ad, units=['rad', 'deg'], dim=2),..
  CL_defParam("Sun beta angle - min", beta_min, id="$betamin", units=['rad', 'deg'], valid='$x>=-90 & $x<=90'),..
  CL_defParam("Sun beta angle - max", beta_max, units=['rad', 'deg'], valid='$x>=-90 & $x<=90 & $x>$betamin'),..
  CL_defParam("Plot factor for: 1=lighting, 2=PRS, 3=drag", opt, accv=1:3)..
);

[alt, sa_type, sa_n_ad, sa_axis_ad, beta_min, beta_max, opt] = CL_inputParam(desc_param);

n = CL_co_sph2car([sa_n_ad(1); sa_n_ad(2); 1]); 
axis = CL_co_sph2car([sa_axis_ad(1); sa_axis_ad(2); 1]); 

if (sa_type == 1); sa_typ_str = "fixed"; else sa_typ_str = "rotating"; end

// option (string)
OPTS = ["lighting", "PRS", "drag"]; 
opt_str = OPTS(opt); 

// print info on axes directions
mprintf("\n"); 
mprintf("Solar array type = %s\n", sa_typ_str); 
mprintf("SA normal in qsw frame = [%.3f; %.3f; %.3f]\n", n(1), n(2), n(3)); 
if (sa_type == 2)
  mprintf("SA rotation axis in qsw frame = [%.3f; %.3f; %.3f]\n", axis(1), axis(2), axis(3)); 
end

N = 201;  // number of "alpha" (arg of latitude, origin = Sun)
P = 200;  // number of "beta"


// -----------------------------------------------------------
// Computation
// -----------------------------------------------------------
function [M] = interv_linspace(a,b,n)
  M = repmat(a,n,1) + repmat((0:n-1)',1,size(a,2)) .* repmat(b-a,n,1) / (n-1); 
endfunction

sma = %CL_eqRad + alt; 

betaa = linspace(beta_min, beta_max, P); 

// numerical ajustements to avoid |beta| = %pi/2
eps = 1.e-10; 
betaa = max(min(betaa, %pi/2-eps), -%pi/2+eps); 

// alpha_ecl: PSO (from Sun) beginning of eclipse
// forced to %pi if eclipses not taken into account
alpha_ecl = %pi * ones(betaa); 
if (opt_str == "lighting" | opt_str == "PRS")
  alpha_ecl = %pi - CL_gm_betaEclipse(sma, betaa); 
end

Alpha = interv_linspace(-alpha_ecl, alpha_ecl, N);  // NxP
Betaa = repmat(betaa, N, 1); // NxP

alpha1 = matrix(Alpha, 1, -1); 
beta1 = matrix(Betaa, 1, -1); 

sun = CL_co_sph2car([zeros(beta1); beta1; ones(beta1)]); 

q_qsw = CL_rot_angles2quat(3, alpha1); 
sun_qsw = CL_rot_rotVect(q_qsw', sun); 

if (sa_type == 1) 
  sa_dir_qsw = repmat(n, 1, N*P); 
else
  q = CL_rot_defRot1Ax(n, axis, sun_qsw, 0); 
  sa_dir_qsw = CL_rot_rotVect(q, n); 
end

if (opt_str == "lighting")
  eff = matrix(max(CL_dot(sa_dir_qsw, sun_qsw), 0), N, P); 
elseif (opt_str == "PRS")
  eff = matrix(abs(CL_dot(sa_dir_qsw, sun_qsw)), N, P); 
elseif (opt_str == "drag")
  eff = matrix(abs(CL_dot(sa_dir_qsw, [0;1;0])), N, P); 
else
  error("Invalid value for opt_str");
end


// average value over orbit (including period in eclipse)
mean_eff = CL_mean(eff, meth="boole", "r") .* alpha_ecl / %pi; 


// -----------------------------------------------------------
// Plot
// -----------------------------------------------------------

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

plot(betaa * %CL_rad2deg, 100 * mean_eff, "-"); 

a = gca();

a.title.text = msprintf("Mean factor for %s - %s array - Alt = %.1f km", opt_str, sa_typ_str, alt/1000);
a.x_label.text = "Sun beta angle (deg)"; 
a.y_label.text = "Factor (%)";

a.data_bounds = [beta_min*%CL_rad2deg, max(min(100*mean_eff)-1, 0); 
                 beta_max*%CL_rad2deg, min(max(100*mean_eff)+1, 100)]; 
a.tight_limits = "on";

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

CL_g_stdaxes(a);

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