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

// -----------------------------------------------------------
//> Repeating ground tracks:  
//> Gives the semi-major axis of orbits such that the ground tracks 
//> repeat after a certain number of planet revolutions (Q). 
//> 
//> Move the mouse over the plot and click near a symbol to get more 
//> information. 
//
// Author: A. Lamy
// -----------------------------------------------------------

smamin = 6878.e3;
smamax = 7278.e3;
Qmin = 1;
Qmax = 20;
sso = 1; // sun-synchronous
inc = 0; // unused if SSO
ecc = 0; 
model = 2;
iplot = 1; // x axis = sma

desc_param = list(..
   CL_defParam("Semi-major axis - min", smamin, units=['m', 'km'], id='$smamin', valid='$smamin>=0' ),..
   CL_defParam("Semi-major axis - max", smamax, units=['m', 'km'], id='$smamax', valid='$smamax>$smamin' ),..
   CL_defParam("Eccentricity", ecc, valid='$x>=0 & $x <1' ),..
   CL_defParam("Sun-synchronous orbit? (1=yes, 0=no)", sso, id='$sso', accv=[0,1], valid="$sso == 0 | $model == 2"),..
   CL_defParam("Inclination (if not sun-synchronous)", inc, units=['rad', 'deg'], valid="$sso == 1 | ($x >=0 & $x<=180)"),..
   CL_defParam("Q (number of planet revolutions / node) - min", Qmin, id='$Qmin', accv = 1:100),..
   CL_defParam("Q (number of planet revolutions / node) - max", Qmax, id='$Qmax', accv = 1:100, valid='$Qmax >= $Qmin'),..
   CL_defParam("Model (1=kepler, 2=J2)", model, accv = [1,2], id='$model'),..
   CL_defParam("X_axis (1=semi-major axis, 2=altitude)", iplot, accv = [1,2])..
);

[smamin,smamax,ecc,sso,inc,Qmin,Qmax,model,iplot] = CL_inputParam(desc_param);

j2 = %CL_j1jn(2);
if (model == 1); j2 = 0; end

result = CL_op_searchRepeatOrbits(smamin,smamax,Qmin,Qmax,ecc,sso,inc,j2=j2);


// -----------------------------------------------------------
// Utilities
// -----------------------------------------------------------

// trace des positions des traces au sol
function trace_NPQ(f, N, P, Q)
  f_old = gcf();
  scf(f);
  f.immediate_drawing = "off";
  clf();
  Gris = addcolor(0.6*[1,1,1]); 
   
  a=gca();
 
  for (i=-2:Q+2)
    plot2d([i,i], [-2,Q+2], style=Gris);
    plot2d([-2,Q+2], [i,i], style=Gris);
  end
  CL_g_tag(a,2);

  res = zeros(1,4);
  if (Q >= 1)
     res = CL_op_repeatGroundTracks(N,P,Q, option=1); 
     plot(res(:,2), res(:,4),'s');
  end; 
  CL_g_tag(a,1);

  h = CL_g_select(a, "Polyline",1);
  h.mark_size=6;
  h.mark_background=24;
  h.mark_foreground=22;
  
  h = CL_g_select(a, "Polyline",2);
  h.line_style=1;

  mg = 0.4*max(Q,1)/10;
  a.data_bounds=[min(res(:,2))-mg, min(res(:,4))-mg; max(res(:,2))+mg, max(res(:,4))+mg];
  a.tight_limits = "on";
  CL_g_stdaxes(a);
  a.grid_position = "background";
  a.grid = [-1,-1];
  a.sub_ticks = [0,0];

  a.x_label.text = "Number of planet revolutions";
  a.y_label.text = "Positions of ground tracks";
  if (Q > 0); a.title.text = sprintf("N,P,Q = (%d,%d,%d)", N, P, Q);
  else; a.title.text = sprintf("N,P,Q = (?,?,?)"); end

  f.immediate_drawing = "on";
  scf(f_old);

endfunction 

