#include "Sundance.h"
#include "NewtonSolver.h"

namespace Sundance
{
  class Bessel2D : public UserDefinedFunction
  {
  public:
    Bessel2D(int n) : UserDefinedFunction(), n_(n) {;}

    virtual ~Bessel2D() {;}

    virtual double evaluateAtPoint(const Point& x) const ;

  private:
    int n_;
  };
}

double Bessel2D::evaluateAtPoint(const Point& x) const
{
  double r = ::sqrt(x*x);
  return jn(n_, r);
}

int main(int argc, void** argv)
{
  try
    {
			
      Sundance::init(&argc, &argv);

      /* create a simple mesh on the unit line */
      double L=6.0;
      int nx = 80;
      int ny = 80;
      MeshGenerator mesher = new RectangleMesher(-L, L, nx, -L, L, ny);
      Mesh mesh = mesher.getMesh();

      TSFVectorSpace discreteSpace 
        = new SundanceVectorSpace(mesh, new Lagrange(2));

      Expr J0 = new UserFuncExpr(new Bessel2D(0), "J0");
      Expr u0 = new DiscreteFunction(discreteSpace, J0);

      /* write to matlab */
      FieldWriter writer = new MatlabWriter("J0.dat");
      writer.writeField(u0);	

    }
  catch(exception& e)
    {
      Sundance::handleError(e, __FILE__);
    }
  Sundance::finalize();

}