//  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 [Zdot] = CL__3b_RHSDP(t, Z, MU)
// Motion model in the 2-body system (including monodromy matrix)
// Inputs: 
// t: time
// Z: input vector (42x1) 
//  (6: position+velocity + 36 for the monodromy matrix)
// MU: mass ratio (> 0 and <= 1)
//
// Outputs: 
// Zdot: derivative (42x1) 
//  (6: position+velocity + 36 for the monodromy matrix)
//
// Author:
// B. Meyssignac (CNES DCT/SB/MO)
// R. Alacevich  (Stage CNES DCT/SB/MO)

// Declarations:


// Code:

r0 = sqrt((Z(1)+MU)^2+Z(2)^2+Z(3)^2);
r1 = sqrt((Z(1)-1+MU)^2+Z(2)^2+Z(3)^2);

Xdot = [Z(4);
   Z(5);
   Z(6);
   2*Z(5)+Z(1)-(1-MU)*(Z(1)+MU)/(r0^3)-MU*(Z(1)-1+MU)/(r1^3);
   -2*Z(4)+Z(2)-(1-MU)*Z(2)/(r0^3)-MU*Z(2)/(r1^3);
   -(1-MU)*Z(3)/(r0^3)-MU*Z(3)/(r1^3)];

C2 = (1-MU)/r0^3+MU/r1^3;

// derivees partielles de OMEGA:
OMEGAxx = 1-C2+3*(1-MU)*(Z(1)+MU)^2/r0^5+3*MU*(Z(1)-1+MU)^2/r1^5;
OMEGAxy = 3*(1-MU)*(Z(1)+MU)*Z(2)/r0^5+3*MU*(Z(1)-1+MU)*Z(2)/r1^5;
OMEGAxz = 3*(1-MU)*(Z(1)+MU)*Z(3)/r0^5+3*MU*(Z(1)-1+MU)*Z(3)/r1^5;
OMEGAyy = 1-C2+3*(1-MU)*Z(2)^2/r0^5+3*MU*Z(2)^2/r1^5;
OMEGAyz = 3*(1-MU)*Z(2)*Z(3)/r0^5+3*MU*Z(2)*Z(3)/r1^5;
OMEGAzz = -C2+3*(1-MU)*Z(3)^2/r0^5+3*MU*Z(3)^2/r1^5;

Ydot = [Z(10);
   Z(11);
   Z(12);
   OMEGAxx*Z(7)+OMEGAxy*Z(8)+OMEGAxz*Z(9)+2*Z(11);
   OMEGAxy*Z(7)+OMEGAyy*Z(8)+OMEGAyz*Z(9)-2*Z(10);
   OMEGAxz*Z(7)+OMEGAyz*Z(8)+OMEGAzz*Z(9);
   Z(16);
   Z(17);
   Z(18);
   OMEGAxx*Z(13)+OMEGAxy*Z(14)+OMEGAxz*Z(15)+2*Z(17);
   OMEGAxy*Z(13)+OMEGAyy*Z(14)+OMEGAyz*Z(15)-2*Z(16);
   OMEGAxz*Z(13)+OMEGAyz*Z(14)+OMEGAzz*Z(15);
   Z(22);
   Z(23);
   Z(24);
   OMEGAxx*Z(19)+OMEGAxy*Z(20)+OMEGAxz*Z(21)+2*Z(23);
   OMEGAxy*Z(19)+OMEGAyy*Z(20)+OMEGAyz*Z(21)-2*Z(22);
   OMEGAxz*Z(19)+OMEGAyz*Z(20)+OMEGAzz*Z(21);
   Z(28);
   Z(29);
   Z(30);
   OMEGAxx*Z(25)+OMEGAxy*Z(26)+OMEGAxz*Z(27)+2*Z(29);
   OMEGAxy*Z(25)+OMEGAyy*Z(26)+OMEGAyz*Z(27)-2*Z(28);
   OMEGAxz*Z(25)+OMEGAyz*Z(26)+OMEGAzz*Z(27);
   Z(34);
   Z(35);
   Z(36);
   OMEGAxx*Z(31)+OMEGAxy*Z(32)+OMEGAxz*Z(33)+2*Z(35);
   OMEGAxy*Z(31)+OMEGAyy*Z(32)+OMEGAyz*Z(33)-2*Z(34);
   OMEGAxz*Z(31)+OMEGAyz*Z(32)+OMEGAzz*Z(33);
   Z(40);
   Z(41);
   Z(42);
   OMEGAxx*Z(37)+OMEGAxy*Z(38)+OMEGAxz*Z(39)+2*Z(41);
   OMEGAxy*Z(37)+OMEGAyy*Z(38)+OMEGAyz*Z(39)-2*Z(40);
   OMEGAxz*Z(37)+OMEGAyz*Z(38)+OMEGAzz*Z(39)]; 

Zdot = [Xdot; Ydot];

endfunction