// event handler
function myHandler(win, xpix, ypix, but)

  f = get_figure_handle(win);
  if isempty(f); return; end

  scf(f);
  sca(f.children(1)); // axes 
  
  ud = f.user_data; 

  // INIT variables %CL...
  %CL_rad2deg = ud.CLDATA.rad2deg; 
  %CL_eqRad = ud.CLDATA.eqRad; 

  result = ud.result;
  f2 = ud.f2; //second window;
  iplot = ud.iplot; 
  fct_trace_NPQ = ud.fct_trace_NPQ;

  // points in pixels (x,y)
  if (iplot == 1)
    x = result(:,1)/1000; // sma
  else
    x = (result(:,1) - %CL_eqRad)/1000; // alt
  end
  y = result(:,6);
  [result_xpix, result_ypix, rect] = xchange(x,y,"f2i"); 

  // find closest solution (x = result(1,:)/1000, y=result(:,6)) 
  [dist,index] = min(sqrt((result_xpix-xpix).^2+(result_ypix-ypix).^2)); 

  N = result(index,4);
  P = result(index,5);
  Q = result(index,6);

  if (but < 0) 
    alt = (result(:,1)-%CL_eqRad)/1000;
    s1 = sprintf("N,P,Q: %d %d %d", N, P, Q);
    if (iplot == 1)
      s2 = sprintf("Sma: %.3f km", (result(index,1))/1000);
    else
      s2 = sprintf("Alt: %.3f km", (result(index,1)-%CL_eqRad)/1000);
    end
    s3 = sprintf("Duration: %.3f days", result(index,9));
    xinfo(sprintf("%-20s %-20s %-20s", s1, s2, s3));
  
  elseif (but == 0 | but == 3)
    fct_trace_NPQ(f2,N,P,Q);
    printf("\n");
    printf("N,P,Q = %d %d %d\n", N,P,Q);
    printf("Semi major axis (km) = %.6f\n", result(index,1)/1000);
    printf("Eccentricity   = %6f\n", result(index,2));
    printf("Inclination (deg) = %.6f\n", result(index,3) * %CL_rad2deg);
    printf("Duration of repeat cycle (days) = %.6f\n", result(index,9));
    printf("Number of orbits per cycle = %d\n", result(index,7));
    printf("Nodal period (s) = %.6f\n", result(index,8));
    printf("Mesh size (deg) = %.6f\n", result(index,10)*%CL_rad2deg);
  end

endfunction 


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

// figure for ground tracks positions
f2 = scf();
f2.event_handler_enable="off";
f2.axes_size = [500,400];
f2.figure_position = [0,0];
f2.visible = "on"; 
trace_NPQ(f2, 0,0,0);

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

if (iplot == 1) 
  xaxis = result(:,1); 
  xlabel = "Semi major axis [km]";
  xaxis_min = smamin; 
  xaxis_max = smamax; 
else
  xaxis = result(:,1) - %CL_eqRad; 
  xlabel = "Altitude (semi major axis minus equ. radius) [km]";
  xaxis_min = smamin - %CL_eqRad; 
  xaxis_max = smamax - %CL_eqRad; 
end

plot(xaxis/1000, result(:,6),'d');  
[xmin, xmax] = CL_graduate(min(xaxis), max(xaxis));

a = gca();
a.data_bounds=[max(xmin,xaxis_min)/1000,Qmin-1;min(xmax,xaxis_max)/1000,Qmax+1];
a.tight_limits = "on";
CL_g_stdaxes(a);

a.grid_position = "background";
a.x_label.text = xlabel;
a.y_label.text = 'Number of planet revolutions (Q)';
if (sso == 0)
  a.title.text = 'Repeat orbits (Ecc = ' + string(ecc) + ', Inc = ' + string(inc*%CL_rad2deg) + ' deg)';
else
  a.title.text = 'Sun-synchronous repeat orbits (Ecc = ' + string(ecc) + ')';
end

h = CL_g_select(a, "Polyline");
h.mark_size=6;
h.mark_background=18;
h.mark_foreground=2;

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


ud = f.user_data;
ud.eventHandler = myHandler; 
ud.fct_trace_NPQ = trace_NPQ; 
ud.result = result;
ud.f2 = f2;
ud.iplot = iplot; 

// initialisation des variables "globales" de celestLab utiles
// au event Handler
ud.CLDATA = struct("rad2deg", %CL_rad2deg, "eqRad", %CL_eqRad);

set(f, 'user_data', ud);

scf(f);
show_window(f);

//seteventhandler('myHandler'); 
f.event_handler_enable = 'on';