//  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 [N] = CL_cw_Nmatrix(alt,delta_t, er,mu)
// Clohessy-Wiltshire N transformation Matrix
//
// Calling Sequence
// N = CL_cw_Nmatrix(alt,delta_t [,er,mu])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the Clohessy-Wiltshire <b>N</b> transformation matrix.</p>
// <p>The local orbital reference frame tied to the target is LVLH (See <link linkend="CL_fr_lvlhMat">CL_fr_lvlhMat</link>).</p>
// <p><inlinemediaobject><imageobject><imagedata fileref="cw.gif"/></imageobject></inlinemediaobject></p>
// <p></p>
// <p><b>Notes:</b></p>
// <p>In the above formula: </p>
// <p>- theta is the phase angle of the target (= "reference") satellite = omega * delta_t, where 
// omega is the target (= "reference") mean motion. </p>
// <p>- delta gamma = difference of acceleration between the chaser and the target (gamma chaser minus gamma target)
// relative to LVLH local frame (of the target). </p>
// </listitem>
// </itemizedlist>
//
// Parameters
// alt:  Target altitude (semi-major axis minus equatorial radius) [m] (1xN or 1x1)
// delta_t: Time step [s] (1xN or 1x1)
// er: (optional) Equatorial radius [m] (default is %CL_eqRad)
// mu: (optional) Gravitational constant [m^3/s^2] (default is %CL_mu)
// N: Clohessy-Wiltshire "N" transformation matrix (6x6xN)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_cw_Mmatrix
// CL_cw_propagate
//
// Bibliography
// 1) Mecanique spatiale, CNES - Cepadues 1995, Tome II, chapter 16 (Le rendez-vous orbital)
//
// Examples
// [N]=CL_cw_Nmatrix(350000,100)
//

// Declarations:

// Code:
if (~exists("er", "local")); er = CL__dataGetEnv("eqRad"); end
if (~exists("mu", "local")); mu = CL__dataGetEnv("mu"); end

// check/resize inputs (n = number of columns)
[alt, delta_t, n] = CL__checkInputs(alt, 1, delta_t, 1); 

N = hypermat([6, 3, n]); 

omega = CL_kp_params('mm', alt+er, mu); // mean motion
teta = omega .* delta_t; 
c = cos(teta); 
s = sin(teta); 

N(1,1,:) = (1 ./ omega.^2) .* ( 4*(1-c) - (3/2) * teta.^2); 
N(3,1,:) = (2 ./ omega.^2) .* (s-teta); 
N(4,1,:) = (4*s - 3*teta) ./ omega; 
N(6,1,:) = (2 ./ omega) .* (c-1); 
N(2,2,:) = (1-c) ./ omega.^2; 
N(5,2,:) = s ./ omega; 
N(1,3,:) = (2 ./ omega.^2) .* (teta-s); 
N(3,3,:) = (1-c) ./ omega.^2; 
N(4,3,:) = (2 ./ omega) .* (1-c); 
N(6,3,:) = s ./ omega; 

if (n == 1); N = N(:,:,1); end

endfunction