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

function [varargout] = CL_op_locTime(cjd,type_par,par,type_par_out, ut1_tref,tt_tref)
// Mean/true local time, longitude, right ascension (ECI/ECF, Earth only)
//
// Calling Sequence
// [par_out] = CL_op_locTime(cjd, type_par, par, type_par_out [,ut1_tref,tt_tref])
// [par_out1,..,par_outN] = CL_op_locTime(cjd, type_par, par, [type1,..,typeN] [,ut1_tref,tt_tref])
// [par_out] = CL_op_locTime(cjd, type_par, par, 'all' [,ut1_tref,tt_tref])
//
// Description
// <itemizedlist><listitem>
// <p>Performs various conversions between quantities defined in the Earth equatorial plane. </p>
// <p>These quantities are: longitude, right ascension, true or mean local time. 
// Additional quantities are computed: sidereal time and right ascension of (true or mean) Sun direction.</p> 
// <p>Quick definitions:</p>
// <p>- The right ascension (<b>ra</b>) is the angle in 
// the equator from the inertial reference axis (in the equatorial plane) chosen.</p>
// <p>- The longitude (<b>lon</b>) is the angle in the 
// equator from the Greenwich meridian.</p>
// <p>- The sidereal time (<b>sidt</b>) is the angle 
// between the inertial axis and the Greenwich meridian.</p>
// <p>- The true solar angle (<b>tsa</b>) is the angle 
// from the projection of the Sun direction onto the equatorial plane.</p>
// <p>- The mean solar angle (<b>msa</b>) is the "mean" 
// value of the true solar angle, that is the angle from a fictitious "mean" direction 
// of the Sun that is rotating uniformly in the equatorial plane (i.e. unaffected by the equation of time). </p>
// <p></p></listitem>
// <listitem>
// <p>The relations between the different parameters are as follows:</p>
// <p>ra = lon + sidt</p>
// <p>tsa = ra - stra</p>
// <p>msa = ra - smra</p>
// <p>smra = sidt - 2*pi*decimal(t) + pi</p>
// <p>tlh = tsa * 12/pi + 12</p>
// <p>mlh = msa * 12/pi + 12</p>
// <p>where: </p> 
// <p>- sidt: sidereal time. </p>
// <p>- stra: right ascension of the "true" Sun direction. </p>
// <p>- smra: right ascension of the "mean" Sun direction. </p>
// <p>- t is the number of days from some origin at 00h00 (UT); </p>
// <p>- decimal(t) = t - floor(t) (decimal part). </p>
// <p> </p>
// <p></p></listitem>
// <listitem>
// <p>The type and number of output arguments depend on "type_par_out":</p>
// <p>- If type_par_out is a (row) vector of strings, the number of output arguments must be equal 
// to the size of type_par_out (and they are returned in the requested order). </p>
// <p>- If type_par_out is 'all', only one output argument is returned in a structure. The field names are: 
// ra, lon, tsa, msa, tlh, mlh, sidt, stra, smra </p>
// <p></p></listitem>
// <listitem>
// <p><b>Important notes: </b></p> 
// <p>- The function is only valid for the Earth. </p> 
// <p>- The inertial axis (origin for right ascensions or sidereal time) is the 
// x-axis of the ECI reference frame. (See <link linkend="Reference frames">Reference frames</link> for more details) </p> 
// <p> </p> 
// <p>Remark: Some output parameters do not depend on all the inputs. Examples: </p> 
// <p>lon = CL_op_locTime(cjd, 'lon', 0, 'lon') => cjd should not be necessary (but it is used in practice!)  </p> 
// <p>sidt = CL_op_locTime(cjd, 'lon', 0, 'sidt') => the result does not depend on 'lon'!)  </p> 
// </listitem>
// </itemizedlist>
//
// Parameters
// cjd: Modified (1950.0) Julian date (Time scale: TREF) (1xN or 1x1)
// type_par: (string) Name of the input parameter: 'lon', 'ra', 'tsa', 'tlh', 'msa' or 'mlh' (1x1)
// par: Value of input parameter. Radians for lon, ra, msa, tsa and hours for mlh and tlh. (1xN or 1x1).  
// type_par_out: (string) Name(s) of output parameter(s): 'lon', 'ra', 'tsa', 'tlh', 'msa', 'mlh', 'sidt', 'stra', 'smra', or a vector of any of these, or 'all'. (1xM)  
// ut1_tref: (optional) UT1-TREF [seconds]. Default is %CL_UT1_TREF (1xN or 1x1)
// tt_tref: (optional) TT-TREF [seconds]. Default is %CL_TT_TREF. (1xN or 1x1)   
// par_out: Value(s) of output parameter(s). See description for more details. 
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_mod_siderealTime
// CL_eph_sun
//
// Examples
// // Right ascension to local time:
// cjd = CL_dat_cal2cjd(2009,09,07,10,29,30.5); // September 7th 2009, 10:29:30.5 TREF
// ra = CL_deg2rad(10);
// mlh = CL_op_locTime(cjd, 'ra', ra, 'mlh') // => mean local time in hours
// tlh = CL_op_locTime(cjd, 'ra', ra, 'tlh') // => true local time in hours
//
// // Local time to right ascension:
// ra = CL_op_locTime(cjd, 'mlh', mlh, 'ra') // using mean local time in hours
// ra = CL_op_locTime(cjd, 'tlh', tlh, 'ra') // using true local time in hours
//
// // Vector of right ascensions:
// ra = CL_deg2rad(0:90:390);
// mlh = CL_op_locTime(cjd, 'ra', ra, 'mlh')
// tlh = CL_op_locTime(cjd, 'ra', ra, 'tlh')
//
// // List of output parameters:
// [tlh,mlh,sidt] = CL_op_locTime(cjd,'ra',ra,['tlh','mlh','sidt'])
//
// // All output parameters in a structure: 
// par_out = CL_op_locTime(cjd, 'ra', ra, 'all');
// par_out.tlh


