gravitationallens.h

Go to the documentation of this file.

/*

  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_GRAVITATIONALLENS_H

#define GRALE_GRAVITATIONALLENS_H

#include "graleconfig.h"
#include "vector2d.h"
#include <serut/serializationinterface.h>
#include <string>

namespace grale
{

class GRALE_IMPORTEXPORT GravitationalLensParams : public errut::ErrorBase
{
public:
        GravitationalLensParams()                                                               { }
        virtual ~GravitationalLensParams()                                                      { }

        virtual bool write(serut::SerializationInterface &si) const                                             { setErrorString("Not implemented"); return false; }

        virtual bool read(serut::SerializationInterface &si)                                                    { setErrorString("Not implemented"); return false; }

        virtual GravitationalLensParams *createCopy() const = 0;
};

class GRALE_IMPORTEXPORT GravitationalLens : public errut::ErrorBase
{
public:
        enum LensType 
        {
                Gaussian,               
                MultiplePlummers,       
                Plummer,                
                Pointmass,              
                SIS,                    
                NSIE,                   
                NSIS,                   
                SIE,                    
                Square,                 
                MultipleSquares,        
                MultipleGaussians,      
                MassSheet,              
                Composite,              
                MassDisk,               
                Profile,                
                PolynomialMassProfile,  
                MultipleWendland,       
                DeflectionGrid,         
                NFW,                    
                EllipticNFW,            
                Sersic,                 
                EllipticSersic,         
                TestLensType,
                MaxLensType
        };
protected:
        GravitationalLens(LensType t);
public:
        ~GravitationalLens();

        LensType getLensType() const                                                            { return m_lensType; }
        
        virtual bool init(double D_d, const GravitationalLensParams *pLensParams = 0);

        bool traceTheta(double D_s, double D_ds, Vector2D<double> theta, Vector2D<double> *pBeta) const;

        virtual bool getAlphaVector(Vector2D<double> theta, Vector2D<double> *pAlpha) const = 0;

        virtual bool getAlphaVectorDerivatives(Vector2D<double> theta, double &axx, 
                                               double &ayy, double &axy) const;

        static double getInverseMagnification(double D_s, double D_ds, double axx, 
                                              double ayy, double axy)                           { double f = D_ds/D_s; return (1.0-f*axx)*(1.0-f*ayy)-f*f*axy*axy; }

        static void getShearInfo(double D_s, double D_ds, double axx, double ayy, 
                                 double axy, double *pShearAngle, double *pShearSize)           { double f = D_ds/D_s; double g1 = 0.5*(axx-ayy); double g2 = axy; Vector2D<double> gVector(g1, g2); double g = gVector.getLength(); if (g2 >= 0) *pShearAngle = 0.5*ACOS(g1/g); else *pShearAngle = -0.5*ACOS(g1/g); *pShearSize = g*f; }

        static void getShearInfo(double gamma1, double gamma2, 
                                 double *pShearAngle, double *pShearSize)                       { Vector2D<double> gVector(gamma1, gamma2); double g = gVector.getLength(); if (gamma2 >= 0) *pShearAngle = 0.5*ACOS(gamma1/g); else *pShearAngle = -0.5*ACOS(gamma1/g); *pShearSize = g; }


        virtual double getInverseMagnification(double D_s, double D_ds, Vector2D<double> theta) const;

        virtual double getSurfaceMassDensity(Vector2D<double> theta) const = 0;
        
        double getLensDistance() const                                                          { return m_Dd; }

        void setLensDistance(double D_d)                                                        { m_Dd = D_d; }

        bool getTimeDelay(double z_d, double D_s, double D_ds, Vector2D<double> theta, 
                          Vector2D<double> beta, double *pTimeDelay) const;
        virtual bool getProjectedPotential(double D_s, double D_ds, Vector2D<double> theta, 
                                           double *pPotentialValue) const;
        
        const GravitationalLensParams *getLensParameters() const                                { return m_pParameters; }

        virtual void setDerivativeAngularDistanceScale(double distanceScale)                    { m_derivDistScale = ABS(distanceScale); }
        
        bool write(serut::SerializationInterface &si) const;

        static bool read(serut::SerializationInterface &si, GravitationalLens **pLens, std::string &errorString);

        bool save(const std::string &fileName) const;

        static bool load(const std::string &fileName, GravitationalLens **pLens, std::string &errorString);

        GravitationalLens *createCopy() const;

        // TODO: experimental
        // OpenCL stuff
        
        virtual bool getSuggestedScales(double *pDeflectionScale, double *pPotentialScale) const;
        virtual bool getCLParameterCounts(int *pNumIntParams, int *pNumFloatParams) const;
        virtual bool getCLParameters(double deflectionScale, double potentialScale, int *pIntParams, float *pFloatParams) const;
        
        std::string getCLLensProgram(std::string &subRoutineName) const;

        // like this:
        // LensQuantities subRoutineName(float2 coord, __global const int *pIntParams, __global const float *pFloatParams)
        //
        // typedef struct
        // {
        //      float alphaX;
        //      float alphaY;
        //      float potential;
        //      float axx;
        //      float ayy;
        //      float axy;
        // } LensQuantities;
        virtual std::string getCLProgram(std::string &subRoutineName) const;

        virtual int getCLSubLenses() const                                                              { return -1; }
protected:
        virtual bool processParameters(const GravitationalLensParams *params) = 0;

        std::string getCLLensQuantitiesStructure() const;
private:
        bool m_init;
        double m_Dd;
        LensType m_lensType;
        GravitationalLensParams *m_pParameters;
        double m_derivDistScale;
};

} // end namespace

#endif // GRALE_GRAVITATIONALLENS_H