//  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 [ra,rp,ha,hp,va,vp,nmoy,period,vinf,anvinf,dinf,vlibp,vliba,l_con,c_con,k,caire,smb,anm,anex,temps,r,x,y,pen,vit] = CL__kp_keplerParams(sma,ecc,v,orb_type,er,mu)


// Declarations:


// Code:

if find(orb_type~=1 & orb_type~=2)~=[] then CL__error('invalid orbit type'); end

//initialisation
ra = %nan * ones(sma)
rp = %nan * ones(sma)
ha = %nan * ones(sma)
va = %nan * ones(sma)
vp = %nan * ones(sma)
period = %nan * ones(sma)
vinf = %nan * ones(sma)
anvinf = %nan * ones(sma)
dinf = %nan * ones(sma)
vlibp = %nan * ones(sma)
vliba = %nan * ones(sma)
l_con = %nan * ones(sma)
k = %nan * ones(sma)
caire = %nan * ones(sma)
smb = %nan * ones(sma)
vit = %nan * ones(sma)

//calculations

//hyperbolic parameters
ih = find(orb_type==2)
  //periapsis radius
  rp(ih) = sma(ih).*(ecc(ih)-1)
  //velocity at periapsis
  vp(ih) = sqrt( mu*( (1.0 ./ sma(ih)) + (2.0 ./ (sma(ih).*(ecc(ih)-1)))) )
  //velocity at infinity
  vinf(ih) = sqrt(mu./sma(ih))
  //true anomaly at infinity
  anvinf(ih) = acos(-1.0 ./ ecc(ih))
  //distance from asymptote to focus
  dinf(ih) = sma(ih).*ecc(ih).*sin(anvinf(ih))
  //parameter l of conic
  l_con(ih) = sma(ih).*(ecc(ih).^2 - 1)
  //energy
  k(ih) = mu./(2.0.*sma(ih))
  //area rule constant
  caire(ih) = sqrt( mu.*sma(ih).*(ecc(ih).^2 - 1))
  //semi minor axis
  smb(ih) = sma(ih).*sqrt(ecc(ih).^2 - 1)

//elliptic parameters
ie = find(orb_type==1)
  //apoapsis radius
  ra(ie) = sma(ie).*(1+ecc(ie))
  //apoapsis altitude
  ha(ie) = ra(ie) - er
  //velocity at apoapsis
  va(ie) = sqrt( (mu./sma(ie)).*((1-ecc(ie))./(1+ecc(ie))) )
  //radius at periapsis
  rp(ie) = sma(ie).*(1-ecc(ie))
  //velocity at periapsis
  vp(ie) = sqrt( (mu./sma(ie)).*((1+ecc(ie))./(1-ecc(ie))) )
  //period
  period(ie) = CL_kp_params('per',sma(ie),mu)
  //liberation velocities
  vlibp(ie) = sqrt( (2*mu)./rp(ie))
  vliba(ie) = sqrt( (2*mu)./ra(ie))
  //parameter l of conic
  l_con(ie) = sma(ie).*(1-ecc(ie).^2)
  //energy
  k(ie) = -mu./(2.0.*sma(ie))
  //area rule constant
  caire(ie) = sqrt( mu.*sma(ie).*(1-ecc(ie).^2))
  //semi minor axis
  smb(ie) = sma(ie).*sqrt(1-ecc(ie).^2)

//common parameters to both type of orbits (elliptic and hyperbolic)
//periapsis altitude
hp = rp-er
//mean motion
nmoy  = CL_kp_params('mm',sma,mu)
//parameter c of conic
c_con = sma.*ecc

//mean and eccentric anomaly
anex = CL_kp_v2E(ecc,v)
anm = CL_kp_E2M(ecc,anex)

//temps
temps = anm./nmoy
//r vector
r = l_con./(1+ecc.*cos(v))
//x cartesian coordinate
x = r.*cos(v)
//y cartesian coordinate
y = r.*sin(v)
//slope of velocity
pen = atan( ecc.*sin(v) , 1+ecc.*cos(v) )
//modulus of velocity for elliptic orbits
vit(ie) = sqrt(mu .* (2.0 ./ r(ie) - 1.0 ./ sma(ie)))
//modulus of velocity for hyperbolic orbits
vit(ih) = sqrt(mu .* (2.0 ./ r(ih) + 1.0 ./ sma(ih)))

endfunction