#include <cmath>
#include <vector>
#include <iostream>

using namespace std;

/*
  bundle two coordinates into one position
*/

class position
{
private:
  double r,p;

public:
  position(double rad=0, double phi=0) : r(rad), p(phi) {}
  double x() {return r*cos(p);}
  double y() {return r*sin(p);}
  position &operator=(const position &pos)
  {
    if (this==&pos) return *this;
    r=pos.r;
    p=pos.p;
    return *this;
  }
};

class orbit
{
private:
  static const double eps=1.0e-6;
  static const double tpi=6.28318531;
  double e;

public:
  orbit(double ex=0) : e(ex)
  {
    if (e>1||e<0)
      {
	cerr << "not an ellitical orbit, exiting" << endl;
	exit(-1);
      }
  }

  position getPosition(double t=0) // solves Kepler's eq.
  {
    double p0=tpi*t,p1;
    for (;;)
      {
	p1=tpi*t+e*sin(p0);
	if (fabs(p1-p0)<eps) break;
	p0=p1;
      }
    double r=1-e*cos(p1);
    double c=acos((cos(p1)-e)/(1-e*cos(p1)));
    if (t>0.5) c=tpi-c;
    position pos(r,c);
    return pos;
  }
};

/*
  main routine: initialise orbit, then compute planet's position for
  n times between two perihelion transits
*/

int main ()
{
  const int n=128;
  orbit orb(0.75);

  for (int i=0;i<n;i++)
    {
      double t=i/double(n-1);
      position p=orb.getPosition(t);
      cout << t << " " << p.x() << " " << p.y() << endl;
    }
  return 0;
}

