-
Andrew Strong authoredEroBubbles_for_gitlab.cc EroBubbles_for_gitlab.h bubbles_model_for_gitlab.cc latest erobubbles model with apec emissivities
Andrew Strong authoredEroBubbles_for_gitlab.cc EroBubbles_for_gitlab.h bubbles_model_for_gitlab.cc latest erobubbles model with apec emissivities
EroBubbles_for_gitlab.cc 14.68 KiB
using namespace std;
#include<cstdio>
#include<string>
#include<fstream>
#include<iostream>
#include<cmath>
#include<valarray> // AWS20221031
#include"EroBubbles.h"
#include "healpix_map.h"
#include "healpix_map_fitsio.h"
#include "arr.h"
#include "fitshandle.h"
#include "pointing.h"
double apec_emissivity( double density, double T, double Emin,double Emax, int debug); // temporary until put in class
///////////////////////////////////////////////////////////////
int EroBubbles::init(double r_minN_, double r_maxN_,
double r_minS_, double r_maxS_,
double xcN_, double ycN_, double zcN_,
double xcS_, double ycS_, double zcS_,
string RLH_,
double ds_, double s_max_,
double R0_)
{
cout<<"EroBubbles.init(r_min_,r_max_)"<<endl;
r_minN=r_minN_;
r_maxN=r_maxN_;
r_minS=r_minS_;
r_maxS=r_maxS_;
xcN=xcN_;
ycN=ycN_;
zcN=zcN_;
xcS=xcS_;
ycS=ycS_;
zcS=zcS_;
RLH=RLH_;
ds = ds_;
s_max = s_max_;
R0=R0_;
return 0;
};
////////////////////////////////////////////////////////////////////
double EroBubbles::emissivity(double x, double y, double z, double E)
{
int debug=0;
if(debug==1)
cout<<"EroBubbles.emissivity(x,y,z,E)"<<endl;
double emissivity_;
double dx, dy, dz;
if(z>=0.)
{
dx=x-xcN;
dy=y-ycN;
dz=z-zcN;
}
if(z<0.)
{
dx=x-xcS;
dy=y-ycS;
dz=z-zcS;
}
double r=sqrt(dx*dx+dy*dy+dz*dz);
emissivity_=0;
double emissivity_nominal;
emissivity_nominal = 1.0e-29 ; // erg cm-3 sr-1 s-1 band 022 from notes
double density = 1e-2;
density = 1e-3; // AWS20221102
double T=5e6;
double Emin=0.2;
double Emax=0.7;
Emin=0.2; //AWS20221101
Emax=0.7; //AWS20221101
int use_apec=1;
debug=0; // debug emissivity here top level was 1
double apec_emissivity_;
// apec_emissivity_=apec_emissivity(density, T, Emin, Emax, debug);
if(use_apec==0)
emissivity_nominal= thermal_emissivity( density, T, Emin, Emax, debug);
if(use_apec==1)
emissivity_nominal= apec_emissivity( density, T, Emin, Emax, debug);
if(z>=0. && r>r_minN && r<r_maxN)emissivity_ = emissivity_nominal ;
if(z< 0. && r>r_minS && r<r_maxS)emissivity_ = emissivity_nominal ;
if(debug==1)
cout<<"r="<<r<<" new emissivity_="<<emissivity_<<endl;
return emissivity_ ;
};
/////////////////////////////////////////////////////////////
double EroBubbles::intensity (double l, double b, double E)
{
int debug=0;
if(debug==1)
cout<<"EroBubbles.intensity(l,b,E)"<<endl;
double intensity_;
double x,y,z;
double s; // line-of-sight distance
double dtr=acos(-1.0)/180.; // degrees to radians pi/180
// RLH: right or left-handed system
// RHS z=x^y y decreases with l
// LHS z=y^x y increases with l
intensity_=0.;
s=0.;
for(s=0; s<=s_max;s+=ds)
{
x = R0 - s * cos(b*dtr) * cos(l*dtr);
if(RLH=="R")
y = - s * cos(b*dtr) * sin(l*dtr);
if(RLH=="L")
y = s * cos(b*dtr) * sin(l*dtr);
z = s * sin(b*dtr);
intensity_+=emissivity(x,y,z,E);
debug=0;
if(debug==1)
cout<<"l="<<l<<" b="<<b<<" s="<<s<<" RLH="<<RLH<<" x="<<x<<" y= "<<y<<" z="<<z
<<" E="<<E<<" intensity_= "<<intensity_<<endl;
}
intensity_*=ds;
double kpctocm;
kpctocm=3.08e21;
intensity_*=kpctocm;
double arcmin2_in_sr=1./pow(57.3*60.,2);
// intensity_*=arcmin2_in_sr;
debug=0;
if(debug==1)
cout<<"l="<<l<<" b="<<b<<" s="<<s<<" RLH="<<RLH<<" x="<<x<<" y= "<<y<<" z="<<z
<<" E="<<E<<" intensity_= "<<intensity_<<" cm-2 sr-1 s-1"<<endl;
return intensity_ ;
};
//////////////////////////////////////////////////////////////////
int EroBubbles::healpix_skymap(int order, string outfile, double E, int debug)
{
cout<<endl<< "=== healpix_skymap"<<endl<<endl;
PDT datatype = PLANCK_FLOAT64; // HealPix defines it this way
datatype = PLANCK_FLOAT32;
arr<double>data; // for HealPix array class
int hdu;
int ncolnum,colnum;
//====================
ncolnum=1;
cout<<"writing healpix skymap to "<<outfile<<endl;
fitshandle out;
out.create("!"+outfile); // ! to overwrite
arr<string> colname;
colname.alloc(ncolnum);
for(int j=0;j<ncolnum;j++) colname[j]="bubbles";
Healpix_Map<double> map_RING_out(order,RING);
prepare_Healpix_fitsmap(out,map_RING_out, datatype, colname);
out.set_key(string("EXTNAME"),string("SKYMAP"));// extension HDU name
colnum=1;
cout<< "map_RING_out: ";
cout<<" Npix = "<<map_RING_out.Npix();
cout<<" Nside = "<<map_RING_out.Nside();
cout<<" Order = "<<map_RING_out.Order();
cout<<" Scheme= "<<map_RING_out.Scheme()<<endl; // 0 = RING, 1 = NESTED
data.alloc(map_RING_out.Npix());
pointing pointing_;
double l,b;
double rtd=180./acos(-1.); // radians to degrees 180/pi
E=3; // should be in init
for (int ipix=0; ipix<map_RING_out.Npix(); ipix++)
{
pointing_ = map_RING_out.pix2ang(ipix);
l = pointing_.phi *rtd;
b = 90.-pointing_.theta*rtd;
data[ipix] = intensity(l,b,E);
if(debug==1)
cout<<"ipix="<<ipix<<" theta="<<pointing_.theta<<" phi="<<pointing_.phi
<<" l="<<l<<" b="<<b<<" data[ipix]= "<<data[ipix]<<endl;
}
out.write_column(colnum,data); // fitshandle.h
out.close();
// for future use
map_RING_out.Set(data,RING); // data is zero size after call!
for (int ipix=0; ipix<map_RING_out.Npix(); ipix++)
{
pointing_ = map_RING_out.pix2ang(ipix);
l = pointing_.phi *rtd;
b = 90.-pointing_.theta*rtd;
if(debug==1)
cout<<"ipix="<< ipix<<" theta="<<pointing_.theta<<" phi="<<pointing_.phi
<<" l="<<l<<" b="<<b<<" map_RING_out[ipix]= "<<map_RING_out[ipix]<<endl;
}
//////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////
cout<<endl<<" ==== healpix_skymap complete"<<endl;
return 0;
};
///////////////////////////////////////////////////////////////////////////
double EroBubbles:: thermal_emission( double em, double T, double E)
{
// thermal formula:
//https://wiki.mpe.mpg.de/eRosita/erosnr_project_Fermi_bubbles?action=AttachFile&do=view&target=erosita_garching_2020_FermiBubbles.pdf
// Longair eq 6.47
// thermal formula:
//https://wiki.mpe.mpg.de/eRosita/erosnr_project_Fermi_bubbles?action=AttachFile&do=view&target=erosita_garching_2020_FermiBubbles.pdf
// Longair eq 6.47
double g=1.2; // Gaunt factor
double k=1.38e-16 ; // erg/K
double keV_to_erg = 1.602e-9;
//
double thermal_emission_= .004 * em * pow(T,-0.5) * g * exp(-(E*keV_to_erg)/(k*T)) ;// cf pdf above. erg cm-2 sr-1 s-1 keV-1
cout<<"thermal_emission_="<<thermal_emission_<<" erg cm-2 sr-1 s-1 keV-1"<<endl;
return thermal_emission_;
};
//////////////////////////////////////////////////////////////////////
double EroBubbles:: thermal_emissivity( double density, double T, double Emin,double Emax, int debug)
{
// thermal formula:
//https://wiki.mpe.mpg.de/eRosita/erosnr_project_Fermi_bubbles?action=AttachFile&do=view&target=erosita_garching_2020_FermiBubbles.pdf
// Longair eq 6.47
// thermal formula:
//https://wiki.mpe.mpg.de/eRosita/erosnr_project_Fermi_bubbles?action=AttachFile&do=view&target=erosita_garching_2020_FermiBubbles.pdf
// Longair eq 6.47
double g=1.2; // Gaunt factor
double k=1.38e-16 ; // erg/K
double keV_to_erg = 1.602e-9;
//
// double thermal_emission_= .004 * em * pow(T,-0.5) * g * exp(-(E*keV_to_erg)/(k*T)) ;// cf pdf above. erg cm-2 sr-1 s-1 keV-1
// cout<<"thermal_emission_="<<thermal_emission_<<" erg cm-2 sr-1 s-1 keV-1"<<endl;
// want this for calculation: emissivity_nominal = 1.0e-29 ; // erg cm-3 sr-1 s-1 band 022 from notes
// em in cm-6 pc
// density^2 in cm-6
// density=1e-2 em=1e-4/3e18 = 1e-22 T=1e7 pow(T,.5)=1e-4, * .004 about right
double pc_to_cm= 3.08e18;
double em = pow(density,2) / pc_to_cm; // em for 1 cm
double E=(Emin+Emax)/2.;
double dE=Emax-Emin;
double thermal_emissivity_= .004 * em * pow(T,-0.5) * g * exp(-(E*keV_to_erg)/(k*T)) ;// cf pdf above. erg cm-2 sr-1 s-1 keV-1
if(debug==1) cout<<"thermal_emissivity_="<<thermal_emissivity_<<" erg cm-3 sr-1 s-1 keV-1"<<endl;
thermal_emissivity_ /= (E * keV_to_erg);
if(debug==1) cout<<"thermal_emissivity_="<<thermal_emissivity_<<" counts cm-3 sr-1 s-1 keV-1"<<endl;
thermal_emissivity_ *= dE;
if(debug==1) cout<<"thermal_emissivity_="<<thermal_emissivity_<<" counts cm-3 sr-1 s-1 "<<endl;
return thermal_emissivity_;
};
//////////////////////////////////////////////////////////////////////////////////
// copied from healpix_spectra.cc 20221031
double apec_emission( double em, double T, double E)
{
// thermal formula:
//https://wiki.mpe.mpg.de/eRosita/erosnr_project_Fermi_bubbles?action=AttachFile&do=view&target=erosita_garching_2020_FermiBubbles.pdf
// Longair eq 6.47
double g=1.2; // Gaunt factor
double k=1.38e-16 ; // erg/K
double keV_to_erg = 1.602e-9;
//
double thermal_emission_= .004 * em * pow(T,-0.5) * g * exp(-(E*keV_to_erg)/(k*T)) ;// cf pdf above. erg cm-2 sr-1 s-1 keV-1
// test.py:
// energies = np.linspace(0.1, 5., 100) # Set up the energy grid
int debug=0; // added to copy AWS20221101
double E1,E2,dE;
E1=0.1;
E2=5.0;
valarray<double> apec_values;
int n_values;
n_values=100;
dE=(E2-E1)/n_values;
apec_values.resize(n_values);
apec_values={1.,2.,3.};
int i;
apec_values={
20.098434984330172,13.513800017054615,7.2605190309170995,6.265155768700881,3.415934874319297,2.8441181827005724,2.0266453286658224,1.3495033990686836,1.473295435632978,1.1689807744063774,1.1793142433854058,4.125035543114648,1.072060096619616,11.301010266914872,5.445650296844631,6.818022915893874,3.0880991159422404,1.3562764501191384,1.0525801864184765,1.590112434960223,0.8985661943799107,0.37893758826527524,0.43732476662266145,0.14835402675438417,0.13756979793308732,0.5228649284767455,0.06788818481564048,0.06502793726412887,0.12639937381439864,0.04681827517378769,0.08810081634871758,0.03464016858985263,0.04160740098681964,0.03916258313054317,0.02387169125703757,0.10361007225473656,0.11361138885477981,0.017174590070385142,0.01700271883867702,0.015517489344494276,0.01139989107656388,0.013910308697404903,0.015659000939298368,0.009139725951549242,0.00965226018809285,0.0071298164690791555,0.006497145360879474,0.01776515246052341,0.01624512821815538,0.004122611717192725,0.0036429515083825884,0.0032697413305504868,0.002845266143416478,0.002536844326253808,0.002507274033784574,0.002699265210896621,0.0023135039500451807,0.001589935917095757,0.0016903226227150227,0.0015372892775794651,0.001165993359260682,0.00244198674274096,0.0009268697876071981,0.0008097545405433635,0.0007276756001221772,0.0006447669878802399,0.0005772577382388693,0.0005157969937586036,0.0004701425440123626,0.00041739017209329145,0.00036616433349341375,0.00037533276939374964,0.00030517738899954884,0.0002596303254276886,0.00023528035093067998,0.000228059703690082,0.0002620882993850366,0.00019384598734946607,0.00014992196497134774,0.00013281062456515594,0.00011954772868416258,0.00010643213645462733,9.503156136659602e-05,8.453537558638502e-05,7.62638177320671e-05,6.822574602503365e-05,6.0832340691193637e-05,5.420491127964766e-05,4.894958607893324e-05,4.673120799356092e-05,3.8927180352296054e-05,3.525124055495572e-05,3.141319219453202e-05,2.8811426290607224e-05,2.5235375613370988e-05,2.287111631484649e-05,2.039619319129798e-05,1.795782226995683e-05,1.6131672513989147e-05,1.453629606063521e-05
};
if(debug==1)
{cout<<"raw apec_values="; for(i=0;i<n_values;i++)cout<<apec_values[i]<<" "; cout<<endl;}
i=(E-E1)/dE + .0001;
if(debug==1)
cout<<"apec E="<<E<< " i="<<i<< " raw apec_values[i]="<<apec_values[i]<<endl;
double apec_emission_;
// apec units are 1.e-14 photons cm^3 per energy bin
// em is ne*np cm-6 pc
double pc_to_cm = 3.08e18;
double pi = acos(-1);
apec_emission_= 1.e-14 * E * keV_to_erg * apec_values[i] / dE * em * pc_to_cm / (4.*pi) ;
if(debug==1)
cout<<"apec_emission_="<<apec_emission_<<" erg cm-2 sr-1 s-1 keV-1"<<endl;
return apec_emission_;
};
//////////////////////////////////////////////////////////////////////
double apec_emissivity( double density, double T, double Emin,double Emax, int debug)
{
// thermal formula:
//https://wiki.mpe.mpg.de/eRosita/erosnr_project_Fermi_bubbles?action=AttachFile&do=view&target=erosita_garching_2020_FermiBubbles.pdf
// Longair eq 6.47
// thermal formula:
//https://wiki.mpe.mpg.de/eRosita/erosnr_project_Fermi_bubbles?action=AttachFile&do=view&target=erosita_garching_2020_FermiBubbles.pdf
// Longair eq 6.47
double g=1.2; // Gaunt factor
double k=1.38e-16 ; // erg/K
double keV_to_erg = 1.602e-9;
//
// double thermal_emission_= .004 * em * pow(T,-0.5) * g * exp(-(E*keV_to_erg)/(k*T)) ;// cf pdf above. erg cm-2 sr-1 s-1 keV-1
// cout<<"thermal_emission_="<<thermal_emission_<<" erg cm-2 sr-1 s-1 keV-1"<<endl;
// want this for calculation: emissivity_nominal = 1.0e-29 ; // erg cm-3 sr-1 s-1 band 022 from notes
// em in cm-6 pc
// density^2 in cm-6
// density=1e-2 em=1e-4/3e18 = 1e-22 T=1e7 pow(T,.5)=1e-4, * .004 about right
double pc_to_cm= 3.08e18;
double em = pow(density,2) / pc_to_cm; // em for 1 cm so cm-2 -> cm-3
double E=(Emin+Emax)/2.;
double dE=Emax-Emin;
double thermal_emissivity_= .004 * em * pow(T,-0.5) * g * exp(-(E*keV_to_erg)/(k*T)) ;// cf pdf above. erg cm-2 sr-1 s-1 keV-1
if(debug==1) cout<<"thermal_emissivity_="<<thermal_emissivity_<<" erg cm-3 sr-1 s-1 keV-1"<<endl;
thermal_emissivity_ /= (E * keV_to_erg);
if(debug==1) cout<<"thermal_emissivity_="<<thermal_emissivity_<<" counts cm-3 sr-1 s-1 keV-1"<<endl;
thermal_emissivity_ *= dE;
if(debug==1) cout<<"thermal_emissivity_="<<thermal_emissivity_<<" counts cm-3 sr-1 s-1 "<<endl;
double apec_emissivity_;
apec_emissivity_= apec_emission(em, T, E);// erg cm-3 sr-1 s-1 keV-1
if(debug==1)cout<<"apec_emissivity_="<<apec_emissivity_<<" erg cm-3 sr-1 s-1 keV-1 "<<endl;
apec_emissivity_ /= (E * keV_to_erg);
if(debug==1) cout<<"apec_emissivity_="<<apec_emissivity_<<" counts cm-3 sr-1 s-1 keV-1 "<<endl;
apec_emissivity_ *= dE;
if(debug==1) cout<<"apec_emissivity_="<<apec_emissivity_<<" counts cm-3 sr-1 s-1 "<<endl;
return apec_emissivity_;
};