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

// -----------------------------------------------------------
//> Frame conversion utility: 
//> Converts position and velocity from one frame to another.
//> 
//> The managed (Earth centered) frames are: 
//>   TIRS         -> Terrestrial Intermediate Reference Frame
//>   ITRS         -> International Terrestrial Reference Frame
//>   PEF          -> Pseudo Earth Fixed
//>   CIRS         -> Celestial Intermediate Reference Frame
//>   Veis         -> Veis 1950
//>   TOD          -> True Of Date
//>   TEME         -> True Equator Mean Equinox
//>   GCRS or ICRS -> Geocentric (or International) Celestial Reference Frame
//>   EME2000      -> Earth Mean Equator at J2000
//>   MOD          -> Mean Of Date
//>   EOD          -> Ecliptic Of Date
//>   ECI          -> Earth Centered (Earth) Inertial
//>   ECF          -> Earth Centered (Earth) Fixed
//>
//> In addition: FRAME1/FRAME2 means: all coordinates given in FRAME1 
//> but velocity is relative to FRAME2. 
//
// Author: A. Lamy
// -----------------------------------------------------------


// ------------------------------
// utility functions
// ------------------------------

// Tranform input string into 2 parts 
// (separator = "/") + remove spaces
function [frames] = decodeframe(sframe) 
   n = strindex(sframe, "/");
   if (n == [] | size(n,"*") > 1)
     frame1 = stripblanks(sframe, %t); 
     frame2 = frame1; 
   else
     frame1 = stripblanks(part(sframe,1:n-1),%t);
     frame2 = stripblanks(part(sframe,n+1:length(sframe)),%t);
   end
   frames = [frame1, frame2];  
endfunction

// check string is valid: 
// frame1 OR frame1/frame2
function [ok] = isvalidframe(sframe)  
   [frames] = decodeframe(sframe) 
   ok  = ( find(frames(1) == FRAME_NAMES) <> [] & ..
           find(frames(2) == FRAME_NAMES) <> []); 
endfunction

// check that the export name (for results) is valid
// The global variable "exportname" should be empty ([])
function [ok] = isvalidexportname(exportname)  
   if (exportname == ""); ok=%t; return; end
   
   cmd = msprintf("global %s;", exportname); 
   err = execstr(cmd, "errcatch"); 
   if (err <> 0); ok=%f; return; end
     
   cmd = msprintf("ok = (%s == []);", exportname); 
   err = execstr(cmd, "errcatch"); 
   if (err <> 0); ok=%f; return; end
endfunction


// print the results
function print_results()  

  // recompute calendar date 
  cal = CL_dat_cjd2cal(cjd); 

  mprintf("\n"); 
  mprintf("Date (TREF)           = %s (%.8f)\n", CL_dat_cal2str(cal), cjd); 
  mprintf("TT - TREF (s)         = %.15g\n", tt_tref); 
  mprintf("UT1 - TREF (s)        = %.15g\n", ut1_tref); 
  mprintf("Pole: xp, yp (arcsec) = %.15g %.15g\n", ..
        CL_unitConvert(xyp(1), "rad", "arcsec"), ..
        CL_unitConvert(xyp(2), "rad", "arcsec")); 
  mprintf("\n"); 

  mprintf("Frame: %s", frames1(1)); 
  if (frames1(1) <> frames1(2))
    mprintf(" / velocity relative to %s", frames1(2)); 
  end
  mprintf("\n"); 

  mprintf("Position - X (km)     = %.15g\n", pos1(1)/1000); 
  mprintf("Position - Y (km)     = %.15g\n", pos1(2)/1000); 
  mprintf("Position - Z (km)     = %.15g\n", pos1(3)/1000); 
  mprintf("Velocity - X (km/s)   = %.15g\n", vel1(1)/1000); 
  mprintf("Velocity - Y (km/s)   = %.15g\n", vel1(2)/1000); 
  mprintf("Velocity - Z (km/s)   = %.15g\n", vel1(3)/1000); 
  mprintf("\n"); 

  mprintf("Frame: %s", frames2(1)); 
  if (frames2(1) <> frames2(2))
    mprintf(" / velocity relative to %s", frames2(2)); 
  end
  mprintf("\n"); 

  mprintf("Position - X (km)     = %.15g\n", pos2(1)/1000); 
  mprintf("Position - Y (km)     = %.15g\n", pos2(2)/1000); 
  mprintf("Position - Z (km)     = %.15g\n", pos2(3)/1000); 
  mprintf("Velocity - X (km/s)   = %.15g\n", vel2(1)/1000); 
  mprintf("Velocity - Y (km/s)   = %.15g\n", vel2(2)/1000); 
  mprintf("Velocity - Z (km/s)   = %.15g\n", vel2(3)/1000); 
  mprintf("\n"); 

