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

// -----------------------------------------------------------
//> Ground stations visibility:  
//> 
//> Main view: ground tracks, portions in view of ground stations, 
//>    longest track between visibility passes (black curve). 
//> Secondary view: visibility passes. 
//> Notes: 
//> - The orbit is propagated considering the J2 secular effects only. 
//> - There is a small unconsistency regarding ground stations: 
//>   coordinates are geodetic in the computation, but spherical 
//>   in the graph. 
//> - Visibility passes of less than 1 mn are discarded.
//> 
//> Use right button to create / delete ground stations. 
//
// Author: A. Lamy
// -----------------------------------------------------------


t0 = 21915;
pas = 60.0 / 86400.0;
duree = 1; // day

par0 = [%CL_eqRad+700.e3; 0.; 98.2*%CL_deg2rad; 0; 0; 0];
//par0 = [42164.e3; 0.2; 63.4*%CL_deg2rad; 3*%pi/2; -%pi/2; %pi];

elevmin = 5*%CL_deg2rad;
distmax = 1.e9; 

desc_param = list(..
   CL_defParam("Initial time (1950.0)", t0),..
   CL_defParam("Semi-major axis", par0(1), units = ["m", "km"], valid='$x > 0'),..
   CL_defParam("Eccentricity", par0(2), valid='$x >= 0 & $x < 1'),..
   CL_defParam("Inclination", par0(3), units = ["rad", "deg"], valid='$x >= 0 & $x <= 180'),..
   CL_defParam("Argument of perigee", par0(4), units = ["rad", "deg"]),..
   CL_defParam("Longitude of asc. node", par0(5), units = ["rad", "deg"]),..
   CL_defParam("Mean anomaly", par0(6), units = ["rad", "deg"]),..
   CL_defParam("Simulation duration", duree, units = ["day"], valid='$x > 0'),..
   CL_defParam("Simulation time step", pas, units = ["day", "s"], valid='$x > 0'),..
   CL_defParam("Minimum elevation", elevmin, units=["rad", "deg"], valid='$x>=-90 & $x<=90')..
);

[t0,par0(1),par0(2),par0(3),par0(4),par0(5),par0(6),duree,pas,elevmin] = CL_inputParam(desc_param);

kep0 = par0; 
kep0(5) = kep0(5) + CL_mod_siderealTime(t0); // right ascending of ascending node

// Orbit calculation
t1 = t0 + duree; 
t = t0 : pas : t1;
kep = CL_ex_secularJ2(t0, kep0, t); // propagation

pos_ECI = CL_oe_kep2car(kep); // position inertial frame
pos_ECF = CL_fr_convert("ECI", "ECF", t, pos_ECI);


tabcol = [5,3,4,6,7,12,14,18,29,32];

data = struct("t0", t0, "t1", t1, "step", pas, "kep0", kep0, "elevmin", elevmin, "distmax", distmax, "t", t, "pos_ECF", pos_ECF, "tabcol", tabcol); 


// -----------------------------------------------------------
// Event handler
// -----------------------------------------------------------
function myHandler(win, x, y, ibut)


// ===========================================================
// COMPUTATIONS
// ===========================================================

// -----------------------------------------------------------
// Ephemeris between t1 and t2 
// (interior points identical to that of data.t)
// -----------------------------------------------------------
function [pos_ECF] = compute_ephem_ter(data, t1, t2)

  I = find(data.t > t1 & data.t < t2); 
  t = [t1, data.t(I), t2]; 
  kep = CL_ex_secularJ2(data.t0, data.kep0, t);
  pos_ECI = CL_oe_kep2car(kep);
  pos_ECF = CL_fr_convert("ECI", "ECF", t, pos_ECI);

endfunction


