#include"Distribution.h"
#include"Skymap.h"        //AWS20090107
#include<vector>          //AWS20100817
#include<valarray>        //AWS20100817

class SourcePopulation
{
 public:
  // parameters of population
  // luminosity in photons s-1 in given energy band

  double pi, rtd, dtr,  kpc_to_cm;

  char title[200]; //AWS20140102  was 100

  char galdef_ID [200]; //AWS20200624
  char psfile_tag[200]; //AWS20200624


  int luminosity_function_type;// luminosity-function: 1=delta-function, 2=power-law
  double L;  // luminosity if delta-function

  //  double density0,alpha_R,alpha_z;      // parameters of exponential disk
  double density0,alpha_R,beta_R,zscale;    // parameters of source distribution
  double alpha_z;     
 
  double oversample;//  source oversampling to avoid statistical fluctuations   AWS20060118


  double L_min,L_max,alpha_L;           // parameters of power-law luminosity function

  int    spectral_model;                // 1 = 2 breaks fixed,  2 = 2 breaks with dispersion AWS20100824

  double spectrum_g,spectrum_norm_E;    // spectral index and normalization  energy of sources
  double spectrum_g_0 ;//       spectrum with two breaks                           AWS20051219
  double spectrum_br0 ;
  double spectrum_g_1 ;
  double spectrum_br1 ;
  double spectrum_g_2 ;


  double spectrum_g_0_sigma ;//       spectrum with two breaks: Gaussian distribution AWS20100823
  double spectrum_br0_sigma ;
  double spectrum_g_1_sigma ;
  double spectrum_br1_sigma ;
  double spectrum_g_2_sigma ;

  int   n_E_bin;                        //     number of energy bins
  double *E_bin_low,*E_bin_high;        // boundaries of energy bins
  double  E_ref_low, E_ref_high;        // range for which luminosity is defined


  int   n_E;                           // number of energies for skymap spectra
  double *E;                           //           energies for skymap spectra

  double NL(double L);
  double density(double R, double z); // number per kpc^3

  // derived quantities
  long  n_sample;
  double *sample_R, *sample_x, *sample_y, *sample_z;
  double *sample_theta,*sample_d,*sample_l,*sample_b,*sample_flux;
  double *sample_flux100,*sample_flux1000,*sample_flux10000;               //AWS20120615 integral fluxes
  double *sample_flux_spectrum_E_ref_low,*sample_flux_spectrum_E_ref_high; //AWS20111001
  double sample_d_min,sample_d_max,sample_flux_min,sample_flux_max;
  double sample_total_flux,sample_total_L;
  double sample_selected_total_flux, sample_selected_total_flux_sublimit, sample_selected_total_flux_soplimit;
  double sample_selected_sourcecnt , sample_selected_sourcecnt_sublimit, sample_selected_sourcecnt_soplimit;
  double sample_selected_d_min,         sample_selected_d_max,         sample_selected_d_av;
  double sample_selected_d_min_sublimit,sample_selected_d_max_sublimit,sample_selected_d_av_sublimit;
  double sample_selected_d_min_soplimit,sample_selected_d_max_soplimit,sample_selected_d_av_soplimit;
  double  *sample_l_min,*sample_b_min, *sample_l_max,*sample_b_max, *sample_l_av,*sample_b_av;
  

  double **sample_flux_spectra, **sample_binned_flux_spectra; //AWS20100817 
  double *sample_spectrum_g_0,*sample_spectrum_br0,*sample_spectrum_g_1,*sample_spectrum_br1,*sample_spectrum_g_2; //AWS20100823 parameters per source

  double   latprof_l_min, latprof_l_max, latprof_dlat ; // for profiles
  double  longprof_b_min,longprof_b_max,longprof_dlong;

  double long_min1,long_max1; // for logN-logS
  double long_min2,long_max2;
  double  lat_min1, lat_max1;
  double  lat_min2, lat_max2;

  double *longprof_intensity,*longprof_intensity_sublimit ;
  double *longprof_sourcecnt,*longprof_sourcecnt_sublimit ;
  int     longprof_nlong;
  double * latprof_intensity,* latprof_intensity_sublimit ;
  double * latprof_sourcecnt,* latprof_sourcecnt_sublimit ;
  int      latprof_nlat; 
  
