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

// -----------------------------------------------------------
//> Intersection of coplanar orbits  
//
// Author: A. Lamy
// -----------------------------------------------------------

// Inputs  
sma1 = 10; 
ecc1 = 0.8; 
pom1 = -%pi/4; 

sma2 = 10; 
ecc2 = 1.1; 
pom2 = %pi/4; 

rmax = 15; 

desc_param = list(..
   CL_defParam("Orbit1: semi-major axis", sma1, valid='$x>0'),..
   CL_defParam("Orbit1: eccentricity", ecc1, valid='$x>=0 & abs($x-1) > 1.e-6'),..
   CL_defParam("Orbit1: argument of periapsis", pom1, units=['rad', 'deg']),..
   CL_defParam("Orbit2: semi-major axis", sma2, valid='$x>0'),..
   CL_defParam("Orbit2: eccentricity", ecc2, valid='$x>=0 & abs($x-1) > 1.e-6'),..
   CL_defParam("Orbit2: argument of periapsis", pom2, units=['rad', 'deg']),..
   CL_defParam("Maximum radius for hyperbola plots", rmax, valid='$x>0')..
);

[sma1, ecc1, pom1, sma2, ecc2, pom2, rmax] = CL_inputParam(desc_param);

// -----------------------------------------------
// Utility functions
// -----------------------------------------------
// generate [x,y] for orbit plot
function [x, y] = genxy(sma, ecc, pom, rmax)
 
  anv_max = %pi; 
  if (ecc > 1)
	anv_max = abs(CL_gm_intersectCoplanOrb(sma, ecc, 0, rmax, 0, 0));
  end
  
  anv = linspace(-anv_max, anv_max, 201); 
  res = CL_kp_characteristics(sma, ecc, anv); 
  r = res.r; 
  pso = pom + anv; 
  
  x = r .* cos(pso); 
  y = r .* sin(pso); 

endfunction

// generate [x,y] of orbit intersection
function [x, y, pso, r] = genxyinters(sma1, ecc1, pom1, sma2, ecc2, pom2)

  [psoa, psob, nbsol] = CL_gm_intersectCoplanOrb(sma1, ecc1, pom1, sma2, ecc2, pom2);
  
  pso = [psoa, psob]; 
  res = CL_kp_characteristics(sma1, ecc1, CL_rMod(pso - pom1, -%pi, %pi)); 
  r = res.r; 
    
  x = r .* cos(pso); 
  y = r .* sin(pso); 
  
endfunction
  
  
// -----------------------------------------------
// Computation
// -----------------------------------------------
[xi, yi, psoi, ri] = genxyinters(sma1, ecc1, pom1, sma2, ecc2, pom2); 

rmax = max(1.3 * max(ri), rmax); 
[x1, y1] = genxy(sma1, ecc1, pom1, rmax); 
[x2, y2] = genxy(sma2, ecc2, pom2, rmax); 

// Display info
mprintf("\n----------------\n");
I = find(~isnan(psoi)); 
mprintf("Number of intersections = %d\n", length(I));
if (length(I) >= 1); mprintf("Intersection 1: PSO = %.2f deg\n", psoi(I(1)) * %CL_rad2deg); end
if (length(I) >= 2); mprintf("Intersection 2: PSO = %.2f deg\n", psoi(I(2)) * %CL_rad2deg); end


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

f = scf();
f.visible = "on"; 
f.immediate_drawing = "off"; 
 
plot(x1, y1, "b", "thickness", 2); 
plot(x2, y2, "r", "thickness", 2); 

plot(xi(1), yi(1), "ko", "thickness", 3); 
plot(xi(2), yi(2), "ko", "thickness", 3); 

// adjust plot (data bounds in particular)
a = gca(); 
a.isoview = "on"; 
a.tight_limits = "on"; 

db = a.data_bounds; 
L = 1.05 * max(db(2,:)-db(1,:), "c"); 
db = [mean(db, "r") - L/2; mean(db, "r") + L/2]; 
a.data_bounds = db;

CL_g_stdaxes(a); 

f.immediate_drawing = "on";