diff options
-rw-r--r-- | mass_profile/beta.hpp | 4 | ||||
-rw-r--r-- | mass_profile/calc_lx.cpp | 16 | ||||
-rw-r--r-- | mass_profile/calc_lx_beta.cpp | 79 | ||||
-rw-r--r-- | mass_profile/calc_lx_dbeta.cpp | 103 | ||||
-rw-r--r-- | mass_profile/chisq.hpp | 26 | ||||
-rw-r--r-- | mass_profile/dbeta.hpp | 10 | ||||
-rw-r--r-- | mass_profile/dump_fit_qdp.cpp | 14 | ||||
-rw-r--r-- | mass_profile/fit_beta_sbp.cpp | 40 | ||||
-rw-r--r-- | mass_profile/fit_dbeta_sbp.cpp | 55 | ||||
-rw-r--r-- | mass_profile/fit_direct_beta.cpp | 10 | ||||
-rw-r--r-- | mass_profile/fit_lt_bpl.cpp | 42 | ||||
-rw-r--r-- | mass_profile/fit_lt_pl.cpp | 22 | ||||
-rw-r--r-- | mass_profile/fit_mt_bpl.cpp | 44 | ||||
-rw-r--r-- | mass_profile/fit_mt_pl.cpp | 18 | ||||
-rw-r--r-- | mass_profile/fit_nfw_mass.cpp | 22 | ||||
-rw-r--r-- | mass_profile/fit_nfw_sbp.cpp | 34 | ||||
-rw-r--r-- | mass_profile/fit_wang2012_model.cpp | 20 | ||||
-rw-r--r-- | mass_profile/nfw_ne.hpp | 20 | ||||
-rw-r--r-- | mass_profile/projector.hpp | 8 | ||||
-rw-r--r-- | mass_profile/vchisq.hpp | 32 |
20 files changed, 300 insertions, 319 deletions
diff --git a/mass_profile/beta.hpp b/mass_profile/beta.hpp index 4e6264c..6ae70ac 100644 --- a/mass_profile/beta.hpp +++ b/mass_profile/beta.hpp @@ -13,7 +13,7 @@ namespace opt_utilities { this->push_param_info(param_info<std::vector<T>,std::string>("n0",1,0,1E99)); this->push_param_info(param_info<std::vector<T>,std::string>("beta",.66,0,1E99)); - this->push_param_info(param_info<std::vector<T>,std::string>("rc",100,0,1E99)); + this->push_param_info(param_info<std::vector<T>,std::string>("rc",100,0,1E99)); } public: @@ -30,7 +30,7 @@ namespace opt_utilities T rc=p[2]; std::vector<T> result(x.size()-1); - for(int i=1;i<x.size();++i) + for(size_t i=1;i<x.size();++i) { T xi=(x[i]+x[i-1])/2; T yi=0; diff --git a/mass_profile/calc_lx.cpp b/mass_profile/calc_lx.cpp index e5cff62..b1137ff 100644 --- a/mass_profile/calc_lx.cpp +++ b/mass_profile/calc_lx.cpp @@ -61,7 +61,7 @@ int main(int argc,char* argv[]) //perform a 1-beta fitting//// ////////////////////////////// fitter<double,double,vector<double>,double,string> f; - + f.set_statistic(chisq<double,double,vector<double>,double,string>()); f.set_opt_method(powell_method<double,vector<double> >()); f.set_model(beta1d<double>()); @@ -74,7 +74,7 @@ int main(int argc,char* argv[]) double rmax=f.get_data_set().get_data(f.get_data_set().size()-1).get_x(); ofstream lx_fit_result("lx_fit_result.qdp"); lx_fit_result<<"read terr 1 2\nskip single\n"; - for(int i=0;i<f.get_data_set().size();++i) + for(size_t i=0;i<f.get_data_set().size();++i) { lx_fit_result<<f.get_data_set().get_data(i).get_x()<<"\t"<< -abs(f.get_data_set().get_data(i).get_x_lower_err())<<"\t"<< @@ -89,13 +89,13 @@ int main(int argc,char* argv[]) { lx_fit_result<<i<<"\t0\t0\t"<<f.eval_model(i,f.get_all_params())<<"\t0\t0"<<endl; } - - for(int i=0;i<f.get_num_params();++i) + + for(size_t i=0;i<f.get_num_params();++i) { cerr<<f.get_param_info(i).get_name()<<"="<< f.get_param_info(i).get_value()<<endl; } - + vector<double> p=f.get_all_params(); double r500kpc=atof(argv[4]); @@ -131,7 +131,7 @@ int main(int argc,char* argv[]) sum_flux+=cnt*ratio; } double lx=sum_flux*4*pi*d_l_cm*d_l_cm; - + double l_mean=0; double l2_mean=0; double cnt=0; @@ -140,7 +140,7 @@ int main(int argc,char* argv[]) cerr<<"."; opt_utilities::default_data_set<double,double> ds(dl.get_data_set()); opt_utilities::default_data_set<double,double> ds1; - for(int i=0;i<ds.size();++i) + for(size_t i=0;i<ds.size();++i) { double yc=ds.get_data(i).get_y(); double ye=(std::abs(ds.get_data(i).get_y_lower_err())+std::abs(ds.get_data(i).get_y_lower_err()))/2; @@ -159,7 +159,7 @@ int main(int argc,char* argv[]) //cout<<f.get_param_value("beta")<<endl; double sum_cnt=0; double sum_flux=0; - + for(double r=0;r<r500pixel;r+=dr) { double r1=r<.2*r500pixel?.2*r500pixel:r; diff --git a/mass_profile/calc_lx_beta.cpp b/mass_profile/calc_lx_beta.cpp index 5bb2570..aae8b2b 100644 --- a/mass_profile/calc_lx_beta.cpp +++ b/mass_profile/calc_lx_beta.cpp @@ -30,20 +30,8 @@ const double Mpc=kpc*1000; double beta_func(double r, double n0,double rc,double beta) { - - return abs(n0)*pow(1+r*r/rc/rc,-3./2.*abs(beta)); -} - - //calculate critical density from z, under following cosmological constants -static double calc_critical_density(double z, - const double H0=2.3E-18, - const double Omega_m=.27) -{ - const double G=6.673E-8;//cm^3 g^-1 s^2 - const double E=std::sqrt(Omega_m*(1+z)*(1+z)*(1+z)+1-Omega_m); - const double H=H0*E; - return 3*H*H/8/pi/G; + return abs(n0)*pow(1+r*r/rc/rc,-3./2.*abs(beta)); } @@ -61,7 +49,7 @@ public: { return spl.get_value(x); } - + //we need this function, when this object is performing a clone of itself spline_func_obj* do_clone()const { @@ -74,7 +62,7 @@ public: { spl.push_point(x,y); } - + //before getting the intepolated value, the spline should be initialzied by calling this function void gen_spline() { @@ -93,9 +81,9 @@ int main(int argc,char* argv[]) ifstream cfg_file(argv[1]); assert(cfg_file.is_open()); cfg_map cfg=parse_cfg_file(cfg_file); - + //check the existence of following parameters - + const double z=cfg.z; //initialize the radius list, sbp list and sbp error list @@ -151,7 +139,7 @@ int main(int argc,char* argv[]) { rmin=cfg.rmin_kpc*kpc/cm_per_pixel; } - + cerr<<"rmin="<<rmin<<endl; std::list<double> radii_tmp,sbps_tmp,sbpe_tmp; radii_tmp.resize(radii.size()); @@ -184,7 +172,7 @@ int main(int argc,char* argv[]) for(ifstream ifs(cfg.cfunc_file.c_str());;) { assert(ifs.is_open()); - double x,y,y1,y2; + double x,y; ifs>>x>>y; if(!ifs.good()) { @@ -199,7 +187,7 @@ int main(int argc,char* argv[]) cf.add_point(x,y);//change with source } cf.gen_spline(); - + //read temperature profile data spline_func_obj Tprof; int tcnt=0; @@ -224,7 +212,7 @@ int main(int argc,char* argv[]) Tprof.gen_spline(); - + default_data_set<std::vector<double>,std::vector<double> > ds; ds.add_data(data<std::vector<double>,std::vector<double> >(radii,sbps,sbpe,sbpe,radii,radii)); @@ -247,11 +235,12 @@ int main(int argc,char* argv[]) //optimization method f.set_opt_method(powell_method<double,std::vector<double> >()); //initialize the initial values + /* double n0=0; - //double beta=atof(arg_map["beta"].c_str()); double beta=0; double rc=0; - double bkg=0; + */ + double bkg_level=0; for(std::map<std::string,std::vector<double> >::iterator i=cfg.param_map.begin(); i!=cfg.param_map.end();++i) @@ -280,12 +269,14 @@ int main(int argc,char* argv[]) f.fit(); f.fit(); std::vector<double> p=f.get_all_params(); + /* n0=f.get_param_value("n0"); rc=f.get_param_value("rc"); beta=f.get_param_value("beta"); + */ //output the datasets and fitting results ofstream param_output("lx_beta_param.txt"); - for(int i=0;i<f.get_num_params();++i) + for(size_t i=0;i<f.get_num_params();++i) { if(f.get_param_info(i).get_name()=="rc") { @@ -322,7 +313,7 @@ int main(int argc,char* argv[]) ofs_sbp<<"color 3 on 3"<<endl; ofs_sbp<<"color 4 on 4"<<endl; ofs_sbp<<"color 5 on 5"<<endl; - + ofs_sbp<<"win 1"<<endl; ofs_sbp<<"yplot 1 2 3 4 5"<<endl; ofs_sbp<<"loc 0 0 1 1"<<endl; @@ -356,7 +347,7 @@ int main(int argc,char* argv[]) ofs_sbp<<"color 4 on 3"<<endl; ofs_sbp<<"color 5 on 4"<<endl; //ofs_sbp<<"ma 1 on 2"<<endl; - + ofs_sbp<<"win 1"<<endl; ofs_sbp<<"yplot 1 2 3 4"<<endl; ofs_sbp<<"loc 0 0 1 1"<<endl; @@ -377,7 +368,7 @@ int main(int argc,char* argv[]) } // cout<<sbps_all.size()<<"\t"<<sbps.size()<<"\t"<<sbps_inner_cut_size<<endl; - for(int i=1;i<sbps_all.size();++i) + for(size_t i=1;i<sbps_all.size();++i) { double x=(radii_all[i]+radii_all[i-1])/2; double y=sbps_all[i-1]; @@ -395,7 +386,7 @@ int main(int argc,char* argv[]) } } ofs_sbp<<"no no no"<<endl; - for(int i=sbps_inner_cut_size;i<sbps_all.size();++i) + for(size_t i=sbps_inner_cut_size;i<sbps_all.size();++i) { double x=(radii_all[i]+radii_all[i-1])/2; double ym=mv[i-1]; @@ -403,8 +394,8 @@ int main(int argc,char* argv[]) } ofs_sbp<<"no no no"<<endl; //bkg level - double bkg_level=abs(f.get_param_value("bkg")); - for(int i=0;i<sbps_all.size();++i) + bkg_level=abs(f.get_param_value("bkg")); + for(size_t i=0;i<sbps_all.size();++i) { double x=(radii_all[i]+radii_all[i-1])/2; ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<bkg_level<<"\t0"<<endl; @@ -420,23 +411,23 @@ int main(int argc,char* argv[]) } //resid ofs_sbp<<"no no no"<<endl; - for(int i=1;i<sbps.size();++i) + for(size_t i=1;i<sbps.size();++i) { double x=(radii[i]+radii[i-1])/2; - double y=sbps[i-1]; - double ye=sbpe[i-1]; + //double y=sbps[i-1]; + //double ye=sbpe[i-1]; double ym=mv[i-1]; ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<(ym-sbps[i-1])/sbpe[i-1]<<"\t"<<1<<endl; } //zero level of resid ofs_sbp<<"no no no"<<endl; - for(int i=1;i<sbps.size();++i) + for(size_t i=1;i<sbps.size();++i) { double x=(radii[i]+radii[i-1])/2; - double y=sbps[i-1]; - double ye=sbpe[i-1]; - double ym=mv[i-1]; + //double y=sbps[i-1]; + //double ye=sbpe[i-1]; + //double ym=mv[i-1]; ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<0<<"\t"<<0<<endl; } @@ -457,13 +448,13 @@ int main(int argc,char* argv[]) p.back()=0; radii.clear(); double rout=atof(argv[2])*kpc; - + for(double r=0;r<rout;r+=1*kpc)//step size=1kpc { double r_pix=r/cm_per_pixel; radii.push_back(r_pix); } - + double Da=cm_per_pixel/(.492/3600./180.*pi); double Dl=Da*(1+z)*(1+z); cout<<"dl="<<Dl/kpc<<endl; @@ -481,19 +472,19 @@ int main(int argc,char* argv[]) break; } //cerr<<x<<"\t"<<y<<endl; - + cf_bolo_erg.add_point(x,y);//change with source } cf_bolo_erg.gen_spline(); projector<double>& pj=dynamic_cast<projector<double>&>(f.get_model()); pj.attach_cfunc(cf_bolo_erg); - - - + + + mv=f.eval_model_raw(radii,p); double flux_erg=0; - for(int i=0;i<radii.size()-1;++i) + for(size_t i=0;i<radii.size()-1;++i) { double S=pi*(radii[i+1]+radii[i])*(radii[i+1]-radii[i]); flux_erg+=S*mv[i]; diff --git a/mass_profile/calc_lx_dbeta.cpp b/mass_profile/calc_lx_dbeta.cpp index f105869..f8f2c2a 100644 --- a/mass_profile/calc_lx_dbeta.cpp +++ b/mass_profile/calc_lx_dbeta.cpp @@ -26,20 +26,8 @@ double dbeta_func(double r, double n01,double rc1,double beta1, double n02,double rc2,double beta2) { - - return abs(n01)*pow(1+r*r/rc1/rc1,-3./2.*abs(beta1))+abs(n02)*pow(1+r*r/rc2/rc2,-3./2.*abs(beta2)); -} - - //calculate critical density from z, under following cosmological constants -static double calc_critical_density(double z, - const double H0=2.3E-18, - const double Omega_m=.27) -{ - const double G=6.673E-8;//cm^3 g^-1 s^2 - const double E=std::sqrt(Omega_m*(1+z)*(1+z)*(1+z)+1-Omega_m); - const double H=H0*E; - return 3*H*H/8/pi/G; + return abs(n01)*pow(1+r*r/rc1/rc1,-3./2.*abs(beta1))+abs(n02)*pow(1+r*r/rc2/rc2,-3./2.*abs(beta2)); } @@ -57,7 +45,7 @@ public: { return spl.get_value(x); } - + //we need this function, when this object is performing a clone of itself spline_func_obj* do_clone()const { @@ -70,7 +58,7 @@ public: { spl.push_point(x,y); } - + //before getting the intepolated value, the spline should be initialzied by calling this function void gen_spline() { @@ -89,9 +77,9 @@ int main(int argc,char* argv[]) ifstream cfg_file(argv[1]); assert(cfg_file.is_open()); cfg_map cfg=parse_cfg_file(cfg_file); - + //check the existence of following parameters - + const double z=cfg.z; //initialize the radius list, sbp list and sbp error list @@ -147,7 +135,7 @@ int main(int argc,char* argv[]) { rmin=cfg.rmin_kpc*kpc/cm_per_pixel; } - + cerr<<"rmin="<<rmin<<endl; std::list<double> radii_tmp,sbps_tmp,sbpe_tmp; radii_tmp.resize(radii.size()); @@ -180,7 +168,7 @@ int main(int argc,char* argv[]) for(ifstream ifs(cfg.cfunc_file.c_str());;) { assert(ifs.is_open()); - double x,y,y1,y2; + double x,y; ifs>>x>>y; if(!ifs.good()) { @@ -195,7 +183,7 @@ int main(int argc,char* argv[]) cf.add_point(x,y);//change with source } cf.gen_spline(); - + //read temperature profile data spline_func_obj Tprof; int tcnt=0; @@ -220,10 +208,10 @@ int main(int argc,char* argv[]) Tprof.gen_spline(); - + default_data_set<std::vector<double>,std::vector<double> > ds; ds.add_data(data<std::vector<double>,std::vector<double> >(radii,sbps,sbpe,sbpe,radii,radii)); - + //initial fitter fitter<vector<double>,vector<double>,vector<double>,double> f; f.load_data(ds); @@ -255,7 +243,7 @@ int main(int argc,char* argv[]) //attach the cooling function a.attach_cfunc(cf); a.set_cm_per_pixel(cm_per_pixel); - + f.set_model(a); //chi^2 statistic vchisq<double> c; @@ -265,12 +253,15 @@ int main(int argc,char* argv[]) //optimization method f.set_opt_method(powell_method<double,std::vector<double> >()); //initialize the initial values - double n01=0; + /* + double =0; double rc1=0; double n02=0; double rc2=0; double beta=0; - double bkg=0; + */ + double bkg_level=0; + if(tie_beta) { f.set_param_value("beta",.7); @@ -311,7 +302,7 @@ int main(int argc,char* argv[]) } - + //perform the fitting, first freeze beta1, beta2, rc1, and rc2 if(tie_beta) { @@ -328,9 +319,9 @@ int main(int argc,char* argv[]) freeze_param<std::vector<double>,std::vector<double>,std::vector<double>,std::string>("rc2") ); } - + f.fit(); - + f.clear_param_modifier(); //then perform the fitting, freeze beta1 and beta2 @@ -341,9 +332,11 @@ int main(int argc,char* argv[]) //finally thaw all parameters f.fit(); + + /* double beta1=0; double beta2=0; - + n01=f.get_param_value("n01"); rc1=f.get_param_value("rc1"); n02=f.get_param_value("n02"); @@ -359,10 +352,12 @@ int main(int argc,char* argv[]) beta1=f.get_param_value("beta1"); beta2=f.get_param_value("beta2"); } + */ + //output the params ofstream param_output("lx_dbeta_param.txt"); //output the datasets and fitting results - for(int i=0;i<f.get_num_params();++i) + for(size_t i=0;i<f.get_num_params();++i) { if(f.get_param_info(i).get_name()=="rc1") { @@ -379,19 +374,19 @@ int main(int argc,char* argv[]) } cerr<<"chi square="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl; param_output<<"chi square="<<f.get_statistic_value()/(radii.size()-f.get_model().get_num_free_params())<<endl; - + //c.verbose(false); //f.set_statistic(c); //f.fit(); std::vector<double> p=f.get_all_params(); f.clear_param_modifier(); std::vector<double> mv=f.eval_model(radii,p); - - + + ofstream ofs_sbp("lx_sbp_fit.qdp"); ofs_sbp<<"read serr 2"<<endl; ofs_sbp<<"skip single"<<endl; - + ofs_sbp<<"line on 2"<<endl; ofs_sbp<<"line on 3"<<endl; ofs_sbp<<"line on 4"<<endl; @@ -401,7 +396,7 @@ int main(int argc,char* argv[]) ofs_sbp<<"ls 2 on 3"<<endl; ofs_sbp<<"ls 2 on 4"<<endl; ofs_sbp<<"ls 2 on 5"<<endl; - + ofs_sbp<<"!LAB POS Y 4.00"<<endl; @@ -423,27 +418,27 @@ int main(int argc,char* argv[]) ofs_sbp<<"log x"<<endl; ofs_sbp<<"log y off"<<endl; ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[sbps.size()-2]+radii[sbps.size()-1])/2*cm_per_pixel/kpc<<endl; - for(int i=1;i<sbps.size();++i) + for(size_t i=1;i<sbps.size();++i) { double x=(radii[i]+radii[i-1])/2; double y=sbps[i-1]; double ye=sbpe[i-1]; - double ym=mv[i-1]; + //double ym=mv[i-1]; ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<y<<"\t"<<ye<<endl; } ofs_sbp<<"no no no"<<endl; - for(int i=1;i<sbps.size();++i) + for(size_t i=1;i<sbps.size();++i) { double x=(radii[i]+radii[i-1])/2; - double y=sbps[i-1]; - double ye=sbpe[i-1]; + //double y=sbps[i-1]; + //double ye=sbpe[i-1]; double ym=mv[i-1]; ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<ym<<"\t"<<0<<endl; } //bkg ofs_sbp<<"no no no"<<endl; - double bkg_level=abs(f.get_param_value("bkg")); - for(int i=0;i<sbps.size();++i) + bkg_level=abs(f.get_param_value("bkg")); + for(size_t i=0;i<sbps.size();++i) { double x=(radii[i]+radii[i-1])/2; ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<bkg_level<<"\t0"<<endl; @@ -466,22 +461,22 @@ int main(int argc,char* argv[]) } //resid ofs_sbp<<"no no no"<<endl; - for(int i=1;i<sbps.size();++i) + for(size_t i=1;i<sbps.size();++i) { double x=(radii[i]+radii[i-1])/2; - double y=sbps[i-1]; - double ye=sbpe[i-1]; + //double y=sbps[i-1]; + //double ye=sbpe[i-1]; double ym=mv[i-1]; ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<(ym-sbps[i-1])/sbpe[i-1]<<"\t"<<1<<endl; } //zero level in resid map ofs_sbp<<"no no no"<<endl; - for(int i=1;i<sbps.size();++i) + for(size_t i=1;i<sbps.size();++i) { double x=(radii[i]+radii[i-1])/2; - double y=sbps[i-1]; - double ye=sbpe[i-1]; - double ym=mv[i-1]; + //double y=sbps[i-1]; + //double ye=sbpe[i-1]; + //double ym=mv[i-1]; ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<0<<"\t"<<0<<endl; } @@ -502,7 +497,7 @@ int main(int argc,char* argv[]) */ p.back()=0; - + radii.clear(); double rout=atof(argv[2])*kpc; for(double r=0;r<rout;r+=1*kpc)//step size=1kpc @@ -527,17 +522,17 @@ int main(int argc,char* argv[]) break; } //cerr<<x<<"\t"<<y<<endl; - + cf_bolo_erg.add_point(x,y);//change with source } cf_bolo_erg.gen_spline(); - + projector<double>& pj=dynamic_cast<projector<double>&>(f.get_model()); pj.attach_cfunc(cf_bolo_erg); - + mv=f.eval_model_raw(radii,p); double flux_erg=0; - for(int i=0;i<radii.size()-1;++i) + for(size_t i=0;i<radii.size()-1;++i) { double S=pi*(radii[i+1]+radii[i])*(radii[i+1]-radii[i]); flux_erg+=S*mv[i]; diff --git a/mass_profile/chisq.hpp b/mass_profile/chisq.hpp index 575f1c0..bdfadcb 100644 --- a/mass_profile/chisq.hpp +++ b/mass_profile/chisq.hpp @@ -46,7 +46,7 @@ namespace opt_utilities bool verb;
bool limit_bound;
int n;
-
+
statistic<Ty,Tx,Tp,Ts,Tstr>* do_clone()const
{
// return const_cast<statistic<Ty,Tx,Tp>*>(this);
@@ -76,15 +76,15 @@ namespace opt_utilities chisq()
:verb(true),limit_bound(false)
{}
-
-
+
+
Ty do_eval(const Tp& p)
{
if(limit_bound)
{
Tp p1=this->get_fitter().get_model().reform_param(p);
- for(int i=0;i<p1.size();++i)
+ for(size_t i=0;i<p1.size();++i)
{
if(p1[i]>this->get_fitter().get_param_info(i).get_upper_limit()||
p1[i]<this->get_fitter().get_param_info(i).get_lower_limit())
@@ -111,7 +111,7 @@ namespace opt_utilities vye2.resize(this->get_data_set().size());
my.resize(this->get_data_set().size());
}
-
+
}
for(int i=(this->get_data_set()).size()-1;i>=0;--i)
@@ -132,7 +132,7 @@ namespace opt_utilities Ty y_model=this->eval_model(this->get_data_set().get_data(i).get_x(),p);
Ty y_obs=this->get_data_set().get_data(i).get_y();
Ty y_err;
-
+
Ty errx=0;
if(errx1<errx2)
{
@@ -170,7 +170,7 @@ namespace opt_utilities Ty chi=(y_obs-y_model)/std::sqrt(y_err*y_err+errx*errx);
result+=chi*chi;
-
+
if(verb&&n%display_interval==0)
{
vx.at(i)=this->get_data_set().get_data(i).get_x();
@@ -178,7 +178,7 @@ namespace opt_utilities vye1.at(i)=std::abs(this->get_data_set().get_data(i).get_y_lower_err());
vye2.at(i)=std::abs(this->get_data_set().get_data(i).get_y_upper_err());
my.at(i)=y_model;
-
+
xmin=std::min(vx.at(i),xmin);
ymin=std::min(vy.at(i),ymin-vye1[i]);
xmax=std::max(vx.at(i),xmax);
@@ -192,7 +192,7 @@ namespace opt_utilities if(n%display_interval==0)
{
cerr<<result<<"\t";
- for(int i=0;i<(int)get_size(p);++i)
+ for(size_t i=0;i<get_size(p);++i)
{
cerr<<get_element(p,i)<<",";
}
@@ -205,15 +205,13 @@ namespace opt_utilities }
}
-
+
return result;
}
};
-
-
+
+
}
#endif
//EOF
-
-
diff --git a/mass_profile/dbeta.hpp b/mass_profile/dbeta.hpp index 8c3a7c5..7246b44 100644 --- a/mass_profile/dbeta.hpp +++ b/mass_profile/dbeta.hpp @@ -24,7 +24,7 @@ namespace opt_utilities this->push_param_info(param_info<std::vector<T>,std::string>("n02",1)); this->push_param_info(param_info<std::vector<T>,std::string>("beta2",.67)); this->push_param_info(param_info<std::vector<T>,std::string>("rc2",110)); - + } public: @@ -44,10 +44,10 @@ namespace opt_utilities T beta2=p[4]; T rc2=p[5]; - + std::vector<T> result(x.size()-1); - for(int i=1;i<x.size();++i) + for(size_t i=1;i<x.size();++i) { T xi=(x[i]+x[i-1])/2; T yi=0; @@ -70,7 +70,7 @@ namespace opt_utilities this->push_param_info(param_info<std::vector<T>,std::string>("n02",1)); this->push_param_info(param_info<std::vector<T>,std::string>("rc2",110)); this->push_param_info(param_info<std::vector<T>,std::string>("beta",.67)); - + } public: @@ -92,7 +92,7 @@ namespace opt_utilities T beta2=beta; std::vector<T> result(x.size()-1); - for(int i=1;i<x.size();++i) + for(size_t i=1;i<x.size();++i) { T xi=(x[i]+x[i-1])/2; T yi=0; diff --git a/mass_profile/dump_fit_qdp.cpp b/mass_profile/dump_fit_qdp.cpp index 556edea..85307d1 100644 --- a/mass_profile/dump_fit_qdp.cpp +++ b/mass_profile/dump_fit_qdp.cpp @@ -10,13 +10,13 @@ namespace opt_utilities os<<"la x \"radius (kpc)\""<<std::endl; os<<"la y \"surface brightness (cts s\\u-1\\d pixel\\u-2\\d)\""<<std::endl; os<<"li on 2"<<std::endl; - for(int i=1;i<r.size();++i) + for(size_t i=1;i<r.size();++i) { os<<(r[i]+r[i-1])/2*cm_per_pixel/kpc<<"\t"<<(r[i]-r[i-1])/2*cm_per_pixel/kpc<<"\t"<<y[i-1]<<"\t"<<ye[i-1]<<std::endl; } os<<"no no no"<<std::endl; std::vector<double> p=f.get_all_params(); - for(int i=1;i<r.size();++i) + for(size_t i=1;i<r.size();++i) { double x=(r[i]+r[i-1])/2; os<<x*cm_per_pixel/kpc<<"\t"<<0<<"\t"<<f.eval_model_raw(x,p)<<"\t"<<0<<std::endl; @@ -30,26 +30,26 @@ namespace opt_utilities os<<"la x \"radius (kpc)\""<<std::endl; os<<"la y \"density (cm\\u-3\\d)\""<<std::endl; os<<"li on 2"<<std::endl; - - for(int i=1;i<r.size();++i) + + for(size_t i=1;i<r.size();++i) { double x=(r[i]+r[i-1])/2; double y=sbps[i-1]; double ye=sbpe[i-1]; os<<x*cm_per_pixel/kpc<<"\t0\t"<<y<<"\t"<<ye<<std::endl; } - + os<<"no no no"<<std::endl; std::vector<double> p=f.get_all_params(); std::vector<double> mv=f.eval_model_raw(r,p); - for(int i=1;i<r.size();++i) + for(size_t i=1;i<r.size();++i) { double x=(r[i]+r[i-1])/2; double y=mv[i-1]; double ye=0; os<<x*cm_per_pixel/kpc<<"\t0\t"<<y<<"\t"<<ye<<std::endl; } - + } }; diff --git a/mass_profile/fit_beta_sbp.cpp b/mass_profile/fit_beta_sbp.cpp index 41d6192..9db48cc 100644 --- a/mass_profile/fit_beta_sbp.cpp +++ b/mass_profile/fit_beta_sbp.cpp @@ -184,7 +184,7 @@ int main(int argc,char* argv[]) for(ifstream ifs(cfg.cfunc_file.c_str());;) { assert(ifs.is_open()); - double x,y,y1,y2; + double x,y; ifs>>x>>y; if(!ifs.good()) { @@ -251,7 +251,7 @@ int main(int argc,char* argv[]) //double beta=atof(arg_map["beta"].c_str()); double beta=0; double rc=0; - double bkg=0; + double bkg_level=0; for(std::map<std::string,std::vector<double> >::iterator i=cfg.param_map.begin(); i!=cfg.param_map.end();++i) @@ -285,7 +285,7 @@ int main(int argc,char* argv[]) beta=f.get_param_value("beta"); //output the datasets and fitting results ofstream param_output("beta_param.txt"); - for(int i=0;i<f.get_num_params();++i) + for(size_t i=0;i<f.get_num_params();++i) { if(f.get_param_info(i).get_name()=="rc") { @@ -377,7 +377,7 @@ int main(int argc,char* argv[]) } // cout<<sbps_all.size()<<"\t"<<sbps.size()<<"\t"<<sbps_inner_cut_size<<endl; - for(int i=1;i<sbps_all.size();++i) + for(size_t i=1;i<sbps_all.size();++i) { double x=(radii_all[i]+radii_all[i-1])/2; double y=sbps_all[i-1]; @@ -395,7 +395,7 @@ int main(int argc,char* argv[]) } } ofs_sbp<<"no no no"<<endl; - for(int i=sbps_inner_cut_size;i<sbps_all.size();++i) + for(size_t i=sbps_inner_cut_size;i<sbps_all.size();++i) { double x=(radii_all[i]+radii_all[i-1])/2; double ym=mv[i-1]; @@ -403,8 +403,8 @@ int main(int argc,char* argv[]) } ofs_sbp<<"no no no"<<endl; //bkg level - double bkg_level=abs(f.get_param_value("bkg")); - for(int i=0;i<sbps_all.size();++i) + bkg_level=abs(f.get_param_value("bkg")); + for(size_t i=0;i<sbps_all.size();++i) { double x=(radii_all[i]+radii_all[i-1])/2; ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<bkg_level<<"\t0"<<endl; @@ -420,23 +420,23 @@ int main(int argc,char* argv[]) } //resid ofs_sbp<<"no no no"<<endl; - for(int i=1;i<sbps.size();++i) + for(size_t i=1;i<sbps.size();++i) { double x=(radii[i]+radii[i-1])/2; - double y=sbps[i-1]; - double ye=sbpe[i-1]; + //double y=sbps[i-1]; + //double ye=sbpe[i-1]; double ym=mv[i-1]; ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<(ym-sbps[i-1])/sbpe[i-1]<<"\t"<<1<<endl; } //zero level of resid ofs_sbp<<"no no no"<<endl; - for(int i=1;i<sbps.size();++i) + for(size_t i=1;i<sbps.size();++i) { double x=(radii[i]+radii[i-1])/2; - double y=sbps[i-1]; - double ye=sbpe[i-1]; - double ym=mv[i-1]; + //double y=sbps[i-1]; + //double ye=sbpe[i-1]; + //double ym=mv[i-1]; ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<0<<"\t"<<0<<endl; } @@ -455,16 +455,16 @@ int main(int argc,char* argv[]) } */ - double lower,upper; + //double lower,upper; double dr=1; //calculate the mass profile - const double G=6.673E-8;//cm^3 g^-1 s^-2 + //const double G=6.673E-8;//cm^3 g^-1 s^-2 // Molecular weight per electron // Reference: Ettori et al. 2013, Space Sci. Rev., 177, 119-154; Eq.(9) below static const double mu=1.155; static const double mp=1.67262158E-24;//g static const double M_sun=1.98892E33;//g - static const double k=1.38E-16; + //static const double k=1.38E-16; ofstream ofs_mass("mass_int.qdp"); ofstream ofs_mass_dat("mass_int.dat"); @@ -494,14 +494,14 @@ int main(int argc,char* argv[]) double T_keV=Tprof(r); double T1_keV=Tprof(r1); - double T_K=T_keV*11604505.9; - double T1_K=T1_keV*11604505.9; + //double T_K=T_keV*11604505.9; + //double T1_K=T1_keV*11604505.9; double dlnT=log(T1_keV/T_keV); double dlnr=log(r+dr)-log(r); double dlnn=log(ne1/ne); - double r_kpc=r_cm/kpc; + //double r_kpc=r_cm/kpc; double r_Mpc=r_cm/Mpc; //double M=-r_cm*T_K*k/G/mu/mp*(dlnT/dlnr+dlnn/dlnr); //ref:http://adsabs.harvard.edu/abs/2012MNRAS.422.3503W diff --git a/mass_profile/fit_dbeta_sbp.cpp b/mass_profile/fit_dbeta_sbp.cpp index 833059c..7fb1c7d 100644 --- a/mass_profile/fit_dbeta_sbp.cpp +++ b/mass_profile/fit_dbeta_sbp.cpp @@ -180,7 +180,7 @@ int main(int argc,char* argv[]) for(ifstream ifs(cfg.cfunc_file.c_str());;) { assert(ifs.is_open()); - double x,y,y1,y2; + double x,y; ifs>>x>>y; if(!ifs.good()) { @@ -270,7 +270,7 @@ int main(int argc,char* argv[]) double n02=0; double rc2=0; double beta=0; - double bkg=0; + double bkg_level=0; if(tie_beta) { f.set_param_value("beta",.7); @@ -362,7 +362,7 @@ int main(int argc,char* argv[]) //output the params ofstream param_output("dbeta_param.txt"); //output the datasets and fitting results - for(int i=0;i<f.get_num_params();++i) + for(size_t i=0;i<f.get_num_params();++i) { if(f.get_param_info(i).get_name()=="rc1") { @@ -423,27 +423,27 @@ int main(int argc,char* argv[]) ofs_sbp<<"log x"<<endl; ofs_sbp<<"log y off"<<endl; ofs_sbp<<"r x "<<(radii[1]+radii[0])/2*cm_per_pixel/kpc<<" "<<(radii[sbps.size()-2]+radii[sbps.size()-1])/2*cm_per_pixel/kpc<<endl; - for(int i=1;i<sbps.size();++i) + for(size_t i=1;i<sbps.size();++i) { double x=(radii[i]+radii[i-1])/2; double y=sbps[i-1]; double ye=sbpe[i-1]; - double ym=mv[i-1]; + //double ym=mv[i-1]; ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<y<<"\t"<<ye<<endl; } ofs_sbp<<"no no no"<<endl; - for(int i=1;i<sbps.size();++i) + for(size_t i=1;i<sbps.size();++i) { double x=(radii[i]+radii[i-1])/2; - double y=sbps[i-1]; - double ye=sbpe[i-1]; + //double y=sbps[i-1]; + //double ye=sbpe[i-1]; double ym=mv[i-1]; ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<ym<<"\t"<<0<<endl; } //bkg ofs_sbp<<"no no no"<<endl; - double bkg_level=abs(f.get_param_value("bkg")); - for(int i=0;i<sbps.size();++i) + bkg_level=abs(f.get_param_value("bkg")); + for(size_t i=0;i<sbps.size();++i) { double x=(radii[i]+radii[i-1])/2; ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<bkg_level<<"\t0"<<endl; @@ -466,22 +466,22 @@ int main(int argc,char* argv[]) } //resid ofs_sbp<<"no no no"<<endl; - for(int i=1;i<sbps.size();++i) + for(size_t i=1;i<sbps.size();++i) { double x=(radii[i]+radii[i-1])/2; - double y=sbps[i-1]; - double ye=sbpe[i-1]; + //double y=sbps[i-1]; + //double ye=sbpe[i-1]; double ym=mv[i-1]; ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<(ym-sbps[i-1])/sbpe[i-1]<<"\t"<<1<<endl; } //zero level in resid map ofs_sbp<<"no no no"<<endl; - for(int i=1;i<sbps.size();++i) + for(size_t i=1;i<sbps.size();++i) { double x=(radii[i]+radii[i-1])/2; - double y=sbps[i-1]; - double ye=sbpe[i-1]; - double ym=mv[i-1]; + //double y=sbps[i-1]; + //double ye=sbpe[i-1]; + //double ym=mv[i-1]; ofs_sbp<<x*cm_per_pixel/kpc<<"\t"<<0<<"\t"<<0<<endl; } @@ -501,16 +501,16 @@ int main(int argc,char* argv[]) } */ - double lower,upper; + //double lower,upper; double dr=1; //calculate the mass profile - const double G=6.673E-8;//cm^3 g^-1 s^-2 + //const double G=6.673E-8;//cm^3 g^-1 s^-2 // Molecular weight per electron // Reference: Ettori et al. 2013, Space Sci. Rev., 177, 119-154; Eq.(9) below static const double mu=1.155; static const double mp=1.67262158E-24;//g static const double M_sun=1.98892E33;//g - static const double k=1.38E-16; + //static const double k=1.38E-16; ofstream ofs_mass("mass_int.qdp"); ofstream ofs_mass_dat("mass_int.dat"); @@ -536,26 +536,23 @@ int main(int argc,char* argv[]) gas_mass+=dmgas; ofs_gas_mass<<r*cm_per_pixel/kpc<<"\t"<<gas_mass<<endl; - double ne_beta1=dbeta_func(r,n01,rc1,beta1, - 0,rc2,beta2); + double ne_beta1=dbeta_func(r,n01,rc1,beta1, 0,rc2,beta2); - double ne_beta2=dbeta_func(r,0,rc1,beta1, - n02,rc2,beta2); + double ne_beta2=dbeta_func(r,0,rc1,beta1, n02,rc2,beta2); - double ne1=dbeta_func(r1,n01,rc1,beta1, - n02,rc2,beta2);//cm^3 + double ne1=dbeta_func(r1,n01,rc1,beta1, n02,rc2,beta2);//cm^3 double T_keV=Tprof(r); double T1_keV=Tprof(r1); - double T_K=T_keV*11604505.9; - double T1_K=T1_keV*11604505.9; + //double T_K=T_keV*11604505.9; + //double T1_K=T1_keV*11604505.9; double dlnT=log(T1_keV/T_keV); double dlnr=log(r+dr)-log(r); double dlnn=log(ne1/ne); - double r_kpc=r_cm/kpc; + //double r_kpc=r_cm/kpc; double r_Mpc=r_cm/Mpc; //double M=-r_cm*T_K*k/G/mu/mp*(dlnT/dlnr+dlnn/dlnr); //ref:http://adsabs.harvard.edu/abs/2012MNRAS.422.3503W diff --git a/mass_profile/fit_direct_beta.cpp b/mass_profile/fit_direct_beta.cpp index cc1d7c9..ce9c63c 100644 --- a/mass_profile/fit_direct_beta.cpp +++ b/mass_profile/fit_direct_beta.cpp @@ -22,7 +22,7 @@ int main(int argc,char* argv[]) } fitter<double,double,vector<double>,double,string> f; - + f.set_statistic(chisq<double,double,vector<double>,double,string>()); f.set_opt_method(powell_method<double,vector<double> >()); f.set_model(beta1d<double>()); @@ -31,11 +31,11 @@ int main(int argc,char* argv[]) ifs>>dl; f.load_data(dl.get_data_set()); f.fit(); - + double rmin=f.get_data_set().get_data(0).get_x(); double rmax=f.get_data_set().get_data(f.get_data_set().size()-1).get_x(); cout<<"read terr 1 2\nskip single\n"; - for(int i=0;i<f.get_data_set().size();++i) + for(size_t i=0;i<f.get_data_set().size();++i) { cout<<f.get_data_set().get_data(i).get_x()<<"\t"<< -abs(f.get_data_set().get_data(i).get_x_lower_err())<<"\t"<< @@ -44,7 +44,7 @@ int main(int argc,char* argv[]) -abs(f.get_data_set().get_data(i).get_y_lower_err())<<"\t"<< abs(f.get_data_set().get_data(i).get_y_upper_err())<<endl; - + } cout<<"no no no\n"; @@ -53,7 +53,7 @@ int main(int argc,char* argv[]) cout<<i<<"\t0\t0\t"<<f.eval_model(i,f.get_all_params())<<"\t0\t0"<<endl; } - for(int i=0;i<f.get_num_params();++i) + for(size_t i=0;i<f.get_num_params();++i) { cerr<<f.get_param_info(i).get_name()<<"="<< f.get_param_info(i).get_value()<<endl; diff --git a/mass_profile/fit_lt_bpl.cpp b/mass_profile/fit_lt_bpl.cpp index 5637dca..0c06b22 100644 --- a/mass_profile/fit_lt_bpl.cpp +++ b/mass_profile/fit_lt_bpl.cpp @@ -45,7 +45,7 @@ double shuffle_data(double xc,double xl,double xu) double lxc=log(xc); double lxl=log(xc-xl)-log(xc); double lxu=log(xc+xu)-log(xc); - + if(std_norm_rand()>0) { double result=std::exp(lxc-std::abs(std_norm_rand()*lxl)); @@ -108,7 +108,7 @@ int main(int argc,char* argv[]) } line+=" "; istringstream iss(line); - + if(line[0]=='#') { if(!is_first_nonono) @@ -190,16 +190,16 @@ int main(int argc,char* argv[]) double gamma0l=k0l; double gamma0u=k0u; - + double ampl0l=exp(b0l)*pow(Tb,gamma0l); double ampl0u=exp(b0u)*pow(Tb,gamma0u);; - + ofs_result<<"no no no"<<endl; fitter<double,double,vector<double>,double,std::string> fit; fit.set_opt_method(powell_method<double,vector<double> >()); - + fit.set_statistic(logchisq<double,double,vector<double>,double,std::string>()); //fit.set_statistic(leastsq<double,double,vector<double>,double,std::string>()); fit.set_model(bpl1d<double>()); @@ -208,7 +208,7 @@ int main(int argc,char* argv[]) cerr<<"k0l="<<k0l<<endl; cerr<<"k0u="<<k0u<<endl; cerr<<"Ampl0="<<(ampl0l+ampl0u)/2<<endl; - + fit.set_param_value("bpx",Tb); fit.set_param_value("bpy",(ampl0l+ampl0u)/2); fit.set_param_value("gamma1",gamma0l); @@ -225,7 +225,7 @@ int main(int argc,char* argv[]) ofs_result<<i<<"\t0\t0\t"<<fit.eval_model_raw(i,p)/yunit<<"\t0\t0\n"; } - + std::vector<double> mean_p(p.size()); std::vector<double> mean_p2(p.size()); int cnt=0; @@ -238,7 +238,7 @@ int main(int argc,char* argv[]) double sxl=0; double syl=0; double sxyl=0; - + double sxxu=0; double s1u=0; double sxu=0; @@ -246,7 +246,7 @@ int main(int argc,char* argv[]) double sxyu=0; opt_utilities::default_data_set<double,double> ds1; - for(int i=0;i<ds.size();++i) + for(size_t i=0;i<ds.size();++i) { double x=ds.get_data(i).get_x(); double y=ds.get_data(i).get_y(); @@ -254,14 +254,14 @@ int main(int argc,char* argv[]) double xu=ds.get_data(i).get_x_upper_err(); double yl=ds.get_data(i).get_y_lower_err(); double yu=ds.get_data(i).get_y_upper_err(); - + double new_x=shuffle_data(x, xl, xu); double new_y=shuffle_data(y, yl, yu); - + ds1.add_data(data<double,double>(new_x,new_y, yl/y*new_y, yu/y*new_y, @@ -292,40 +292,40 @@ int main(int argc,char* argv[]) double Mbl=sxxl*syl-sxl*sxyl; double k0l=Mal/Ml; double b0l=Mbl/Ml; - + double Mu=sxxu*s1u-sxu*sxu; double Mau=sxyu*s1u-syu*sxu; double Mbu=sxxu*syu-sxu*sxyu; double k0u=Mau/Mu; double b0u=Mbu/Mu; - + double gamma0l=k0l; double gamma0u=k0u; - + double ampl0l=exp(b0l)*pow(Tb,gamma0l); double ampl0u=exp(b0u)*pow(Tb,gamma0u);; - + fit.set_param_value("bpx",Tb); fit.set_param_value("bpy",(ampl0l+ampl0u)/2); fit.set_param_value("gamma1",gamma0l); fit.set_param_value("gamma2",gamma0u); - - + + fit.load_data(ds1); - + fit.fit(); vector<double> p=fit.fit(); - for(int i=0;i<p.size();++i) + for(size_t i=0;i<p.size();++i) { mean_p[i]+=p[i]; mean_p2[i]+=p[i]*p[i]; } //cerr<<fit.get_param_value("gamma1")<<"\t"<<fit.get_param_value("gamma2")<<endl; - + } vector<double> std_p(p.size()); cerr<<endl; - for(int i=0;i<mean_p.size();++i) + for(size_t i=0;i<mean_p.size();++i) { mean_p[i]/=cnt; mean_p2[i]/=cnt; diff --git a/mass_profile/fit_lt_pl.cpp b/mass_profile/fit_lt_pl.cpp index e1686e5..fa37249 100644 --- a/mass_profile/fit_lt_pl.cpp +++ b/mass_profile/fit_lt_pl.cpp @@ -28,13 +28,13 @@ double std_norm_rand() double x=0; double u=0; double v=0; - + do { u=rand()/(double)RAND_MAX; rand(); v=rand()/(double)RAND_MAX; - + x=std::sqrt(-log(u))*cos(2*pi*v); }while(isnan(x)); return x; @@ -65,7 +65,7 @@ int main(int argc,char* argv[]) cerr<<"Usage:"<<argv[0]<<" <a 5 column file with T -Terr +Terr L Lerr> <T lower limit>"<<endl; return -1; } - double T_lower_limit(atof(argv[2])); + double T_lower_limit(atof(argv[2])); ifstream ifs_data(argv[1]); default_data_set<double,double> ds; ofstream ofs_result("l-t_result.qdp"); @@ -91,15 +91,15 @@ int main(int argc,char* argv[]) double L,Lerr; std::string line; getline(ifs_data,line); - - + + if(!ifs_data.good()) { break; } line+=" "; istringstream iss(line); - + if(line[0]=='#') { if(!is_first_nonono) @@ -163,7 +163,7 @@ int main(int argc,char* argv[]) double Mb=sxx*sy-sx*sxy; double k0=Ma/M; double b0=Mb/M; - + ofs_result<<"no no no"<<endl; fitter<double,double,vector<double>,double,std::string> fit; fit.set_opt_method(powell_method<double,vector<double> >()); @@ -185,7 +185,7 @@ int main(int argc,char* argv[]) ofs_result<<i<<"\t0\t0\t"<<exp(fit.eval_model_raw(log(i),p))/yunit<<"\t0\n"; } - + double mean_A=0; double mean_A2=0; double mean_g=0; @@ -196,7 +196,7 @@ int main(int argc,char* argv[]) ++cnt; cerr<<"."; opt_utilities::default_data_set<double,double> ds1; - for(int i=0;i<ds.size();++i) + for(size_t i=0;i<ds.size();++i) { double new_x=shuffle_data(ds.get_data(i).get_x(), ds.get_data(i).get_x_lower_err(), @@ -212,7 +212,7 @@ int main(int argc,char* argv[]) //cerr<<new_x<<"\t"<<new_y<<endl; } fit.load_data(ds1); - + fit.fit(); double k=fit.get_param_value("k"); double b=fit.get_param_value("b"); @@ -234,6 +234,6 @@ int main(int argc,char* argv[]) std::cerr<<"L=L0*T^gamma"<<endl; std::cout<<"L0= "<<exp(p[1])<<"+/-"<<std_A<<endl; std::cout<<"gamma= "<<p[0]<<"+/-"<<std_g<<endl; - std::cout<<"Num of sources:"<<ds.size()<<endl; + std::cout<<"Num of sources:"<<ds.size()<<endl; } diff --git a/mass_profile/fit_mt_bpl.cpp b/mass_profile/fit_mt_bpl.cpp index 7cf22b2..e956f4f 100644 --- a/mass_profile/fit_mt_bpl.cpp +++ b/mass_profile/fit_mt_bpl.cpp @@ -45,7 +45,7 @@ double shuffle_data(double xc,double xl,double xu) double lxc=log(xc); double lxl=log(xc-xl)-log(xc); double lxu=log(xc+xu)-log(xc); - + if(std_norm_rand()>0) { double result=std::exp(lxc-std::abs(std_norm_rand()*lxl)); @@ -107,7 +107,7 @@ int main(int argc,char* argv[]) } line+=" "; istringstream iss(line); - + if(line[0]=='#') { if(!is_first_nonono) @@ -145,7 +145,7 @@ int main(int argc,char* argv[]) cerr<<line<<endl; continue; } -#endif +#endif if(std::abs(Mu)+std::abs(Ml)<M*.1) { @@ -204,16 +204,16 @@ int main(int argc,char* argv[]) double gamma0l=k0l; double gamma0u=k0u; - + double ampl0l=exp(b0l)*pow(Tb,gamma0l); double ampl0u=exp(b0u)*pow(Tb,gamma0u);; - + ofs_result<<"no no no"<<endl; fitter<double,double,vector<double>,double,std::string> fit; fit.set_opt_method(powell_method<double,vector<double> >()); - + fit.set_statistic(logchisq<double,double,vector<double>,double,std::string>()); //fit.set_statistic(leastsq<double,double,vector<double>,double,std::string>()); fit.set_model(bpl1d<double>()); @@ -222,7 +222,7 @@ int main(int argc,char* argv[]) cerr<<"k0l="<<k0l<<endl; cerr<<"k0u="<<k0u<<endl; cerr<<"Ampl0="<<(ampl0l+ampl0u)/2<<endl; - + fit.set_param_value("bpx",Tb); fit.set_param_value("bpy",(ampl0l+ampl0u)/2); fit.set_param_value("gamma1",gamma0l); @@ -239,7 +239,7 @@ int main(int argc,char* argv[]) ofs_result<<i<<"\t0\t0\t"<<fit.eval_model_raw(i,p)<<"\t0\t0\n"; } - + std::vector<double> mean_p(p.size()); std::vector<double> mean_p2(p.size()); int cnt=0; @@ -252,7 +252,7 @@ int main(int argc,char* argv[]) double sxl=0; double syl=0; double sxyl=0; - + double sxxu=0; double s1u=0; double sxu=0; @@ -260,7 +260,7 @@ int main(int argc,char* argv[]) double sxyu=0; opt_utilities::default_data_set<double,double> ds1; - for(int i=0;i<ds.size();++i) + for(size_t i=0;i<ds.size();++i) { double x=ds.get_data(i).get_x(); double y=ds.get_data(i).get_y(); @@ -268,14 +268,14 @@ int main(int argc,char* argv[]) double xu=ds.get_data(i).get_x_upper_err(); double yl=ds.get_data(i).get_y_lower_err(); double yu=ds.get_data(i).get_y_upper_err(); - + double new_x=shuffle_data(x, xl, xu); double new_y=shuffle_data(y, yl, yu); - + ds1.add_data(data<double,double>(new_x,new_y, yl/y*new_y, yu/y*new_y, @@ -306,40 +306,40 @@ int main(int argc,char* argv[]) double Mbl=sxxl*syl-sxl*sxyl; double k0l=Mal/Ml; double b0l=Mbl/Ml; - + double Mu=sxxu*s1u-sxu*sxu; double Mau=sxyu*s1u-syu*sxu; double Mbu=sxxu*syu-sxu*sxyu; double k0u=Mau/Mu; double b0u=Mbu/Mu; - + double gamma0l=k0l; double gamma0u=k0u; - + double ampl0l=exp(b0l)*pow(Tb,gamma0l); double ampl0u=exp(b0u)*pow(Tb,gamma0u);; - + fit.set_param_value("bpx",Tb); fit.set_param_value("bpy",(ampl0l+ampl0u)/2); fit.set_param_value("gamma1",gamma0l); fit.set_param_value("gamma2",gamma0u); - - + + fit.load_data(ds1); - + fit.fit(); vector<double> p=fit.fit(); - for(int i=0;i<p.size();++i) + for(size_t i=0;i<p.size();++i) { mean_p[i]+=p[i]; mean_p2[i]+=p[i]*p[i]; } //cerr<<fit.get_param_value("gamma1")<<"\t"<<fit.get_param_value("gamma2")<<endl; - + } vector<double> std_p(p.size()); cerr<<endl; - for(int i=0;i<mean_p.size();++i) + for(size_t i=0;i<mean_p.size();++i) { mean_p[i]/=cnt; mean_p2[i]/=cnt; diff --git a/mass_profile/fit_mt_pl.cpp b/mass_profile/fit_mt_pl.cpp index a2246a7..bfccaf4 100644 --- a/mass_profile/fit_mt_pl.cpp +++ b/mass_profile/fit_mt_pl.cpp @@ -92,7 +92,7 @@ int main(int argc,char* argv[]) } line+=" "; istringstream iss(line); - + if(line[0]=='#') { if(!is_first_nonono) @@ -129,7 +129,7 @@ int main(int argc,char* argv[]) cerr<<line<<endl; continue; } -#endif +#endif if(std::abs(Mu)+std::abs(Ml)<M*.1) { double k=M*.1/(std::abs(Mu)+std::abs(Ml)); @@ -170,7 +170,7 @@ int main(int argc,char* argv[]) double Mb=sxx*sy-sx*sxy; double k0=Ma/M; double b0=Mb/M; - + ofs_result<<"no no no"<<endl; fitter<double,double,vector<double>,double,std::string> fit; fit.set_opt_method(powell_method<double,vector<double> >()); @@ -203,12 +203,12 @@ int main(int argc,char* argv[]) ofs_resid<<"skip single"<<endl; ofs_resid<<"ma 3 on 1"<<endl; ofs_resid<<"log x"<<endl; - for(int i=0;i<ds.size();++i) + for(size_t i=0;i<ds.size();++i) { double x=ds.get_data(i).get_x(); double y=ds.get_data(i).get_y(); - double xe1=-ds.get_data(i).get_x_lower_err()*0; - double xe2=ds.get_data(i).get_x_upper_err()*0; + //double xe1=-ds.get_data(i).get_x_lower_err()*0; + //double xe2=ds.get_data(i).get_x_upper_err()*0; double ye1=-ds.get_data(i).get_y_lower_err(); double ye2=ds.get_data(i).get_y_upper_err(); ofs_resid<<exp(x)<<"\t"<<0<<"\t"<<0<<"\t"<<y-fit.eval_model_raw(x,p)<<"\t"<<ye1<<"\t"<<ye2<<"\t"<<"0\t0\t0"<<endl; @@ -223,7 +223,7 @@ int main(int argc,char* argv[]) ++cnt; cerr<<"."; opt_utilities::default_data_set<double,double> ds1; - for(int i=0;i<ds.size();++i) + for(size_t i=0;i<ds.size();++i) { double new_x=shuffle_data(ds.get_data(i).get_x(), ds.get_data(i).get_x_lower_err(), @@ -238,7 +238,7 @@ int main(int argc,char* argv[]) ds.get_data(i).get_y_upper_err())); } fit.load_data(ds1); - + fit.fit(); double k=fit.get_param_value("k"); double b=fit.get_param_value("b"); @@ -260,6 +260,6 @@ int main(int argc,char* argv[]) std::cerr<<"M=M0*T^gamma"<<endl; std::cout<<"M0= "<<exp(p[1])<<"+/-"<<std_A<<endl; std::cout<<"gamma= "<<p[0]<<"+/-"<<std_g<<endl; - std::cout<<"Num of sources:"<<ds.size()<<endl; + std::cout<<"Num of sources:"<<ds.size()<<endl; } diff --git a/mass_profile/fit_nfw_mass.cpp b/mass_profile/fit_nfw_mass.cpp index 48064be..fae74ef 100644 --- a/mass_profile/fit_nfw_mass.cpp +++ b/mass_profile/fit_nfw_mass.cpp @@ -1,7 +1,7 @@ /* Fitting nfw mass profile model Author: Junhua Gu - Last modification 20120721 + Last modification 20120721 */ #include "nfw.hpp" @@ -26,7 +26,7 @@ static double calc_critical_density(double z, { const double G=6.673E-8;//cm^3 g^-1 s^2 const double E=std::sqrt(Omega_m*(1+z)*(1+z)*(1+z)+1-Omega_m); - const double H=H0*E; + const double H=H0*E; return 3*H*H/8/pi/G; } @@ -53,14 +53,14 @@ int main(int argc,char* argv[]) ifstream ifs(argv[1]); //cout<<"read serr 2"<<endl; ofstream ofs_fit_result("nfw_fit_result.qdp"); - + ofs_fit_result<<"read serr 1 2"<<endl; ofs_fit_result<<"skip single"<<endl; ofs_fit_result<<"line off"<<endl; ofs_fit_result<<"li on 2"<<endl; ofs_fit_result<<"li on 4"<<endl; ofs_fit_result<<"ls 2 on 4"<<endl; - + ofs_fit_result<<"win 1"<<endl; ofs_fit_result<<"yplot 1 2"<<endl; ofs_fit_result<<"loc 0 0 1 1"<<endl; @@ -108,7 +108,7 @@ int main(int argc,char* argv[]) vector<double> p=fit.fit(); //output parameters ofstream ofs_param("nfw_param.txt"); - for(int i=0;i<fit.get_num_params();++i) + for(size_t i=0;i<fit.get_num_params();++i) { cout<<fit.get_param_info(i).get_name()<<"\t"<<fit.get_param_info(i).get_value()<<endl; ofs_param<<fit.get_param_info(i).get_name()<<"\t"<<fit.get_param_info(i).get_value()<<endl; @@ -117,11 +117,11 @@ int main(int argc,char* argv[]) ofs_param<<"reduced chi^2="<<fit.get_statistic_value()<<endl; ofstream ofs_model("nfw_dump.qdp"); ofstream ofs_overdensity("overdensity.qdp"); - const double G=6.673E-8;//cm^3 g^-1 s^-2 - static const double mu=1.4074; - static const double mp=1.67262158E-24;//g + //const double G=6.673E-8;//cm^3 g^-1 s^-2 + //static const double mu=1.4074; + //static const double mp=1.67262158E-24;//g static const double M_sun=1.98892E33;//g - static const double k=1.38E-16; + //static const double k=1.38E-16; double xmax=0; for(double x=std::max(rmin_kpc,ds.get_data(0).get_x());;x+=1) @@ -141,7 +141,7 @@ int main(int argc,char* argv[]) } } ofs_fit_result<<"no no no"<<endl; - for(int i=0;i<ds.size();++i) + for(size_t i=0;i<ds.size();++i) { data<double,double> d=ds.get_data(i); double x=d.get_x(); @@ -155,5 +155,5 @@ int main(int argc,char* argv[]) { ofs_fit_result<<x<<"\t"<<0<<"\t"<<0<<"\t"<<0<<endl; } - + } diff --git a/mass_profile/fit_nfw_sbp.cpp b/mass_profile/fit_nfw_sbp.cpp index 6343dea..a7d0537 100644 --- a/mass_profile/fit_nfw_sbp.cpp +++ b/mass_profile/fit_nfw_sbp.cpp @@ -46,7 +46,7 @@ public: */ return spl.get_value(x); } - + //we need this function, when this object is performing a clone of itself spline_func_obj* do_clone()const { @@ -59,7 +59,7 @@ public: { spl.push_point(x,y); } - + //before getting the intepolated value, the spline should be initialzied by calling this function void gen_spline() { @@ -137,9 +137,9 @@ int main(int argc,char* argv[]) this will be convenient to calculate the volume of each spherical shell, and can naturally ensure the annuli are adjacent with each other, with out any gaps. - + */ - + std::vector<double> radii;//to store radius std::vector<double> sbps;//to store the surface brightness value std::vector<double> sbpe;//to store the sbp error @@ -148,8 +148,8 @@ int main(int argc,char* argv[]) About the format of the radius file: the radius file contains only radius, separated by space or line feed (i.e., the <ENTER> key). the unit should be pixel - - The number of radius can be larger than the number of annuli+1, the exceeded radius can be used + + The number of radius can be larger than the number of annuli+1, the exceeded radius can be used to calculate the influence of outer shells. */ int ncut=0; @@ -208,7 +208,7 @@ int main(int argc,char* argv[]) for(ifstream ifs(arg_map["cfunc_file"].c_str());;) { assert(ifs.is_open()); - double x,y,y1,y2; + double x,y; ifs>>x>>y; if(!ifs.good()) { @@ -250,7 +250,7 @@ int main(int argc,char* argv[]) nfw.set_cm_per_pixel(cm_per_pixel); //define the temperature profile model spline_func_obj tf; - + for(ifstream ifs_tfunc(arg_map["T_file"].c_str());;) { assert(ifs_tfunc.is_open()); @@ -286,8 +286,8 @@ int main(int argc,char* argv[]) f.set_param_value("n0",n0); f.set_param_value("rho0",rho0); f.set_param_value("rs",rs); - - + + f.set_param_value("bkg",bkg); cout<<f.get_data_set().size()<<endl; @@ -302,8 +302,8 @@ int main(int argc,char* argv[]) #endif //output the parameters ofstream param_output("nfw_param.txt"); - - for(int i=0;i<f.get_num_params();++i) + + for(size_t i=0;i<f.get_num_params();++i) { cout<<f.get_param_info(i).get_name()<<"\t"<<abs(f.get_param_info(i).get_value())<<endl; param_output<<f.get_param_info(i).get_name()<<"\t"<<abs(f.get_param_info(i).get_value())<<endl; @@ -326,7 +326,7 @@ int main(int argc,char* argv[]) //output the surface brightness profile ofs_sbp<<"read serr 2"<<endl; ofs_sbp<<"skip single"<<endl; - for(int i=1;i<sbps.size();++i) + for(size_t i=1;i<sbps.size();++i) { double x=(radii[i]+radii[i-1])/2; double y=sbps[i-1]; @@ -338,7 +338,7 @@ int main(int argc,char* argv[]) //output the electron density mv=nfw.eval(radii,p); ofstream ofs_rho("rho_fit.qdp"); - for(int i=1;i<sbps.size();++i) + for(size_t i=1;i<sbps.size();++i) { double x=(radii[i]+radii[i-1])/2; double ym=mv[i-1]; @@ -354,7 +354,7 @@ int main(int argc,char* argv[]) } //calculate the overdensity profile ofstream ofs_overdensity("overdensity.qdp"); - + std::vector<double> radius_list; std::vector<double> delta_list; @@ -380,8 +380,8 @@ int main(int argc,char* argv[]) */ ofs_overdensity<<r/kpc<<"\t"<<delta<<endl; } - - for(int i=0;i<radius_list.size()-1;++i) + + for(size_t i=0;i<radius_list.size()-1;++i) { double r=radius_list[i]; if(delta_list[i]>=200&&delta_list[i+1]<200) diff --git a/mass_profile/fit_wang2012_model.cpp b/mass_profile/fit_wang2012_model.cpp index 3fe14b6..f15a2b1 100644 --- a/mass_profile/fit_wang2012_model.cpp +++ b/mass_profile/fit_wang2012_model.cpp @@ -2,7 +2,7 @@ Fitting Jy Wang's temperature profile model Author: Jingying Wang Last modification 20120819 - + */ #include "wang2012_model.hpp" @@ -34,7 +34,7 @@ int main(int argc,char* argv[]) { cm_per_pixel=atof(argv[3]); } - + //define the fitter fitter<double,double,vector<double>,double,std::string> fit; //define the data set @@ -88,7 +88,7 @@ int main(int argc,char* argv[]) fit.set_statistic(chisq_object); //fit.set_statistic(chisq<double,double,vector<double>,double,std::string>()); fit.set_model(wang2012_model<double>()); - + if(argc>=3&&std::string(argv[2])!="NONE") { std::vector<std::string> freeze_list; @@ -123,7 +123,7 @@ int main(int argc,char* argv[]) { freeze_param<double,double,std::vector<double>,std::string> fp(freeze_list[0]); fit.set_param_modifier(fp); - for(int i=1;i<freeze_list.size();++i) + for(size_t i=1;i<freeze_list.size();++i) { dynamic_cast<freeze_param<double,double,std::vector<double>,std::string>&>(fit.get_param_modifier())+=freeze_param<double,double,std::vector<double>,std::string>(freeze_list[i]); } @@ -147,11 +147,11 @@ int main(int argc,char* argv[]) } #endif //output parameters - for(int i=0;i<fit.get_num_params();++i) + for(size_t i=0;i<fit.get_num_params();++i) { std::string pname=fit.get_param_info(i).get_name(); std::string pstatus=fit.get_model().report_param_status(pname); - + if(pstatus==""||pstatus=="thawed") { pstatus="T"; @@ -168,20 +168,20 @@ int main(int argc,char* argv[]) } #endif } - + //cout<<"T0"<<"\t"<<fit.get_param_value("T0")<<endl; //cout<<"T1"<<"\t"<<fit.get_param_value("T1")<<endl; //cout<<"xt"<<"\t"<<fit.get_param_value("xt")<<endl; //cout<<"eta"<<"\t"<<fit.get_param_value("eta")<<endl; //dump the data for checking ofstream ofs_model("wang2012_dump.qdp"); - - + + min_r=0; for(double x=min_r;x<3000;x+=10) { double model_value=fit.eval_model_raw(x,p); - + ofs_model<<x<<"\t"<<model_value<<endl; if(cm_per_pixel>0) { diff --git a/mass_profile/nfw_ne.hpp b/mass_profile/nfw_ne.hpp index 5842bb8..8f23820 100644 --- a/mass_profile/nfw_ne.hpp +++ b/mass_profile/nfw_ne.hpp @@ -36,7 +36,7 @@ namespace opt_utilities { return rho0; } - + return nfw_mass_enclosed(r,rho0,rs)/(4.*pi/3*r*r*r); } @@ -46,7 +46,7 @@ namespace opt_utilities const double Omega_m=.27) { const double E=std::sqrt(Omega_m*(1+z)*(1+z)*(1+z)+1-Omega_m); - const double H=H0*E; + const double H=H0*E; return 3*H*H/8/pi/G; } @@ -66,7 +66,7 @@ namespace opt_utilities nfw_ne() :pTfunc(0),cm_per_pixel(1) { - + this->push_param_info(param_info<std::vector<T>,std::string>("rho0",1));//in mp this->push_param_info(param_info<std::vector<T>,std::string>("rs",100)); this->push_param_info(param_info<std::vector<T>,std::string>("n0",.01)); @@ -89,7 +89,7 @@ namespace opt_utilities this->push_param_info(param_info<std::vector<T>,std::string>("rs",rhs.get_param_info("rs").get_value())); this->push_param_info(param_info<std::vector<T>,std::string>("n0",rhs.get_param_info("n0").get_value())); } - + //assignment operator nfw_ne& operator=(const nfw_ne& rhs) { @@ -146,19 +146,19 @@ namespace opt_utilities */ std::vector<T> do_eval(const std::vector<T> & r, const std::vector<T>& p) - { + { assert(pTfunc); //const T kT_erg=k*5; T rho0=std::abs(p[0])*mp; T rs=std::abs(p[1]); T n0=std::abs(p[2]); T rs_cm=rs*cm_per_pixel; - + std::vector<T> yvec(r.size()); const T kT_erg0=pTfunc->eval((r.at(0)+r.at(1))/2)*k; //calculate the integration #pragma omp parallel for schedule(dynamic) - for(int i=0;i<r.size();++i) + for(size_t i=0;i<r.size();++i) { T r_cm=r[i]*cm_per_pixel; T kT_erg=pTfunc->eval(r[i])*k; @@ -170,10 +170,10 @@ namespace opt_utilities //std::cout<<r_cm/1e20<<"\t"<<nfw_mass_enclosed(r_cm,rho0,rs_cm)/1e45<<std::endl; //std::cout<<r_cm/1e20<<"\t"<<G*nfw_mass_enclosed(r_cm,rho0,rs_cm)*mu*mp/kT_erg/r_cm/r_cm<<std::endl; } - + std::vector<T> ydxvec(r.size()-1); #pragma omp parallel for schedule(dynamic) - for(int i=1;i<r.size();++i) + for(size_t i=1;i<r.size();++i) { T dr=r[i]-r[i-1]; T dr_cm=dr*cm_per_pixel; @@ -183,7 +183,7 @@ namespace opt_utilities //construct the result std::vector<T> result(r.size()-1); #pragma omp parallel for schedule(dynamic) - for(int i=0;i<r.size()-1;++i) + for(size_t i=0;i<r.size()-1;++i) { T y=-ydxvec.at(i); T kT_erg=pTfunc->eval(r[i])*k; diff --git a/mass_profile/projector.hpp b/mass_profile/projector.hpp index a9d5f2a..6f2f192 100644 --- a/mass_profile/projector.hpp +++ b/mass_profile/projector.hpp @@ -99,7 +99,7 @@ namespace opt_utilities void attach_model(const model<std::vector<T>,std::vector<T>,std::vector<T> >& m) { this->clear_param_info(); - for(int i=0;i<m.get_num_params();++i) + for(size_t i=0;i<m.get_num_params();++i) { this->push_param_info(m.get_param_info(i)); } @@ -167,7 +167,7 @@ namespace opt_utilities } } std::vector<T> p2(p1.size()-1); - for(int i=0;i<p1.size()-1;++i) + for(size_t i=0;i<p1.size()-1;++i) { p2.at(i)=p1[i]; } @@ -183,10 +183,10 @@ namespace opt_utilities std::vector<T> unprojected(pmodel->eval(x,p)); std::vector<T> projected(unprojected.size()); - for(int nrad=0;nrad<x.size()-1;++nrad) + for(size_t nrad=0;nrad<x.size()-1;++nrad) { double v1=0; - for(int nsph=nrad;nsph<x.size()-1;++nsph) + for(size_t nsph=nrad;nsph<x.size()-1;++nsph) { double v=calc_v(x,nsph,nrad)*cm_per_pixel*cm_per_pixel*cm_per_pixel; v1+=v; diff --git a/mass_profile/vchisq.hpp b/mass_profile/vchisq.hpp index 0780ce8..c3d38b8 100644 --- a/mass_profile/vchisq.hpp +++ b/mass_profile/vchisq.hpp @@ -28,7 +28,7 @@ namespace opt_utilities int n;
bool limit_bound;
typedef std::vector<T> Tp;
-
+
vchisq<T>* do_clone()const
{
return new vchisq<T>(*this);
@@ -38,7 +38,7 @@ namespace opt_utilities {
return "chi^2 statistic";
}
-
+
public:
void verbose(bool v)
{
@@ -49,7 +49,7 @@ namespace opt_utilities {
limit_bound=true;
}
-
+
void clear_limit()
{
limit_bound=false;
@@ -58,8 +58,8 @@ namespace opt_utilities vchisq()
:verb(false),limit_bound(false)
{}
-
-
+
+
T do_eval(const std::vector<T>& p)
{
@@ -71,7 +71,7 @@ namespace opt_utilities }
}
T result(0);
-
+
std::vector<float> vx;
std::vector<float> vy;
std::vector<float> vye;
@@ -80,7 +80,7 @@ namespace opt_utilities if(verb)
{
n++;
-
+
if(n%100==0)
{
vx.resize(this->get_data_set().get_data(0).get_y().size());
@@ -88,24 +88,24 @@ namespace opt_utilities vye.resize(this->get_data_set().get_data(0).get_y().size());
my.resize(this->get_data_set().get_data(0).get_y().size());
}
-
+
}
for(int i=(this->get_data_set()).size()-1;i>=0;--i)
{
const std::vector<double> y_model(this->eval_model(this->get_data_set().get_data(i).get_x(),p));
const std::vector<double>& y=this->get_data_set().get_data(i).get_y();
const std::vector<double>& ye=this->get_data_set().get_data(i).get_y_lower_err();
- for(int j=0;j<y.size();++j)
+ for(size_t j=0;j<y.size();++j)
{
double chi=(y_model[j]-y[j])/ye[j];
result+=chi*chi;
}
-
-
+
+
if(verb&&n%100==0)
{
-
- for(int j=0;j<y.size();++j)
+
+ for(size_t j=0;j<y.size();++j)
{
vx.at(j)=((this->get_data_set().get_data(i).get_x().at(j)+this->get_data_set().get_data(i).get_x().at(j+1))/2.);
vy.at(j)=(y[j]);
@@ -121,9 +121,9 @@ namespace opt_utilities my[j]=log10(my[j]);
}
-
+
}
-
+
}
if(verb)
{
@@ -148,7 +148,7 @@ namespace opt_utilities pr.plot_line(vx,my);
}
}
-
+
return result;
}
};
|