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

// -----------------------------------------------------------
//> Restricted 3-body problem:
//> (Normalized) acceleration or potential relative to the 
//> rotating frame for a point mass (with negligible mass) at rest 
//> in a 2-body system. 
//> 
//> The normalized values are such that: 
//> - distance between the 2 bodies = 1 (length unit = L)
//> - rotation period of the line joining the 2 bodies = 2*pi (time unit = T) 
//> It follows that: 
//> - sum of the gravitational constants of the 2 bodies = 1 (L^3/T^2)
//> - rotation rate of the line joining the 2 bodies = 1 (rad/T) 
//> 
//> The unnormalized values can be obtained from the normalized ones by: 
//> acc = acc_norm * mu / D^2
//> pot = pot_norm * mu / D
//> where: 
//> mu = sum of the gravitational constants of the 2 bodies
//> D = (constant) distance between the 2 bodies
//> 
//> The results depend on the ratio mu (body2) / mu (body1) where 
//> body1 is the body with the biggest mass. 
//> 
//> The frame is defined by: 
//> - center:  center of mass of the 2 main bodies. 
//> - X axis: directed from body1 (heaviest body) to body2. 
//> The 2 bodies are then fixed in this frame. 
//
// Author: A. Lamy
// -----------------------------------------------------------

// distance between body1 and body2 (here: normalized) 
D = 1; 

// sum of gravitational constants (here: normalized)
mu = 1; 

// mass ratio: m2/m1 (> 0 and <=1) 
ratio = 0.2; 

// plot option (see below for meaning)
iplot = 1; 

desc = list(..
   CL_defParam("Body mass ratio (>0 and <=1)", ratio, valid = "$x > 0 & $x <= 1"),..
   CL_defParam("Plot: 1=acceleration, 2=potential", iplot, accv = [1,2])..
); 

[ratio, iplot] = CL_inputParam(desc);


// -----------------------------------------------------------
// computation
// -----------------------------------------------------------

// The formulas that follow are also valid for unnormalized data 
// (i.e. "real" values for mu and D) 

// gravitational constants (mu1+mu2 = mu)
mu1 = mu / (1 + ratio); 
mu2 = mu * ratio / (1 + ratio); 

// angular rate (2-body problem)
omega = sqrt((mu1+mu2) / D^3); 

// positions of bodies (2 body problem)
p1 =  -(D * mu2 / (mu1 + mu2)) * [1;0;0]; 
p2 = (D * mu1 / (mu1 + mu2)) *  [1;0;0]; 

n = 150; 
x = linspace(-D, D, n) * 1.5; 
y = linspace(-D, D, n) * 1.5; 

[X, Y] = ndgrid(x, y); 

tab_x = X(:)'; 
tab_y = Y(:)'; 

pos = [tab_x; tab_y; zeros(tab_x)]; 

N = size(pos, 2); 
pos1 = repmat(p1, 1, N); 
pos2 = repmat(p2, 1, N); 

// acceleration and potentiel: 
// gravitational effects (bodies 1 and 2) + centrigugal force
 
acc = CL_fo_centralAcc(pos - pos1, mu1) + ..
      CL_fo_centralAcc(pos - pos2, mu2) + ..
      CL_fo_centrifugAcc(pos, [0; 0; omega]); 

pot = CL_fo_centralPot(pos - pos1, mu1) + ..
      CL_fo_centralPot(pos - pos2, mu2) + ..
      CL_fo_centrifugPot(pos, [0; 0; omega]); 

// text used in figure titles
txt = msprintf("Mass ratio = %.5g", ratio); 


// -----------------------------------------------------------
// plot function for acceleration
// -----------------------------------------------------------
function plot_acc(x, y, acc, txt)
  f = scf(); 
  f.visible = "on"; 
  f.immediate_drawing = "off"; 

  nbcol = 128;
  f.color_map = 0.4 + 0.6*jetcolormap(nbcol);

  z = CL_norm(acc); 

  zmin = 0; 
  zmax = 1.5; 

  z = min(z, zmax); 
  Z = matrix(z, size(X)); 

  colorbar(zmin,zmax,colminmax=[1,nbcol],fmt="%.2f");

  Sgrayplot(x, y, Z, colminmax=[1,nbcol], zminmax=[zmin, zmax], colout = [0,0]);

  levels = [0.06, 0.2, 0.4, 0.6]; 
  contour2d(x, y, Z, levels, style=color("black")*ones(levels));

  h = CL_g_select(gce(), "Text"); 
  CL_g_delete(h); 

  // soleil
  plot(p1(1), p1(2), "ko"); 
  h = CL_g_select(gce(), "Polyline"); 
  h.mark_size = 5; 

  plot(p2(1), p2(2), "ko");  
  h = CL_g_select(gce(), "Polyline"); 
  h.mark_size = 5; 

  a = gca(); 
  a.title.text = "(Normalized) Acceleration in rotating frame - " + txt; 
  a.x_label.text = "(Normalized) X"; 
  a.y_label.text = "(Normalized) Y"; 

  a.tight_limits = "on"; 

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

endfunction


// -----------------------------------------------------------
// plot function for potential
// -----------------------------------------------------------
function plot_pot(x, y, pot, txt)

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

  nbcol = 128;
  f.color_map = 0.4 + 0.6*jetcolormap(nbcol);

  z = pot; 

  zmin = 1.3; 
  zmax = 2; 
  z = min(z, zmax); 

  Z = matrix(z, size(X)); 

  colorbar(zmin,zmax,colminmax=[1,nbcol],fmt="%.2f");

  Sgrayplot(x, y, Z, colminmax=[1,nbcol], zminmax=[zmin, zmax]);

  levels = 1.4:0.1:1.8; 
  contour2d(x, y, Z, levels, style=color("black")*ones(levels));

  h = CL_g_select(gce(), "Text"); 
  CL_g_delete(h); 

  // soleil
  plot(p1(1), p1(2), "ko"); 
  h = CL_g_select(gce(), "Polyline"); 
  h.mark_size = 5; 

  plot(p2(1), p2(2), "ko");  
  h = CL_g_select(gce(), "Polyline"); 
  h.mark_size = 5; 

  a = gca(); 
  a.title.text = "(Normalized) Potential in rotating frame - " + txt; 
  a.x_label.text = "(Normalized) X"; 
  a.y_label.text = "(Normalized) Y"; 

  a.tight_limits = "on"; 

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

endfunction


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

if (iplot == 1)
  plot_acc(x, y, acc);
else
  plot_pot(x, y, pot);
end