// -----------------------------------------------------------
// visibility calculation for one ground station
// -----------------------------------------------------------
function [res] = compute_visi(data, station)

  tvisi = CL_ev_stationVisibility(data.t, data.pos_ECF, station, data.elevmin);

  arcs_ter = list();
  for k = 1:size(tvisi,2)
    t1 = tvisi(1,k);
    t2 = tvisi(2,k);
    arcs_ter(k) = compute_ephem_ter(data, t1, t2); 
  end

  // results : 
  res = struct("id", 0, "sta", station, "tvisi", tvisi, "arcs_ter", arcs_ter); 

endfunction


// -----------------------------------------------------------
// visibility calculation for all ground stations
// -----------------------------------------------------------
function [res_all] = compute_visi_all(data, list_res)

  res_all= [];
  if isempty(list_res); return; end

  tvisi = [];
  for k = 1:length(list_res) 
    tvisi = CL_intervUnion(tvisi, list_res(k).tvisi); 
  end

  // longest track
  tdeb = [data.t0, tvisi(2,:)];
  tfin = [tvisi(1,:), data.t1];

  [dt, k] = max(tfin-tdeb);
  t1 = tdeb(k);
  t2 = tfin(k);
  pos_ECF = compute_ephem_ter(data, t1, t2);

  res_all = struct("tvisi", tvisi, "t1_lg", t1, "t2_lg", t2, "pos_ECF_lg", pos_ECF); 

endfunction


// -----------------------------------------------------------
// print statistics
// -----------------------------------------------------------
function print_stat(data, list_res, res_all)

  if isempty(list_res); return; end
  if isempty(res_all); return; end
  if isempty(res_all.tvisi); return; end

  duree = (res_all.tvisi(2,:)-res_all.tvisi(1,:));  
  mean_dur = sum(duree/(data.t1-data.t0)); 
  max_hole = res_all.t2_lg - res_all.t1_lg; 

  printf("\n-----------------------\n");
  printf("Ground stations: ");
  for k=1:length(list_res); printf("GS%d ", gsort(list_res(k).id, 'c', 'i')); end 
  printf("\n");

  printf("Mean visibility duration (mn/day) : %.2f\n", mean_dur*1440); 
  printf("Max (combined) pass visibility duration (mn) : %.2f\n", max(duree)*1440);
  printf("Max duration without visibility (h) : %.2f\n", max_hole * 24); 

endfunction


// ===========================================================
// PLOT
// ===========================================================

// -----------------------------------------------------------
// Plot visibility passes
// -----------------------------------------------------------
function plot_visi_pass(f, data, res, onoff)

  scf(f); 
  f.immediate_drawing="off";
  a=gca();

  icol = data.tabcol(modulo(res.id-1, 10)+1); 
  tag = 10000 + res.id; 

  for k=1:1; 

    h = CL_g_select(a, 'Compound', tag);
    CL_g_delete(h);

    if (onoff==0); continue; end
    if isempty(res.tvisi); continue; end

    duree = res.tvisi(2,:) - res.tvisi(1,:); 

    for i=1:size(res.tvisi,2)
      x = ([res.tvisi(1,i),res.tvisi(2,i)]-data.t0)*24; 
      y = [0, duree(i)]*1440; 
      plot2d(x, y, style=icol); 
    end
    CL_g_tag(a, tag);

    a.data_bounds = [0, 0; (data.t1-data.t0)*24, max(duree*1440)*1.05]; 
    a.tight_limits = "on"; 
    CL_g_stdaxes(a);
    a.grid_position = "background";

    h = CL_g_select(a, "Polyline", tag);
    h.thickness = 3; 
  
  end

  // legend
  // removed because of scilab bug (depending on version)
  // h = CL_g_select(a, "Compound"); 
  // if (length(h) > 0)
  //   [str,I] = unique(string(floor(modulo(h.user_data,10000)))); 
  //   h = h(I);
  //   str = matrix(str, 1, -1); // row vector
  //   if (str <> []) str = "GS" + str; end
  //   hl = legend(h,str);
  //   hl.legend_location = "by_coordinates";
  //   hl.position = [0.8, 0.2]; 
  //   hl.line_mode = "off";
  // end

  f.immediate_drawing="on";

