#include "Vector.hh"
#include "Matrix.hh"
#include "lorenz.hh"
#include "lorenz1.hh"
#include "lorenz4.hh"
#include "PeriodicOrbit.hh"
#include <iostream>

double getAngle(Matrix& po, double T){
  assert(po.getm()==3);
  int n = po.getn();
  Matrix A(9, n);
  ComputeLinearization(po, A, lorenzJ);
  A.scale(T/(2*PI)); //forgot to do this initially!
  Matrix Y(9, n);
  FindFundSolnB(A, Y);
  Matrix m(3,3);
  Matrix mm;
  int i, j;
  for(i=0; i < 3; i++)
    for(j=0; j < 3; j++)
      m(i,j) = (i==j)?1:0;
  double factor = 1;
  for(i=0; i < n; i++){
    mm.shadow(Y.getColumn(i), 3, 3);
    mult(mm, m);
    m.scale(1/factor);
  } 
  Matrix V(3,3);
  lorenz(po.getColumn(0), V.getColumn(0));
  V(0,1) = 1; V(1,1) = 0; V(2,1) = 0;
  V(0,2) = 0; V(1,2) = 1; V(2,2) = 0;
  Matrix iV(3,3);
  invert(V, iV);
  mult(iV, m);
  mult(m, V);
  m = V;
  Matrix ev(2,2);
  Vector er(2), ei(2);
  eig(m.getShadow(1,2,1,2), ev, er, ei);
  //cout<<er(0)<<endl<<er(1)<<endl<<endl;
  Vector v1(2), v2(2);
  v1 = ev.getColumn(0);
  v2 = ev.getColumn(1);
  double cosangle = v1(0)*v2(0)+v1(1)*v2(1);
  cosangle /= (v1.norm()*v2.norm());
  return acos(cosangle);
}

//returns angle of A^n B^n
double getAngle(int n){
  char s[100];
  int i;
  for(i=0; i<n; i++){
    s[i] = 'A';
    s[n+i] = 'B';
  }
  s[2*n] = '\0';
  
  Matrix *po;
  char dir[100] = "DATA3/";
  double T; 
  double err;
  po = loadANDreturnOrbit(dir, n, T, err);
  shift2FirstPoint(*po, s);
  return getAngle(*po, T);
}
  
void printAngles(ostream& ostr){
  int i;
  for(i=1; i<=25; i++){
    double angle;
    angle= getAngle(i);
    ostr<<"n = "<<i<<"  angle = "<< angle/PI <<endl;
  }
}