//  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 [s] = CL_stumpS(z)
// Stumpff function S(z)
//
// Calling Sequence
// s = CL_stumpS(z)
//
// Description
// <itemizedlist><listitem>
// <p>Evaluates the Stumpff function S(z) according to the following equation :</p>
// <p><inlinemediaobject><imageobject><imagedata fileref="stumpS.gif"/></imageobject></inlinemediaobject></p>
// </listitem>
// </itemizedlist>
//
// Parameters
// z: input argument (1xN)
// s: value of S(z) (1xN)
//
// Authors
// CNES - DCT/SB
//
// Bibliography
// 1) Orbital Mechanics for engineering students, H D Curtis, Chapter 3 and Appendix D.4
// 2) Modern astrodynamics Fundamentals and perturbation methods, V Bond and M Allman, Chapter 6.2
//
// See also
// CL_stumpC
//
// Examples
// z=0:1:500;
// s=CL_stumpS(z);
// plot(z,s);
//

// Declarations:


// Code:

s = %nan * ones(z);   // CL_stumpC(%nan) = %nan

I1 = find( z > 0);
I2 = find( z < 0);
I3 = find( z == 0); // CL_stumpS(0) = 1/6

s(I1) = (sqrt(z(I1)) - sin(sqrt(z(I1))))./(sqrt(z(I1))).^3;
s(I2) = (sinh(sqrt(-z(I2))) - sqrt(-z(I2)))./(sqrt(-z(I2))).^3;
s(I3) = 1/6;


endfunction