endfunction 


// -----------------------------------------------------------
// Shows ground tracks in visibility
// -----------------------------------------------------------
function plot_track_visi(f, data, res, onoff)

  scf(f);
  f.immediate_drawing="off";
  a = gca();

  icol = data.tabcol(modulo(res.id-1, 10)+1); 
  tag = 20000 + res.id; 

  h = CL_g_select(a, 'Compound', tag);
  CL_g_delete(h);

  for k=1:1

    if (onoff==0); continue; end
    if isempty(res.tvisi); continue; end

    // "%nan" added between arcs
    pos_ECF = [%nan; %nan; %nan]; 
    for i=1:length(res.arcs_ter)
      pos_ECF = [pos_ECF, res.arcs_ter(i), [%nan; %nan; %nan]]; 
    end
    CL_plot_ephem(pos_ECF,color_id=10); 
    CL_g_tag(a, tag);

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

  end

  CL_g_stdaxes(a);
  f.immediate_drawing="on";

endfunction 

// -----------------------------------------------------------
// Shows ground station 
// -----------------------------------------------------------
function plot_station(f, data, res, onoff)

  scf(f);
  f.immediate_drawing="off";
  a = gca();

  icol = data.tabcol(modulo(res.id-1, 10)+1); 
  tag = 30000 + res.id; 

  h = CL_g_select(a, 'Compound', tag);
  CL_g_delete(h);
  h = CL_g_select(a, 'Text', tag);
  CL_g_delete(h);

  for k=1:1
  
    if (onoff==0); continue; end

    plot(res.sta(1)*%CL_rad2deg, res.sta(2)*%CL_rad2deg, "sk");
    xstring(res.sta(1)*%CL_rad2deg, res.sta(2)*%CL_rad2deg, ...
            " GS" + string(res.id)); 
    CL_g_tag(a, tag);

    h = CL_g_select(a, "Polyline", tag);
    h.mark_mode = "on"; 
    h.line_mode = "off"; 
    h.mark_background = icol; 
    h.mark_foreground = 1; 
    h.mark_size = 7; 

    h = CL_g_select(a, "Text", tag);
    h.font_size=2;
    h.font_style=8;

  end

  CL_g_stdaxes(a);
  a.grid_position = "background";
  f.immediate_drawing="on";

endfunction 

// -----------------------------------------------------------
// Shows longest ground tracks 
// -----------------------------------------------------------
function plot_longest_track(f, data, res_all)

  scf(f);
  f.immediate_drawing="off";
  a = gca();

  tag = 40000; 
  h = CL_g_select(a, 'Compound', tag);
  CL_g_delete(h);

  for k=1:1
  
    if isempty(res_all); continue; end
    if isempty(res_all.pos_ECF_lg); continue; end

    CL_plot_ephem(res_all.pos_ECF_lg,color_id=1); 
    CL_g_tag(a, tag);

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

  end

  CL_g_stdaxes(a);
  a.grid_position = "background";
  f.immediate_drawing="on";

endfunction 

