################################################################################
#
# Orbital library for gnuplot
#   - Version     :  3.0 beta
      HSG_libver  = "3.0 beta";
#   - Author      : Helton da Silva Gaspar
#   - Mail        : helton.s.gaspar@ufsc.br
#   - Date        : 04 / 08 / 2023
      HSG_libdate = "2023-08-04";
#
################################################################################

################################################################################
# 
# DESCRIPTION
# 
#   This library gathers useful functions for plotting orbits. It contains 
#   four main packages:
#     1. rv2** whose functions convert fixed cartesian coordinates to 
#              keplerian elements 
#     2. eo2** whose functions convert keplerian elements to fixed cartesian 
#              coordinates
#     3. polar and perifocal frame functions
#     4. time/position functions
#     5. constants
#
# DEFINITIONS AND UNIFORMITY
#    p     stands for the semilatus rectum
#    a     stands for the semimajor axis ( a==p case e==1 )
#    e     stands for the eccentricity
#    I     stands for the inclination
#    O     stands for the longitude of the ascending node
#    w     stands for the argument of periapsis
#    f     stands for the true anomaly
#    M     stands for the mean anomaly
#    E     stands for the eccetric anomaly
#    n     stands for the mean motion
#    T     stands for the period
#    t     stands for the flight time since periapsis
#    C3    stands for the specific mechanical energy
#    mu    stands for the mass parameter (G*mass)
#    gamma stands of the flight path angle
#
# FUNCTIONS WITHIN THE rv2** PACKAGE
#   In order to keep uniformity, all the package functions must have the 
#   following parameters as argument: mu, x, y, z, vx, vy and vz
#     rv2**(mu,x,y,z,vx,vy,vz)
#   despite not all of them are necessary for some of the functions.
#   
#   SYNOPSIS rv2**(mu,x,y,z,vx,vy,vz), where ** may be one of the following:
#               ** =  p,a,e,I,O,w,f,n,T,C3,gamma
#
# FUNCTIONS WITHIN THE eo2** PACKAGE
#   In order to keep uniformity, all the package functions must have the 
#   following parameters as argument: mu, a, e, I, O, w and  f.
#     rv2**(mu,a,e,I,O,w,f)
#   despite not all of them are necessary for some of the functions.
#   
#   SYNOPSIS eo2**(mu,a,e,I,O,w,f), where ** may be one of the following:
#               ** = x,y,z,vx,vy,vz,n,T,t
#   
# FUNCTIONS WITHIN THE POLAR PACKAGE
#   - polp(a,e)      returns semilatus rectum     given a and e
#   - polr(a,e,f)    returns range from focus     given a, e and f
#   - poldrdf(a,e,f) returns the derivative dr/dr given a, e and f
#   
# FUNCTIONS OF THE PERIFOCAL FRAME PACKAGE
#   - x_pq(a,e,w,f)     returns the x coordinate within the perifocal frame
#   - y_pq(a,e,w,f)     returns the y coordinate within the perifocal frame
#   - vr(mu,a,e,f)      returns the radial      component of the velocity
#   - vt(mu,a,e,f)      returns the transversal component of the velocity
#   - vx_pq(mu,a,e,w,f) returns the x           component of the velocity
#   - vy_pq(mu,a,e,w,f) returns the y           component of the velocity
#   - gamma(e,f)        returns the flight path angle
#   
# FUNCTIONS WITHIN THE TIME/POSITION PACKAGE
#   - f2E(e,f)    returns eccetric anomaly            given e and f
#   - E2M(e,E)    returns mean     anomaly            given e and E
#   - f2M(e,f)    returns mean     anomaly            given e and f
#   - M2E(e,M)    returns eccetric anomaly            given e and M
#   - E2f(e,E)    returns true     anomaly            given e and E
#   - M2f(e,M)    returns true     anomay             given e and M
#   - t2E(n,e,t)  returns eccetric anomaly            given n, e and t
#   - t2f(n,e,t)  returns true     anomaly            given n, e and t
#   - M2t(n,M)    returns flight time since periapsis given n and M
#   - f2t(n,e,f)  returns flight time since periapsis given n, e and f
#   - t2M(n,t)    returns flight time since periapsis given n and t
# 
################################################################################



################################################################################
#
## Package of constants

  liborb_nrmax = 100;                       #maximum number of iterations for Newton Raphson algorithm
  liborb_nrdE  = 0.0;                       #residual value from Newton Raphson iterations
  liborb_zero  = 1.000000000000000000e-08;  #limit value to consider somme quantities as null
  liborb_inf   = 1.23456789876543210e1000;  #this value stands for infinity value
  liborb_deg   = 1.745329251994329547e-02;  #1 degree in radians
  liborb_gon   = 1.570796326794896697e-02;  #1 gradian in radians
  liborb_hpi   = 1.570796326794896558e+00;  #half pi
  liborb_dpi   = 6.283185307179586232e+00;  #double pi
  liborb_au    = 1.495978707000000000e+11;  #m
  liborb_G     = 6.674080000000000335e-11;  #m^3/kg s^2
  liborb_musun = 1.327541252800000000e+20;  #m^3/s^2
  liborb_mumer = 2.203900081014675781e+13;  #m^3/s^2
  liborb_muven = 3.249606376090141250e+14;  #m^3/s^2
  liborb_muear = 3.987256375298450625e+14;  #m^3/s^2
  liborb_mulua = 4.904338898288770508e+12;  #m^3/s^2
  liborb_mumar = 4.284176672342150781e+13;  #m^3/s^2
  liborb_mujup = 1.267525686099165120e+17;  #m^3/s^2
  liborb_musat = 3.795254329314348000e+16;  #m^3/s^2
  liborb_muura = 5.796369087341472000e+15;  #m^3/s^2
  liborb_munep = 6.838681433981858000e+15;  #m^3/s^2
  liborb_ioang = 1.0;                       #input and output angles (1.0 for radians | liborb_deg for degrees)
  liborb_curang= "rad";




