11 #include "utilities_slsim.h"
17 double mag_to_jansky_AB(
double);
18 double jansky_to_mag_AB(
double flux);
21 double mag_to_flux_AB(
double);
23 double flux_to_mag_AB(
double flux);
26 T flux_to_mag(T flux,T zeropoint){
27 if(flux <=0)
return 100;
28 return -2.5 * log10(flux) + zeropoint;
32 T mag_to_flux(T m,T zeropoint){
33 if(m == 100)
return 0;
34 return pow(10,-0.4*(m - zeropoint));
68 mag_zero_point(s.mag_zero_point){ }
71 if(
this == &s)
return *
this;
75 sb_limit = s.sb_limit;
77 mag_zero_point = s.mag_zero_point;
88 virtual PosType getTotalFlux()
const = 0;
89 virtual void printSource() = 0;
92 double SBlimit_magarcsec(
double limit) {
return mag_to_flux(limit,mag_zero_point)*pow(180*60*60/PI,2);}
102 void setZ(PosType my_z){
zsource = my_z;}
122 void setMagZeroPoint(
float zeropoint){mag_zero_point=zeropoint;}
123 double getMagZeroPoint(){
return mag_zero_point;}
125 PosType
changeFilter(std::string filter_in, std::string filter_out, std::string sed);
126 PosType
integrateFilter(std::vector<PosType> wavel_fil, std::vector<PosType> fil);
127 PosType
integrateFilterSED(std::vector<PosType> wavel_fil, std::vector<PosType> fil, std::vector<PosType> wavel_sed, std::vector<PosType> sed);
137 double total_flux = 0;
139 for(
size_t i=0 ; i<N ; ++i){
140 for(
size_t j=0 ; j<N ; ++j){
147 total_flux *= res*res;
149 std::cout <<
"------- check flux ----------" << std::endl;
150 std::cout <<
" total flux in pixels : " << total_flux << std::endl;
151 std::cout <<
" getTotalFlux() : " << getTotalFlux() << std::endl;
152 std::cout <<
" fractional difference : " << (total_flux - getTotalFlux()) / getTotalFlux() << std::endl;
153 std::cout <<
" magnitude from flux in pixels : " << flux_to_mag(total_flux) << std::endl;
154 std::cout <<
" mag from getTotalFlux() : " << flux_to_mag(getTotalFlux()) << std::endl;
155 std::cout <<
"-----------------------------" << std::endl;
157 if(abs( (total_flux - getTotalFlux()) / getTotalFlux() ) > 0.1 ){
164 long getID()
const {
return id;}
165 void setID(
long i){
id=i;}
180 double flux_to_mag(
double flux)
const {
181 if(flux <=0)
return 100;
182 return -2.5 * log10(flux) + mag_zero_point;
191 double mag_zero_point;
197 class SourceColored :
public Source{
199 SourceColored(PosType magnitude,PosType r,
Point_2d x,PosType z
200 ,PosType SBlimit,
double zero_point)
201 :
Source(r,x,z,SBlimit,zero_point)
203 zeropoints[Band::NoBand] = zero_point;
204 current_band = Band::NoBand;
206 mag_map[current_band]=magnitude;
207 flux_total = mag_to_flux(mag,zero_point);
211 SourceColored(
const SourceColored &s):
Source(s){
214 current_band = s.current_band;
215 sed_type = s.sed_type;
216 flux_total = s.flux_total;
219 SourceColored & operator= (
const SourceColored &s){
220 if(
this == &s)
return *
this;
222 Source::operator=(s);
225 current_band = s.current_band;
226 sed_type = s.sed_type;
227 flux_total = s.flux_total;
233 SourceColored & copy_color(
const SourceColored &s){
234 if(
this == &s)
return *
this;
236 Source::operator=(s);
239 current_band = s.current_band;
240 sed_type = s.sed_type;
241 flux_total = s.flux_total;
246 virtual ~SourceColored(){};
248 PosType getMag()
const { assert(current_band != Band::NoBand);
return mag;}
249 PosType getMag(
Band band)
const {
return (mag_map.size() > 0) ? mag_map.at(band) : mag;}
250 Band getBand()
const{
return current_band;}
251 float getSEDtype()
const {
return sed_type;}
252 double getMagZeroPoint(
Band band)
const {
return zeropoints.at(band);}
254 void setSEDtype(
float sed) {sed_type = sed;}
255 static void setMagZeroPoint(
Band band,
double zeropoint){zeropoints[band]=zeropoint;}
256 static void setMagZeroPoints(std::map<Band,PosType> &zero_points){zeropoints=zero_points;}
258 void setActiveBand(
Band band);
260 inline PosType getTotalFlux()
const {
return flux_total;}
261 void setMag(
float magnitude,
Band band,
double zeropoint){
262 mag_map[band]=magnitude;
263 zeropoints[band]=zeropoint;
264 if(band==current_band) flux_total = mag_to_flux(mag,zeropoint);
268 virtual void rotate(PosType theta) = 0;
271 void shiftmag(
float delta_mag){
273 flux_total = mag_to_flux(mag,zeropoints[current_band]);
274 for(
auto &b : mag_map){
275 b.second += delta_mag;
281 static std::map<Band,PosType> zeropoints;
288 std::map<Band,PosType> mag_map;
300 SourcePixelled(PosType my_z, PosType* center,
int Npixels, PosType resolution, PosType* arr_val,PosType zero_point);
307 inline PosType getTotalFlux()
const {
return flux;}
308 inline PosType getRadius()
const {
return source_r;}
309 inline PosType* getEll(){
return ell;};
310 inline PosType getQuad(
int i,
int j){
return quad[i][j];};
311 inline PosType getSize(){
return size;};
312 inline PosType* getCentroid(){
return centroid;};
313 inline PosType getMag(){
return flux_to_mag(flux);};
319 void calcTotalFlux();
328 std::valarray<PosType> values;
345 SourceShapelets(PosType my_z, PosType my_mag, PosType my_scale, std::valarray<PosType> my_coeff, PosType* my_center, PosType my_ang,PosType zero_point);
346 SourceShapelets(PosType my_z, PosType my_mag, std::string shap_file, PosType *my_center, PosType my_ang, PosType zero_point);
347 SourceShapelets(std::string shap_file, PosType my_ang, PosType zero_point);
356 cos_sin[0] = s.cos_sin[0];
357 cos_sin[1] = s.cos_sin[1];
358 coeff_flux = s.coeff_flux;
363 if(
this == &s)
return *
this;
365 SourceColored::operator=(s);
370 cos_sin[0] = s.cos_sin[0];
371 cos_sin[1] = s.cos_sin[1];
372 coeff_flux = s.coeff_flux;
380 cos_sin[0] = cos(ang);
381 cos_sin[1] = sin(ang);
395 for(
int i=0 ; i < n ; ++i){
397 double theta = 2*PI*ran();
398 double ff = f*( t - t2 * ran() );
400 Point_2d x( r*cos(theta) , r*sin(theta) );
402 grains.emplace_back(ff,s,x);
406 double sub_structure_sb(
double *x){
409 for(
auto &g : grains){
418 PepperCorn(
double f,
double sig,
Point_2d &x):
419 x(x),sigma(sig),fo(f){};
421 double operator()(
double *y){
422 return fo * exp(-( (x[0]-y[0])*(x[0]-y[0]) - (x[1]-y[1])*(x[1]-y[1]) ) /sigma/sigma);
435 void Hermite(std::vector<PosType> &hg,
int N, PosType x)
const;
437 void NormalizeFlux();
438 std::valarray<PosType> coeff;
445 std::vector<SourceShapelets::PepperCorn> grains;
456 ,PosType radius_in_radians
474 ,PosType radius_in_radians
477 ):SourceColored(magnitude,radius_in_radians,position,z,SBlimit,zero_point)
485 void printSource(){};
487 void rotate(PosType t){};
494 class SourceGaussian :
public Source{
504 Source(0,position,z,SBlimit,zero_point),source_gauss_r2(r_size*r_size)
514 PosType source_gauss_r2;
519 PosType getTotalFlux()
const {
return 2*PI*source_gauss_r2;}
529 PosType getTotalFlux()
const {std::cout <<
"No total flux in SourceBLR yet" << std::endl; exit(1);}
531 virtual inline PosType getRadius()
const {
return source_r_out;}
544 float source_opening_angle;
552 inline PosType getDlDs()
const{
return DlDs;}
564 PosType getTotalFlux()
const {std::cout <<
"No total flux in SourceBLRDisk yet" << std::endl; exit(1);}
574 PosType getTotalFlux()
const {std::cout <<
"No total flux in SourceBLRSph1 yet" << std::endl; exit(1);}
584 PosType getTotalFlux()
const {std::cout <<
"No total flux in SourceBLRSph2 yet" << std::endl; exit(1);}
605 PosType getNorm() {
return pow(10,log_phi)*norm;};
625 PosType ave_colors[4];
626 PosType color_dev[4];
628 std::string kcorr_file, colors_file;
639 LF_kernel(PosType alpha,PosType beta,PosType mstar)
640 : alpha(alpha),beta(beta),mstar(mstar){};
646 double operator () (
double mag) {
647 return 1.0/( pow(10,0.4*(alpha+1)*(mag-mstar)) + pow(10,0.4*(beta+1)*(mag-mstar)) );
658 double operator()(
double x,
double y)
661 double z[2] = {x, y};