monopolebasisfunction.h

/*

  This file is a part of GRALE, a library to facilitate the simulation
  and inversion of gravitational lenses.

  Copyright (C) 2008-2012 Jori Liesenborgs

  Contact: jori.liesenborgs@gmail.com
  
  This program 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 2 of the License, or
  (at your option) any later version.
  
  This program 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 this program; if not, write to the Free Software
  Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  
*/

#ifndef GRALE_MONOPOLEBASISFUNCTION_H

#define GRALE_MONOPOLEBASISFUNCTION_H

#include "graleconfig.h"
#include "real2dderivablefunction.h"

namespace grale
{

class MonopoleBasisFunction : public Real2DDerivableFunction
{
public:
        MonopoleBasisFunction(double zeroPoint, double lengthScale)
        { 
                m_zeroPoint = zeroPoint; 
                m_lengthScale = lengthScale; 
        
                m_a2 = -1.0/(m_zeroPoint*m_zeroPoint);
                m_a0 = 1.0;

                double m1 = m_zeroPoint-1;
                double factor = (m_zeroPoint*m_zeroPoint)/4.0*(1.0/(m1*m1*m1));
                
                m_b1 = factor*(-6.0);
                m_b0 = factor*6.0*(1.0+m_zeroPoint);
                m_b_1 = factor*(-6.0)*m_zeroPoint;
                
        }
        
        ~MonopoleBasisFunction()
        {
        }

        double operator()(Vector2D<double> point) const
        {
                Vector2D<double> p = point/m_lengthScale;
                double x = p.getLength();
                double val = 0;

                if (x < m_zeroPoint)
                {
                        double x2 = x*x;
                        
                        val = m_a2*x2 + m_a0;
                }
                else if (x < 1.0)
                {
                        val = m_b1*x + m_b0 + m_b_1/x;
                }

                return val;
        }

        Vector2D<double> getGradient(Vector2D<double> point) const
        {
                Vector2D<double> p = point/m_lengthScale;
                double x = p.getLength();
                double val = 0;

                if (x < m_zeroPoint)
                {
                        val = 2.0*m_a2*x;
                }
                else if (x < 1.0)
                {
                        double x2 = x*x;
                        
                        val = m_b1-m_b_1/x2;
                }

                if (val != 0)
                        val /= (x*m_lengthScale);
                return val*p;
        }

        double getLengthScale() const
        {
                return m_lengthScale;
        }

        double getZeroPoint() const
        {
                return m_zeroPoint;
        }
private:
        double m_lengthScale;
        double m_zeroPoint;
        double m_a2, m_a0;
        double m_b1, m_b0, m_b_1;
};
        
} // end namespace

#endif // GRALE_MONOPOLEBASISFUNCTION_H