################################################################################
#
## math sub-pack

  #Magnitude and squared magnitude of vector U
  liborb_mag2(Ux,Uy,Uz)      = Ux*Ux + Uy*Uy + Uz*Uz;
  liborb_mag (Ux,Uy,Uz)      = sqrt( liborb_mag2(Ux,Uy,Uz) );

  #Wrap around x within the domain xmin <= x <= xmax
  liborb_wrap(x,xmin,xmax)   = x<xmin\
                             ? liborb_wrap(x+xmax-xmin,xmin,xmax)\
                             : x>xmax\
                               ? liborb_wrap(x-xmax+xmin,xmin,xmax)\
                               : x;

  # linear interpolation of the points (0,a) (1,b)
  lin(a,b,u)                 = (1-u)*a+u*b;

  # Euler angles rotations
  liborb_Rx(x,y,z,a,b,c)     = x*( cos(a*liborb_ioang)*cos(c*liborb_ioang)-sin(a*liborb_ioang)*cos(b*liborb_ioang)*sin(c*liborb_ioang) ) - y*( cos(a*liborb_ioang)*sin(c*liborb_ioang) + sin(a*liborb_ioang)*cos(b*liborb_ioang)*cos(c*liborb_ioang) ) + z*sin(a*liborb_ioang)*sin(b*liborb_ioang);
  liborb_Ry(x,y,z,a,b,c)     = x*( sin(a*liborb_ioang)*cos(c*liborb_ioang)+cos(a*liborb_ioang)*cos(b*liborb_ioang)*sin(c*liborb_ioang) ) - y*( sin(a*liborb_ioang)*sin(c*liborb_ioang) - cos(a*liborb_ioang)*cos(b*liborb_ioang)*cos(c*liborb_ioang) ) - z*cos(a*liborb_ioang)*sin(b*liborb_ioang);
  liborb_Rz(x,y,z,a,b,c)     = x*  sin(b*liborb_ioang)*sin(c*liborb_ioang)                                                               + y*  sin(b*liborb_ioang)*cos(c*liborb_ioang)                                                                 + z*cos(b*liborb_ioang);



################################################################################
#
## conversion from "fixed" cartesian coordinates to keplerian elements

  #vector specific angular momentum
  liborb_hx(x,y,z,vx,vy,vz) = y*vz-z*vy;
  liborb_hy(x,y,z,vx,vy,vz) = z*vx-x*vz;
  liborb_hz(x,y,z,vx,vy,vz) = x*vy-y*vx;
  liborb_h2(x,y,z,vx,vy,vz) = liborb_mag2( y*vz-z*vy, z*vx-x*vz, x*vy-y*vx );

  #eccentricity vector
  liborb_ex(mu,x,y,z,vx,vy,vz) = x*(liborb_mag2(vx,vy,vz)*1.0/mu-1.0/liborb_mag(x,y,z)) - vx*(x*vx+y*vy+z*vz)*1.0/mu;
  liborb_ey(mu,x,y,z,vx,vy,vz) = y*(liborb_mag2(vx,vy,vz)*1.0/mu-1.0/liborb_mag(x,y,z)) - vy*(x*vx+y*vy+z*vz)*1.0/mu;
  liborb_ez(mu,x,y,z,vx,vy,vz) = z*(liborb_mag2(vx,vy,vz)*1.0/mu-1.0/liborb_mag(x,y,z)) - vz*(x*vx+y*vy+z*vz)*1.0/mu;

  liborb_rv2C3(mu,x,y,z,vx,vy,vz) = 0.5*liborb_mag2(vx,vy,vz) - mu*1.0/liborb_mag(x,y,z);
  liborb_rv2p(mu,x,y,z,vx,vy,vz) = liborb_h2(x,y,z,vx,vy,vz)*1.0/mu;
  liborb_rv2a(mu,x,y,z,vx,vy,vz) = ( aux1=liborb_rv2C3(mu,x,y,z,vx,vy,vz) )**2 < liborb_zero**2\
                                 ? liborb_rv2p(mu,x,y,z,vx,vy,vz)\
                                 : abs(0.5*mu/aux1);

  liborb_rv2e(mu,x,y,z,vx,vy,vz) = liborb_rv2C3(mu,x,y,z,vx,vy,vz) == 0.0\
                                 ? 1.0\
                                 : 0*( aux1=(liborb_mag2(vx,vy,vz)*1.0/mu-1.0/liborb_mag(x,y,z)) )\
                                    *( aux2=(x*vx+y*vy+z*vz)*1.0/mu )\
                                 + liborb_mag( x*aux1-vx*aux2,\
                                               y*aux1-vy*aux2,\
                                               z*aux1-vz*aux2);

  liborb_rv2I(mu,x,y,z,vx,vy,vz) = 0*( aux1=liborb_mag(liborb_hx(x,y,z,vx,vy,vz),liborb_hy(x,y,z,vx,vy,vz),0.0) )\
                                   + atan2(aux1,liborb_hz(x,y,z,vx,vy,vz))/liborb_ioang;
  liborb_rv2O(mu,x,y,z,vx,vy,vz) = (abs(z)/liborb_mag(x,y,z)+abs(vz)/liborb_mag(vx,vy,vz))<liborb_zero\
                                 ? 0.0\
                                 : liborb_wrap( atan2(+liborb_hx(x,y,z,vx,vy,vz),-liborb_hy(x,y,z,vx,vy,vz)),0.0,liborb_dpi )/liborb_ioang;
  liborb_rv2f(mu,x,y,z,vx,vy,vz) = (aux1=liborb_rv2e(mu,x,y,z,vx,vy,vz))>liborb_zero\
                                 ? acos( (liborb_rv2p(mu,x,y,z,vx,vy,vz)/liborb_mag(x,y,z)-1.0)/aux1)*((x*vx+y*vy+z*vz)<0.0?-1:+1)/liborb_ioang\
                                 : 0*(aux2=(abs(z)/liborb_mag(x,y,z)+abs(vz)/liborb_mag(vx,vy,vz)))\
                                    *(aux3=aux2<liborb_zero?1.0:-liborb_hy(x,y,z,vx,vy,vz))\
                                    *(aux4=aux2<liborb_zero?0.0:+liborb_hx(x,y,z,vx,vy,vz))\
                                    *(aux5=-liborb_hz(x,y,z,vx,vy,vz)*aux4)\
                                    *(aux6=+liborb_hz(x,y,z,vx,vy,vz)*aux3)\
                                    *(aux7=liborb_mag2(liborb_hx(x,y,z,vx,vy,vz),liborb_hy(x,y,z,vx,vy,vz),0.0))\
                                 +  atan2((x*aux5+y*aux6+z*aux7)*1.0/liborb_mag(aux5,aux6,aux7),(x*aux3+y*aux4)*1.0/liborb_mag(aux3,aux4,0.0))/liborb_ioang;
  liborb_rv2w(mu,x,y,z,vx,vy,vz) = (aux1=liborb_rv2e(mu,x,y,z,vx,vy,vz))==0\
                                 ? 0.0\
                                 : 0*(aux2=(abs(z)/liborb_mag(x,y,z)+abs(vz)/liborb_mag(vx,vy,vz)))\
                                    *(aux3=aux2==0?1.0:-liborb_hy(x,y,z,vx,vy,vz))\
                                    *(aux4=aux2==0?0.0:+liborb_hx(x,y,z,vx,vy,vz))\
                                    *(aux5=-liborb_hz(x,y,z,vx,vy,vz)*aux4)\
                                    *(aux6=+liborb_hz(x,y,z,vx,vy,vz)*aux3)\
                                    *(aux7=+liborb_hx(x,y,z,vx,vy,vz)*aux4-liborb_hy(x,y,z,vx,vy,vz)*aux3)\
                                 + liborb_wrap( atan2( liborb_ex(mu,x,y,z,vx,vy,vz)*aux5 + liborb_ey(mu,x,y,z,vx,vy,vz)*aux6 + liborb_ez(mu,x,y,z,vx,vy,vz)*aux7\
                                                     , liborb_ex(mu,x,y,z,vx,vy,vz)*aux3 + liborb_ey(mu,x,y,z,vx,vy,vz)*aux4)\
                                              ,0.0,liborb_dpi)/liborb_ioang;
  liborb_rv2n(mu,x,y,z,vx,vy,vz) = 0*(aux1=liborb_rv2e(mu,x,y,z,vx,vy,vz))\
                                    *(aux2=aux1==1.0?1.0:abs(aux1**2-1))\
                                 +  mu**2*(aux2/liborb_h2(x,y,z,vx,vy,vz))**(1.5)/liborb_ioang;
  liborb_rv2T(mu,x,y,z,vx,vy,vz) = liborb_dpi/liborb_rv2n(mu,x,y,z,vx,vy,vz)*liborb_ioang;
  liborb_rv2gamma(mu,x,y,z,vx,vy,vz) = asin((x*vx+y*vy+z*vz)*1.0/liborb_mag(x,y,z)/liborb_mag(vx,vy,vz))/liborb_ioang;



