// GenFem - A high-level finite element library
// Copyright (C) 2010-2026 Eric Bechet
//
// See the LICENSE file for license information and contributions.
// Please report all bugs and problems to <bechet@cadxfem.org>.

#ifndef _GEN_SIMPLE_FUNCTION_UTR_H_
#define _GEN_SIMPLE_FUNCTION_UTR_H_


#include "genSimpleFunction.h"
#include "genSimpleFunctionExtended.h"
#include <boost/math/special_functions/bessel.hpp>
#include <cmath>

template<class T1,class T2> class genSimpleFunctionUtr : public genSimpleFunction<T1>
{
public :
  typedef typename TensorialTraits<T1>::Scalar Scalar;
  typedef typename TensorialTraits<T1>::ScalarType ScalarType;
  typedef typename TensorialTraits<T1>::ValType ValType;
  typedef typename TensorialTraits<T1>::GradType GradType;
  typedef typename TensorialTraits<T1>::HessType HessType;
    
  std::vector<std::vector<Scalar> > Coefficients;
  genSimpleFunction<T2> *polar;
  double kp_int, ks_int;
  int comp;
public :
  genSimpleFunctionUtr(std::vector<std::vector<Scalar> >& Coefficients_,double &kp_,double &ks_, int comp_ = -1) : Coefficients(Coefficients_) , kp_int(kp_), ks_int(ks_) , comp(comp_) 
  {
        polar = new genSimpleFunctionCylindricalCoordSyst<T2>();
  }
  virtual ValType operator()(double x, double y, double z) const {
    
    std::vector<double> polar_coor = vectorize((*polar)(x,y,z));
    double r = polar_coor[0];
    double theta = polar_coor[1];
    Scalar u_3r, u_3theta, ux, uy;    
    for(int n=0; n< Coefficients.size(); n++){
      u_3r += (Coefficients[n][4]*kp_int*( -boost::math::cyl_bessel_j(n+1,kp_int*r) + n/(kp_int*r)*boost::math::cyl_bessel_j(n,kp_int*r) )+
	      n/r*Coefficients[n][5]*boost::math::cyl_bessel_j(n,ks_int*r) )*cos(n*theta);
      
      u_3theta += (-n/r*Coefficients[n][4]*boost::math::cyl_bessel_j(n, kp_int*r) - 
		Coefficients[n][5]*ks_int*(-boost::math::cyl_bessel_j(n+1,ks_int*r) + n/(ks_int*r)*boost::math::cyl_bessel_j(n,ks_int*r) ) )*sin(n*theta);
    }
    
    ux = u_3r*cos(theta) - u_3theta*sin(theta);
    uy = u_3r*sin(theta) + u_3theta*cos(theta);
    std::vector<Scalar> res;
    res.push_back(ux);
    res.push_back(uy);
    if(comp<0)return ValType(res);
    else return ValType(res[comp]);

  };
  
  virtual GradType grad(double x, double y, double z) const { return GradType(); };
  virtual HessType hess(double x, double y, double z) const { return HessType(); };
  ~genSimpleFunctionUtr(){ delete polar;}
};

#endif