#ifndef __MYPO896__
#define __MYPO896__
#include "Matrix.hh"
#include "Vector.hh"
#include "SparseMatrix.hh"

class PeriodicOrbit;
typedef void (* vfTYPE)(const Vector&, Vector&);
typedef void (* jTYPE)(const Vector&, Matrix&);

//utility routines
void ComputeResidual(const Matrix& X, double w, const Matrix& dX,
		     Matrix& Res, vfTYPE vfield);
void ComputeResidual(const Matrix& X, double w, 
		     Matrix& Res, vfTYPE vfield);
void ComputeLinearization(const Matrix& X, Matrix& A, jTYPE jacobian);
void ComputeGradDirn(double w, const Matrix& dX, const Matrix& Res, 
		     const Matrix& A, Matrix& gradX, double& gradw);
void ComputeGradDirn(const Matrix& X, double w,
		     Matrix& gradX, double& gradw,
		     vfTYPE vfield, jTYPE jacobian);
double ComputeProjectionFactor(const Matrix& A, const Matrix& B);
void FindFundSolnA(const Matrix& A, Matrix& Y);
void FindFundSolnB(const Matrix& A, Matrix& Y);
void ApplyFundSolnA(const Matrix& A, const Matrix& g, Matrix& gg);
void ApplyFundSolnB(const Matrix& A, const Matrix& g, Matrix& gg);
void ApplyFundSolnC(const Matrix& A, const Matrix& g, Matrix& gg);
void ApplyFundSolnD(const Matrix& A, const Matrix& g, Matrix& gg);
void printstuff();
void printSTUFF();
//trust region method
void ComputeTrustStep(PeriodicOrbit& po, double& delta,
		      Matrix& step, double& wstep);

class PeriodicOrbit{
public:
  PeriodicOrbit(const Matrix& po, double t, vfTYPE vf, jTYPE jac);//po and t
  // are the initial approximations. 
  ~PeriodicOrbit();
  void addStep(const Matrix& step, double wstep);
  double residualError();
  double relativeResidualError();
  double residualErrorAfterStep(const Matrix& step, double wstep);
  double predictedResidualError(const Matrix& step, double wstep);
  double dotStepAndGrad(const Matrix& step, double wstep);//computes
  //dot product of (step,wstep) with the gradient of the residual 
  void getPeriodicOrbit(Matrix& M, double& Period);
  void smooth(double alpha=0.9, double beta=1.8);
  void getPowellStep(Matrix& step, double& wstep);
  void getNewtonStep(Matrix& step, double& wstep, int flag=0);
  void getLinearResidual(Matrix& step, double wstep, Matrix& linRes);
  void setSintvl(int intvl=64);
  void setResidual(Matrix& res);
  int getd();
  int getn();
private:
  vfTYPE vfield;
  jTYPE jacobian;
  Matrix pOrbit; //matrix holding the periodic orbit
  double T; //period
  double tinyFactor;
private:
  //cache dX, A, Res, GradX, Gradw here
  Matrix dX;
  Matrix A;
  Matrix Res;
  Matrix GradX;
  double Gradw;
  void updateCache();
private:
  //working variables and functions for Newton's method
  Matrix f;
  Matrix g;
  Matrix Y;
  int *spts;
  int numOfSpts;
  Matrix *ff;
  Matrix *gg;
  Matrix *YY;
  Matrix *vf;
  SparseMatrix *LPMatrix; 
  void allocateYYggffvf();
  void computeY();
  void computef();
  void computeg();
  void computeYY();
  void computegg();
  void computeff();
  void computevf();
  void computeLPMatrix();
  void computeShootingStep(Matrix& sstep, double& wstep);
  void propagatey(Vector& yi, int i, double wstep,
		  int k, Vector& yk);
  void convertSStep2NStep(Matrix& sstep, double wstep, Matrix& nstep);
};


#endif