  function g = gSep(t,plist)
%
% Pre:  t is a 2-vector and plist encodes the orbits
%       (A1,P1,phi1) and (A2,P2,phi2).
%
% Post: The gradient of sep(t,plist).

   A1 = plist(1); P1 = plist(2); phi1 = plist(3); 
   A2 = plist(4); P2 = plist(5); phi2 = plist(6); 
  

   alfa1  = (P1-A1)/2; beta1  = (P1+A1)/2; gamma1 = sqrt(P1*A1);
   alfa2  = (P2-A2)/2; beta2  = (P2+A2)/2; gamma2 = sqrt(P2*A2);
   s1     = sin(t(1)); c1     = cos(t(1));
   s2     = sin(t(2)); c2     = cos(t(2));
   cphi1  = cos(phi1); sphi1  = sin(phi1);
   cphi2  = cos(phi2); sphi2  = sin(phi2);
   Rot1   = [cphi1 sphi1; -sphi1 cphi1];
   Rot2   = [cphi2 sphi2; -sphi2 cphi2];
   P1     = Rot1*[alfa1+beta1*c1;gamma1*s1];
   P2     = Rot2*[alfa2+beta2*c2;gamma2*s2];
   dP1    = Rot1*[-beta1*s1;gamma1*c1];
   dP2    = Rot2*[-beta2*s2;gamma2*c2]; 
   g = [-dP1';dP2']*(P2-P1);
   
   