################################################################################
#
## polar functions

  liborb_polp(a,e)      = e==1.0\
                        ? a\
                        : a*abs(1.0-e*e);
  liborb_polr(a,e,f)    = (aux1=1.0+e*cos(f*liborb_ioang))==0\
                        ? liborb_inf\
                        : liborb_polp(a,e)*1.0/aux1;
  liborb_poldrdf(a,e,f) = liborb_polp(a,e)*e*sin(f*liborb_ioang)/(1.0+e*cos(f*liborb_ioang))**2;



################################################################################
#
## perifocal frame

  liborb_xpq(a,e,w,f)     = liborb_polr(a,e,f)*cos((w+f)*liborb_ioang);
  liborb_ypq(a,e,w,f)     = liborb_polr(a,e,f)*sin((w+f)*liborb_ioang);
  liborb_vr(mu,a,e,f)     = sqrt(mu/liborb_polp(a,e))*e*sin(f*liborb_ioang);
  liborb_vt(mu,a,e,f)     = sqrt(mu/liborb_polp(a,e))*(1.0+e*cos(f*liborb_ioang));
  liborb_vxpq(mu,a,e,w,f) = liborb_vr(mu,a,e,f)*cos((w+f)*liborb_ioang) - liborb_vt(mu,a,e,f)*sin((w+f)*liborb_ioang);
  liborb_vypq(mu,a,e,w,f) = liborb_vr(mu,a,e,f)*sin((w+f)*liborb_ioang) + liborb_vt(mu,a,e,f)*cos((w+f)*liborb_ioang);
  liborb_fpgamma(e,f)     = atan2(e*sin(f*liborb_ioang),1.0+e*cos(f*liborb_ioang))/liborb_ioang;
  liborb_gamma(e,f)       = atan2(e*sin(f*liborb_ioang),1.0+e*cos(f*liborb_ioang))/liborb_ioang;



################################################################################
#
## conversion from keplerian elements to "fixed" cartesian coordinates

  liborb_eo2x (mu,a,e,I,O,w,f) = liborb_Rx( liborb_xpq(   a,e,0.0,f), liborb_ypq(   a,e,0.0,f),0.0,O,I,w);
  liborb_eo2y (mu,a,e,I,O,w,f) = liborb_Ry( liborb_xpq(   a,e,0.0,f), liborb_ypq(   a,e,0.0,f),0.0,O,I,w);
  liborb_eo2z (mu,a,e,I,O,w,f) = liborb_Rz( liborb_xpq(   a,e,0.0,f), liborb_ypq(   a,e,0.0,f),0.0,O,I,w);
  liborb_eo2vx(mu,a,e,I,O,w,f) = liborb_Rx(liborb_vxpq(mu,a,e,0.0,f),liborb_vypq(mu,a,e,0.0,f),0.0,O,I,w);
  liborb_eo2vy(mu,a,e,I,O,w,f) = liborb_Ry(liborb_vxpq(mu,a,e,0.0,f),liborb_vypq(mu,a,e,0.0,f),0.0,O,I,w);
  liborb_eo2vz(mu,a,e,I,O,w,f) = liborb_Rz(liborb_vxpq(mu,a,e,0.0,f),liborb_vypq(mu,a,e,0.0,f),0.0,O,I,w);
  liborb_eo2n (mu,a,e,I,O,w,f) = sqrt( mu/a**3 )/liborb_ioang;
  liborb_eo2T (mu,a,e,I,O,w,f) = liborb_dpi/liborb_eo2n(mu,a,e,I,O,w,f);



################################################################################
#
## time and position functions

  liborb_f2E(e,f)              = e>1.0\
                               ? 2.0*atanh( sqrt((e-1.0)/(e+1.0))*tan(0.5*f*liborb_ioang) )/liborb_ioang\
                               : 2.0*atan2( sqrt(1.0-e)*sin(0.5*f*liborb_ioang),sqrt(1.0+e)*cos(0.5*f*liborb_ioang) )/liborb_ioang;
  liborb_E2M(e,E)              = e>1.0\
                               ? (e*sinh(E*liborb_ioang)-E*liborb_ioang)/liborb_ioang\
                               : (E*liborb_ioang-e*sin(E*liborb_ioang))/liborb_ioang;
  liborb_f2M(e,f)              = e==1.0\
                               ? 0.5*(aux1=tan(0.5*f*liborb_ioang))*( 1.0 + 0.33333333*aux1*aux1 )/liborb_ioang\
                               : liborb_E2M(e,liborb_f2E(e,f));
  liborb_M2t(n,M)              = M*1.0/n;
  liborb_f2t(n,e,f)            = liborb_M2t(n,liborb_f2M(e,f));
  liborb_eo2t(mu,a,e,I,O,w,f)  = liborb_M2t(liborb_eo2n(mu,a,e,I,O,w,f),liborb_f2M(e,f));

  liborb_nr(e,M,E0,itr)        = 0*(dE=e<1.0\
                               ? (E0-e*sin(E0)-M)/(1.0-e*cos(E0))\
                               : (e*sinh(E0)-E0-M)/(e*cosh(E0)-1.0))\
                               + itr>0\
                                 ? dE**2>liborb_zero**2\
                                   ? liborb_nr(e,M,E0-dE,itr-1)\
                                   : E0-dE\
                                 : 0*(liborb_nrdE=dE)\
                                   + E0-dE;

  liborb_t2M(n,t)              = n*t;
  liborb_M2E(e,M)              = e==1.0\
                               ? 1/0\
                               : liborb_nr(e,M*liborb_ioang,M*liborb_ioang,liborb_nrmax)/liborb_ioang;
  M2Err(e,M)                   = e==1.0\
                               ? 1/0\
                               : (M*liborb_ioang-liborb_E2M(e,liborb_nr(e,M*liborb_ioang,M*liborb_ioang,liborb_nrmax)))/liborb_ioang;
  liborb_E2f(e,E)              = e>1.0\
                               ? 2.0*atan2( sqrt(1.0+e)*sinh(0.5*E*liborb_ioang) , sqrt(e-1.0)*cosh(0.5*E*liborb_ioang) )/liborb_ioang\
                               : 2.0*atan2( sqrt(1.0+e)*sin( 0.5*E*liborb_ioang) , sqrt(1.0-e)*cos( 0.5*E*liborb_ioang) )/liborb_ioang;
  liborb_M2f(e,M)              = e==1.0\
                               ? 2.0*atan( (aux1=3.0*M*liborb_ioang+sqrt(9.0*M*M*liborb_ioang*liborb_ioang+1))**(1.0/3)-aux1**(-1.0/3) )/liborb_ioang\
                               : liborb_E2f(e,liborb_M2E(e,M));
  liborb_t2E(n,e,t)            = liborb_M2E(e,liborb_t2M(n,t));
  liborb_t2f(n,e,t)            = liborb_M2f(e,liborb_t2M(n,t));