  double flux_detection_limit;
  double flux_detection_limit_sensitivity_file; //AWS20200713

  long i_flux_min,i_flux_max;



  double *dlnN_dlnS;
  double *dlnN_dlnS_deg;

  double *dlnN_dlnS_int;
  double *dlnN_dlnS_int_deg;

  double *FS,*FS_int;

  double lnS_min,lnS_max,dlnS;
  double dlnN_dlnS_total_N, dlnN_dlnS_total_S;
  int   n_dlnN_dlnS;

 
  double *dlnN_dlnS_soplimit;  // N(S) above flux limit  AWS20170123
  double *dlnN_dlnS_sublimit;  // N(S) below flux limit  AWS20170123

  double *dlnN_dlnS_int_soplimit;                      //AWS20170123
  double *dlnN_dlnS_int_sublimit;                      //AWS20170123

  double *FS_soplimit,*FS_int_soplimit;                //AWS20170123
  double *FS_sublimit,*FS_int_sublimit;                //AWS20170123

  double R0 ;// Sun-Galactic centre in kpc

// skymaps of source counts N(S)
  Distribution skymap_sourcecnt;
  Distribution skymap_sourcecnt_sublimit;
  Distribution skymap_sourcecnt_soplimit;

// skymaps integrated over energy ranges
  Distribution skymap_intensity;
  Distribution skymap_intensity_sublimit;
  Distribution skymap_intensity_soplimit;

// skymaps for full energy spectrum
  Distribution skymap_intensity_spectrum;
  Distribution skymap_intensity_spectrum_sublimit;
  Distribution skymap_intensity_spectrum_soplimit;

// ---  healpix versions of the same skymaps ---
  int healpix_skymap_order; // order of healpix maps

// skymaps of source counts N(S)
  Skymap  <double> healpix_skymap_sourcecnt;
  Skymap  <double> healpix_skymap_sourcecnt_sublimit;
  Skymap  <double> healpix_skymap_sourcecnt_soplimit;

// skymaps integrated over energy ranges
  Skymap  <double> healpix_skymap_intensity;
  Skymap  <double> healpix_skymap_intensity_sublimit;
  Skymap  <double> healpix_skymap_intensity_soplimit;

// skymaps for full energy spectrum
  Skymap  <double> healpix_skymap_intensity_spectrum;
  Skymap  <double> healpix_skymap_intensity_spectrum_sublimit;
  Skymap  <double> healpix_skymap_intensity_spectrum_soplimit;

// Fermi-LAT source sensitivity maps                       AWS20170112
  Healpix_Map <double> healpix_skymap_Fermi_sensitivity; //AWS20170112
  char* Fermi_sensitivity_file;                          //AWS20170112

  int verbose;

  ofstream *txt_stream;

  // control idl plotting
  int IDL_plot_control;                                  //AWS20200623
  int IDL_psym;                                          //AWS20200623

  int init();
  int print(int options=1);
  int gen_sample(unsigned seed); 
  int analyse_sample();
  int gen_sample_EGRET_catalogue(int options);
  int gen_sample_SPI_catalogue  (char *directory, char *filename, int options);
  int gen_sample_Fermi_catalogue(char *directory, char *filename, int options);//AWS20081215
  
  int gen_Fermi_catalogue_from_sample(char *directory, char *filename, int options);//AWS20110829

  int read_sample_MSPs               (char *directory, char *filename, int options);//AWS20120521


  SourcePopulation operator=(SourcePopulation); // copy constructor

  double integratepowbr(double A0,double g0,double g1,double g2, double br0,double br1, double E1,double E2);
  double     shapepowbr(double A0,double g0,double g1,double g2, double br0,double br1, double E);

  double flux_cutoff_powerlaw(double E                , double g, double Ecutoff,double Emin_ref, double Emax_ref, double flux_ref);
  double flux_cutoff_powerlaw(double Emin, double Emax, double g, double Ecutoff,double Emin_ref, double Emax_ref, double flux_ref);

 };