//  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 [M] = CL_cw_Mmatrix(alt,delta_t, er,mu)
// Clohessy-Wiltshire M transformation matrix
//
// Calling Sequence
// M = CL_cw_Mmatrix(alt,delta_t [,er,mu])
//
// Description
// <itemizedlist><listitem>
// <p>Computes the Clohessy-Wiltshire <b>M</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>
// </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)
// M: Clohessy-Wiltshire "M" transformation matrix (6x6xN)
//
// Authors
// CNES - DCT/SB
//
// See also
// CL_cw_Nmatrix
// CL_cw_propagate
//
// Bibliography
// 1) Mecanique spatiale, CNES - Cepadues 1995, Tome II, chapter 16
//
// Examples
// alt = [350.e3 550.e3];
// delta_t = [100 150];
// [M]=CL_cw_Mmatrix(alt,delta_t)
//

// 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); 

M = hypermat([6, 6, n]); 

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

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

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

endfunction