endfunction


// export to global variable "export_name"
function export_results()
  
  res = struct(); 
  res.tt_tref = tt_tref;
  res.ut1_tref = ut1_tref;
  res.t_cjd = cjd;
  res.t_cal = scal;	  
  res.xp = xyp(1);
  res.yp = xyp(2);
  res.frame1_coord = frames1(1);
  res.frame1_vel = frames1(2);
  res.pos1 = pos1;
  res.vel1 = vel1;
  res.frame2_coord = frames2(1);
  res.frame2_vel = frames2(2);
  res.pos2 = pos2;
  res.vel2 = vel2;

  cmd = msprintf("global %s; %s = res;", export_name, export_name); 
  status = execstr(cmd, "errcatch"); 
  if (status == 0)
    mprintf("-> Results exported to global variable: %s\n\n", export_name); 
  end
  
endfunction


// frame names
FRAME_NAMES = [ "ECF", "TIRS", "ITRS", "PEF", ..
                "ECI", "CIRS", "Veis", "TOD", "TEME", ..
			          "GCRS", "EME2000", "MOD", "EOD", "ICRS"]; 


// -----------------------------------------------------------
// Initializations
// -----------------------------------------------------------

pos1 = [7.e6; 0; 0];
vel1 = [0; 7.e3; 0];
scal = []; // calendar date: auto

sframe1 = "GCRS"; 
sframe2 = "ITRS"; 

xyp = [0, 0]; // pole (rad)
tt_tref = %CL_TT_TREF; 
ut1_tref = %CL_UT1_TREF; 


// name of variable for export
export_name = ""; 

desc = list(..
   CL_defParam("Date: (calendar format, TREF time scale)", scal, typ='cal'),..
   CL_defParam("Position / frame1 - X", pos1(1), units=['m','km']),..
   CL_defParam("Position / frame1 - Y", pos1(2), units=['m','km']),..
   CL_defParam("Position / frame1 - Z", pos1(3), units=['m','km']),..
   CL_defParam("Velocity / frame1 - X", vel1(1), units=['m/s','km/s']),..
   CL_defParam("Velocity / frame1 - Y", vel1(2), units=['m/s','km/s']),..
   CL_defParam("Velocity / frame1 - Z", vel1(3), units=['m/s','km/s']),..
   CL_defParam("TT - TREF", tt_tref, units=['s']),..
   CL_defParam("UT1 - TREF", ut1_tref, units=['s']),..
   CL_defParam("Pole axis: xp, yp", xyp, units=['rad', 'arcsec'], dim=2),..
   CL_defParam("Frame1 (=> see help for list)", sframe1, typ="s", valid='isvalidframe($x)'),..
   CL_defParam("Frame2 (=> see help for list)", sframe2, typ="s", valid='isvalidframe($x)'),..
   CL_defParam("Export results to (empty) global variable...", export_name, typ="s", valid='isvalidexportname($x)')..
);

[scal, pos1(1), pos1(2), pos1(3), vel1(1), vel1(2), vel1(3), .. 
 tt_tref, ut1_tref, xyp, sframe1, sframe2, export_name] = CL_inputParam(desc);

// -----------------------------------------------------------
// Conversion
// -----------------------------------------------------------

// calendar -> cjd
cal = CL_dat_str2cal(scal); 
cjd = CL_dat_cal2cjd(cal); 

// convert frames
frames1 = decodeframe(sframe1); 
frames2 = decodeframe(sframe2); 

// intermediate variables
P = pos1; 
V = vel1; 

// change coordinates only
if (frames1(1) <> frames1(2))
  mat = CL_fr_convertMat(frames1(1), frames1(2), cjd, tt_tref=tt_tref, ut1_tref=ut1_tref, xp=xyp(1), yp=xyp(2));
  P = mat * P;
  V = mat * V; 
end

[P, V] = CL_fr_convert(frames1(2), frames2(2), cjd, P, V, tt_tref=tt_tref, ut1_tref=ut1_tref, xp=xyp(1), yp=xyp(2));

// change coordinates only
if (frames2(1) <> frames2(2))
  mat = CL_fr_convertMat(frames2(2), frames2(1), cjd, tt_tref=tt_tref, ut1_tref=ut1_tref, xp=xyp(1), yp=xyp(2));
  P = mat * P;
  V = mat * V; 
end

pos2 = P; 
vel2 = V; 


// -----------------------------------------------------------
// Results
// -----------------------------------------------------------

print_results(); 

if (export_name <> "")
  export_results(); 
end