CshapeFunction_global_brasilian.h

Go to the documentation of this file.

/**********************************************************************
 * Copyright (C) 2012, The YaDICs Project Developers. 
 * See the COPYRIGHT file at the top-level directory of this distribution ./COPYRIGHT.
 * See ./COPYING file for copying and redistribution conditions.
 * 
 * This file is part of YaDICs.
 * YaDICs is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 * 
 * YaDICs is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with YaDICs.  If not, see <http://www.gnu.org/licenses/>.
 * 
 * Information about how to use the software are provided at http://yadic.univ-lille1.fr/
 **********************************************************************/


#ifndef CSHAPEFUNCTION_GLOBAL_BRASILIAN
#define CSHAPEFUNCTION_GLOBAL_BRASILIAN

#include "CshapeFunction_global_homogeneous.h"

template<typename T,typename Timg> 
class CshapeFunction_global_brasilian : public CshapeFunction_global_homogeneous<T,Timg>
{
  
  public:

    std::vector<int> m_modB;//Required Brasilian mode
    double m_nu, m_r;//Poisson, radius
    std::vector<double> m_f;//force
    double m_k0, m_k1;//parameters function of m_nu, m_f and m_r
    
    CshapeFunction_global_brasilian() : CshapeFunction_global_homogeneous<T,Timg>()
    {// constructor

      this->class_name = "CshapeFunction_global_brasilian";
      
    }

    virtual void init(CParameterNetCDF &fp)
    {//fill m_modB vector with the indice of required Brasilian mode from parameter.nc
    
      //look for the Homogeneous Deformation set of required modes.
      CshapeFunction_global_homogeneous<T,Timg>::init(fp);
      
      //identify the Brasilian mode withins this->m_mods[] and store it within m_modB[]
      m_modB.clear();
      for (int i=0; i<this->m_mods.size(); i++){if (this->m_mods[i].compare(0,1,"B") == 0){m_modB.push_back(i);}}
      if (this->m_verbose){std::cout<<"\tmodes B : [ ";for (int i=0; i<m_modB.size(); i++){std::cout<<m_modB[i]<<" ";}std::cout<<"]\n";}
      
      std::string att_name("Poisson");
      fp.loadAttribute(att_name,m_nu);
      if (this->m_verbose){std::cout<<"\tPoisson factor : "<<m_nu<<"\n";}
      
      att_name = "radius";
      fp.loadAttribute(att_name,m_r);
      if (this->m_verbose){std::cout<<"\tCircle radius : "<<m_r<<"\n";}
      
      att_name = "force";
      fp.loadAttribute(att_name,m_f);
      if (this->m_verbose){std::cout<<"\tforce : "<<m_f[frame_number]<<"\n";}
      
      m_k0 = (m_f[frame_number]*m_r*(1+m_nu))/(2*this->m_pi);
      m_k1 = (3-m_nu)/(1+m_nu);
      
    }
 
    virtual void exec(const CImg<Timg> &box, const std::vector<T> &x, CImgList<T> &N)
    {// assign CImgList<T> N[dimensions][modes] for required Brasilian and HD modes at coordonate (x,y,z)
      
      //fill HD modes
      CshapeFunction_global_homogeneous<T,Timg>::exec(box, x, N);
      
      double a, b, xx, yy;
      
      //90° rotation of the reference system (Z* = Z x i)
      xx=-x[1];
      yy=x[0];
      
      a = m_r + xx;
      b = m_r - xx;
      
      if ((a==0 or b==0) and yy==0)
      {// two points are singular : (-m_r,0) and (m_r,0)
        if (this->m_verbose){printf("\nwarning singular points\n");}
      }
      else
      {
        //fill Brasilian mode
        for (int mode=0; mode<m_modB.size(); mode++)
        {
          N[1][m_modB[mode]] = m_k0*( m_k1*(a/m_r - std::log(std::sqrt(std::pow(a,2)+std::pow(yy,2))/std::sqrt(std::pow(b,2)+std::pow(yy,2)))) + std::log(std::sqrt(std::pow(a,2)+std::pow(yy,2))/std::sqrt(std::pow(b,2)+std::pow(yy,2))) + (std::pow(a,2)-std::pow(yy,2))/(std::pow(a,2)+std::pow(yy,2)) + (std::pow(b,2)-std::pow(yy,2))/(std::pow(b,2)+std::pow(yy,2)) - xx/m_r);

          N[0][m_modB[mode]] = m_k0*( m_k1*(yy/m_r - std::atan2(yy,b) - std::atan2(yy,a)) - std::atan2(yy,b) + (2*yy*a)/(std::pow(a,2)+std::pow(yy,2)) - (2*yy*b)/(std::pow(b,2)+std::pow(yy,2)) - yy/m_r - std::atan2(yy,a));
        }
        
      }
      
    }
       
};

#endif