// ===========================================================
// MANAGEMENT
// ===========================================================

  function [id] = nextid(result)
    tabid = []; 
    for k=1:length(result); tabid = [tabid, result(k).id]; end
    tabid = gsort(tabid, 'c', 'i'); 
    id = length(result)+1; 
    for k=1:length(tabid); if (k <> tabid(k)); id = k; break; end; end
  endfunction 

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

  ud = f.user_data;

  // INIT variables %CL...
  %CL_rad2deg = ud.CLDATA.rad2deg; 
  %CL_deg2rad = ud.CLDATA.deg2rad; 

  data = ud.data ; 
  result = ud.list_res; 
  f2 = ud.f2; 
  
  // pointer location in pixels (x,y) => (xu, yu) in degrees
  scf(f); // needed for the conversion
  [xu, yu, rect] = xchange(x,y,"i2f"); 

  if (x <= rect(1) | x >= rect(1)+rect(3) | y <= rect(2) | y >= rect(2)+rect(4))
     f.info_message = "Out of bounds!";
     return;
  end
  
  lon = xu * %CL_deg2rad;
  lat = yu * %CL_deg2rad;
  sta = [ lon; lat; 0]; 

  num=0; 
  distmax = %inf; 
  for k=1:length(result)
    dist = CL_vectAngle(CL_co_ell2car(sta), CL_co_ell2car(result(k).sta)); 
    if (dist < distmax & dist < 5*%CL_deg2rad); num=k; end
  end

  if (ibut == -1) 
    s3 = "[Right click => Add]"; 
    if (num > 0)
      s3 = "[GS" + string(result(num).id) + ": Right click => Delete]"; 
      xu = result(num).sta(1) * %CL_rad2deg; 
      yu = result(num).sta(2) * %CL_rad2deg; 
    end 
    s1 = sprintf("Long: %.3f", xu); 
    s2 = sprintf("Lat: %.3f", yu); 

    f.info_message = sprintf("%-15s %-15s %s", s1, s2, s3);
  end

  if (ibut == 5)
    // initialisation des variables utiles aux calculs
    CL_init(init_glob = %f); 

    if (num == 0) 
      res = compute_visi(data, sta); 
      res.id = nextid(result); 
      f.info_message = "Creating GS" + string(res.id) + "..."; 
      num = length(result)+1; 
      result(num) = res; 
    else
      res = result(num);
      f.info_message = "Deleting GS" + string(res.id) + "..."; 
      result(num) = null();
      num = 0;
    end
    plot_visi_pass(f2, data, res, num); 

    scf(f);
    plot_track_visi(f, data, res, num); 
    plot_station(f, data, res, num); 
    res_all = compute_visi_all(data, result); 
    plot_longest_track(f, data, res_all); 
    print_stat(data, result, res_all); 
    f.info_message = ""; 

    ud.list_res = result; 
    set(f, 'user_data', ud);

  end

endfunction 


// ===========================================================
// MAIN
// ===========================================================

// -----------------------------------------------------------
// visibility passes view
// -----------------------------------------------------------
f2 = scf();
f2.event_handler_enable="off";
f2.immediate_drawing="off";
f2.axes_size = [400,300];
f2.figure_position = [0,0];
f2.visible = "on";

a=gca();
a.x_label.text = 'Time from beginning of simulation [hours]';
a.y_label.text = 'Time from beginning of pass [mn]';
a.title.text = 'Visibility passes';
CL_g_stdaxes();

// margin for x axis (text)
a.margins(4) = 0.15;

f2.immediate_drawing="on";
f2.visible = "on";


// -----------------------------------------------------------
// ground track view
// -----------------------------------------------------------
f=scf();
f.visible = "on"; 
f.immediate_drawing="off";

f.event_handler_enable="off";
Col = addcolor([0.6,0.7,0.8]);

CL_plot_earthMap(color_id=Col); 
a=gca();
CL_g_tag(a,0);
h=CL_g_select(a, "Polyline", 0);
h.thickness=2;

CL_plot_ephem(pos_ECF, color_id=10); 
CL_g_tag(a, 0); 

a.x_label.text = 'Longitude [deg]';
a.y_label.text = 'Latitude [deg]';
a.title.text = 'Ground tracks';

CL_g_stdaxes(a);

f.immediate_drawing="on";


// -----------------------------------------------------------
// other initializations
// -----------------------------------------------------------

// set handler and store results
ud = f.user_data;
ud.eventHandler = myHandler; 

ud.data = data; 
ud.list_res = list(); 
ud.f2 = f2; 

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

set(f, 'user_data', ud);

scf(f);
show_window(f);

//seteventhandler('myHandler');  // uncomment this line if demo is to be used on its own
f.event_handler_enable = 'on';