################################################################################
#
## Help "functions"

  more_orbital_anomalies=sprintf("\
  Anomalies function are:\n\
      f2E(e,f)    returns eccentric anomaly            given e and f\n\
      E2M(e,E)    returns mean      anomaly            given e and E\n\
      f2M(e,f)    returns mean      anomaly            given e and f\n\
      M2E(e,M)    returns eccentric anomaly            given e and M\n\
      E2f(e,E)    returns true      anomaly            given e and E\n\
      M2f(e,M)    returns true      anomay             given e and M\n\
      t2E(n,e,t)  returns eccentric anomaly            given n, e and t\n\
      t2f(n,e,t)  returns true      anomaly            given n, e and t\n\
      M2t(n,M)    returns flight time since periapsis given n and M\n\
      f2t(n,e,f)  returns flight time since periapsis given n, e and f\n\
      t2M(n,t)    returns flight time since periapsis given n and t\n\n\
  Please, type:\n\
     print more_orbital_nr\n\
  to learn more about accuracy in Newton-Raphson method\n");

  liborb_newtonraphson=sprintf("\
  Newton-Raphson iteration methods are employed to solve Kepler equations:\n\
     M = E - e*sin(E)\n\
     M = e*sinh(E) - E\n\
  for eccentric anomaly E, for eliptic and hiperbolic cases, respectively.\n\
  Convergence may fail depending upon the value of the eccentricity. The closer\n\
  to 1 is the eccentricity the harder is the convergence. Iterations are carried\n\
  out until:\n\
     | E_n - E_{n-1} | < liborb_zero (%.0e by default)\n\
  or\n\
     number of iterations = liborb_nrmax (%d by default)\n\n\
  When recursion limit liborb_nrmax is reached, the \"error\" is saved as:\n\
     liborb_nrdE = E_n - E_{n-1}\n\n\
  The library also provides the function:\n\
     M2Err(e,M)\n\
  which allows one to infer the residual error of Newton-Raphson iteration. Such error is the\n\
  difference of the resulting E converted back to M' to the original value of M.\n", liborb_zero, liborb_nrmax );



  more_orbital_definitions=sprintf("\
  Definitions and uniformity:\n\
   - p  semilatus rectum\n\
   - a  semimajor axis. In particular case parabolic orbits \"a\" stands for the\n\
        semilaltus rectum while for the circular ones, the circumference radius.\n\
   - e  eccentricity\n\
   - I  inclination of the orbital plane\n\
   - O  longitude of the ascending node\n\
   - w  argument of periapsis\n\
   - f  true anomaly\n\
   - M  mean anomaly\n\
   - E  eccetric anomaly\n\
   - n  mean motion\n\
   - T  period\n\
   - t  flight time since periapsis\n\
   - C3 specific mechanical energy\n\
   - mu mass parameter (G*mass)\n");



  more_orbital_constants=sprintf("\
  Builtin constants are:\n\
      liborb_zero   is the threshold value considered as zero within some\n\
                    functions such as Newton-Raphson iteration or the planar\n\
                    orbits constrained to |z+vz|<liborb_zero (%.0e by default)\n\
      liborb_nrmax  is the maximum number of iterations performed in Newton-Raphson\n\
                    method (%d by default).\n\
      liborb_nrdE   the resulting error of Newton-Raphson when it doesn't converge\n\
                    after %d iterations\n\
      liborb_au     is the astronomical unit (%.6e m by default)\n\
      liborb_G      is the Gravitational constant (%.6e m^3/kg s^2 by default)\n\
      liborb_muxxx  is G times mass for xxx planet where:\n\
                      xxx = (sun, mer, ven, ear, moon, mar, jup, sat, ura, nep, plu)\n\
      liborb_dpi    is (d)ouble pi\n\
      liborb_hpi    is (h)alf pi\n\
      liborb_deg    is 1 (deg)ree in radians\n\
  ",liborb_zero, liborb_nrmax, liborb_nrmax, liborb_au, liborb_G )


 liborb_true_anomaly=sprintf("  The following functions convert true anomaly into:\n\
   - f2E(e,f)   || true_anomaly_to_eccentric_anomaly(e,f)\n\
   - f2M(e,f)   || true_anomaly_to_mean_anomaly(e,f)\n\
   - f2t(n,e,f) || true_anomaly_to_time(n,e,f)");
 
 liborb_mean_anomaly=sprintf("  The following functions convert mean anomaly into:\n\
   - M2E(e,M)   || mean_anomaly_to_eccentric_anomaly(e,M)\n\
   - M2f(e,M)   || mean_anomaly_to_true_anomaly(e,M)\n\
   - M2t(n,M)   || mean_anomaly_to_time(n,M)");
  
 liborb_eccentric_anomaly=sprintf("  The following functions convert eccentric anomaly into:\n\
   - E2f(e,E)   || eccentric_anomaly_to_true_anomaly(e,E)\n\
   - E2M(e,E)   || eccentric_anomaly_to_mean_anomaly(e,E)");
 
 liborb_time=sprintf("  The following functions convert time into:\n\
   - t2M(n,t)   || time_to_mean_anomaly(n,t)\n\
   - t2f(n,e,t) || time_to_true_anomaly(n,e,t)\n\
   - t2E(n,e,t) || time_to_eccentric_anomaly(n,e,t)");
  
 liborb_anomalies_funcs=sprintf("\n#\n# Set of functions to deal with time from periapses and anomalies.\n\
%s\n\n%s\n\n%s\n\n%s",\
    liborb_true_anomaly,\
    liborb_mean_anomaly,\
    liborb_eccentric_anomaly,\
    liborb_time);
 
 liborb_anomalies=sprintf("%s\n\n\
  Functions arguments are:\n\
   - f true anomaly         - e eccentricity\n\
   - M mean anomaly         - t time since periapsis\n\
   - E eccentric anomaly    - n mean motion\n\n\
  Notice the short name functions have long name aliases!\n\
  All the angles in radians by default (see more: print liborb_unities)\n\
  You may also want to know more about Newton-Raphson iteration. If so, please\n\
     type: print liborb_newtonraphson\n", liborb_anomalies_funcs);



 liborb_perifocal_funcs=sprintf("#\n# Set of functions to plot 2D orbits within the perifocal frame (pq).\n\
  The following functions give the position coordinates:\n\
   - x_pq(a,e,w,f)\n\
   - y_pq(a,e,w,f)\n\n\
  The following functions give the velocity coordinates:\n\
   - vx_pq(mu,a,e,w,f)\n\
   - vy_pq(mu,a,e,w,f)\n\n\
  The following functions give the polar components of the velocity vector:\n\
   - vr(mu,a,e,f) #radial\n\
   - vt(mu,a,e,f) #transversal\n\n\
  The following one returns the flight path angle gamma\n\
   - liborb_gamma(e,f) || flight_path_angle(e,f)");
   
 liborb_perifocal=sprintf("%s\n\n\
  Functions arguments are:\n\
   - a  semimajor axis            - f  true anomaly\n\
   - e  eccentricity              - mu parameter of mass\n\
   - w  argument of the periapsis\n\n\
  Notice that (w) is measured counterclockwise since the x_pq axis.\n", liborb_perifocal_funcs);
  
 liborb_orbital_to_funcs=sprintf("#\n# Set of functions to deal with orbital elements\n\
  The following functions convert the orbital elements into cartesian position\n\
   - oe2x(mu,a,e,I,O,w,f)  || orbital_to_x(mu,a,e,I,O,w,f)\n\
   - oe2y(mu,a,e,I,O,w,f)  || orbital_to_y(mu,a,e,I,O,w,f)\n\
   - oe2z(mu,a,e,I,O,w,f)  || orbital_to_z(mu,a,e,I,O,w,f)\n\
     *Notice that (mu) is not necessary!\n\n\
  The following functions convert the orbital elements into cartesian velocities\n\
   - oe2vx(mu,a,e,I,O,w,f) || orbital_to_vx(mu,a,e,I,O,w,f)\n\
   - oe2vy(mu,a,e,I,O,w,f) || orbital_to_vy(mu,a,e,I,O,w,f)\n\
   - oe2vz(mu,a,e,I,O,w,f) || orbital_to_vz(mu,a,e,I,O,w,f)\n\n\
  The following functions convert the orbital elements into:\n\
   - oe2n(mu,a,e,I,O,w,f)  || orbital_to_mean_motion(mu,a,e,I,O,w,f)\n\
   - oe2T(mu,a,e,I,O,w,f)  || orbital_to_period(mu,a,e,I,O,w,f)\n\
   - oe2t(mu,a,e,I,O,w,f)  || orbital_to_time(mu,a,e,I,O,w,f)\n\
     *Notice that only (mu) and (a) are necessary for the two first, and to the\n\
     last one (f) in addition"); 
   
 liborb_orbital_to=sprintf("%s\n\n\
  Functions arguments are:\n\
   - mu parameter of mass    - O  long. of the ascending node\n\
   - a  semimajor axis       - w  argument of the periapsis\n\
   - e  eccentricity         - f  true anomaly\n\
   - I  inclination\n\n\
  Notice the short name functions have long name aliases!\n\
  Not all the parameters are necessary for all of the functions but the author choose\n\
  to keep the uniformity of the arguments.", liborb_orbital_to_funcs);

  
  
 liborb_cartesian_to_funcs=sprintf("#\n# Set of functions to deal with cartesian coordinates\n\
  The following functions convert the cartesian coordinates into orbital elements\n\
   - rv2p(mu,x,y,z,vx,vy,vz)     || cartesian_to_semilatus_rectum(mu,x,y,z,vx,vy,vz)\n\
   - rv2a(mu,x,y,z,vx,vy,vz)     || cartesian_to_semimajor_axis(mu,x,y,z,vx,vy,vz)\n\
   - rv2e(mu,x,y,z,vx,vy,vz)     || cartesian_to_eccentricity(mu,x,y,z,vx,vy,vz)\n\
   - rv2I(mu,x,y,z,vx,vy,vz)     || cartesian_to_inclination(mu,x,y,z,vx,vy,vz)\n\
   - rv2O(mu,x,y,z,vx,vy,vz)     || cartesian_to_ascending_node(mu,x,y,z,vx,vy,vz)\n\
   - rv2w(mu,x,y,z,vx,vy,vz)     || cartesian_to_argument_periapsis(mu,x,y,z,vx,vy,vz)\n\
   - rv2f(mu,x,y,z,vx,vy,vz)     || cartesian_to_true_anomaly(mu,x,y,z,vx,vy,vz)\n\
   - rv2n(mu,x,y,z,vx,vy,vz)     || cartesian_to_mean_motion(mu,x,y,z,vx,vy,vz)\n\
   - rv2T(mu,x,y,z,vx,vy,vz)     || cartesian_to_period(mu,x,y,z,vx,vy,vz)\n\
   - rv2gamma(mu,x,y,z,vx,vy,vz) || cartesian_to_flight_path_angle(mu,x,y,z,vx,vy,vz)");
  
 liborb_cartesian_to=sprintf("%s\n\n\
  Notice the short name functions have long name aliases!\n\
  In short names (rv) are mnemonics to position (r) and (v)elocity vectors.\n\
  Not all the parameters are necessary for all of the functions but the author choose\n\
  to keep the uniformity of the arguments.\n\n",liborb_cartesian_to_funcs );


 liborb_constants=sprintf("#\n# Set of constants:\n\
   - liborb_dpi (d)ouble PI\n\
   - liborb_hpi (h)alf   PI\n\
   - liborb_deg 1 degree in radians\n\
   - liborb_au  1 astronomical unit in [m]\n\
   - liborb_G   Universal gravitational constant [m^3/kg/s^2]\n\n\
  The following ones are parameters of mass in [m^3/s^2]\n\
   - liborb_musun\n\
   - liborb_mumer\n\
   - liborb_muven\n\
   - liborb_muear\n\
   - liborb_mulua\n\
   - liborb_mumar\n\
   - liborb_mujup\n\
   - liborb_musat\n\
   - liborb_muura\n\
   - liborb_munep\n");
 liborb_constant=liborb_constants;
 liborb_all=sprintf("%s\n\n%s\n\n%s\n\n%s\n\n%s\n\n",\
 liborb_constants,\
 liborb_perifocal_funcs,\
 liborb_anomalies_funcs,\
 liborb_orbital_to_funcs,\
 liborb_cartesian_to_funcs);



 liborb_angles=sprintf("All the angles are assumed to be in radians, the input and the output ones.\n\
  If you prefer to deal with angles in degrees, please call once the following macro:\n\
     type: @liborb_angles_deg\n\n\
  If you prefer to deal with angles in gradians, please call once the following macro:\n\
     type: @liborb_angles_gon\n\n\
  To query which is the current unity of angle:\n\
     type: print liborb_curang\n\n\
  If you want to choose any other unity for angles, please set the following:\n\
     liborb_curang = \"Your specification name\";  #only for your convenience;\n\
     liborb_ioang  = X_into_rad;\n\
  where X_into_rad is the conversion of the unity X into radians. BE CAREFUL!\n\n\
  Evidently, you may restore the package default option with angles in radians via:\n\
     type: @liborb_angles_rad\n\n\
  Be aware that the mean motion must agree with the unity angles per unity time!");



 liborb_unities=sprintf("#\n# You may use choose any set of unities you want for mass, length and time.\n\
  However, all the parameters must obey the same set. Let us assume that:\n\
    - [um] is the chosen unity of mass\n\
    - [ul] is the chosen unity of length\n\
    - [ut] is the chosen unity of time\n\
  Then, the unity for the parameter of mass (mu) must be:\n\
    - [mu] = [ul]^3/[ut]^2,\n\
  the unity for the velocities must be:\n\
    -[vel] = [ul]/[ut]\n\
  and the unity for the specific angular momentum will be:\n\
    -[h] = [ul]^2/[ut]\n\n\
  %s", liborb_angles);



 liborb_all=sprintf("%s\n\n%s\n\n%s\n\n%s\n\n%s\n\n",\
 liborb_constants,\
 liborb_perifocal_funcs,\
 liborb_anomalies_funcs,\
 liborb_orbital_to_funcs,\
 liborb_cartesian_to_funcs);
 
 liborb_more=sprintf("\n Orbital library %s (%s) for gnuplot \n\n\
 HSGlib_orbital library gathers a set of functions that makes gnuplot powerful\n\
 to deal with the ordinary tasks of the orbital mechanics. Plotting and\n\
 conversions at hand.\n\n\
 Please, type:\n\
  print liborb_all          #to see ALL the available functions\n\
  print liborb_anomalies    #to see the anomalies functions\n\
  print liborb_constants    #to see the the builtin constants\n\
  print liborb_perifocal    #to see the Perifocal frame functions\n\
  print liborb_orbital_to   #to see the conversions functions from orbital elements\n\
  print liborb_cartesian_to #to see the conversions functions from cartesian coordinates\n\
  print liborb_angles       #to see the how to change angle unities\n\
  \n\
  ",\
  HSG_libver, HSG_libdate);



