uniform_lens.h

/*
 * uniform_lens.h
 *
 *  Created on: Jan 21, 2013
 *      Author: D. Leier
 */
 
#ifndef UNIFORM_LENS_H_
#define UNIFORM_LENS_H_
 
#include "base_analens.h"
 
class LensHaloUniform: public LensHalo{
public:
   LensHaloUniform(double zlens  
                  ,double Sigma  
                  ,Point_2d &Shear  
                  ,COSMOLOGY &cosmo
                  );
 
  ~LensHaloUniform();
 
  float getSigma() const {return perturb_modes[0];};
  Point_2d getShear() const{
    Point_2d p(perturb_modes[1],perturb_modes[2]);
    return p;
  }
  
  void resetGamma(const Point_2d &Shear  
                  ){
    perturb_modes[1] = Shear[0];
    perturb_modes[2] = Shear[1];
  }
  
  void resetKappa(double Sigma  
                  ){
    perturb_modes[0] = Sigma;;
  }
  
  void force_halo(PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi,PosType const *xcm,bool subtract_point=false,PosType screening = 1.0);
  
  LensHaloUniform(const LensHaloUniform &h):LensHalo(h){
    perturb_modes[0] = h.perturb_modes[0];
    perturb_modes[1] = h.perturb_modes[1];
    perturb_modes[2] = h.perturb_modes[2];
  }
  
  LensHaloUniform & operator=(const LensHaloUniform &h){
    if(this == &h) return *this;
    LensHalo::operator=(h);
    perturb_modes[0] = h.perturb_modes[0];
    perturb_modes[1] = h.perturb_modes[1];
    perturb_modes[2] = h.perturb_modes[2];
    return *this;
 }
 
protected:
 
  PosType perturb_modes[3];  
  PosType lens_expand(PosType *mod,PosType const *x,PosType *alpha,KappaType *gamma,KappaType *phi);
};
 
class LensHaloKappaDisk: public LensHalo{
public:
   LensHaloKappaDisk(double zlens  
                     ,double Sigma_in  
                     ,double Rmax   
                     ,const Point_2d ang_center  
                     ,COSMOLOGY &cosmo
                     )
  :LensHalo(zlens,cosmo),Sigma(Sigma_in)
  {
  
    if(Rmax <= 0){
      std::cerr << "LensHaloKappaDisk: Rmax cannot be <= 0" << std::endl;
      throw std::invalid_argument("Rmax<=0");
    }
    LensHalo::setMass(PI*Rmax*Rmax*Sigma);
    LensHalo::Rmax = Rmax;
    LensHalo::setRsize(Rmax);
    LensHalo::set_rsize(Rmax);
    LensHalo::setTheta(ang_center);
   }
 
  ~LensHaloKappaDisk(){};
 
  float getSigma() const {return Sigma;};
    
  LensHaloKappaDisk(const LensHaloKappaDisk &h):LensHalo(h){
    Sigma = h.Sigma;
  }
  
  LensHaloKappaDisk & operator=(const LensHaloKappaDisk &h){
    if(this == &h) return *this;
    LensHalo::operator=(h);
    Sigma = h.Sigma;
    return *this;
 }
 
  void force_halo(PosType *alpha,KappaType *kappa,KappaType *gamma,KappaType *phi,PosType const *xcm,bool subtract_point=false,PosType screening = 1.0){
 
      float tmp;
      PosType tmp_mass;
      
      PosType r=sqrt(xcm[0]*xcm[0] + xcm[1]*xcm[1]);
 
      if(r < Rmax){
        alpha[0] -= Sigma*xcm[0];
        alpha[1] -= Sigma*xcm[1];
        *kappa += Sigma;
      }else{
        tmp_mass = mass/PI/r/r; //Sigma * Rmax*Rmax/r/r; 
        alpha[0] -= tmp_mass*xcm[0];
        alpha[1] -= tmp_mass*xcm[1];
 
        tmp = 2*tmp_mass/r/r;
        gamma[0] -= 0.5*(xcm[0]*xcm[0] - xcm[1]*xcm[1])*tmp;
        gamma[1] -= xcm[0]*xcm[1]*tmp;
      }
    
    return;
  }
protected:
  double Sigma; // surface density
};
#endif /* UNIFORM_LENS_H_ */