// VERSION 1.1  (March 22, 2004)
//
// Given a satellite in orbit around a much more massive body at the
// origin, with DIST1 (DIST2) as the nearest (furthest) distance of
// the orbit and ANG0 as the angle at t=0, this will return the angle
// of the body at time t=CLK with an error of no more than MAXERR.
//
// Presumes DIST1 along +X, DIST2 along -X, and angles (all in
// radians) measured ccw from +X in the X,Y plane.
//
#macro OrbitAngle( clk, ang0, period, dist1, dist2, maxerr )
   #local _ee   = sqrt(1-dist1*dist1/dist2/dist2);
   #local _mang = ang0+clk/period*2*pi;
   #local _err  = maxerr*2;
   #local _ang  = _mang;
   #while ( _err > maxerr )
      #local _ang1 = _mang + _ee * sin(_ang);
      #local _err  = abs(_ang1 - _ang);
      #local _ang  = _ang1;
   #end
   _ang
#end

//
// Given the same inputs as for OrbitAngle above, this returns the
// vector of the satellite from the massive body.
//
#macro OrbitLocation( clk, ang0, period, dist1, dist2, maxerr )
   #if (dist1>dist2)
      #local   dist_ = dist1;
      #declare dist1 = dist2;
      #declare dist2 = dist_;
   #end
   #local _a   = (dist1+dist2)/2;
   #local _b   = (dist2-dist1)/2; #local _b = sqrt(_a*_a-_b*_b);
   #local _e   = sqrt(_a*_a-_b*_b)/_a;
   #if (_e>1)
      #error concat( "Sorry, but eccentricity of ellipse is over 1: ",str(_e,9,6),"\n")
   #end
   #local _ang = OrbitAngle( clk, ang0, period, dist1, dist2, maxerr );
   #local _rad = _a * (1-_e*_e) / (1+_e*cos(_ang));
   #local _xx  = _rad*sin(_ang);
   #local _yy  = _rad*cos(_ang);
#debug concat("OrbLoc[A=",str(_a,8,6),",E=",str(_e,8,6),",ANG=",str(_ang,8,6),",RAD=",str(_rad,8,6),"] = <",str(_xx,8,6),",",str(_yy,8,6),",0>\n")
   <_xx,_yy,0>
#end

//
// Given positional vectors POS1 and POS2 for orbiting bodies (with
// MASS1, MASS2) this macro will alter POS1 and POS2 to place the 
// barycenter for pair at POS0 (but retain their relative postions).
//
// Ref: http://encyclopedia.thefreedictionary.com/barycenter
//
#macro OrbitBaryCenter( pos0, pos1, pos2, mass1, mass2 )
   #local _cntr = pos0;
   #local _vect = pos2 - pos1;
   #local _bary = _vect * mass2 / ( mass1 + mass2 );
   #declare pos1 = _cntr - _bary;
   #declare pos2 = pos1 + _vect;
#end