// Declarations:

// Code:
if (argn(2) < 4)
  CL__error("Not enough input arguments");
end

if (~exists("ut1_tref", "local")); ut1_tref = CL__dataGetEnv("UT1_TREF"); end
if (~exists("tt_tref", "local")); tt_tref = CL__dataGetEnv("TT_TREF"); end

if (typeof(type_par) <> "string" | typeof(type_par_out) <> "string")
   CL__error("String expected for type_par and type_par_out"); 
end

if (size(type_par,"*") <> 1); CL__error("Wrong dimension for type_par"); end;

// Check/resize input arguments
// NB: ut1_tref and tt_tref checked but not resized !
[cjd, par] = CL__checkInputs(cjd,1, par,1, ut1_tref,1, tt_tref,1);

// Number of output arguments (depends on type_par_out)
if ( (argn(1) <> 1 & type_par_out == "all") | size(type_par_out,"*") <> argn(1))
   CL__error("Wrong number of output arguments"); 
end

// Validity of strings
type_par_names = ["ra", "lon", "tsa", "msa", "tlh", "mlh"];
if (find(type_par == type_par_names) == []); CL__error("Unrecognized string in type_par"); end;

type_par_out_names = ["ra", "lon", "tsa", "msa", "tlh", "mlh", "sidt", "stra", "smra", "all"];
if (setdiff(type_par_out,type_par_out_names) <> []); CL__error("Unrecognized string in type_par_out"); end;


// -------------------------------------------------------------
res = struct('ra', 0, 'lon', 0, 'tsa', 0, 'msa', 0, ...
             'tlh', 0, 'mlh', 0, 'sidt', 0, 'stra', 0, 'smra', 0);

// Sidereal time (ECI to ECF)
cjd_ut1 = cjd + ut1_tref / 86400;
res.sidt = CL_mod_siderealTime(cjd, ut1_tref);

// Sun cartesian coordinates in ECI frame
pos_sun = CL_eph_sun(cjd, tt_tref);

// Sun right ascension in ECI frame
pos_sun_sph = CL_co_car2sph(pos_sun); 
res.stra = pos_sun_sph(1,:); 

// Right ascension of "mean" Sun in G50
res.smra = res.sidt + %pi - (cjd_ut1 - floor(cjd_ut1))*2*%pi;  // function of ut1 !! 

// -------------------------------------------------------------
// Computes right ascension and derive all the other quantities. 

select type_par
 case "ra"  ; res.ra = par; 
 case "lon" ; res.ra = res.sidt + par;  
 case "tsa" ; res.ra = res.stra + par; 
 case "msa" ; res.ra = res.smra + par; 
 case "tlh" ; res.ra = res.stra + (par - 12) * %pi/12; 
 case "mlh" ; res.ra = res.smra + (par - 12) * %pi/12; 
end

// Computes everything from ra
res.lon  = pmodulo(res.ra - res.sidt, 2*%pi); 
res.tsa  = pmodulo(res.ra - res.stra, 2*%pi); 
res.msa  = pmodulo(res.ra - res.smra, 2*%pi); 
res.tlh  = pmodulo(12 + res.tsa * 12/%pi, 24); 
res.mlh  = pmodulo(12 + res.msa * 12/%pi, 24); 

// rest: apply modulo
res.ra   = pmodulo(res.ra, 2*%pi);
res.sidt = pmodulo(res.sidt, 2*%pi);
res.stra = pmodulo(res.stra, 2*%pi);
res.smra = pmodulo(res.smra, 2*%pi);

// Outputs
if (type_par_out == "all")
  varargout(1) = res; 

else 
  for k = 1:size(type_par_out, "*")
     varargout(k) = res(type_par_out(k));
  end
end

endfunction