//  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 tilt angle is the angle between the orbit plane and the SA
//> normal, positive if the angle between the SA normal and
//> the orbit's angular momentum is less than 90 degrees.  
//> 
//> - Fixed solar array: 
//> The orientation of the SA normal is defined by the tilt angle and
//> by the rotation angle (alpha) around w (origin = q axis). Alpha and
//> tilt are the spherical coordinates of the SA normal in qsw. 
//> 
//> - Rotating solar array: 
//> The solar array is rotating around the normal of the orbit plane
//> and is oriented towards the Sun as much as possible. 
//> 
//> 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;  // altitude (m)
sa_type = 1; 
sa_az = 0 * %CL_deg2rad; 
sa_tilts = (0:5:30) * %CL_deg2rad; 
beta_min = 0 * %CL_deg2rad; 
beta_max = 30 * %CL_deg2rad; 

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 rotation angle from q, around w (if fixed)", sa_az, units=['rad','deg']),..
  CL_defParam("Values of solar array tilt", sa_tilts, units=['rad','deg'], valid='$x >= -90 & $x <= 90', dim=[1, 20]),..
  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')..
);

[alt, sa_type, sa_az, sa_tilts, beta_min, beta_max] = CL_inputParam(desc_param);

// -----------------------------------------------------------
// Computation
// -----------------------------------------------------------
function [eff]= efficiency_rotating(betaa, tilt, alpha_ecl)
  // Mean value of max(cos(ang), 0)
  // cos(ang) = cos(beta - tilt) 
  eff = max(cos(betaa - tilt), 0) .* (alpha_ecl/%pi); 
endfunction

function [eff]= efficiency_fixed(betaa, az, tilt, alpha_ecl)
  // Mean value of max(cos(ang), 0)
  // cos(ang) = cos(beta)*cos(tilt)*cos(alpha+az) + sin(beta)*sin(tilt)

  // alpha1, alpha2: such that cos(ang) has same sign in [alpha1, alpha2]
  alpha1 = -az - %pi * ones(betaa); 
  alpha2 = -az + %pi * ones(betaa); 
  
  cond_sign = abs(sin(betaa) * sin(tilt)) < abs(cos(betaa) * cos(tilt)); 
  I = find(cond_sign); 
  if (I <> [])
    x = acos(-tan(betaa(I)) * tan(tilt)); 
    alpha1(I) = -az - x; 
    alpha2(I) = -az + x; 
  end
  
  // check value at middle of interval
  // if value is negative and sign changes => take complementary interval
  // if value is negative and sign does not change => 0 length interval

  I = find(cos(betaa-tilt) < 0 & cond_sign); 
  // does not happen if beta and tilt are in [-%pi/2, %pi/2]
  if (I <> [])
    [alpha1(I), alpha2(I)] = (alpha2(I), alpha1(I)+2*%pi); 
  end
  I = find(cos(betaa-tilt) < 0 & ~cond_sign); 
  if (I <> [])
    alpha1(I) = -az; 
    alpha2(I) = -az; 
  end

  // alpha1_ecl, alpha2_ecl: in shadow and inside [alpha1, alpha2]
  alpha_mid = (alpha1 + alpha2) / 2; 
  alpha0_ecl = CL_rMod(%pi, alpha_mid-%pi, alpha_mid+%pi);
  alpha1_ecl = alpha0_ecl - (%pi - alpha_ecl); 
  alpha2_ecl = alpha0_ecl + (%pi - alpha_ecl); 
  alpha1_ecl = max(alpha1, min(alpha2, alpha1_ecl)); 
  alpha2_ecl = max(alpha1, min(alpha2, alpha2_ecl)); 
  
  eff = (cos(betaa)*cos(tilt).*(sin(alpha2+az) - sin(alpha1+az) - sin(alpha2_ecl+az) + sin(alpha1_ecl+az)) + ..
         sin(betaa)*sin(tilt).*(alpha2 - alpha1 - alpha2_ecl + alpha1_ecl)) / (2*%pi); 

endfunction


// -----------------------------------------------------------

sma = alt + %CL_eqRad; // semi major axis
betaa = linspace(beta_min, beta_max, 200); 
alpha_ecl = %pi - CL_gm_betaEclipse(sma, betaa); 

// result = eff(%)  nb_tilts x N 
result = []; 

for tilt = sa_tilts
  if (sa_type == 1)
    eff = efficiency_fixed(betaa, sa_az, tilt, alpha_ecl); 
  else
    eff = efficiency_rotating(betaa, tilt, alpha_ecl); 
  end
  result = [result; 100 * eff]; 
end


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

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

nb = length(sa_tilts);
f.color_map = jetcolormap(min(max(nb,3),100));
Noir = addcolor([0,0,0]);

for k = 1:length(sa_tilts)
  plot(betaa * %CL_rad2deg, result(k,:));
  h = CL_g_select(gce(), "Polyline");
  h.thickness = 2;
  h.foreground = k; 
end

a = gca();

if (sa_type == 1); str = "fixed"; else; str = "rotating"; end
a.title.text = msprintf("Mean (%s) solar array efficiency - Alt = %.1f km", str, alt/1000);
a.x_label.text = "Sun beta angle (deg)"; 
a.y_label.text = "Efficiency (%)";

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

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

CL_g_legend(a, header = "SA tilt (deg)", str=string(sa_tilts * %CL_rad2deg));

CL_g_stdaxes(a);

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