// 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_CASE1_H_
#define _GEN_SIMPLE_FUNCTION_CASE1_H_


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

template<class T1,class T2> class genSimpleFunctionCase1 : 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, ks;
  int comp;
public :
  genSimpleFunctionCase1(std::vector<std::vector<Scalar> >& Coefficients_,double &kp_,double &ks_, int comp_ = -1) : Coefficients(Coefficients_) , kp(kp_), ks(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_1r, u_1theta, u_2r, u_2theta, ux, uy;    
    for(int n=0; n< Coefficients.size(); n++){
      if(n==0){
      u_1r += -Coefficients[n][0]*kp*boost::math::cyl_hankel_1(n+1,kp*r)*cos(n*theta);
      
      u_2r += -Coefficients[n][2]*kp*boost::math::cyl_hankel_2(n+1,kp*r)*cos(n*theta);
      
      u_1theta += Coefficients[n][1]*ks*boost::math::cyl_hankel_1(n+1,ks*r)*sin(n*theta);
      
      u_2theta += Coefficients[n][3]*ks*boost::math::cyl_hankel_2(n+1,ks*r)*sin(n*theta);
      }
      else{
      u_1r += (Coefficients[n][0]*kp*( -boost::math::cyl_hankel_1(n+1,kp*r) + n/(kp*r)*boost::math::cyl_hankel_1(n,kp*r) )+
	      n/r*Coefficients[n][1]*boost::math::cyl_hankel_1(n,ks*r) )*cos(n*theta);
      
      u_2r += (Coefficients[n][2]*kp*( boost::math::cyl_hankel_2(n-1,kp*r) - n/(kp*r)*boost::math::cyl_hankel_2(n,kp*r) )+
	      n/r*Coefficients[n][3]*boost::math::cyl_hankel_2(n,ks*r) )*cos(n*theta);
      
      u_1theta += (-n/r*Coefficients[n][0]*boost::math::cyl_hankel_1(n, kp*r) - 
		Coefficients[n][1]*ks*(-boost::math::cyl_hankel_1(n+1,ks*r) + n/(ks*r)*boost::math::cyl_hankel_1(n,ks*r) ) )*sin(n*theta);
		
      u_2theta += (-n/r*Coefficients[n][2]*boost::math::cyl_hankel_2(n, kp*r) - 
		Coefficients[n][3]*ks*(boost::math::cyl_hankel_2(n-1,ks*r) - n/(ks*r)*boost::math::cyl_hankel_2(n,ks*r) ) )*sin(n*theta);
      }
    }
    
    ux = (u_1r + u_2r)*cos(theta) - (u_1theta + u_2theta)*sin(theta);
    uy = (u_1r + u_2r)*sin(theta) + (u_1theta + u_2theta)*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(); };
  ~genSimpleFunctionCase1(){ delete polar;}
};

#endif // _GEN_SIMPLE_FUNCTION_CASE1_H_