################################################################################
#
## Macros

 liborb_angles_deg="liborb_ioang=liborb_deg; liborb_curang='deg'; print 'liborbital: Input and output angles set to degrees!\n'"
 liborb_angles_gon="liborb_ioang=liborb_gon; liborb_curang='gon'; print 'liborbital: Input and output angles set to gradians!\n'"
 liborb_angles_rad="liborb_ioang=1.0;        liborb_curang='rad'; print 'liborbital: Input and output angles set to radians!\n'"



################################################################################
#
## Package Header

  print sprintf("\n        HSG_lib orbital");
  print sprintf("        Version %-20slast modified %-10s\n", HSG_libver, HSG_libdate);
  print sprintf("        Author:              Helton da Silva Gaspar");
  print sprintf("        Further info:        type print liborb_more");
  print sprintf("        liborbital home:     http://helton.paginas.ufsc.br/");
  print sprintf("        Author's contact:    helton.s.gaspar@ufsc.br");
  print sprintf("");
  print sprintf("        Please, take care not to overwrite this package function names.");
  print sprintf("        Avoid define your functions with prefix \"liborb_\" ");
  print sprintf("\n");
  print sprintf("print liborb_more;    #for further info");

  
  
################################################################################
#
## Library default aliases

  deg        = liborb_deg;      # 1 degree in radians
  dpi        = liborb_dpi;      # (d)ouble pi
  hpi        = liborb_hpi;      # (h)alf pi
  au         = liborb_au;       # astronomical unity [m]
  G          = liborb_G;        # Universal gravitational constant [m^3/kg/s^2]
  
  #Parameters of mass 
  musun      = liborb_musun;
  mumer      = liborb_mumer;
  muven      = liborb_muven;
  muear      = liborb_muear;
  mumoon     = liborb_mulua;
  mumar      = liborb_mumar;
  mujup      = liborb_mujup;
  musat      = liborb_musat;
  muura      = liborb_muura;
  munep      = liborb_munep;


 liborb_constants=sprintf("#\n# Set of constants:\n\
   - dpi (d)ouble PI\n\
   - hpi (h)alf   PI\n\
   - deg 1 degree in radians\n\
   - au  1 astronomical unit in [m]\n\
   - G   Universal gravitational constant [m^3/kg/s^2]\n\n\
  The following ones are parameters of mass in [m^3/s^2]\n\
   - mu_sun\n\
   - mu_mercury\n\
   - mu_venus\n\
   - mu_earth\n\
   - mu_moon\n\
   - mu_mars\n\
   - mu_jupiter\n\
   - mu_saturn\n\
   - mu_uranus\n\
   - mu_neptune\n");
 liborb_constant=liborb_constants;
 
  #Quadratic and standard magnitudes
  mag2(Ux,Uy,Uz)            = liborb_mag2(Ux,Uy,Uz);
  mag (Ux,Uy,Uz)            = liborb_mag (Ux,Uy,Uz);
  
  #wrapping x into the domain xmin <= x <= xmax
  wrap(x,xmin,xmax)         = liborb_wrap(x,xmin,xmax);

  #Specific angular momentum vector. Components x, y and z and its quadratic magnitude
  hx(x,y,z,vx,vy,vz)        = liborb_hx(x,y,z,vx,vy,vz);
  hy(x,y,z,vx,vy,vz)        = liborb_hy(x,y,z,vx,vy,vz);
  hz(x,y,z,vx,vy,vz)        = liborb_hz(x,y,z,vx,vy,vz);
  h2(x,y,z,vx,vy,vz)        = liborb_h2(x,y,z,vx,vy,vz);

  #Vector eccentricity
  ex(mu,x,y,z,vx,vy,vz)     = liborb_ex(mu,x,y,z,vx,vy,vz);
  ey(mu,x,y,z,vx,vy,vz)     = liborb_ey(mu,x,y,z,vx,vy,vz);
  ez(mu,x,y,z,vx,vy,vz)     = liborb_ez(mu,x,y,z,vx,vy,vz);
  
  #Conversion from state vectors of position (r) and velocity (v) to orbital elements
  rv2p(mu,x,y,z,vx,vy,vz)   = liborb_rv2p(mu,x,y,z,vx,vy,vz);
  rv2a(mu,x,y,z,vx,vy,vz)   = liborb_rv2a(mu,x,y,z,vx,vy,vz);
  rv2e(mu,x,y,z,vx,vy,vz)   = liborb_rv2e(mu,x,y,z,vx,vy,vz);
  rv2I(mu,x,y,z,vx,vy,vz)   = liborb_rv2I(mu,x,y,z,vx,vy,vz);
  rv2i(mu,x,y,z,vx,vy,vz)   = liborb_rv2I(mu,x,y,z,vx,vy,vz);
  rv2O(mu,x,y,z,vx,vy,vz)   = liborb_rv2O(mu,x,y,z,vx,vy,vz);
  rv2w(mu,x,y,z,vx,vy,vz)   = liborb_rv2w(mu,x,y,z,vx,vy,vz);
  rv2f(mu,x,y,z,vx,vy,vz)   = liborb_rv2f(mu,x,y,z,vx,vy,vz);
  rv2n(mu,x,y,z,vx,vy,vz)   = liborb_rv2n(mu,x,y,z,vx,vy,vz);
  rv2T(mu,x,y,z,vx,vy,vz)   = liborb_rv2T(mu,x,y,z,vx,vy,vz);
  rv2Omg(mu,x,y,z,vx,vy,vz) = liborb_rv2O(mu,x,y,z,vx,vy,vz);
  rv2omg(mu,x,y,z,vx,vy,vz) = liborb_rv2w(mu,x,y,z,vx,vy,vz);
  rv2gamma(mu,x,y,z,vx,vy,vz) = liborb_rv2gamma(mu,x,y,z,vx,vy,vz);
  
  #Polar functions
  polp(a,e)      = liborb_polp(a,e);
  polr(a,e,f)    = liborb_polr(a,e,f);
  poldrdf(a,e,f) = liborb_poldrdf(a,e,f);
  
  #Perifocal frame
  x_pq(a,e,w,f)          = liborb_xpq(a,e,w,f);
  y_pq(a,e,w,f)          = liborb_ypq(a,e,w,f);
  vr(mu,a,e,f)           = liborb_vr(mu,a,e,f);
  vt(mu,a,e,f)           = liborb_vt(mu,a,e,f);
  vx_pq(mu,a,e,w,f)      = liborb_vxpq(mu,a,e,w,f);
  vy_pq(mu,a,e,w,f)      = liborb_vypq(mu,a,e,w,f);
  gamma(e,f)             = liborb_gamma(e,f);
  flight_path_angle(e,f) = liborb_gamma(e,f);

  #Conversion from orbital elements into cartesian coordinates
  oe2x (mu,a,e,I,O,w,f)                  = liborb_eo2x (mu,a,e,I,O,w,f);
  oe2y (mu,a,e,I,O,w,f)                  = liborb_eo2y (mu,a,e,I,O,w,f);
  oe2z (mu,a,e,I,O,w,f)                  = liborb_eo2z (mu,a,e,I,O,w,f);
  oe2vx(mu,a,e,I,O,w,f)                  = liborb_eo2vx(mu,a,e,I,O,w,f);
  oe2vy(mu,a,e,I,O,w,f)                  = liborb_eo2vy(mu,a,e,I,O,w,f);
  oe2vz(mu,a,e,I,O,w,f)                  = liborb_eo2vz(mu,a,e,I,O,w,f);
  oe2n (mu,a,e,I,O,w,f)                  = liborb_eo2n (mu,a,e,I,O,w,f);
  oe2T (mu,a,e,I,O,w,f)                  = liborb_eo2T (mu,a,e,I,O,w,f);
  oe2t (mu,a,e,I,O,w,f)                  = liborb_eo2t (mu,a,e,I,O,w,f);

  #Time-angle and anomalies correlations
  f2E(e,f)                               = liborb_f2E(e,f);
  f2M(e,f)                               = liborb_f2M(e,f);
  f2t(n,e,f)                             = liborb_f2t(n,e,f);
  E2M(e,E)                               = liborb_E2M(e,E);
  E2f(e,E)                               = liborb_E2f(e,E);
  M2t(n,M)                               = liborb_M2t(n,M);
  M2E(e,M)                               = liborb_M2E(e,M);
  M2f(e,M)                               = liborb_M2f(e,M);
  t2M(n,t)                               = liborb_t2M(n,t);
  t2E(n,e,t)                             = liborb_t2E(n,e,t);
  t2f(n,e,t)                             = liborb_t2f(n,e,t);



