// 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_UINC_H_
#define _GEN_SIMPLE_FUNCTION_UINC_H_


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

template<class T1,class T2> class genSimpleFunctionUinc : 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;
    
  genSimpleFunction<T2> *polar;
  double kp;
  std::vector<double>::size_type size;
  int comp;
public :
  genSimpleFunctionUinc( std::vector<double>::size_type size_, double &kp_, int comp_ =-1) : kp(kp_), size(size_) , 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_r, u_theta, ux, uy;  
    std::complex<double> i(0,1);
    for(int n=0; n< (int)size; n++){
      if(n==0){
	    u_r += pow(i,n)*kp*( -boost::math::cyl_bessel_j(n+1,kp*r) + n/(kp*r)*boost::math::cyl_bessel_j(n,kp*r) )*cos(n*theta);
      
	    u_theta += -n/r*pow(i,n)*boost::math::cyl_bessel_j(n,kp*r)*sin(n*theta);
      }
      else{
            u_r += 2.0*pow(i,n)*kp*( -boost::math::cyl_bessel_j(n+1,kp*r) + n/(kp*r)*boost::math::cyl_bessel_j(n,kp*r) )*cos(n*theta);
      
	    u_theta += -2.0*n/r*pow(i,n)*boost::math::cyl_bessel_j(n,kp*r)*sin(n*theta);
      }
      
    }
    
    
    
    ux = u_r*cos(theta) - u_theta*sin(theta);
    uy = u_r*sin(theta) + u_theta*cos(theta);
    // */
     /*
    Scalar ux, uy;
    double ampl = 0.00628318530717959;
    ux = -ampl*exp(std::complex<double>(0, -1)*(kp*x + M_PI/2) );
    uy = 0;
     */
    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(); };
  ~genSimpleFunctionUinc(){ delete polar;}
};



#endif // _GEN_SIMPLE_FUNCTION_UINC_H_