//  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_locTimeG50(cjd,type_par,par,type_par_out, ut1_utc)
// Conversion between longitude, right ascension and local time - DEPRECATED
//
// Calling Sequence
// [par_out] = CL_op_locTimeG50(cjd, type_par, par, type_par_out [,ut1_utc])
// [res1,..,resN] = CL_op_locTimeG50(cjd, type_par, par, [type1,..,typeN] [,ut1_utc])
// [par_out] = CL_op_locTimeG50(cjd, type_par, par, 'all' [,ut1_utc])
//
// Description
// <itemizedlist><listitem>
// <p>This function is deprecated. </p>
// <p>Replacement function: <link linkend="CL_op_locTime">CL_op_locTime</link></p>
// <p></p></listitem>
// <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 direction. </p>
// <p>- smra: right ascension of the "mean" planet to Sun direction. </p>
// <p>- t is the number of days from some origin at midnight (UT); </p>
// <p>- decimal(t) = t - floor(t) (decimal part). </p>
// <p> </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 Gamma 50 (Veis) reference frame. </p> 
// <p> </p> 
// <p>Remark: Some output parameters do not depend on all the inputs. Examples: </p> 
// <p>lon = CL_op_locTimeG50(cjd, 'lon', 0, 'lon') => cjd should not be needed in theory(but it is used in practice!)  </p> 
// <p>sidt = CL_op_locTimeG50(cjd, 'lon', 0, 'sidt') => the result does not depend on 'lon'!)  </p> 
// <p></p></listitem>
// <listitem>
// <p><b>Deprecated use:</b></p> 
// <p>- The use of tlt and mlt is deprecated. </p>
// <p>- Use respectively tsa and msa instead. </p> 
// </listitem>
// </itemizedlist>
//
// Parameters
// cjd: Time from 1950.0 in the UTC time scale (1xN or 1x1)
// type_par: (string) Name of the input parameter. It can be 'lon' (longitude), 'ra' (right ascension), 'tsa' (true solar angle), 'tlh' (true local time expressed in hours), 'msa' (mean solar angle), 'mlh' (mean local time expressed in hours). 
// par: Value of the input parameter expressed in degrees for lon, ra, msa, tsa and hours for mlh and tlh. (1xN or 1x1).  
// type_par_out: (string) Name(s) of the output parameter(s). It can be 'lon', 'ra', 'tsa', 'tlh', 'msa', 'mlh' (see type_par), or 'sidt' (sidereal time), 'stra' (Sun true right ascension) 'smra' (Sun mean right ascension), or a vector of any of these, or 'all'. (1xM)  
// ut1_utc: Offset between UT1 and UTC (s) = Time of some event in UT1 minus time of the same event in UTC.   
// varargout: Value(s) of the output parameters. As many arguments are needed as there are names in 'type_par_out', except if type_par_out is 'all' in which case varargout is a structure (with fields names equal to the names in type_par_out). Units are radians or hours depending on the parameter. 
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_mod_sidTimeG50
// CL_mod_moonSunG50
//
// Examples
// // Right ascension to local time:
// t = CL_dat_cal2cjd(2009,09,07,10,29,30.5); // September 7th 2009, 10:29:30.5 UTC
// ra   = CL_deg2rad(10);
// mlh = CL_op_locTimeG50(t, 'ra', ra, 'mlh') // => mean local time in hours
// tlh = CL_op_locTimeG50(t, 'ra', ra, 'tlh') // => true local time in hours
//
// // Local time to right ascension:
// ra = CL_op_locTimeG50(t, 'mlh', mlh, 'ra') // using mean local time in hours
// ra = CL_op_locTimeG50(t, 'tlh', tlh, 'ra') // using true local time in hours
//
// // Local time to longitude:
// lon = CL_op_locTimeG50(t, 'mlh', mlh, 'lon') // using mean local time in hours
// lon = CL_op_locTimeG50(t, 'tlh', tlh, 'lon') // using true local time in hours
//
// // Vector of right ascensions :
// ra = CL_deg2rad(0:90:390);
// mlh = CL_op_locTimeG50(t, 'ra', ra, 'mlh')
// tlh = CL_op_locTimeG50(t, 'ra', ra, 'tlh')
//
// // List of output parameters :
// [tlh,mlh,sidt] = CL_op_locTimeG50(t,'ra',ra,['tlh','mlh','sidt'])
//
// // All output parameters in a structure: 
// par_out = CL_op_locTimeG50(t, 'ra', ra, 'all');
// par_out.tlh


// Declarations:

// Code:

CL__warnDeprecated(); // deprecated function

[lhs,rhs] = argn(); // number of input (rhs) / output (lhs) arguments

if (~exists('ut1_utc','local'))
   ut1_utc = 0;   // OK if scalar! 
else 
   rhs = rhs - 1; // ut1_utc not counted in what follows. 
end

// error checking
if (rhs <> 4) 
   CL__error('Wrong number of input arguments'); 
end

// only 1 row accepted
if (size(cjd,1) <> 1 | size(type_par,1) <> 1 | size(par,1) <> 1 | size(type_par_out,1) <> 1)
   CL__error('Wrong dimension(s) of arguments'); 
end

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

if ~( (lhs == 1 & type_par_out == 'all') | size(type_par_out,2) == lhs)
   CL__error('Wrong number of output arguments'); 
end

Ncjd = size(cjd,2);
Npar = size(par,2);
N = max(Ncjd, Npar);

if ~((Ncjd == 1 | Ncjd == N) & (Npar == 1 | Npar == N))
   CL__error('Wrong size of input arguments'); 
end

// force same sizes for input arguments
cjd = cjd .* ones(1,N);
par = par .* ones(1,N); 

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

// calculation of sidereal time
cjd_ut1 = cjd + ut1_utc/86400.; // time in ut1
res.sidt = CL_mod_sidTimeG50(cjd_ut1);

// Sun cartesian coordinates in G50
[u_sun, r_sun] = CL_mod_moonSunG50(cjd,'s'); // u_sun : unit vector, r_sun: norm

// transformation from cartesian to spherical coordinates 
// (=> right ascension, declination, radius)
u_sun_sph = CL_co_car2sph(u_sun); 
res.stra = u_sun_sph(1,:); // Sun right ascension

// 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; 
 else CL__error('Unknown input parameter name');
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);

// for compatibility ...
res.tlt = res.tsa; 
res.mlt = res.msa; 

// outputs
if (type_par_out == 'all')
  varargout(1) = res; 

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

endfunction