################################################################################
#
## Long name aliases

  #Parameters of mass
  mu_sun     = liborb_musun;
  mu_mercury = liborb_mumer;
  mu_venus   = liborb_muven;
  mu_earth   = liborb_muear;
  mu_moon    = liborb_mulua;
  mu_mars    = liborb_mumar;
  mu_jupiter = liborb_mujup;
  mu_saturn  = liborb_musat;
  mu_uranus  = liborb_muura;
  mu_neptune = liborb_munep;
  
  #Polar functions
  polar_p(a,e)          = liborb_polp(a,e);
  polar_r(a,e,f)        = liborb_polr(a,e,f);
  polar_drdf(a,e,f)     = liborb_poldrdf(a,e,f);

  #Conversion from state vectors of position (r) and velocity (v) to orbital elements
  rvtop(mu,x,y,z,vx,vy,vz)       = liborb_rv2p(mu,x,y,z,vx,vy,vz);
  rvtoa(mu,x,y,z,vx,vy,vz)       = liborb_rv2a(mu,x,y,z,vx,vy,vz);
  rvtoe(mu,x,y,z,vx,vy,vz)       = liborb_rv2e(mu,x,y,z,vx,vy,vz);
  rvtoI(mu,x,y,z,vx,vy,vz)       = liborb_rv2I(mu,x,y,z,vx,vy,vz);
  rvtoi(mu,x,y,z,vx,vy,vz)       = liborb_rv2I(mu,x,y,z,vx,vy,vz);
  rvtoO(mu,x,y,z,vx,vy,vz)       = liborb_rv2O(mu,x,y,z,vx,vy,vz);
  rvtow(mu,x,y,z,vx,vy,vz)       = liborb_rv2w(mu,x,y,z,vx,vy,vz);
  rvtof(mu,x,y,z,vx,vy,vz)       = liborb_rv2f(mu,x,y,z,vx,vy,vz);
  rvton(mu,x,y,z,vx,vy,vz)       = liborb_rv2n(mu,x,y,z,vx,vy,vz);
  rvtoT(mu,x,y,z,vx,vy,vz)       = liborb_rv2T(mu,x,y,z,vx,vy,vz);
  rvtogamma(mu,x,y,z,vx,vy,vz)   = liborb_rv2gamma(mu,x,y,z,vx,vy,vz);

  rv2Omg(mu,x,y,z,vx,vy,vz)      = liborb_rv2O(mu,x,y,z,vx,vy,vz);
  rv2omg(mu,x,y,z,vx,vy,vz)      = liborb_rv2w(mu,x,y,z,vx,vy,vz);

  cartesian_to_semilatus_rectum(mu,x,y,z,vx,vy,vz)   = liborb_rv2p(mu,x,y,z,vx,vy,vz);
  cartesian_to_semilatus(mu,x,y,z,vx,vy,vz)          = liborb_rv2p(mu,x,y,z,vx,vy,vz);
  cartesian_to_semimajor_axis(mu,x,y,z,vx,vy,vz)     = liborb_rv2a(mu,x,y,z,vx,vy,vz);
  cartesian_to_semimajor(mu,x,y,z,vx,vy,vz)          = liborb_rv2a(mu,x,y,z,vx,vy,vz);
  cartesian_to_eccentricity(mu,x,y,z,vx,vy,vz)       = liborb_rv2e(mu,x,y,z,vx,vy,vz);
  cartesian_to_inclination(mu,x,y,z,vx,vy,vz)        = liborb_rv2I(mu,x,y,z,vx,vy,vz);
  cartesian_to_ascending_node(mu,x,y,z,vx,vy,vz)     = liborb_rv2O(mu,x,y,z,vx,vy,vz);
  cartesian_to_node(mu,x,y,z,vx,vy,vz)               = liborb_rv2O(mu,x,y,z,vx,vy,vz);
  cartesian_to_argument_periapsis(mu,x,y,z,vx,vy,vz) = liborb_rv2w(mu,x,y,z,vx,vy,vz);
  cartesian_to_periapsis(mu,x,y,z,vx,vy,vz)          = liborb_rv2w(mu,x,y,z,vx,vy,vz);
  cartesian_to_true_anomaly(mu,x,y,z,vx,vy,vz)       = liborb_rv2f(mu,x,y,z,vx,vy,vz);
  cartesian_to_mean_motion(mu,x,y,z,vx,vy,vz)        = liborb_rv2n(mu,x,y,z,vx,vy,vz);
  cartesian_to_period(mu,x,y,z,vx,vy,vz)             = liborb_rv2T(mu,x,y,z,vx,vy,vz);
  cartesian_to_flight_path_angle(mu,x,y,z,vx,vy,vz)  = liborb_rv2gamma(mu,x,y,z,vx,vy,vz);
  
  cartesian_to_p(mu,x,y,z,vx,vy,vz)   = liborb_rv2p(mu,x,y,z,vx,vy,vz);
  cartesian_to_a(mu,x,y,z,vx,vy,vz)   = liborb_rv2a(mu,x,y,z,vx,vy,vz);
  cartesian_to_e(mu,x,y,z,vx,vy,vz)   = liborb_rv2e(mu,x,y,z,vx,vy,vz);
  cartesian_to_I(mu,x,y,z,vx,vy,vz)   = liborb_rv2I(mu,x,y,z,vx,vy,vz);
  cartesian_to_i(mu,x,y,z,vx,vy,vz)   = liborb_rv2I(mu,x,y,z,vx,vy,vz);
  cartesian_to_O(mu,x,y,z,vx,vy,vz)   = liborb_rv2O(mu,x,y,z,vx,vy,vz);
  cartesian_to_Omg(mu,x,y,z,vx,vy,vz) = liborb_rv2O(mu,x,y,z,vx,vy,vz);
  cartesian_to_w(mu,x,y,z,vx,vy,vz)   = liborb_rv2w(mu,x,y,z,vx,vy,vz);
  cartesian_to_omg(mu,x,y,z,vx,vy,vz) = liborb_rv2w(mu,x,y,z,vx,vy,vz);
  cartesian_to_f(mu,x,y,z,vx,vy,vz)   = liborb_rv2f(mu,x,y,z,vx,vy,vz);
  cartesian_to_n(mu,x,y,z,vx,vy,vz)   = liborb_rv2n(mu,x,y,z,vx,vy,vz);
  cartesian_to_T(mu,x,y,z,vx,vy,vz)   = liborb_rv2T(mu,x,y,z,vx,vy,vz);

  #Conversion from orbital elements into cartesian coordinates
  orbital_to_x (mu,a,e,I,O,w,f)          = liborb_eo2x (mu,a,e,I,O,w,f);
  orbital_to_y (mu,a,e,I,O,w,f)          = liborb_eo2y (mu,a,e,I,O,w,f);
  orbital_to_z (mu,a,e,I,O,w,f)          = liborb_eo2z (mu,a,e,I,O,w,f);
  orbital_to_vx(mu,a,e,I,O,w,f)          = liborb_eo2vx(mu,a,e,I,O,w,f);
  orbital_to_vy(mu,a,e,I,O,w,f)          = liborb_eo2vy(mu,a,e,I,O,w,f);
  orbital_to_vz(mu,a,e,I,O,w,f)          = liborb_eo2vz(mu,a,e,I,O,w,f);
  orbital_to_n (mu,a,e,I,O,w,f)          = liborb_eo2n (mu,a,e,I,O,w,f);
  orbital_to_T (mu,a,e,I,O,w,f)          = liborb_eo2T (mu,a,e,I,O,w,f);
  orbital_to_t (mu,a,e,I,O,w,f)          = liborb_eo2t(mu,a,e,I,O,w,f);
  orbital_to_mean_motion(mu,a,e,I,O,w,f) = liborb_eo2n (mu,a,e,I,O,w,f);
  orbital_to_period(mu,a,e,I,O,w,f)      = liborb_eo2T (mu,a,e,I,O,w,f);
  orbital_to_time(mu,a,e,I,O,w,f)        = liborb_eo2t(mu,a,e,I,O,w,f);

  #Time-angle and anomalies correlations
  true_anomaly_to_eccentric_anomaly(e,f) = liborb_f2E(e,f);
  true_anomaly_to_mean_anomaly(e,f)      = liborb_f2M(e,f);
  true_anomaly_to_time(n,e,f)            = liborb_f2t(n,e,f);
  eccentric_anomaly_to_mean_anomaly(e,E) = liborb_E2M(e,E);
  eccentric_anomaly_to_true_anomaly(e,E) = liborb_E2f(e,E);
  mean_anomaly_to_time(n,M)              = liborb_M2t(n,M);
  mean_anomaly_to_eccentric_anomaly(e,M) = liborb_M2E(e,M);
  mean_anomaly_to_true_anomaly(e,M)      = liborb_M2f(e,M);
  time_to_mean_anomaly(n,t)              = liborb_t2M(n,t);
  time_to_eccentric_anomaly(n,e,t)       = liborb_t2E(n,e,t);
  time_to_true_anomaly(n,e,t)            = liborb_t2f(n,e,t);