using namespace std;
#include<iostream>
#include<fstream>
#include"SourcePopulation.h"
#include"EGRET_catalogue.h"
#include"CAT.h"
#include"Sources.h" //AWS20081215
#include"pointing.h"// healpix class//AWS20090108
#include "fitsio.h" //AWS20120521

#include<valarray> //AWS20090902

#include <gsl/gsl_rng.h>            //AWS20100824
#include <gsl/gsl_randist.h>        //AWS20100824
#include <gsl/gsl_sf.h>             //AWS20100827

#include"slalib.h"                  //AWS20120612

#include "healpix_map.h"            //AWS20170112   
#include "healpix_map_fitsio.h"     //AWS20170112
#include "pointing.h"               //AWS20170112

int poisson(double  x);
int poisson(unsigned int initialize_seed);
double power_law_sampler(double  x_min, double x_max, double g);
int power_law_sampler(unsigned int initialize_seed);

////////////////////////////////////////////////////////////

int  SourcePopulation::init()
{

  cout<< "SourcePopulation::init"<<endl;


  pi        =  acos(-1.);
  rtd       =  180./pi;  // radians to degrees
  dtr       =   pi/180.; // degrees to radians
  kpc_to_cm =  3.08e21;

  R0 = 8.5 ;// Sun-Galactic centre in kpc
 
  strcpy(title,"no title");

  IDL_plot_control =  1; //AWS20200623
  IDL_psym         = -6; //AWS20200623

  strcpy(galdef_ID, "sourcepop_needs_galdef_ID" ); //AWS20200701
  strcpy(psfile_tag,"sourcepop_needs_psfile_tag"); //AWS20200701


  //AWS20170112  read and test sensitivity map
  {     
   cout<<"  reading and testing Fermi sensitivity skymap"<<endl;
   string mapname;
          mapname = "../FITS/"+ string(Fermi_sensitivity_file);

    cout<<"  reading  healpix skymap Fermi sensitivity from file "<<mapname<<endl;
    int colnum=1, hdunum=2;
    read_Healpix_map_from_fits(mapname,healpix_skymap_Fermi_sensitivity,colnum,hdunum);
    int pix=100;
    cout<<"value="<<healpix_skymap_Fermi_sensitivity[pix]<<endl;
   
  
   
    for (double b=-90;b<  90; b+=10) 
    for (double l=0  ;l<=360; l+=10) 
    {
     double theta=(90-b)*dtr;
     double phi  =    l *dtr;
     pointing pointing_=pointing(theta,phi);
     cout<<" Fermi sensitivity  l="<< l<<" b="<<b<<" theta="<<theta<<" phi="<<phi;
     cout<<" interpolated value="<<healpix_skymap_Fermi_sensitivity.interpolated_value(pointing_)<<endl;
    }
  }

  return 0;
}



//////////////////////////////////////////////////////

int  SourcePopulation::print(int options) // default options=1 in header
{
  int i;

  *txt_stream<< "SourcePopulation::print, options = "<<options<<endl;

 if(options==1)
 {



  *txt_stream<<" ***** Source Population parameters *****"<<endl<<endl;
  *txt_stream<< "title : "<< title        <<endl;
  *txt_stream<<"luminosity_function_type="<<luminosity_function_type<<endl;
  if(luminosity_function_type==1)
  *txt_stream<<"delta function L="<<L<<endl;
  if(luminosity_function_type==2)
  *txt_stream<<"power-law L_min="<<L_min <<" L_max="<<L_max <<" alpha_L="<<alpha_L <<endl;

  *txt_stream<<"alpha_R="<<alpha_R<<endl;
  *txt_stream<<" beta_R="<< beta_R<<endl;
  *txt_stream<<"zscale ="<<zscale <<endl;
  *txt_stream<<"alpha_z="<<alpha_z<<endl;
  *txt_stream<<"density0="<<density0<<endl;
  *txt_stream<<"oversample="<<oversample<<endl;


  *txt_stream<<"spectral_model="<<spectral_model<<endl;
  *txt_stream<<"spectrum_g_0="<<spectrum_g_0<<endl;
  *txt_stream<<"spectrum_br0="<<spectrum_br0<<endl;
  *txt_stream<<"spectrum_g_1="<<spectrum_g_1<<endl;
  *txt_stream<<"spectrum_br1="<<spectrum_br1<<endl;
  *txt_stream<<"spectrum_g_2="<<spectrum_g_2<<endl;

  *txt_stream<<"spectrum_g_0_sigma="<<spectrum_g_0_sigma<<endl;
  *txt_stream<<"spectrum_br0_sigma="<<spectrum_br0_sigma<<endl;
  *txt_stream<<"spectrum_g_1_sigma="<<spectrum_g_1_sigma<<endl;
  *txt_stream<<"spectrum_br1_sigma="<<spectrum_br1_sigma<<endl;
  *txt_stream<<"spectrum_g_2_sigma="<<spectrum_g_2_sigma<<endl;

  *txt_stream<<"user-defined source flux detection limit="<<flux_detection_limit<<endl;

  *txt_stream<<"longitude range 1 "<<long_min1<<" to "<<long_max1<<endl;
  *txt_stream<<"longitude range 2 "<<long_min2<<" to "<<long_max2<<endl;
  *txt_stream<<" latitude range 1 "<< lat_min1<<" to "<< lat_max1<<endl;
  *txt_stream<<" latitude range 2 "<< lat_min2<<" to "<< lat_max2<<endl;

  *txt_stream<<"number of spectral energies ="<<n_E<<endl;
  *txt_stream<<"spectral energies:" ; for (i=0;i<n_E;i++) *txt_stream<<E[i]<<" "; *txt_stream<<endl;

  *txt_stream<<"healpix skymap order="<<healpix_skymap_order<<endl;

  *txt_stream<<"Fermi sensitivity file="<< Fermi_sensitivity_file<<endl; //AWS201701113

  *txt_stream<<"IDL_plot_control="<<IDL_plot_control<<endl; //AWS20200623
  *txt_stream<<"IDL_psym="        <<IDL_psym<<endl;         //AWS20200623

  *txt_stream<<"galdef_ID=" <<galdef_ID<<endl;              //AWS20200701
  *txt_stream<<"psfile_tag="<<psfile_tag<<endl;             //AWS20200701


 
}



 if(options==2)
 {
   *txt_stream<<"source sample listing"<<endl<<endl;

  for(int i=0;i<n_sample;i++)

	{
	  *txt_stream<<"source "<<i<<":  R=" <<sample_R[i]<< " z="<<sample_z[i]<<" theta="<<sample_theta[i]
                                      <<" l="<<sample_l[i] <<" b="<<sample_b[i]
		     <<" distance="<<sample_d[i]    <<" L="<<L <<" flux="<<sample_flux[i];

          if(spectral_model==2||spectral_model==3)
           *txt_stream<<" g_0="<<sample_spectrum_g_0[i]<< " br0="<<sample_spectrum_br0[i]
                      <<" g_1="<<sample_spectrum_g_1[i]<< " br1="<<sample_spectrum_br1[i]<<" g_2="<<sample_spectrum_g_2[i];

          *txt_stream << endl <<" flux spectrum= " ;              for (int iE=0;iE<n_E    ;iE++){*txt_stream<< sample_flux_spectra       [i][iE]<<" ";}
          *txt_stream << endl <<" binned flux spectrum= " ;       for (int iE=0;iE<n_E_bin;iE++){*txt_stream<< sample_binned_flux_spectra[i][iE]<<" ";}

          *txt_stream      <<endl<<endl;
	}
  }

 if(options==3)
 {
   *txt_stream<<"source sample listing"<<endl<<endl;
   *txt_stream<<"l    b   flux>100 MeV index break index break index"<<endl<<endl;

  for(int i=0;i<n_sample;i++)

	{
	  *txt_stream<<  sample_l[i] <<" "<<sample_b[i]	<<" "<<sample_flux[i];

          if(spectral_model==2||spectral_model==3)
           *txt_stream<<"     "<<sample_spectrum_g_0[i]<< "     "<<sample_spectrum_br0[i]
                      <<"     "<<sample_spectrum_g_1[i]<< "     "<<sample_spectrum_br1[i]<<"     "<<sample_spectrum_g_2[i];

          *txt_stream      <<endl      ;
	}
  }



  *txt_stream<<"SourcePopulation:: print complete"<<endl<<endl;

  return 0;
}
///////////////////////////////////////////////////
double  SourcePopulation::density(double R, double z)
{
  
  double density_;
  //  density_=density0*exp(-abs(z)/alpha_z)*exp(-(R-R0)/alpha_R);

  //zscale=alpha_z;
  //beta_R=1./alpha_R;
  density_=density0*exp(-abs(z)/zscale) * pow(R/R0,alpha_R) * exp(-beta_R*(R-R0)/R0); //AWS20051202 for consistency with galprop

  return density_;

  return 0;
}
///////////////////////////////////////////////////
int    SourcePopulation::gen_sample(unsigned  seed)
{
  *txt_stream<< "SourcePopulation::gen_sample"<<endl;
  cout       << "SourcePopulation::gen_sample"<<endl;

  // general purpose case for any distribution
  double Rmax, zmax, dR, dz;
  double R1,R2,R_,z1,z2,z_,sample_dd;
  double theta, dtheta;
  double sinl, cosl;

  int nR, nz, iR, iz;
  int m_sample,m_sample_total;
  int i,ii;
  double density_,number_total,number_delta ;
  double solid_angle,solid_angle_long,solid_angle_lat,solid_angle_long_lat;
  double l,b,b1,b2;
  int il,ib;


 

  Rmax= 20.;
  zmax= 0.3;
  dR=0.201; // avoid having R0 on grid
  dz=0.01;


  nR=int(  Rmax/dR);
  nz=int(2*zmax/dz);
  
 

  
  number_total=0.;

  for (iR=0;iR<nR;iR++)
  for (iz=0;iz<nz;iz++)
    {
      R1  =         iR*dR; R2=R1+dR;
      z1  = -zmax + iz*dz; z2=z1+dz;
      R_=(R1+R2)/2;
      z_=(z1+z2)/2;
      density_=density(R_,z_);
      density_ *= oversample;                       //AWS20060118
      number_delta=density_ *pi*(R2*R2 - R1*R1) *dz;
      number_total+=number_delta;
    }


  *txt_stream<<"total number of sources ="<<number_total<<endl;

  n_sample    =long(number_total);  
  sample_R    =new double[n_sample];
  sample_z    =new double[n_sample];
  sample_theta=new double[n_sample];
  sample_l    =new double[n_sample];
  sample_b    =new double[n_sample];
  sample_d    =new double[n_sample];
  sample_flux =new double[n_sample];

  sample_flux100  =new double[n_sample]; //AWS20120615
  sample_flux1000 =new double[n_sample]; //AWS20120615
  sample_flux10000=new double[n_sample]; //AWS20120615

  sample_flux_spectrum_E_ref_low  =new double[n_sample]; //AWS20111001
  sample_flux_spectrum_E_ref_high =new double[n_sample]; //AWS20111001

   sample_flux_spectra        = new double*[n_sample];                                                                          //AWS20100817
   sample_binned_flux_spectra = new double*[n_sample];   ;                                                                      //AWS20100817
   for (int i=0;i<n_sample;i++){ sample_flux_spectra[i]=new double[n_E];  sample_binned_flux_spectra[i]=new double[n_E_bin];   }//AWS20100817

   sample_spectrum_g_0        = new double[n_sample];                                                                           //AWS20100823
   sample_spectrum_br0        = new double[n_sample];                                                                           //AWS20100823
   sample_spectrum_g_1        = new double[n_sample];                                                                           //AWS20100823
   sample_spectrum_br1        = new double[n_sample];                                                                           //AWS20100823
   sample_spectrum_g_2        = new double[n_sample];                                                                           //AWS20100823

  // initialize gsl for Gaussian sampling
  // see example in gsl User Manual

  cout<<" initializing gsl random number generator"<<endl;
      const gsl_rng_type * T;
      gsl_rng * r;
         
       /* create a generator chosen by the environment variable GSL_RNG_TYPE */
     
       gsl_rng_env_setup();
       T = gsl_rng_default;
       r = gsl_rng_alloc (T);

  cout<<" initialized  gsl random number generator"<<endl;

   // case of no dispersion in spectra
  if(spectral_model==1)
  {
   for (int i=0;i<n_sample;i++)                                                                                                 //AWS20100823
   {
    sample_spectrum_g_0[i] = spectrum_g_0;
    sample_spectrum_br0[i] = spectrum_br0;
    sample_spectrum_g_1[i] = spectrum_g_1;
    sample_spectrum_br1[i] = spectrum_br1;
    sample_spectrum_g_2[i] = spectrum_g_2;
   }
  }

   // case of uncorrelated dispersion in spectra
  if(spectral_model==2)
  {
    
   for (int i=0;i<n_sample;i++)                                                                                                 //AWS20100823
   {
    sample_spectrum_g_0[i] = spectrum_g_0 +  gsl_ran_gaussian ( r,  spectrum_g_0_sigma);

                             
    sample_spectrum_br0[i] = spectrum_br0 * pow(10.,  gsl_ran_gaussian ( r,  spectrum_br0_sigma)); // sigma for break is log10

    sample_spectrum_g_1[i] = sample_spectrum_g_0[i] + spectrum_g_1 - spectrum_g_0;    // keep size of spectral break
    sample_spectrum_br1[i] = spectrum_br1;                                            // keep second break
    sample_spectrum_g_2[i] = sample_spectrum_g_0[i] + spectrum_g_2 - spectrum_g_0;    // keep size of second break

    if(verbose==1)
    {
    cout<<" sample_spectrum_g_0="<<sample_spectrum_g_0[i]<<" sample_spectrum_br0="<<sample_spectrum_br0[i];
    cout<<" sample_spectrum_g_1="<<sample_spectrum_g_1[i]<<" sample_spectrum_br1="<<sample_spectrum_br1[i];
    cout<<" sample_spectrum_g_2="<<sample_spectrum_g_2[i]<<endl;
    }
   }
  }

   // case of correlated index-break dispersion in spectra
  if(spectral_model==3)
  {
   
    double X,Y,correlation_angle;
    double pi=acos(-1.0);
    correlation_angle=19.9; // degrees, from Alice
    double s=sin(correlation_angle*pi/180.);
    double c=cos(correlation_angle*pi/180.);

 

   for (int i=0;i<n_sample;i++)                                                                                               
   {
    X =   gsl_ran_gaussian ( r,  spectrum_g_0_sigma); // offset in rotated system centred on mean
    Y =   gsl_ran_gaussian ( r,  spectrum_br0_sigma);

    sample_spectrum_g_0[i] =          spectrum_g_0  +       X * c - Y * s;
    sample_spectrum_br0[i] = pow(10., spectrum_br0  +       X * s + Y * c);

    sample_spectrum_g_1[i] = sample_spectrum_g_0[i] + spectrum_g_1 - spectrum_g_0;    // keep size of spectral break
    sample_spectrum_br1[i] = spectrum_br1;                                            // keep second break
    sample_spectrum_g_2[i] = sample_spectrum_g_0[i] + spectrum_g_2 - spectrum_g_0;    // keep size of second break

    if(verbose==1)
    {
    cout<<" sample_spectrum_g_0="<<sample_spectrum_g_0[i]<<" sample_spectrum_br0="<<sample_spectrum_br0[i];
    cout<<" sample_spectrum_g_1="<<sample_spectrum_g_1[i]<<" sample_spectrum_br1="<<sample_spectrum_br1[i];
    cout<<" sample_spectrum_g_2="<<sample_spectrum_g_2[i]<<endl;
    }
   }
  }


  i=0;
  theta=0.;
  dtheta=.1;
  
  // the following are redundant since rand uses same seed in each routine
  srand            (seed); // initialize random number generator
  //  poisson          (seed); // initialize Poisson generator
  //  power_law_sampler(seed); // initialize power-law sampler





  //----------------------------------------------------------------- generate sample

  number_total=0.;
  m_sample_total =0;
  sample_flux_min=1e30;
  sample_flux_max=0.  ;
  sample_d_min   =1e30;
  sample_d_max   =0.  ;

  sample_total_flux=0;
  sample_total_L   =0;

  for (iR=0;iR<nR;iR++)
  for (iz=0;iz<nz;iz++)
    {
      R1  =         iR*dR; R2=R1+dR;
      z1  = -zmax + iz*dz; z2=z1+dz;
      R_=(R1+R2)/2;
      z_=(z1+z2)/2;
      density_=density(R_,z_);    
      density_ *= oversample;                       //AWS20060118                   
      number_delta=density_ *pi*(R2*R2 - R1*R1) *dz;
      number_total+=number_delta;
      m_sample=poisson(number_delta);
      m_sample_total+=m_sample;

      if(verbose==1)
      *txt_stream<<" R_=" <<R_<< " z_="<<z_ <<" number_delta="<<number_delta
           <<" number_total="<<number_total
           << " m_sample="<<m_sample<<" m_sample_total="<<m_sample_total<<endl;

      //    for(ii=0;ii<number_delta-0.5;ii++,i++)

      for(ii=0;ii<m_sample        ;ii++,i++)
        if(i<n_sample)
	{
	  sample_R[i]=R_;
	  sample_z[i]=z_;
 
          theta=2.0*pi*double(rand())/RAND_MAX;                      // angle at GC
          sample_theta[i]=theta;
          sample_d[i]=sqrt(R_*R_ + R0*R0 - 2.0* R_ *R0 *cos(theta) + z_ * z_ ); // distance to source in kpc
          sample_dd  =sqrt(R_*R_ + R0*R0 - 2.0* R_ *R0 *cos(theta)           ); // distance to source in kpc projected onto plane

          sinl= R_ / sample_dd   * sin(theta) ;                              //   sine rule
          cosl=(R0*R0 + sample_dd*sample_dd -  R_*R_)/(2.0 * sample_dd * R0);// cosine rule
	  sample_l[i]=atan2(sinl,cosl)               *rtd        ; // -180/+180 deg

          sample_b[i]=atan (sample_z[i] / sample_dd )*rtd;

	  if(luminosity_function_type==2)L=power_law_sampler(L_min,L_max,alpha_L);
          sample_flux[i]=L/(4.0*pi* pow(sample_d[i]*kpc_to_cm,2));

          if(sample_d[i]   <   sample_d_min) sample_d_min   =sample_d   [i];
          if(sample_d[i]   >   sample_d_max) sample_d_max   =sample_d   [i];
          if(sample_flux[i]<sample_flux_min){sample_flux_min=sample_flux[i]; i_flux_min=i;}
          if(sample_flux[i]>sample_flux_max){sample_flux_max=sample_flux[i]; i_flux_max=i;}


          sample_total_flux+= sample_flux[i];
          sample_total_L   += L;

	  // assign spectra AWS20100817
	  //        double spectrum_norm=integratepowbr(1.0,spectrum_g_0,spectrum_g_1,spectrum_g_2, spectrum_br0,spectrum_br1,E_ref_low,E_ref_high);

	  double spectrum_norm=integratepowbr(1.0,sample_spectrum_g_0[i],sample_spectrum_g_1[i],sample_spectrum_g_2[i],
                                                  sample_spectrum_br0[i],sample_spectrum_br1[i],E_ref_low,E_ref_high);//AWS20100823

          // spectral energies
          for (int iE=0;iE<n_E;iE++)
          {
	    
            
          double spectrum_factor=shapepowbr(1.0,sample_spectrum_g_0[i],sample_spectrum_g_1[i],sample_spectrum_g_2[i],
                                                sample_spectrum_br0[i],sample_spectrum_br1[i],E[iE])
	                                                              //  *pow(E[iE], 2.0 )// NB need *E^2 for compatibility with galprop (l,b) skymaps, applied when required
	                                   / spectrum_norm; //AWS20100823


	  //          *txt_stream<<" E= "<<E[iE]<<" spectrum_norm="<<spectrum_norm<<" spectrum factor with breaks="<<spectrum_factor<<endl;
	  //	  cout<<" br0="<<sample_spectrum_br0[i]<<" "<<" E= "<<E[iE]<<" spectrum_norm="<<spectrum_norm<<" spectrum factor with breaks="<<spectrum_factor<<endl;

           sample_flux_spectra[i][iE]=sample_flux[i]*spectrum_factor;

           if(spectral_model==3)
               sample_flux_spectra[i][iE]
            = flux_cutoff_powerlaw( E[iE], sample_spectrum_g_0[i], sample_spectrum_br0[i], E_ref_low, E_ref_high, sample_flux[i]);//AWS20100830


          }


	 // reference energies AWS20110901

           if(spectral_model==1 || spectral_model==2)
	   {
            double spectrum_factor=shapepowbr(1.0,sample_spectrum_g_0[i],sample_spectrum_g_1[i],sample_spectrum_g_2[i],
                                                sample_spectrum_br0[i],sample_spectrum_br1[i],E_ref_low)
	                                   / spectrum_norm; 

	    sample_flux_spectrum_E_ref_low[i]=sample_flux[i]*spectrum_factor;

                   spectrum_factor=shapepowbr(1.0,sample_spectrum_g_0[i],sample_spectrum_g_1[i],sample_spectrum_g_2[i],
                                                sample_spectrum_br0[i],sample_spectrum_br1[i],E_ref_high)
	                                   / spectrum_norm; 

            sample_flux_spectrum_E_ref_high[i]=sample_flux[i]*spectrum_factor;



	    // >100 MeV, 
	    double E_bin_low_  = 1.e2;
	    double E_bin_high_ = 1.e6;

                   spectrum_factor=integratepowbr(1.0,sample_spectrum_g_0[i],sample_spectrum_g_1[i],sample_spectrum_g_2[i], 
                                                   sample_spectrum_br0[i],sample_spectrum_br1[i],E_bin_low_,E_bin_high_)
                                           /spectrum_norm; //AWS20120615


	     sample_flux100[i]=sample_flux[i]*spectrum_factor;                                           //AWS20120615

             //cout<<"sample_flux100="<<sample_flux100[i]      ;

	    // >1000 MeV, 
	           E_bin_low_  = 1.e3;
	           E_bin_high_ = 1.e6;

                   spectrum_factor=integratepowbr(1.0,sample_spectrum_g_0[i],sample_spectrum_g_1[i],sample_spectrum_g_2[i], 
                                                   sample_spectrum_br0[i],sample_spectrum_br1[i],E_bin_low_,E_bin_high_)
                                           /spectrum_norm; //AWS20120615


	     sample_flux1000[i]=sample_flux[i]*spectrum_factor;                                           //AWS20120615

             //cout<<" sample_flux1000="<<sample_flux1000[i];

	    // >10000 MeV, 
	           E_bin_low_  = 1.e4;
	           E_bin_high_ = 1.e6;

                   spectrum_factor=integratepowbr(1.0,sample_spectrum_g_0[i],sample_spectrum_g_1[i],sample_spectrum_g_2[i], 
                                                   sample_spectrum_br0[i],sample_spectrum_br1[i],E_bin_low_,E_bin_high_)
                                           /spectrum_norm; //AWS20120615


	     sample_flux10000[i]=sample_flux[i]*spectrum_factor;                                           //AWS20120615

             //cout<<" sample_flux10000="<<sample_flux10000[i]<<endl;




           }


           if(spectral_model==3)
	   {
            sample_flux_spectrum_E_ref_low[i]
             = flux_cutoff_powerlaw
                               ( E_ref_low,  sample_spectrum_g_0[i], sample_spectrum_br0[i], E_ref_low, E_ref_high, sample_flux[i]);

            sample_flux_spectrum_E_ref_high[i]
             = flux_cutoff_powerlaw
                               ( E_ref_high, sample_spectrum_g_0[i], sample_spectrum_br0[i], E_ref_low, E_ref_high, sample_flux[i]);
           }


         //  binned energy ranges

         for (int iE=0;iE<n_E_bin;iE++)
         {
    

         double spectrum_factor=integratepowbr(1.0,sample_spectrum_g_0[i],sample_spectrum_g_1[i],sample_spectrum_g_2[i], sample_spectrum_br0[i],sample_spectrum_br1[i],E_bin_low[iE],E_bin_high[iE])	   /spectrum_norm; //AWS20100823


          sample_binned_flux_spectra[i][iE]=sample_flux[i]*spectrum_factor;

           if(spectral_model==3)
	         sample_binned_flux_spectra[i][iE]
             = flux_cutoff_powerlaw( E_bin_low[iE] ,E_bin_high[iE], sample_spectrum_g_0[i], sample_spectrum_br0[i], E_ref_low, E_ref_high, sample_flux[i]);//AWS20100830


	  //      *txt_stream<<" E= "<<E_bin_low[iE]<<" - "<<E_bin_high[iE]<<" source binned spectrum factor="<<spectrum_factor<<endl;
          }

	      if(verbose==1)
	{
	  *txt_stream<<"i="<<i<<" R=" <<sample_R[i]<< " z="<<sample_z[i]<<" theta="<<sample_theta[i]
		      <<" L="<<L
			     <<" distance="<<sample_d[i]   <<" l="<<sample_l[i]    <<" b="<<sample_b[i]      <<" flux="<<sample_flux[i];
          *txt_stream  <<" flux spectrum= " ;              for (int iE=0;iE<n_E    ;iE++){*txt_stream<< sample_flux_spectra      [i][iE]<<" ";}
          *txt_stream  <<" binned flux spectrum= " ;       for (int iE=0;iE<n_E_bin;iE++){*txt_stream<< sample_binned_flux_spectra[i][iE]<<" ";}

          *txt_stream      <<endl;
	}
 



	}// if
    }// for


   *txt_stream<<"oversampling factor for sources     = "<<oversample<<" (following values include this factor ! ) " <<endl;     //AWS2010824
   *txt_stream<<"total number of sources from density= "<<number_total<<endl;
   *txt_stream<<"total number of sources simulated   = "<<m_sample_total<<endl;
   *txt_stream<<"total luminosity of sources = "<<sample_total_L<<endl;
   *txt_stream<<" sample distance min= "<<sample_d_min <<" max= "<<sample_d_max<<endl;                                           //AWS2010824
   *txt_stream<<" flux min= "<<sample_flux_min<<" flux max= "<<sample_flux_max<<endl;
   *txt_stream<<" l of min flux = "<< sample_l[i_flux_min]<<"    l of max flux = "<< sample_l[i_flux_max] <<endl;                //AWS2010824
   *txt_stream<<" b of min flux = "<< sample_b[i_flux_min]<<"    b of max flux = "<< sample_b[i_flux_max] <<endl;
   *txt_stream<<" sample  global         total flux             = "<<sample_total_flux<<endl;

   if(verbose==2){int options=2; print(options);}
   if(verbose==3){int options=3; print(options);}

   *txt_stream<<"SourcePopulation:: gen_sample  complete"<<endl<<endl;

   return 0;

}

/////////////////////////////////////////////////////////

int SourcePopulation::analyse_sample()
{
  int i,ii;
  double density_,number_total,number_delta ;
  double solid_angle,solid_angle_long,solid_angle_lat,solid_angle_long_lat;
  double l,b,b1,b2;
  int il,ib;
  double spectrum_factor;
  double spectrum_norm;
  double long_min1a,long_min2a,long_max1a,long_max2a;

  pointing apointing; // healpix coordinates of a pixel //AWS20090109
  int ipix;           // healpix pixel index            //AWS20090109

  int skymap_method = 2; // 1=scale total intensity (original method), 2 = from source spectra (new method) //AWS20100818

  *txt_stream<< " >> SourcePopulation::analyse_sample"<<endl;
  *txt_stream<< "title :"<< title        <<endl;

   //----------------------------------------------------------  generate differential logN-logS
   lnS_min =-20.0001; // to make printing better
   lnS_max =-00;
   dlnS    =.20;

  n_dlnN_dlnS=  int((lnS_max-lnS_min)/dlnS);
    dlnN_dlnS        =new double [ n_dlnN_dlnS];
    dlnN_dlnS_deg    =new double [ n_dlnN_dlnS];
    dlnN_dlnS_int    =new double [ n_dlnN_dlnS];
    dlnN_dlnS_int_deg=new double [ n_dlnN_dlnS];
    FS               =new double [ n_dlnN_dlnS];
    FS_int           =new double [ n_dlnN_dlnS];

   dlnN_dlnS_soplimit        =new double [ n_dlnN_dlnS]; //AWS20170123
   dlnN_dlnS_sublimit        =new double [ n_dlnN_dlnS]; //AWS20170123
   dlnN_dlnS_int_soplimit    =new double [ n_dlnN_dlnS]; //AWS20170123
   dlnN_dlnS_int_sublimit    =new double [ n_dlnN_dlnS]; //AWS20170123
   FS_soplimit               =new double [ n_dlnN_dlnS]; //AWS20170123
   FS_sublimit               =new double [ n_dlnN_dlnS]; //AWS20170123
   FS_int_soplimit           =new double [ n_dlnN_dlnS]; //AWS20170123
   FS_int_sublimit           =new double [ n_dlnN_dlnS]; //AWS20170123

  sample_l_min=new double[ n_dlnN_dlnS];
  sample_l_max=new double[ n_dlnN_dlnS];
  sample_l_av =new double[ n_dlnN_dlnS];
  sample_b_min=new double[ n_dlnN_dlnS];
  sample_b_max=new double[ n_dlnN_dlnS];
  sample_b_av =new double[ n_dlnN_dlnS];

  longprof_nlong             =   int(360/longprof_dlong);
  longprof_intensity         =new double[longprof_nlong];
  longprof_intensity_sublimit=new double[longprof_nlong];
  longprof_sourcecnt         =new double[longprof_nlong];
  longprof_sourcecnt_sublimit=new double[longprof_nlong];

   latprof_nlat               =   int(180/latprof_dlat );
   latprof_intensity         =new double[ latprof_nlat ];
   latprof_intensity_sublimit=new double[ latprof_nlat ];
   latprof_sourcecnt         =new double[ latprof_nlat ];
   latprof_sourcecnt_sublimit=new double[ latprof_nlat ];

  for (ii=0;ii< n_dlnN_dlnS;ii++)
  {
   dlnN_dlnS      [ii]=0;
   dlnN_dlnS_deg  [ii]=0;

   FS             [ii]=0;

   dlnN_dlnS_soplimit[ii]=0; //AWS20170123
   dlnN_dlnS_sublimit[ii]=0; //AWS20170123

   sample_l_min[ii]= 1e30;
   sample_b_min[ii]= 1e30;
   sample_l_max[ii]=-1e30;
   sample_b_max[ii]=-1e30;

   sample_l_av[ii]=0;
   sample_b_av[ii]=0;
  }

  for (ii=0;ii<longprof_nlong;ii++)
  {
   longprof_intensity         [ii]=0.;
   longprof_intensity_sublimit[ii]=0.;
   longprof_sourcecnt         [ii]=0.;
   longprof_sourcecnt         [ii]=0.;
   }


  for (ii=0;ii< latprof_nlat ;ii++)
  {
    latprof_intensity         [ii]=0.;
    latprof_intensity_sublimit[ii]=0.;
    latprof_sourcecnt         [ii]=0.;
    latprof_sourcecnt         [ii]=0.;
   }


  sample_selected_total_flux         =0.;
  sample_selected_total_flux_sublimit=0.;
  sample_selected_total_flux_soplimit=0.;

  sample_selected_sourcecnt          =0.;
  sample_selected_sourcecnt_sublimit =0.;
  sample_selected_sourcecnt_soplimit =0.;

  sample_selected_d_min              = 1e30;
  sample_selected_d_min_sublimit     = 1e30;
  sample_selected_d_min_soplimit     = 1e30;

  sample_selected_d_max              =-1e30;
  sample_selected_d_max_sublimit     =-1e30;
  sample_selected_d_max_soplimit     =-1e30;


  sample_selected_d_av               = 0.;         
  sample_selected_d_av_sublimit      = 0.;
  sample_selected_d_av_soplimit      = 0.;

  //----------------------------------------------------------------------------------
  // set limits in range -180<l<+180

  // but this is problematic, instead convert to 0-360 and compare with original ranges  AWS20120522

  long_min1a=long_min1;
  long_min2a=long_min2;
  long_max1a=long_max1;
  long_max2a=long_max2;

  if(long_min1a  >= +180.){long_min1a -=360.; long_max1a -=360.;}
  if(long_min2a  >= +180.){long_min2a -=360.; long_max2a -=360.;}

  // more robust in marginal cases like 0.0-180.0, 180.0-360.0 but not correct    AWS20120521
  //  if(long_max1a  >= +180.){long_min1a -=360.; long_max1a -=360.;}           //AWS20120521
  //  if(long_max2a  >= +180.){long_min2a -=360.; long_max2a -=360.;}           //AWS20120521

   for (i=0;i<n_sample;i++)
     {

       ii=int( (log10(sample_flux[i])-lnS_min)/dlnS  );

       if(verbose==3)
        *txt_stream<<"before selection :sample i="<<i<<" l="<<sample_l[i]<<" sample_flux[i]="<<sample_flux[i]
         << " ii="<<ii
	 <<" long_min1 ="<<long_min1 <<" long_max1 ="<<long_max1 
	 <<" long_min1a="<<long_min1a<<" long_max1a="<<long_max1a
	 <<" long_min2 ="<<long_min2 <<" long_max2 ="<<long_max2 
	 <<" long_min2a="<<long_min2a<<" long_max2a="<<long_max2a
         <<endl;

       if(ii>=0 && ii<n_dlnN_dlnS)
       {
         //          -180 <     sample_l     <180
         //             0 <long_min,long_max <360


	 /*
	 if(  (    ( sample_l[i]  >= long_min1a     && sample_l[i]  <= long_max1a  )  
		 ||( sample_l[i]  >= long_min2a     && sample_l[i]  <= long_max2a  ) )
	 */

	 // put l in range 0-360
         double sample_l_0360= sample_l[i]; if ( sample_l_0360 < 0.0)  sample_l_0360+= 360.0; //AWS20120522 

	 if(  (    ( sample_l_0360 >= long_min1  && sample_l_0360<= long_max1   )             //AWS20120522 
		 ||( sample_l_0360 >= long_min2  && sample_l_0360<= long_max2   ) )           //AWS20120522 
               &&
              (    ( sample_b[i]   >=  lat_min1  && sample_b[i]  <=  lat_max1  )  
		 ||( sample_b[i   ]>=  lat_min2  && sample_b[i]  <=  lat_max2  ) )
           )

       {

	 //AWS20170116  replace flux_detection by value in sensitivity map
         {
          double theta=(90-sample_b[i])*dtr;
          double phi  =    sample_l[i] *dtr;
          pointing pointing_=pointing(theta,phi);
          cout<<"source statistics sample i="<<i<<" Fermi sensitivity  l="<<sample_l[i] <<" b="<<sample_b[i]
              <<" theta="<<theta<<" phi="<<phi;
	  //          cout<<" interpolated value="<<healpix_skymap_Fermi_sensitivity.interpolated_value(pointing_)<<endl;
          flux_detection_limit_sensitivity_file = healpix_skymap_Fermi_sensitivity.interpolated_value(pointing_); //AWS20200713 was flux_detection_limit
          cout<<" flux detection limit_sensitivity_file="<<flux_detection_limit_sensitivity_file<<" sample flux="<<sample_flux[i]<<endl;
         }

        dlnN_dlnS  [ii]++;
        FS         [ii]+= sample_flux[i];

        if(sample_flux[i]<flux_detection_limit_sensitivity_file)    //AWS20200713
	{                                          //AWS20170123
         dlnN_dlnS_sublimit[ii]++;                 //AWS20170123
                FS_sublimit[ii]+= sample_flux[i];  //AWS20170123
	}                                          //AWS20170123

        if(sample_flux[i]>flux_detection_limit_sensitivity_file)    //AWS20200713
	{                                          //AWS20170123
         dlnN_dlnS_soplimit[ii]++;                 //AWS20170123
                FS_soplimit[ii]+= sample_flux[i]; //AWS20170123
	}                                          //AWS20170123


	sample_l_av[ii]+=abs(sample_l[i]);
	sample_b_av[ii]+=abs(sample_b[i]);
	
        if (abs(sample_l[i])< sample_l_min[ii]) sample_l_min[ii]= abs(sample_l[i]);
        if (abs(sample_b[i])< sample_b_min[ii]) sample_b_min[ii]= abs(sample_b[i]);
        if (abs(sample_l[i])> sample_l_max[ii]) sample_l_max[ii]= abs(sample_l[i]);
        if (abs(sample_b[i])> sample_b_max[ii]) sample_b_max[ii]= abs(sample_b[i]);



         sample_selected_sourcecnt++;
         sample_selected_total_flux         +=sample_flux[i];
         sample_selected_d_av               +=sample_d   [i];
         if (sample_d[i] < sample_selected_d_min) sample_selected_d_min= sample_d[i];
         if (sample_d[i] > sample_selected_d_max) sample_selected_d_max= sample_d[i];



	 if(sample_flux[i]<flux_detection_limit_sensitivity_file) //AWS20200713
	{
         sample_selected_sourcecnt_sublimit++;
         sample_selected_total_flux_sublimit+=sample_flux[i];
         sample_selected_d_av_sublimit      +=sample_d   [i];

         if (sample_d[i] < sample_selected_d_min_sublimit) sample_selected_d_min_sublimit= sample_d[i];
         if (sample_d[i] > sample_selected_d_max_sublimit) sample_selected_d_max_sublimit= sample_d[i];
	}

	 if(sample_flux[i]>flux_detection_limit_sensitivity_file)  //AWS20200713
	{
         sample_selected_sourcecnt_soplimit++;
         sample_selected_total_flux_soplimit+=sample_flux[i];
         sample_selected_d_av_soplimit      +=sample_d   [i];

         if (sample_d[i] < sample_selected_d_min_soplimit) sample_selected_d_min_soplimit= sample_d[i];
         if (sample_d[i] > sample_selected_d_max_soplimit) sample_selected_d_max_soplimit= sample_d[i];
	}

      if(verbose==3)
	    	  *txt_stream<<"selected sample i="<<i<<" R=" <<sample_R[i]<< " z="<<sample_z[i]<<" theta="<<sample_theta[i]
		      <<" L="<<L
                      <<" distance="<<sample_d[i]   <<" l="<<sample_l[i]    <<" b="<<sample_b[i]      <<" flux="<<sample_flux[i]                  <<endl;

        } // if      

	// longitude profile
	 //      if(sample_b[i]>=longprof_b_min && sample_b[i]<=longprof_b_max)

        if         (( sample_b[i]>=  lat_min1     && sample_b[i]<=  lat_max1  )  
		 ||(  sample_b[i]>=  lat_min2     && sample_b[i]<=  lat_max2  ) )
	  {
	    ii=int((sample_l[i]+180.)/longprof_dlong);
	    if(ii >= 0 && ii < longprof_nlong) 
	    {
             longprof_intensity         [ii] += sample_flux[i];

             if(sample_flux[i]<flux_detection_limit_sensitivity_file) //AWS20200713
             longprof_intensity_sublimit[ii] += sample_flux[i];

             longprof_sourcecnt         [ii] ++               ;

             if(sample_flux[i]<flux_detection_limit_sensitivity_file) //AWS20200713
             longprof_sourcecnt_sublimit[ii] ++               ;

	    }//if ii
	  }//if sample_b
	
      

	//  latitude profile
	//  if(sample_l[i]>= latprof_l_min && sample_l[i]<= latprof_l_max)

	 if   (    ( sample_l[i]  >= long_min1a     && sample_l[i]  <= long_max1a  )  
		 ||( sample_l[i]  >= long_min2a     && sample_l[i]  <= long_max2a  ) )
	  {
	    ii=int((sample_b[i]+ 90.)/ latprof_dlat );
	    if(ii >= 0 && ii <  latprof_nlat ) 
	    {
              latprof_intensity         [ii] += sample_flux[i];

	      if(sample_flux[i]<flux_detection_limit_sensitivity_file) //AWS20200713
              latprof_intensity_sublimit[ii] += sample_flux[i];

              latprof_sourcecnt         [ii] ++               ;

	      if(sample_flux[i]<flux_detection_limit_sensitivity_file) //AWS20200713
              latprof_sourcecnt_sublimit[ii] ++               ;

	    }//if ii
	  }//if sample_l
	
       }//if


     }// i 

   if                                        (sample_selected_sourcecnt >0)
   sample_selected_d_av             /=        sample_selected_sourcecnt         ;
   if                                        (sample_selected_sourcecnt_sublimit >0)
   sample_selected_d_av_sublimit    /=        sample_selected_sourcecnt_sublimit;
   if                                        (sample_selected_sourcecnt_soplimit >0)
   sample_selected_d_av_soplimit    /=        sample_selected_sourcecnt_soplimit;



   //--------------------------------------------------------------------------------

   dlnN_dlnS_total_N=0.;
   dlnN_dlnS_total_S=0.;


  for (ii=0;ii< n_dlnN_dlnS;ii++)
  { 
    if( dlnN_dlnS  [ii] > 0)
    {
     sample_l_av[ii] /= dlnN_dlnS  [ii];
     sample_b_av[ii] /= dlnN_dlnS  [ii];
    }

    solid_angle=  (long_max1-long_min1)*(sin(lat_max1*dtr)-sin(lat_min1*dtr))*rtd
                + (long_max2-long_min2)*(sin(lat_max1*dtr)-sin(lat_min1*dtr))*rtd
                + (long_max1-long_min1)*(sin(lat_max2*dtr)-sin(lat_min2*dtr))*rtd
                + (long_max2-long_min2)*(sin(lat_max2*dtr)-sin(lat_min2*dtr))*rtd;

    dlnN_dlnS_deg[ii]= dlnN_dlnS[ii] / solid_angle;

    dlnN_dlnS_total_N += dlnN_dlnS  [ii];
    dlnN_dlnS_total_S += dlnN_dlnS  [ii] * (pow(10,lnS_min+(ii)*dlnS) + pow(10,lnS_min+(ii+1)*dlnS))/2.0 ;

  }

  // integral N(S)
  for (ii=0;ii< n_dlnN_dlnS;ii++)
  {
    dlnN_dlnS_int    [ii]=0;
    dlnN_dlnS_int_deg[ii]=0;
    FS_int           [ii]=0;

    dlnN_dlnS_int_sublimit    [ii]=0;                         //AWS20170123
           FS_int_sublimit    [ii]=0;                         //AWS20170123

    dlnN_dlnS_int_soplimit    [ii]=0;                         //AWS20170123
           FS_int_soplimit    [ii]=0;                         //AWS20170123

    for(i=ii;i< n_dlnN_dlnS;i++)
    {    
     dlnN_dlnS_int    [ii]+= dlnN_dlnS      [i];
     dlnN_dlnS_int_deg[ii]+= dlnN_dlnS_deg  [i];
     FS_int           [ii]+= FS             [i];

     dlnN_dlnS_int_sublimit    [ii]+= dlnN_dlnS_sublimit [i]; //AWS20170123
            FS_int_sublimit    [ii]+=    FS_int_sublimit [i]; //AWS20170123

     dlnN_dlnS_int_soplimit    [ii]+= dlnN_dlnS_soplimit [i]; //AWS20170123
            FS_int_soplimit    [ii]+=    FS_int_soplimit [i]; //AWS20170123

    }
  }


  //  solid_angle_long=longprof_dlong*(longprof_b_max-longprof_b_min)/(rtd*rtd);

  solid_angle_long=longprof_dlong*dtr * (sin(lat_max1*dtr)-sin(lat_min1*dtr))
                  +longprof_dlong*dtr * (sin(lat_max2*dtr)-sin(lat_min2*dtr));
  *txt_stream<<"longitude profile solid_angle     ="<<solid_angle_long<<" sr"<<endl;

  //  solid_angle_lat = latprof_dlat *( latprof_l_max- latprof_l_min)/(rtd*rtd);
  // *txt_stream<<" latitude profile solid_angle     ="<<solid_angle_lat <<" sr"<<endl;


  for (ii=0;ii<longprof_nlong;ii++)
  {
   longprof_intensity         [ii]/=solid_angle_long;
   longprof_intensity_sublimit[ii]/=solid_angle_long;
  }

  for (ii=0;ii< latprof_nlat ;ii++)
  {
   b=-90.+ (ii+.5)* latprof_dlat;
   // solid_angle_lat = latprof_dlat *( latprof_l_max- latprof_l_min)/(rtd*rtd)*cos(b*dtr);
    b1=-90.+ ii* latprof_dlat;
    b2= b1 +     latprof_dlat;
    solid_angle_lat = (long_max1-long_min1) *dtr  *(sin(b2*dtr) - sin(b1*dtr))
                    + (long_max2-long_min2) *dtr  *(sin(b2*dtr) - sin(b1*dtr));
    latprof_intensity         [ii]/=solid_angle_lat ;
    latprof_intensity_sublimit[ii]/=solid_angle_lat ;
  }



  //-------------------------------------------------------------------- generate skymaps

  // (l,b) maps
   skymap_sourcecnt         .init(longprof_nlong, latprof_nlat,1);
   skymap_sourcecnt_sublimit.init(longprof_nlong, latprof_nlat,1);
   skymap_sourcecnt_soplimit.init(longprof_nlong, latprof_nlat,1);

   skymap_intensity_spectrum         .init(longprof_nlong, latprof_nlat,n_E);
   skymap_intensity_spectrum_sublimit.init(longprof_nlong, latprof_nlat,n_E);
   skymap_intensity_spectrum_soplimit.init(longprof_nlong, latprof_nlat,n_E);

   skymap_intensity         .init(longprof_nlong, latprof_nlat,n_E_bin);
   skymap_intensity_sublimit.init(longprof_nlong, latprof_nlat,n_E_bin);
   skymap_intensity_soplimit.init(longprof_nlong, latprof_nlat,n_E_bin);

   // healpix maps
   healpix_skymap_sourcecnt                    = Skymap<double> (healpix_skymap_order, 1);       //AWS20090107
   healpix_skymap_sourcecnt_sublimit           = Skymap<double> (healpix_skymap_order, 1);       //AWS20090107
   healpix_skymap_sourcecnt_soplimit           = Skymap<double> (healpix_skymap_order, 1);       //AWS20090107


   healpix_skymap_intensity_spectrum           = Skymap<double> (healpix_skymap_order, n_E);     //AWS20090107
   healpix_skymap_intensity_spectrum_sublimit  = Skymap<double> (healpix_skymap_order, n_E);     //AWS20090107
   healpix_skymap_intensity_spectrum_soplimit  = Skymap<double> (healpix_skymap_order, n_E);     //AWS20090107

   healpix_skymap_intensity                    = Skymap<double> (healpix_skymap_order, n_E_bin); //AWS20090107
   healpix_skymap_intensity_sublimit           = Skymap<double> (healpix_skymap_order, n_E_bin); //AWS20090107
   healpix_skymap_intensity_soplimit           = Skymap<double> (healpix_skymap_order, n_E_bin); //AWS20090107

   healpix_skymap_sourcecnt           = 0.;
   healpix_skymap_sourcecnt_sublimit  = 0.;
   healpix_skymap_sourcecnt_soplimit  = 0.;

   healpix_skymap_intensity           = 0.;
   healpix_skymap_intensity_sublimit  = 0.;
   healpix_skymap_intensity_soplimit  = 0.;

   healpix_skymap_intensity_spectrum           = 0.; //AWS20100818
   healpix_skymap_intensity_spectrum_sublimit  = 0.; //AWS20100818
   healpix_skymap_intensity_spectrum_soplimit  = 0.; //AWS20100818


   //   cout<<"testing ang2pix:";
   //   apointing = pointing(0.5,0.1);
   //   cout<<healpix_skymap_sourcecnt.ang2pix(apointing)<<endl; // since Skymap is derived from Healpix_Base


   //---------------------------------------------------------------------------------------------------------


   cout<<" generating skymaps of sources"<<endl;

   for (i=0;i<n_sample;i++)
   {

	 //AWS20170116  replace flux_detection by value in sensitivity map
         {
          double theta=(90-sample_b[i])*dtr;
          double phi  =    sample_l[i] *dtr;
          pointing pointing_=pointing(theta,phi);
          cout<<"skymap generation source sample i="<<i<<" Fermi sensitivity  l="<<sample_l[i] <<" b="<<sample_b[i]
              <<" theta="<<theta<<" phi="<<phi;
	  //          cout<<" interpolated value="<<healpix_skymap_Fermi_sensitivity.interpolated_value(pointing_)<<endl;
          flux_detection_limit_sensitivity_file = healpix_skymap_Fermi_sensitivity.interpolated_value(pointing_); //AWS20200713
          cout<<" flux detection limit_sensitivity_file="<<flux_detection_limit_sensitivity_file<<" sample flux="<<sample_flux[i]<<endl;
         }

       ii=int( (log10(sample_flux[i])-lnS_min)/dlnS  );
       //   *txt_stream<<"sample_flux[i]="<<sample_flux[i]<< " ii="<<ii<<endl;
       if(ii>=0 && ii<n_dlnN_dlnS)
       {
	 // for skymap use -180/+180 degrees as in galplot internally
            if(sample_l[i]>=0.0 )
	    il=int((sample_l[i]     )/longprof_dlong);
            if(sample_l[i]< 0.0 )
	    il=int((sample_l[i]+360.)/longprof_dlong);

	    ib=int((sample_b[i]+ 90.)/ latprof_dlat );

	    if(  il >= 0 && il < longprof_nlong  
	      && ib >= 0 && ib <  latprof_nlat )
	    {
             skymap_intensity.d2[il][ib].s[0]       += sample_flux[i];

             if(sample_flux[i]<flux_detection_limit_sensitivity_file) //AWS20200713
             skymap_intensity_sublimit.d2[il][ib].s[0]       += sample_flux[i];

             if(sample_flux[i]>flux_detection_limit_sensitivity_file) //AWS20200713
             skymap_intensity_soplimit.d2[il][ib].s[0]       += sample_flux[i];

             skymap_sourcecnt.d2[il][ib].s[0]       ++               ;

             if(sample_flux[i]<flux_detection_limit_sensitivity_file) //AWS20200713
             skymap_sourcecnt_sublimit.d2[il][ib].s[0]       ++               ;

             if(sample_flux[i]>flux_detection_limit_sensitivity_file) //AWS20200713
             skymap_sourcecnt_soplimit.d2[il][ib].s[0]       ++               ;


	    }//if il ib


	    // healpix                                                                                          //AWS20090108

	    apointing=pointing( (90.-sample_b[i])*dtr , sample_l[i]*dtr   ); // input is theta, phi => 90-b, l
            ipix=healpix_skymap_sourcecnt.ang2pix(apointing);
	    //           cout<<"SourcePopulation healpix maps: l="<<sample_l[i]<<" b="<<sample_b[i]<<" ipix="<<ipix<<endl;

                                                      healpix_skymap_sourcecnt         [ipix][0]++;
            if(sample_flux[i]<flux_detection_limit_sensitivity_file)   healpix_skymap_sourcecnt_sublimit[ipix][0]++;  
            if(sample_flux[i]>flux_detection_limit_sensitivity_file)   healpix_skymap_sourcecnt_soplimit[ipix][0]++;


           if(skymap_method==1)//AWS20100820
	   {
	                                              healpix_skymap_intensity         [ipix][0]+= sample_flux[i];
            if(sample_flux[i]<flux_detection_limit_sensitivity_file)   healpix_skymap_intensity_sublimit[ipix][0]+= sample_flux[i];
            if(sample_flux[i]>flux_detection_limit_sensitivity_file)   healpix_skymap_intensity_soplimit[ipix][0]+= sample_flux[i];
	   }

	    // using source spectra directly 
            if(skymap_method==2)//AWS20100818
	    {
	      for(int iE=0;iE<n_E;iE++)
	      {
	                                                 healpix_skymap_intensity_spectrum         [ipix][iE]+= sample_flux_spectra[i][iE];
               if(sample_flux[i]<flux_detection_limit_sensitivity_file)   healpix_skymap_intensity_spectrum_sublimit[ipix][iE]+= sample_flux_spectra[i][iE];
               if(sample_flux[i]>flux_detection_limit_sensitivity_file)   healpix_skymap_intensity_spectrum_soplimit[ipix][iE]+= sample_flux_spectra[i][iE];
	      }

	      for(int iE=0;iE<n_E_bin;iE++)
	      {
	                                                 healpix_skymap_intensity                  [ipix][iE]+= sample_binned_flux_spectra[i][iE];
               if(sample_flux[i]<flux_detection_limit_sensitivity_file)   healpix_skymap_intensity_sublimit         [ipix][iE]+= sample_binned_flux_spectra[i][iE];
               if(sample_flux[i]>flux_detection_limit_sensitivity_file)   healpix_skymap_intensity_soplimit         [ipix][iE]+= sample_binned_flux_spectra[i][iE];
	      }

	    }


       }//if ii

       //     cout<<" processed source #"<<i<<endl;


     }// for i


   // check sums of l,b and healpix maps                                                                                            //AWS20090901
   *txt_stream<<" completed adding sources to skymaps"<<endl;                                                                       //AWS20090901
   double skymap_intensity_total         =0.;    for(il=0;il< longprof_nlong;il++)  for(ib=0;ib<  latprof_nlat ;ib++)skymap_intensity_total         +=    skymap_intensity            .d2[il][ib].s[0];
   double skymap_intensity_sublimit_total=0.;    for(il=0;il< longprof_nlong;il++)  for(ib=0;ib<  latprof_nlat ;ib++)skymap_intensity_sublimit_total+=    skymap_intensity_sublimit   .d2[il][ib].s[0];
   double skymap_intensity_soplimit_total=0.;    for(il=0;il< longprof_nlong;il++)  for(ib=0;ib<  latprof_nlat ;ib++)skymap_intensity_soplimit_total+=    skymap_intensity_soplimit   .d2[il][ib].s[0];


   *txt_stream<<"sum  before applying solid angle         skymap_intensity           "<<        skymap_intensity_total         <<endl;//AWS20090901
   *txt_stream<<"sum  before applying solid angle         skymap_intensity_sublimit  "<<        skymap_intensity_sublimit_total<<endl;//AWS20090901
   *txt_stream<<"sum  before applying solid angle         skymap_intensity_soplimit  "<<        skymap_intensity_soplimit_total<<endl;//AWS20090901

   *txt_stream<<"sum  before applying solid angle healpix_skymap_intensity           "<<healpix_skymap_intensity         .sum()<<endl;//AWS20090901
   *txt_stream<<"sum  before applying solid angle healpix_skymap_intensity_sublimit  "<<healpix_skymap_intensity_sublimit.sum()<<endl;//AWS20090901
   *txt_stream<<"sum  before applying solid angle healpix_skymap_intensity_soplimit  "<<healpix_skymap_intensity_soplimit.sum()<<endl;//AWS20090901




   for(il=0;il< longprof_nlong;il++)
   for(ib=0;ib<  latprof_nlat ;ib++)
   {
    b=-90.+ (ib+.5)* latprof_dlat;
    b1=-90.+ ib* latprof_dlat;
    b2= b1 +     latprof_dlat;
    //solid_angle_long_lat = longprof_dlong * latprof_dlat /(rtd*rtd) * cos(b*dtr);*txt_stream<< solid_angle_long_lat;
    solid_angle_long_lat = longprof_dlong / rtd *(sin(b2*dtr) - sin(b1*dtr));
    //*txt_stream<<" "<< solid_angle_long_lat<<endl;
    skymap_intensity         .d2[il][ib].s[0] /=solid_angle_long_lat;
    skymap_intensity_sublimit.d2[il][ib].s[0] /=solid_angle_long_lat;
    skymap_intensity_soplimit.d2[il][ib].s[0] /=solid_angle_long_lat;
   }

   healpix_skymap_intensity          /= healpix_skymap_intensity.solidAngle();  //AWS20090108
   healpix_skymap_intensity_sublimit /= healpix_skymap_intensity.solidAngle();  //AWS20090108
   healpix_skymap_intensity_soplimit /= healpix_skymap_intensity.solidAngle();  //AWS20090108

   healpix_skymap_intensity_spectrum          /= healpix_skymap_intensity.solidAngle();  //AWS20100818
   healpix_skymap_intensity_spectrum_sublimit /= healpix_skymap_intensity.solidAngle();  //AWS20100818
   healpix_skymap_intensity_spectrum_soplimit /= healpix_skymap_intensity.solidAngle();  //AWS20100818



   //   cout<<" completed multiplying  skymaps by solid angle"<<endl;

   if(skymap_method==1) //AWS20100818
     {

   spectrum_norm=integratepowbr(1.0,spectrum_g_0,spectrum_g_1,spectrum_g_2, spectrum_br0,spectrum_br1,E_ref_low,E_ref_high);

   // === skymaps at spectral energies

   for (ii=0;ii<n_E;ii++)
   {
     //   spectrum_factor=pow(E[ii]/E[0],-spectrum_g + 2 ) // *E^2 for compatibility with galprop skymaps
     //     spectrum_factor=pow(E[ii],-spectrum_g + 2 ) // *E^2 for compatibility with galprop skymaps
     //               * (spectrum_g -1)/(pow(E_ref_low    ,-spectrum_g+1)-pow(E_ref_high    ,-spectrum_g+1));

     //    *txt_stream<<" E= "<<E[ii]<<" spectrum factor without breaks="<<spectrum_factor<<endl;

     

      spectrum_factor=shapepowbr(1.0,spectrum_g_0,spectrum_g_1,spectrum_g_2, spectrum_br0,spectrum_br1,E[ii])
                        *pow(E[ii], 2.0 )// *E^2 for compatibility with galprop (l,b) skymaps
               	        / spectrum_norm;

    *txt_stream<<" E= "<<E[ii]<<" spectrum factor with breaks="<<spectrum_factor<<endl;

    for(il=0;il< longprof_nlong;il++)
    for(ib=0;ib<  latprof_nlat ;ib++)
    {

       skymap_intensity_spectrum          .d2[il][ib].s[ii] = skymap_intensity          .d2[il][ib].s[0] * spectrum_factor;
       skymap_intensity_spectrum_sublimit .d2[il][ib].s[ii] = skymap_intensity_sublimit .d2[il][ib].s[0] * spectrum_factor;
       skymap_intensity_spectrum_soplimit .d2[il][ib].s[ii] = skymap_intensity_soplimit .d2[il][ib].s[0] * spectrum_factor;


    } //for il ib


    //  cout<<" completed multiplying  l,b  spectral skymaps by spectral factor"<<endl;

    // healpix  //AWS20090108
    // healpix without *E^2, compatible with galprop healpix skymaps
      spectrum_factor=shapepowbr(1.0,spectrum_g_0,spectrum_g_1,spectrum_g_2, spectrum_br0,spectrum_br1,E[ii])
               	        / spectrum_norm;

    for (ipix=0;ipix<healpix_skymap_intensity_spectrum.Npix();ipix++)
    {
     healpix_skymap_intensity_spectrum         [ipix][ii] = healpix_skymap_intensity         [ipix][0] * spectrum_factor;
     healpix_skymap_intensity_spectrum_sublimit[ipix][ii] = healpix_skymap_intensity_sublimit[ipix][0] * spectrum_factor;
     healpix_skymap_intensity_spectrum_soplimit[ipix][ii] = healpix_skymap_intensity_soplimit[ipix][0] * spectrum_factor;
    }


   } //for ii

  }//skymap_method==1

   //  cout<<" completed multiplying  healpix spectral skymaps by spectral factor"<<endl;



   // all-sky average intensity
   valarray<double> healpix_skymap_intensity_spectrum_sublimit_av(n_E)   ;
   valarray<double> healpix_skymap_intensity_spectrum_soplimit_av(n_E)   ;
                    healpix_skymap_intensity_spectrum_sublimit_av=0.;
                    healpix_skymap_intensity_spectrum_soplimit_av=0.;

   for(ii=0;ii<n_E;ii++) 
   {
    for (ipix=0;ipix<healpix_skymap_intensity_spectrum.Npix();ipix++)
     {
      healpix_skymap_intensity_spectrum_sublimit_av[ii]+=healpix_skymap_intensity_spectrum_sublimit[ipix][ii];
      healpix_skymap_intensity_spectrum_soplimit_av[ii]+=healpix_skymap_intensity_spectrum_soplimit[ipix][ii];
     }
      healpix_skymap_intensity_spectrum_sublimit_av[ii]/= (healpix_skymap_intensity.Npix()*oversample); 
      healpix_skymap_intensity_spectrum_soplimit_av[ii]/= (healpix_skymap_intensity.Npix()*oversample); 
   }

  if(FALSE)//AWS20110826
  {
   for(ii=0;ii<n_E;ii++) *txt_stream<<"healpix all sky average intensity spectrum sublimit at energy "<<E[ii]<<" : "<<  healpix_skymap_intensity_spectrum_sublimit_av[ii]<<" "<<endl;
   for(ii=0;ii<n_E;ii++) *txt_stream<<"healpix all sky average intensity spectrum soplimit at energy "<<E[ii]<<" : "<<  healpix_skymap_intensity_spectrum_soplimit_av[ii]<<" "<<endl;
  } 



   // === skymaps in binned energy ranges


   if(skymap_method==1)
     {
   spectrum_norm=integratepowbr(1.0,spectrum_g_0,spectrum_g_1,spectrum_g_2, spectrum_br0,spectrum_br1,E_ref_low,E_ref_high);//AWS20100818

   //   for (ii=0;ii<n_E_bin;ii++)
   for (ii=1;ii<n_E_bin;ii++) //AWS20100820 to avoid overwriting first energy range at start
   {
     //   spectrum_factor=(pow(E_bin_low[ii],-spectrum_g+1)-pow(E_bin_high[ii],-spectrum_g+1))
     //                   /(pow(E_ref_low    ,-spectrum_g+1)-pow(E_ref_high    ,-spectrum_g+1));
     //    *txt_stream<<" E= "<<E_bin_low[ii]<<" - "<<E_bin_high[ii]<<" spectrum factor="<<spectrum_factor<<endl;

    spectrum_factor=integratepowbr(1.0,spectrum_g_0,spectrum_g_1,spectrum_g_2, spectrum_br0,spectrum_br1,E_bin_low[ii],E_bin_high[ii]) //AWS20060123
                   /spectrum_norm;
    *txt_stream<<" E= "<<E_bin_low[ii]<<" - "<<E_bin_high[ii]<<" spectrum factor="<<spectrum_factor<<endl;

    for(il=0;il< longprof_nlong;il++)
    for(ib=0;ib<  latprof_nlat ;ib++)
    {

       skymap_intensity          .d2[il][ib].s[ii] = skymap_intensity          .d2[il][ib].s[0];
       skymap_intensity_sublimit .d2[il][ib].s[ii] = skymap_intensity_sublimit .d2[il][ib].s[0];;
       skymap_intensity_soplimit .d2[il][ib].s[ii] = skymap_intensity_soplimit .d2[il][ib].s[0];;

       skymap_intensity          .d2[il][ib].s[ii] *= spectrum_factor;
       skymap_intensity_sublimit .d2[il][ib].s[ii] *= spectrum_factor;
       skymap_intensity_soplimit .d2[il][ib].s[ii] *= spectrum_factor;

    } //for il ib

    //  cout<<" completed multiplying  l,b  binned energy range skymaps by spectral factor"<<endl;

    // healpix  
    for (ipix=0;ipix<healpix_skymap_intensity_spectrum.Npix();ipix++)                                             //AWS20090108
    {
     healpix_skymap_intensity         [ipix][ii] = healpix_skymap_intensity          [ipix][0] * spectrum_factor; //AWS20090108
     healpix_skymap_intensity_sublimit[ipix][ii] = healpix_skymap_intensity_sublimit [ipix][0] * spectrum_factor; //AWS20090108
     healpix_skymap_intensity_soplimit[ipix][ii] = healpix_skymap_intensity_soplimit [ipix][0] * spectrum_factor; //AWS20090108
    }


   } //for ii


   ii=0; //AWS20100820 assign the first bin correctly now

    spectrum_factor=integratepowbr(1.0,spectrum_g_0,spectrum_g_1,spectrum_g_2, spectrum_br0,spectrum_br1,E_bin_low[ii],E_bin_high[ii]) //AWS20060123
                   /spectrum_norm;
    *txt_stream<<" E= "<<E_bin_low[ii]<<" - "<<E_bin_high[ii]<<" spectrum factor="<<spectrum_factor<<endl;

    for(il=0;il< longprof_nlong;il++)
    for(ib=0;ib<  latprof_nlat ;ib++)
    {

       skymap_intensity          .d2[il][ib].s[ii] = skymap_intensity          .d2[il][ib].s[0];
       skymap_intensity_sublimit .d2[il][ib].s[ii] = skymap_intensity_sublimit .d2[il][ib].s[0];;
       skymap_intensity_soplimit .d2[il][ib].s[ii] = skymap_intensity_soplimit .d2[il][ib].s[0];;

       skymap_intensity          .d2[il][ib].s[ii] *= spectrum_factor;
       skymap_intensity_sublimit .d2[il][ib].s[ii] *= spectrum_factor;
       skymap_intensity_soplimit .d2[il][ib].s[ii] *= spectrum_factor;

    } //for il ib

    //  cout<<" completed multiplying  l,b  binned energy range skymaps by spectral factor"<<endl;

    // healpix  
    for (ipix=0;ipix<healpix_skymap_intensity_spectrum.Npix();ipix++)                                             //AWS20090108
    {
     healpix_skymap_intensity         [ipix][ii] = healpix_skymap_intensity          [ipix][0] * spectrum_factor; //AWS20090108
     healpix_skymap_intensity_sublimit[ipix][ii] = healpix_skymap_intensity_sublimit [ipix][0] * spectrum_factor; //AWS20090108
     healpix_skymap_intensity_soplimit[ipix][ii] = healpix_skymap_intensity_soplimit [ipix][0] * spectrum_factor; //AWS20090108
    }


  }//skymap_method==1



   //  cout<<" completed multiplying  healpix  binned energy range skymaps by spectral factor"<<endl;

  //-------------------------------------------------------------------- correct all results for oversampling    //AWS20060118
    for (ii=0;ii< n_dlnN_dlnS;ii++)
    {
     dlnN_dlnS        [ii]/= oversample;
     dlnN_dlnS_int    [ii]/= oversample;
     dlnN_dlnS_deg    [ii]/= oversample;
     dlnN_dlnS_int_deg[ii]/= oversample;
     FS               [ii]/= oversample;
     FS_int           [ii]/= oversample;

         dlnN_dlnS_sublimit [ii]/= oversample; //AWS20170123
     dlnN_dlnS_int_sublimit [ii]/= oversample; //AWS20170123
                FS_sublimit [ii]/= oversample; //AWS20170123
            FS_int_sublimit [ii]/= oversample; //AWS20170123

         dlnN_dlnS_soplimit [ii]/= oversample; //AWS20170123
     dlnN_dlnS_int_soplimit [ii]/= oversample; //AWS20170123
                FS_soplimit [ii]/= oversample; //AWS20170123
            FS_int_soplimit [ii]/= oversample; //AWS20170123


    }

     skymap_intensity_spectrum          /= oversample;
     skymap_intensity_spectrum_sublimit /= oversample;
     skymap_intensity_spectrum_soplimit /= oversample;

     skymap_intensity                   /= oversample;
     skymap_intensity_sublimit          /= oversample;
     skymap_intensity_soplimit          /= oversample;
    
     dlnN_dlnS_total_N                  /= oversample;
     dlnN_dlnS_total_S                  /= oversample;

     sample_selected_sourcecnt          /= oversample;
     sample_selected_sourcecnt_soplimit /= oversample;
     sample_selected_sourcecnt_sublimit /= oversample;
     sample_selected_total_flux         /= oversample;
     sample_selected_total_flux_soplimit/= oversample;
     sample_selected_total_flux_sublimit/= oversample;

     //  cout<<" completed dividing  l,b skymaps by oversample"<<endl;

     // healpix
     healpix_skymap_intensity_spectrum          /= oversample;
     healpix_skymap_intensity_spectrum_sublimit /= oversample;
     healpix_skymap_intensity_spectrum_soplimit /= oversample;

     healpix_skymap_intensity                   /= oversample;
     healpix_skymap_intensity_sublimit          /= oversample;
     healpix_skymap_intensity_soplimit          /= oversample;


     //  cout<<" completed dividing  healpix   skymaps by oversample"<<endl;


     // assign energy grid for healpix skymaps 
     valarray<double> skymap_energy_spectrum_grid(n_E);
     for (int iE=0;iE<n_E;iE++){skymap_energy_spectrum_grid[iE]=E[iE];}
     healpix_skymap_intensity_spectrum         .setSpectra(skymap_energy_spectrum_grid);
     healpix_skymap_intensity_spectrum_sublimit.setSpectra(skymap_energy_spectrum_grid);
     healpix_skymap_intensity_spectrum_soplimit.setSpectra(skymap_energy_spectrum_grid);

     valarray<double> skymap_energy_bins_E_low (n_E_bin);
     valarray<double> skymap_energy_bins_E_high(n_E_bin);
     for (int iE_bin=0;iE_bin<n_E_bin;iE_bin++){skymap_energy_bins_E_low [iE_bin]=E_bin_low [iE_bin];}
     for (int iE_bin=0;iE_bin<n_E_bin;iE_bin++){skymap_energy_bins_E_high[iE_bin]=E_bin_high[iE_bin];}

     healpix_skymap_intensity         .setSpectra(skymap_energy_bins_E_low, skymap_energy_bins_E_high);
     healpix_skymap_intensity_sublimit.setSpectra(skymap_energy_bins_E_low, skymap_energy_bins_E_high);
     healpix_skymap_intensity_soplimit.setSpectra(skymap_energy_bins_E_low, skymap_energy_bins_E_high);




     // write to FITS files
     int write_to_fits =1;//0; //AWS20140930 AWS20150612
     if( write_to_fits==1) //AWS20140930

     //         if(TRUE)//AWS20110902 try again AWS20140930
     // if(FALSE)//AWS20110817 was causing CCfitsError

 {
  cout<<" writing population synthesis skymaps to healpix files"<<endl;

  string skymap_filename;             //AWS20140131
  string fits_directory = "../FITS/"; //AWS20140131  not yet specified in class, so define as in Configure.cc for now

  if(skymap_method==1)
    {
     healpix_skymap_intensity         .write("../FITS/galplot_healpix_skymap_intensity_method_1.fits");
     healpix_skymap_intensity_sublimit.write("../FITS/galplot_healpix_skymap_intensity_sublimit_method_1.fits");
     healpix_skymap_intensity_soplimit.write("../FITS/galplot_healpix_skymap_intensity_soplimit_method_1.fits");
 
     healpix_skymap_intensity_spectrum         .write("../FITS/galplot_healpix_skymap_intensity_spectrum_method_1.fits");
     healpix_skymap_intensity_spectrum_sublimit.write("../FITS/galplot_healpix_skymap_intensity_spectrum_sublimit_method_1.fits");
     healpix_skymap_intensity_spectrum_soplimit.write("../FITS/galplot_healpix_skymap_intensity_spectrum_soplimit_method_1.fits");
    }

  if(skymap_method==2)
    {
     cout<<" writing population synthesis method 2 intensity skymaps to healpix files"<<endl;

     

     skymap_filename = fits_directory + string(title) + "_healpix_skymap_intensity_method_2.fits";                  //AWS20140131
                                                          healpix_skymap_intensity         .write(skymap_filename); //AWS20140131


     skymap_filename = fits_directory + string(title) + "_healpix_skymap_intensity_sublimit_method_2.fits";         //AWS20140131
                                                          healpix_skymap_intensity_sublimit.write(skymap_filename); //AWS20140131


     skymap_filename = fits_directory + string(title) + "_healpix_skymap_intensity_soplimit_method_2.fits";         //AWS20140131
                                                          healpix_skymap_intensity_soplimit.write(skymap_filename); //AWS20140131

							  /*                                                          AWS20140131
     healpix_skymap_intensity         .write("../FITS/galplot_healpix_skymap_intensity_method_2.fits");
     healpix_skymap_intensity_sublimit.write("../FITS/galplot_healpix_skymap_intensity_sublimit_method_2.fits");
     healpix_skymap_intensity_soplimit.write("../FITS/galplot_healpix_skymap_intensity_soplimit_method_2.fits");
							  */
     


     cout<<" writing population synthesis intensity spectrum  method 2 skymaps to healpix files"<<endl;

    
     skymap_filename = fits_directory + string(title) + "_healpix_skymap_intensity_spectrum_method_2.fits";                  //AWS20140131
                                                          healpix_skymap_intensity_spectrum         .write(skymap_filename); //AWS20140131


     skymap_filename = fits_directory + string(title) + "_healpix_skymap_intensity_spectrum_sublimit_method_2.fits";         //AWS20140131
                                                          healpix_skymap_intensity_spectrum_sublimit.write(skymap_filename); //AWS20140131


     skymap_filename = fits_directory + string(title) + "_healpix_skymap_intensity_spectrum_soplimit_method_2.fits";         //AWS20140131
                                                          healpix_skymap_intensity_spectrum_soplimit.write(skymap_filename); //AWS20140131

							  /*                                                                   AWS20140131
     healpix_skymap_intensity_spectrum         .write("../FITS/galplot_healpix_skymap_intensity_spectrum_method_2.fits");
     healpix_skymap_intensity_spectrum_sublimit.write("../FITS/galplot_healpix_skymap_intensity_spectrum_sublimit_method_2.fits");
     healpix_skymap_intensity_spectrum_soplimit.write("../FITS/galplot_healpix_skymap_intensity_spectrum_soplimit_method_2.fi

							  */

    }

     cout<<" writing population synthesis source count skymaps to healpix files"<<endl;


     skymap_filename = fits_directory + string(title) + "_healpix_skymap_sourcecnt.fits";                                      //AWS20140131
                                                          healpix_skymap_sourcecnt.write(skymap_filename);                     //AWS20140131

							 /*
     healpix_skymap_sourcecnt         .write("../FITS/galplot_healpix_skymap_sourcecnt.fits");
							 */

 } // if write_to_fits

         cout<<" printing results to file"<<endl;

  //--------------------------------------------------------------------  print results

   *txt_stream<<"source counts for region "<<long_min1 <<" < l < "<<long_max1 <<",  "<<long_min2 <<" < l < "<<long_max2<<" /  "
                                    << lat_min1 <<" < b < "<< lat_max1 <<",  "<< lat_min2 <<" < b < "<< lat_max2<<endl;
   *txt_stream<<"converted to      region "<<long_min1a<<" < l < "<<long_max1a<<",  "<<long_min2a<<" < l < "<<long_max2a<<endl;
                                
   cout<<"writing summary table to results txt for "<< title<<endl;         //AWS20200622

   *txt_stream<<"summary table of N(S) for "<< title <<endl;                //AWS20200622

   for (ii=0;ii< n_dlnN_dlnS;ii++)
   {
       if( dlnN_dlnS[ii]>0)
       *txt_stream <<"ii="<<ii
            <<" lnS="<<       lnS_min+ii*dlnS  <<"/"<<       lnS_min+(ii+1)*dlnS
	    <<" S="  <<pow(10,lnS_min+ii*dlnS) <<"/"<<pow(10,lnS_min+(ii+1)*dlnS)
            <<" N(S)= "<<  dlnN_dlnS    [ii]
            <<" N(>S)= "<< dlnN_dlnS_int[ii]
 //         <<" logN=  " << log10(dlnN_dlnS[ii]+1e-10)
  	    <<" N(S) deg^-2 ="  <<   dlnN_dlnS_deg    [ii]
  	    <<" N(>S) deg^-2 ="  <<  dlnN_dlnS_int_deg[ii]
            <<" F(S) = "<<         FS    [ii]
            <<" F(>S)= "<<         FS_int[ii]
            <<" l_av="<<   sample_l_av [ii]                
            <<" b_av="<<   sample_b_av [ii]  
            <<" l_min="<<  sample_l_min[ii]                
            <<" l_max="         <<  sample_l_max[ii]  
            <<" b_min="<<  sample_b_min[ii]                
            <<" b_max="         <<  sample_b_max[ii]  

            <<endl;
    }

   // output for external plotting routines                                 //AWS20200622

   *txt_stream<<"format for external plotting routines"<<endl;              //AWS20200622

   // idl

   *txt_stream<<"S = ["; 
   for (ii=0;ii< n_dlnN_dlnS;ii++)
       { *txt_stream<<pow(10, (lnS_min+ii*dlnS + lnS_min+(ii+1)*dlnS)/2. );
           if(ii<n_dlnN_dlnS-1)    *txt_stream<<",";
       }
   *txt_stream<<"]";
   *txt_stream<<";  idl format";
   *txt_stream<<endl;

  
   *txt_stream<<"NS = ["; 
   for (ii=0;ii< n_dlnN_dlnS;ii++)
       { *txt_stream<<dlnN_dlnS[ii];
           if(ii<n_dlnN_dlnS-1)    *txt_stream<<",";
       }
   *txt_stream<<"]";
   *txt_stream<<";  idl format";
   *txt_stream<<endl;



   //   *txt_stream<<"device,file='"<<title<<".ps'          ;idl"<<endl;

   if(IDL_plot_control==1)  // start plot
   {
    *txt_stream<<"set_plot,'ps'                         ;idl"<<endl;
    *txt_stream<<"device,file='"<<"sourcepop_NS"<<"_"<<galdef_ID<<"_"<<psfile_tag<<".ps' ;idl"  <<endl;
    *txt_stream<<              "gv sourcepop_NS"<<"_"<<galdef_ID<<"_"<<psfile_tag<<".ps &    "  <<endl;
    *txt_stream<<"!p.font=0                             ;idl"<<endl;
   }
  

   if(IDL_plot_control==1) // start plot
     *txt_stream<<"plot,S,NS,/xlog,/ylog,xrange=[1e-13,1e-5],yrange=[0.1,1e4], xstyle=1,ystyle=1, psym="<<IDL_psym<<", xtitle='S,  cm!e-2!n s!e-1!n',ytitle='N(S)', charsize=1.0,title='"<<galdef_ID<<"_"<<psfile_tag<<"',/nodata ; idl"<<endl;




   if(IDL_plot_control==1||IDL_plot_control==2) // overplot on start plot
     *txt_stream<<"oplot,S,NS,psym="<<IDL_psym<<"   ; idl"<<endl;


   // python format

  *txt_stream<<"title='" <<title<< "' # python format"<<endl;


  *txt_stream<<"S = ["; 
   for (ii=0;ii< n_dlnN_dlnS;ii++)
       { *txt_stream<<pow(10, (lnS_min+ii*dlnS + lnS_min+(ii+1)*dlnS)/2. );
           if(ii<n_dlnN_dlnS-1)    *txt_stream<<",";
       }
   *txt_stream<<"]";
   *txt_stream<<" #  python format";
   *txt_stream<<endl;

  
   *txt_stream<<"NS = ["; 
   for (ii=0;ii< n_dlnN_dlnS;ii++)
       { *txt_stream<<dlnN_dlnS[ii];
           if(ii<n_dlnN_dlnS-1)    *txt_stream<<",";
       }
   *txt_stream<<"]";
   *txt_stream<<" #  python format";
   *txt_stream<<endl;

   // end of format for external plotting routines



   *txt_stream<<"(oversampling corrected) dlnN_dlnS:  total N ="<< dlnN_dlnS_total_N  <<"           total S ="<< dlnN_dlnS_total_S<< endl;

   *txt_stream<<" sample region corrected for oversampling factor ="<<oversample                         <<endl;
   *txt_stream<<" sample region selected total sources            ="<<sample_selected_sourcecnt          <<endl;
   *txt_stream<<" sample region selected total sources above limit="<<sample_selected_sourcecnt_soplimit <<endl;
   *txt_stream<<" sample region selected total sources below limit="<<sample_selected_sourcecnt_sublimit <<endl;
   *txt_stream<<" sample region selected total flux               ="<<sample_selected_total_flux         <<endl;
   *txt_stream<<" sample region selected total flux   above limit ="<<sample_selected_total_flux_soplimit<<endl;
   *txt_stream<<" sample region selected total flux   below limit ="<<sample_selected_total_flux_sublimit<<endl;
   *txt_stream<<" sample region selected average distance         ="<<sample_selected_d_av               <<endl;
   *txt_stream<<" sample region selected average dist below limit ="<<sample_selected_d_av_sublimit      <<endl;
   *txt_stream<<" sample region selected average dist above limit ="<<sample_selected_d_av_soplimit      <<endl;
   *txt_stream<<" sample region selected min     distance         ="<<sample_selected_d_min              <<endl;
   *txt_stream<<" sample region selected min     dist below limit ="<<sample_selected_d_min_sublimit     <<endl;
   *txt_stream<<" sample region selected min     dist above limit ="<<sample_selected_d_min_soplimit     <<endl;
   *txt_stream<<" sample region selected max     distance         ="<<sample_selected_d_max              <<endl;
   *txt_stream<<" sample region selected max     dist below limit ="<<sample_selected_d_max_sublimit     <<endl;
   *txt_stream<<" sample region selected max     dist above limit ="<<sample_selected_d_max_soplimit     <<endl;

   if(FALSE)//AWS20110826
   {
   *txt_stream<<"longitude profile: " << lat_min1<<" < b < "<< lat_max1<<",  "<< lat_min2<<" < b < "<< lat_max2<<endl;

  for (ii=0;ii<longprof_nlong;ii++)
  {
   l=-180.+ (ii+.5)*longprof_dlong;
   *txt_stream<<"l="<<l
       <<" intensity="         << longprof_intensity         [ii]
       <<" below limit="       << longprof_intensity_sublimit[ii]
       <<" source count="      << longprof_sourcecnt         [ii]
       <<" below limit="       << longprof_sourcecnt_sublimit[ii]
       <<" above limit="       << longprof_sourcecnt         [ii]-longprof_sourcecnt_sublimit[ii]
       <<" / "                 <<(longprof_sourcecnt         [ii]-longprof_sourcecnt_sublimit[ii])/(solid_angle_long*rtd*rtd) <<" deg-2"
       <<endl;
  }


  *txt_stream<<" latitude profile: "<<long_min1 <<" < l < "<<long_max1 <<",  "<<long_min2 <<" < l < "<<long_max2<<endl;
  *txt_stream<<" converted to      "<<long_min1a<<" < l < "<<long_max1a<<",  "<<long_min2a<<" < l < "<<long_max2a<<endl;

  for (ii=0;ii< latprof_nlat ;ii++)
  {
   b=-90.+ (ii+.5)* latprof_dlat;
   //   solid_angle_lat = latprof_dlat *( latprof_l_max- latprof_l_min)/(rtd*rtd)/pow(cos(b*dtr),2);
    b1=-90.+ ii* latprof_dlat;
    b2= b1 +     latprof_dlat;
    solid_angle_lat = (long_max1-long_min1) *dtr  *(sin(b2*dtr) - sin(b1*dtr))
                    + (long_max2-long_min2) *dtr  *(sin(b2*dtr) - sin(b1*dtr));
   *txt_stream<<"b="<<b
       <<" intensity="         <<  latprof_intensity         [ii]
       <<" below limit="       <<  latprof_intensity_sublimit[ii]
       <<" source count="      <<  latprof_sourcecnt         [ii]
       <<" below limit="       <<  latprof_sourcecnt_sublimit[ii]
       <<" above limit="       <<  latprof_sourcecnt         [ii]- latprof_sourcecnt_sublimit[ii]
       <<" / "                 <<( latprof_sourcecnt         [ii]- latprof_sourcecnt_sublimit[ii])
                                  /(solid_angle_lat*rtd*rtd) <<" deg-2"
       <<endl;
  }
 }//if(FALSE)

  if(verbose==2)
  {
  *txt_stream<<endl;

  *txt_stream<<"skymap: intensity spectrum"<<endl;
  skymap_intensity_spectrum.print();
  *txt_stream<<"skymap: intensity spectrum   below  detection limit"<<endl;
  skymap_intensity_spectrum_sublimit.print();
  *txt_stream<<"skymap: intensity spectrum   above  detection limit"<<endl;
  skymap_intensity_spectrum_soplimit.print();

  *txt_stream<<"skymap: intensity"<<endl;
  skymap_intensity.print();
  *txt_stream<<"skymap: intensity    below detection limit"<<endl;
  skymap_intensity_sublimit.print();
  *txt_stream<<"skymap: intensity    above detection limit"<<endl;
  skymap_intensity_soplimit.print();

  *txt_stream<<"skymap: source count"<<endl;
  skymap_sourcecnt.print();
  *txt_stream<<"skymap:  source count below detection limit"<<endl;
  skymap_sourcecnt_sublimit.print();
  *txt_stream<<"skymap: source count  above detection limit"<<endl;
  skymap_sourcecnt_soplimit.print();
  }


  *txt_stream<< " << SourcePopulation:analyse_sample complete"<<endl<<endl;

  cout       << " << SourcePopulation:analyse_sample"<<endl;

  return 0;
}

/////////////////////////////////////////////////////////
 int SourcePopulation::gen_sample_EGRET_catalogue(int options)
{

  int i_sample,i_source;
  double R_,z_,theta_,d_;
  EGRET_catalogue egret_catalogue;

  
  egret_catalogue .read(0);
  egret_catalogue. select(options);
  egret_catalogue .print(*txt_stream);

  
  *txt_stream<<"total number of EGRET sources ="<<egret_catalogue .n_source<<endl;

  n_sample    =long(egret_catalogue.n_source);  
  sample_R    =new double[n_sample];
  sample_z    =new double[n_sample];
  sample_theta=new double[n_sample];
  sample_l    =new double[n_sample];
  sample_b    =new double[n_sample];
  sample_d    =new double[n_sample];
  sample_flux =new double[n_sample];

   sample_flux_spectra        = new double*[n_sample];                                                                          //AWS20100818
   sample_binned_flux_spectra = new double*[n_sample];   ;                                                                      //AWS20100818
   for (int i=0;i<n_sample;i++){ sample_flux_spectra[i]=new double[n_E];  sample_binned_flux_spectra[i]=new double[n_E_bin];   }//AWS20100818

  // values from catalogue
  i_sample=0;
  for (i_source=0;i_source<n_sample;i_source++)
    {
      if(egret_catalogue.selected[i_source]==1)
	{
        sample_l[i_sample]= egret_catalogue.l[i_source];
     if(sample_l[i_sample] >180.)
        sample_l[i_sample]-=360.;           // SourcePopulation uses -180<l<+180
     
     sample_b   [i_sample]= egret_catalogue.b[i_source];
     sample_flux[i_sample]= egret_catalogue.F[i_source]*1e-8; // unit of EGRET catalogue fluxes = 1e-8  AWS20051206
    
     if(verbose==4)
       *txt_stream<<"selected EGRET source "<< egret_catalogue.name[i_source]<<"  l="<< sample_l[i_sample]<<" b="<< sample_b[i_sample]<<" flux="<< sample_flux[i_sample]<<endl;

     i_sample++;
	}
    }

  n_sample=i_sample;

  oversample=1.0; // AWS200601018

  // fictitious values to complete SourcePopulation
  R_    = 5.;
  z_    = 0.;
 theta_ = 0.;
 d_     = 5.;



  for (i_sample=0;i_sample<n_sample;i_sample++)
    {
     sample_R    [i_sample]= R_;
     sample_z    [i_sample]= z_;
     sample_theta[i_sample]= theta_;
     sample_d    [i_sample]= d_;
    }


  

  return 0;
}
/////////////////////////////////////////////////////////
 int SourcePopulation::gen_sample_SPI_catalogue(char *directory, char* filename, int options)
{

  int i_sample,i_source;
  double R_,z_,theta_,d_;
  CAT spi_catalogue;

  
  spi_catalogue .read(directory,filename);
  
  //  spi_catalogue .print(*txt_stream);
  spi_catalogue .print();

  
  *txt_stream<<"total number of SPI sources ="<<spi_catalogue .n_sources<<endl;

  n_sample    =long(spi_catalogue.n_sources);  
  sample_R    =new double[n_sample];
  sample_z    =new double[n_sample];
  sample_theta=new double[n_sample];
  sample_l    =new double[n_sample];
  sample_b    =new double[n_sample];
  sample_d    =new double[n_sample];
  sample_flux =new double[n_sample];

   sample_flux_spectra        = new double*[n_sample];                                                                          //AWS20100818
   sample_binned_flux_spectra = new double*[n_sample];   ;                                                                      //AWS20100818
   for (int i=0;i<n_sample;i++){ sample_flux_spectra[i]=new double[n_E];  sample_binned_flux_spectra[i]=new double[n_E_bin];   }//AWS20100818

  // values from catalogue
  i_sample=0;
  for (i_source=0;i_source<n_sample;i_source++)
    {
   
        sample_l[i_sample]= spi_catalogue.longitude[i_source];
     if(sample_l[i_sample] >180.)
        sample_l[i_sample]-=360.;           // SourcePopulation uses -180<l<+180
     
     sample_b   [i_sample]= spi_catalogue.latitude[i_source];
     sample_flux[i_sample]= spi_catalogue.flux    [i_source];
    
     if(verbose==4)
       *txt_stream<<"selected SPI source "<< spi_catalogue.name[i_source]<<"  l="<< sample_l[i_sample]<<" b="<< sample_b[i_sample]<<" flux="<< sample_flux[i_sample]<<endl;

     i_sample++;
	
    }

  n_sample=i_sample;

  oversample=1.0; // AWS200601018

  // fictitious values to complete SourcePopulation
  R_    = 5.;
  z_    = 0.;
 theta_ = 0.;
 d_     = 5.;

  for (i_sample=0;i_sample<n_sample;i_sample++)
    {
     sample_R    [i_sample]= R_;
     sample_z    [i_sample]= z_;
     sample_theta[i_sample]= theta_;
     sample_d    [i_sample]= d_;
    }


  

  return 0;
}

/////////////////////////////////////////////////////////

int SourcePopulation::gen_sample_Fermi_catalogue(char *directory, char *filename,int options) //AWS20081215
{

  // generate a population of sources from an input Fermi catalogue

  int i_sample,i_source;
  double R_,z_,theta_,d_;
 

  Sources sources;
  
  
  cout<<"SourcePopulation::gen_sample_Fermi_catalogue"<<endl;  
  //  cout<<"reading Fermi catalogue from "<<directory<<" "<< filename<<endl;
  string Fermi_cat_file=string(directory)+string(filename);
  cout<<"reading Fermi catalogue from "<<Fermi_cat_file<<endl;
 
  sources.read(Fermi_cat_file);
  cout<<"total number of Fermi sources ="<<sources.size()<<endl;

  if(options==  100)cout<<"Using Fermi source flux >   100 MeV"<<endl;
  if(options== 1000)cout<<"Using Fermi source flux >  1000 MeV"<<endl;
  if(options==10000)cout<<"Using Fermi source flux > 10000 MeV"<<endl;//AWS20120530

  n_sample    =long(sources.size());  
  sample_R    =new double[n_sample];
  sample_z    =new double[n_sample];
  sample_theta=new double[n_sample];
  sample_l    =new double[n_sample];
  sample_b    =new double[n_sample];
  sample_d    =new double[n_sample];
  sample_flux =new double[n_sample];

   sample_flux_spectra        = new double*[n_sample];                                                                          //AWS20100818
   sample_binned_flux_spectra = new double*[n_sample];   ;                                                                      //AWS20100818
   for (int i=0;i<n_sample;i++){ sample_flux_spectra[i]=new double[n_E];  sample_binned_flux_spectra[i]=new double[n_E_bin];   }//AWS20100818

  // values from catalogue
  i_sample=0;
  for (i_source=0;i_source<n_sample;i_source++)
    {
      cout<<"Fermi Catalogue: l="<<sources[i_source].l<<" b="<<sources[i_source].b
          <<" flux100="               <<sources[i_source].flux100
          <<" flux1000="              <<sources[i_source].flux1000  
	  <<" flux10000="             <<sources[i_source].flux10000    //AWS20120530  
          <<" flux_density=prefactor="<<sources[i_source].prefactor    //AWS20140927
          <<" spectral index="        <<sources[i_source].index        //AWS20140927
          <<" pivot energy="          <<sources[i_source].pivot        //AWS20140927
          <<endl;//AWS20110806

        sample_l[i_sample]= sources[i_source].l;
     if(sample_l[i_sample] >180.)
        sample_l[i_sample]-=360.;           // SourcePopulation uses -180<l<+180


	sample_b[i_sample]= sources[i_source].b;

     if(options==  100)                                    //AWS20110808
     sample_flux[i_sample]= sources[i_source].flux100;

     if(options== 1000)                                    //AWS20110808
     sample_flux[i_sample]= sources[i_source].flux1000;    //AWS20110808

     if(options==10000)                                    //AWS20120530
     sample_flux[i_sample]= sources[i_source].flux10000;   //AWS20120530

      i_sample++;
    }



  n_sample=i_sample;

  oversample=1.0; // AWS200601018

  // fictitious values to complete SourcePopulation
  R_    = 5.;
  z_    = 0.;
 theta_ = 0.;
 d_     = 5.;

  

  for (i_sample=0;i_sample<n_sample;i_sample++)
    {
     sample_R    [i_sample]= R_;
     sample_z    [i_sample]= z_;
     sample_theta[i_sample]= theta_;
     sample_d    [i_sample]= d_;
    }
  
  
  
  

  return 0;
}

////////////////////////////////////////////////////////////////////////////////////////////

 int  SourcePopulation::gen_Fermi_catalogue_from_sample(char *directory, char *filename, int options)//AWS20110829

{
  // generate a Fermi catalogue from the population synthesis sample

  Sources         sources; // a catalogue of sources
  Sources::Source source ; // one source, struture defined in Sources



  cout<<">> SourcePopulation::gen_Fermi_catalogue_from_sample"<<endl;
  cout<<" options="<<options<<endl;

  for (int i=1;i<n_sample;i++)
  {
   char source_number_char[20];
   string source_number; sprintf(source_number_char,"%d",i); 
   source.name="popsyn"+string(source_number_char);
   source.l = sample_l[i];
   source.b = sample_b[i];
  

  source.flux100  =sample_flux100  [i]; //AWS20120615
  source.flux1000 =sample_flux1000 [i]; //AWS20120615
  source.flux10000=sample_flux10000[i]; //AWS20120615

  source.fitType=Sources::strToFitType("FIXED");

  if(spectral_model==1)                                        //AWS20120615
  source. spType=Sources::strToFitType("BROKENPOWERLAW");

  if(spectral_model==2)                                        //AWS20120615
  source. spType=Sources::strToFitType("BROKENPOWERLAW");      //AWS20120615

  if(spectral_model==3)                                        //AWS20110901
  source. spType=Sources::strToFitType("POWERLAWEXPCUTOFF");   //AWS20110901


  source.spParameters.resize(5);                               //AWS20120614 was (2)
  source.spParameters[0]=sample_spectrum_g_0[i];               //AWS20110901
  source.spParameters[1]=sample_spectrum_br0[i];               //AWS20110901
  source.spParameters[2]=sample_spectrum_g_1[i];               //AWS20120614
  source.spParameters[3]=sample_spectrum_br1[i];               //AWS20120614
  source.spParameters[4]=sample_spectrum_g_2[i];               //AWS20120614

  source.spParametersUnc.resize(5);                            //AWS20120614 was (2)
  source.spParametersUnc[0]=0.0;
  source.spParametersUnc[1]=0.0;
  source.spParametersUnc[2]=0.0;                               //AWS20120614
  source.spParametersUnc[3]=0.0;                               //AWS20120614
  source.spParametersUnc[4]=0.0;                               //AWS20120614

  source.index=sample_spectrum_g_0[i];                         //AWS20110901
  source.pivot=                         E_ref_low;             //AWS20110901
  source.prefactor=sample_flux_spectrum_E_ref_low[i];          //AWS20110901 since Source outputs this as flux_density OK

  double pi=acos(-1.); double dtr=pi/180.;                                                       //AWS20120612
  slaGaleq(source.l*dtr, source.b*dtr, &source.ra, &source.dec);                                 //AWS20120612

  source.ra/=dtr; source.dec/=dtr;                                                               //AWS20120612

  // cout<<"from slaGaleq: l="<<source.l<<" b="<<source.b<<" ra="<<source.ra<<" dec="<<source.dec<<endl;

  

  sources.add_source(source);

  }

  cout<<" number of sample sources in catalogue for output = "<<sources.size()<<endl;

  //  string Fermi_cat_sample_file=string(directory)+string(filename);

  string Fermi_cat_sample_file=string(filename); //AWS20110830

 int write_Fermi_catalogue =1; //0; //AWS20140930 AWS20150612
 if( write_Fermi_catalogue==1) //AWS20140930
 {
  cout<<"output Fermi catalogue from sample to file  "<<Fermi_cat_sample_file<<endl;
  sources.write(Fermi_cat_sample_file);
 }
  cout<<"<< SourcePopulation::gen_Fermi_catalogue_from_sample"<<endl;
  return 0;
}

///////////////////////////////////////////////////////////////////////////////////

int SourcePopulation::read_sample_MSPs (char *directory, char *filename, int options)//AWS20120521
{
  cout<<">> SourcePopulation::read_sample_MSPs"<<endl;
  int debug=0;

  double flux_factor=0.1;// convert from >100 MeV to >1 GeV, provisional fix

 fitsfile *fptr;
 int status,colnum;
 int hdunum,hdutype,total_hdus;
 long nrows;
 char colname[100];
 
 string MSPs_sample_file;
 MSPs_sample_file=string(directory)+string(filename);

 cout<<"reading from "<<MSPs_sample_file<<endl;
 status=0;
 fits_open_file(&fptr, MSPs_sample_file.c_str(), READONLY, &status);
 fits_report_error(stderr, status);

fits_get_num_hdus(fptr,&total_hdus,&status);
cout<<"total number of header units="<<total_hdus<<endl;
hdunum=1;
fits_movabs_hdu(fptr,hdunum+1,&hdutype,&status );
cout<<"FITS movabs_hdu hdunum="<<hdunum<<" status= "<<status<<" hdutype="<<hdutype<<endl  ;



 strcpy(colname,"FLUX");
 fits_get_colnum  (fptr,1,colname,&colnum,&status); // 1=case-sensitive
 fits_report_error(stderr, status);
 cout<<colname<<" column number="<<colnum<<endl;


 fits_get_num_rows(fptr,          &nrows, &status);
 fits_report_error(stderr, status);
 cout<<"number of rows="<<nrows<<endl;

 n_sample=nrows;


 sample_flux=new double[n_sample];
 double nulval=0.;
 int    anynul=0 ;
 fits_read_col(fptr,TDOUBLE, colnum, 1, 1, nrows, &nulval, sample_flux, &anynul, &status);
 fits_report_error(stderr, status);

 if(debug==1)
 for (int i=0;i<nrows;i++)cout<<"sample_flux= "<<sample_flux[i]<<endl;

 cout<<" applying flux_factor = "<<flux_factor<<endl;
 for (int i=0;i<nrows;i++)sample_flux[i]*=flux_factor;



 strcpy(colname,"LDEG");
 fits_get_colnum  (fptr,1,colname,&colnum,&status); // 1=case-sensitive
 fits_report_error(stderr, status);
 cout<<colname<<" column number="<<colnum<<endl;

 sample_l=new double[n_sample];
 fits_read_col(fptr,TDOUBLE, colnum, 1, 1, nrows, &nulval, sample_l   , &anynul, &status);
 fits_report_error(stderr, status);

 // NB input file has -180<l<180 which is as for generated samples
 

 if(debug==1)
 for (int i=0;i<nrows;i++)cout<<"sample_l= "<<sample_l[i]<<endl;

 strcpy(colname,"BDEG");
 fits_get_colnum  (fptr,1,colname,&colnum,&status); // 1=case-sensitive
 fits_report_error(stderr, status);
 cout<<colname<<" column number="<<colnum<<endl;

 sample_b=new double[n_sample];
 fits_read_col(fptr,TDOUBLE, colnum, 1, 1, nrows, &nulval, sample_b   , &anynul, &status);
 fits_report_error(stderr, status);

 if(debug==1)
 for (int i=0;i<nrows;i++)cout<<"sample_b= "<<sample_b[i]<<endl;



 strcpy(colname,"REARTH");
 fits_get_colnum  (fptr,1,colname,&colnum,&status); // 1=case-sensitive
 fits_report_error(stderr, status);
 cout<<colname<<" column number="<<colnum<<endl;

 sample_d=new double[n_sample];
 fits_read_col(fptr,TDOUBLE, colnum, 1, 1, nrows, &nulval, sample_d   , &anynul, &status);
 fits_report_error(stderr, status);

 if(debug==1)
 for (int i=0;i<nrows;i++)cout<<"sample_d= "<<sample_d[i]<<endl;

 // for completeness, not yet read in
 sample_R    =new double[n_sample];
 sample_z    =new double[n_sample];
 sample_theta=new double[n_sample];

 oversample = 1;

   sample_flux_spectra        = new double*[n_sample];
   sample_binned_flux_spectra = new double*[n_sample];


   //   n_E=1;
   //   n_E_bin=1;
   cout<<"MSPs allocating sample flux spectra for testing only here : n_E="<<n_E<< " n_E_bin="<<n_E_bin<<endl;

  for (int i=0;i<n_sample;i++)
  { sample_flux_spectra[i]=new double[n_E];  sample_binned_flux_spectra[i]=new double[n_E_bin];   }


  // generate spectra

  sample_flux_spectrum_E_ref_low  =new double[n_sample]; 
  sample_flux_spectrum_E_ref_high =new double[n_sample]; 

   sample_flux_spectra        = new double*[n_sample];                                                                         sample_binned_flux_spectra = new double*[n_sample];                                                                       
   for (int i=0;i<n_sample;i++)
     { sample_flux_spectra[i]=new double[n_E];  sample_binned_flux_spectra[i]=new double[n_E_bin]; }

   sample_spectrum_g_0        = new double[n_sample];                                                                          sample_spectrum_br0        = new double[n_sample];                                                                          sample_spectrum_g_1        = new double[n_sample];                                                                          sample_spectrum_br1        = new double[n_sample];                                                                          sample_spectrum_g_2        = new double[n_sample];  

   
   // case of no dispersion in spectra : use for MSPs for now
   spectral_model=1;

   if(spectral_model==1)
  {
   for (int i=0;i<n_sample;i++)                                                                                             
   {
    sample_spectrum_g_0[i] = spectrum_g_0;
    sample_spectrum_br0[i] = spectrum_br0;
    sample_spectrum_g_1[i] = spectrum_g_1;
    sample_spectrum_br1[i] = spectrum_br1;
    sample_spectrum_g_2[i] = spectrum_g_2;
   }
  }

   

   for (int i=0;i<n_sample;i++)
   {
    double spectrum_norm=integratepowbr(1.0,sample_spectrum_g_0[i],sample_spectrum_g_1[i],sample_spectrum_g_2[i],
                                                  sample_spectrum_br0[i],sample_spectrum_br1[i],E_ref_low,E_ref_high);

     // spectral energies
     for (int iE=0;iE<n_E;iE++)
     {
	    
            
      double spectrum_factor=shapepowbr(1.0,sample_spectrum_g_0[i],sample_spectrum_g_1[i],sample_spectrum_g_2[i],
                                            sample_spectrum_br0[i],sample_spectrum_br1[i],E[iE])
	                                                              //  *pow(E[iE], 2.0 )// NB need *E^2 for compatibility with galprop (l,b) skymaps, applied when required
	                                   / spectrum_norm; 


	  //          *txt_stream<<" E= "<<E[iE]<<" spectrum_norm="<<spectrum_norm<<" spectrum factor with breaks="<<spectrum_factor<<endl;
	  //	  cout<<" br0="<<sample_spectrum_br0[i]<<" "<<" E= "<<E[iE]<<" spectrum_norm="<<spectrum_norm<<" spectrum factor with breaks="<<spectrum_factor<<endl;

           sample_flux_spectra[i][iE]=sample_flux[i]*spectrum_factor;

           if(spectral_model==3)
               sample_flux_spectra[i][iE]
            = flux_cutoff_powerlaw( E[iE], sample_spectrum_g_0[i], sample_spectrum_br0[i], E_ref_low, E_ref_high, sample_flux[i]);


          }
  

	 // reference energies

           if(spectral_model==1 || spectral_model==2)
	   {
            double spectrum_factor=shapepowbr(1.0,sample_spectrum_g_0[i],sample_spectrum_g_1[i],sample_spectrum_g_2[i],
                                                sample_spectrum_br0[i],sample_spectrum_br1[i],E_ref_low)
	                                   / spectrum_norm; 

	    sample_flux_spectrum_E_ref_low[i]=sample_flux[i]*spectrum_factor;

                   spectrum_factor=shapepowbr(1.0,sample_spectrum_g_0[i],sample_spectrum_g_1[i],sample_spectrum_g_2[i],
                                                sample_spectrum_br0[i],sample_spectrum_br1[i],E_ref_high)
	                                   / spectrum_norm; 

            sample_flux_spectrum_E_ref_high[i]=sample_flux[i]*spectrum_factor;

           }


           if(spectral_model==3)
	   {
            sample_flux_spectrum_E_ref_low[i]
             = flux_cutoff_powerlaw
                  ( E_ref_low,  sample_spectrum_g_0[i], sample_spectrum_br0[i], E_ref_low, E_ref_high, sample_flux[i]);

            sample_flux_spectrum_E_ref_high[i]
             = flux_cutoff_powerlaw
                  ( E_ref_high, sample_spectrum_g_0[i], sample_spectrum_br0[i], E_ref_low, E_ref_high, sample_flux[i]);
           }


        //  binned energy ranges

         for (int iE=0;iE<n_E_bin;iE++)
         {
    

         double spectrum_factor=integratepowbr(1.0,sample_spectrum_g_0[i],sample_spectrum_g_1[i],sample_spectrum_g_2[i], sample_spectrum_br0[i],sample_spectrum_br1[i],E_bin_low[iE],E_bin_high[iE])	   /spectrum_norm; 


          sample_binned_flux_spectra[i][iE]=sample_flux[i]*spectrum_factor;

           if(spectral_model==3)
	         sample_binned_flux_spectra[i][iE]
             = flux_cutoff_powerlaw( E_bin_low[iE] ,E_bin_high[iE], sample_spectrum_g_0[i], sample_spectrum_br0[i], E_ref_low, E_ref_high, sample_flux[i]);


	  //      *txt_stream<<" E= "<<E_bin_low[iE]<<" - "<<E_bin_high[iE]<<" source binned spectrum factor="<<spectrum_factor<<endl;
          }

	      if(verbose==1)
	{
	  *txt_stream<<"i="<<i<<" R=" <<sample_R[i]<< " z="<<sample_z[i]<<" theta="<<sample_theta[i]
		      <<" L="<<L
			     <<" distance="<<sample_d[i]   <<" l="<<sample_l[i]    <<" b="<<sample_b[i]      <<" flux="<<sample_flux[i];
          *txt_stream  <<" flux spectrum= " ;              for (int iE=0;iE<n_E    ;iE++){*txt_stream<< sample_flux_spectra      [i][iE]<<" ";}
          *txt_stream  <<" binned flux spectrum= " ;       for (int iE=0;iE<n_E_bin;iE++){*txt_stream<< sample_binned_flux_spectra[i][iE]<<" ";}

          *txt_stream      <<endl;
	}
 


     }// sample

 fits_close_file(fptr, &status);
  cout<<"<< SourcePopulation::read_sample_MSPs"<<endl;
 return 0;
}

///////////////////////////////////////////////////////////////////////////////////
 SourcePopulation  SourcePopulation::operator=(SourcePopulation sourcepop)
{
  int i;

  init();
  verbose                 =sourcepop.verbose;
  luminosity_function_type=sourcepop.luminosity_function_type;
  L                       =sourcepop.L;
  density0                =sourcepop.density0;
  oversample              =sourcepop.oversample; //AWS20060118
  alpha_R                 =sourcepop.alpha_R;
   beta_R                 =sourcepop. beta_R;
  zscale                  =sourcepop.zscale ;

  alpha_z                 =sourcepop.alpha_z;
  L_min                   =sourcepop.L_min   ;
  L_max                   =sourcepop.L_max  ;
  alpha_L                 =sourcepop.alpha_L;

  spectral_model          =sourcepop.spectral_model; //AWS20100824
  spectrum_g              =sourcepop.spectrum_g  ;
  spectrum_g_0            =sourcepop.spectrum_g_0;
  spectrum_g_1            =sourcepop.spectrum_g_1;
  spectrum_g_2            =sourcepop.spectrum_g_2;
  spectrum_br0            =sourcepop.spectrum_br0;
  spectrum_br1            =sourcepop.spectrum_br1;

  spectrum_g_0_sigma      =sourcepop.spectrum_g_0_sigma; //AWS20100824
  spectrum_g_1_sigma      =sourcepop.spectrum_g_1_sigma; //AWS20100824
  spectrum_g_2_sigma      =sourcepop.spectrum_g_2_sigma; //AWS20100824
  spectrum_br0_sigma      =sourcepop.spectrum_br0_sigma; //AWS20100824
  spectrum_br1_sigma      =sourcepop.spectrum_br1_sigma; //AWS20100824

  spectrum_norm_E         =sourcepop.spectrum_norm_E;
  n_E_bin                 =sourcepop.n_E_bin        ;
  E_ref_low               =sourcepop.E_ref_low      ;
  E_ref_high              =sourcepop.E_ref_high     ;
  n_E                     =sourcepop.n_E            ;

  long_min1               =sourcepop.long_min1      ;
  long_max1               =sourcepop.long_max1      ;
  long_min2               =sourcepop.long_min2      ;
  long_max2               =sourcepop.long_max2      ;

   lat_min1               =sourcepop. lat_min1      ;
   lat_max1               =sourcepop. lat_max1      ;
   lat_min2               =sourcepop. lat_min2      ;
   lat_max2               =sourcepop. lat_max2      ;

   longprof_dlong         =sourcepop. longprof_dlong;
    latprof_dlat          =sourcepop.  latprof_dlat ;

   healpix_skymap_order   =sourcepop.healpix_skymap_order; //AWS20090107

   flux_detection_limit   =sourcepop. flux_detection_limit;

   E_bin_low =new double[n_E_bin];
   E_bin_high=new double[n_E_bin];
   for (i=0;i<n_E_bin;i++) 
   {     
    E_bin_low [i]=   sourcepop.E_bin_low [i]    ;
    E_bin_high[i]=   sourcepop.E_bin_high[i]    ;
   }

   E         =new double[n_E];
   for (i=0;i<n_E;i++)       E[i]=   sourcepop.E[i]    ;

   txt_stream = sourcepop.txt_stream;

  return *this;
}
///////////////////////////////////////////////////////////////////
 double SourcePopulation::integratepowbr(double A0,double g0,double g1,double g2, double br0,double br1, double E1,double E2)
{

  // integral of power-law with two breaks, from E1 to E2
  // A0 E^-g0      A0 input
  // A1 E^-g1      defined by continuity at br0
  // A2 E^-g2      defined by continuity at br1

  double result;
  double result0,result1,result2;
  double Ea,Eb;
  double A1,A2;



  A1=A0*pow(br0,-g0+g1); //continuity at br0
  A2=A1*pow(br1,-g1+g2); //continuity at br1

  result0=0.;
  result1=0.;
  result2=0.;

  // sketch of ranges
  //      E1-------------E2
  //0<<---------br0-----br1------------>>inf
  //      g0        g1         g2 
 
  //      E1-------------E2
  //0<<---------br0
  //      Ea----Eb


  Ea=    E1;
  Eb=min(E2,br0);
  if(Eb>Ea)      result0=A0*(pow(Ea,-g0+1)-pow(Eb,-g0+1))/(g0-1);
  
  if(verbose==1) //AWS20100817
  *txt_stream<<"SourcePopulation::integratepowbr: Ea="<<Ea<<" Eb="<<Eb<<" result0="<<result0<<endl;

  //      E1-------------E2
  //        br0------------br1
  //        Ea-----------Eb

  Ea=max(E1,br0);
  Eb=min(E2,br1);
  if(Eb>Ea)      result1=A1*(pow(Ea,-g1+1)-pow(Eb,-g1+1))/(g1-1);

  if(verbose==1) //AWS20100817
  *txt_stream<<"SourcePopulation::integratepowbr: Ea="<<Ea<<" Eb="<<Eb<<" result1="<<result1<<endl;


  //      E1-------------E2
  //          br1------------>>inf
  //          Ea---------Eb

  Ea=max(E1,br1);
  Eb=    E2     ;
  if(Eb>Ea)      result2=A2*(pow(Ea,-g2+1)-pow(Eb,-g2+1))/(g2-1);

  if(verbose==1) //AWS20100817
  *txt_stream<<"SourcePopulation::integratepowbr: Ea="<<Ea<<" Eb="<<Eb<<" result2="<<result2<<endl;

  result=result0+result1+result2;

  if(verbose==1) //AWS20100817
  *txt_stream<<"SourcePopulation::integratepowbr: E1="<<E1<<" E2="<<E2<<" result ="<<result <<endl;

  return result;
 
}
/* remove next, duplicate: !
///////////////////////////////////////////////////////////////////
 double SourcePopulation::integratepowbr(double A0,double g0,double g1,double g2, double br0,double br1, double E1,double E2)
{

  // integral of power-law with two breaks, from E1 to E2
  // A0 E^-g0      A0 input
  // A1 E^-g1      defined by continuity at br0
  // A2 E^-g2      defined by continuity at br1

  double result;
  double result0,result1,result2;
  double Ea,Eb;
  double A1,A2;



  A1=A0*pow(br0,-g0+g1); //continuity at br0
  A2=A1*pow(br1,-g1+g2); //continuity at br1

  result0=0.;
  result1=0.;
  result2=0.;

  // sketch of ranges
  //      E1-------------E2
  //0<<---------br0-----br1------------>>inf
  //      g0        g1         g2 
 
  //      E1-------------E2
  //0<<---------br0
  //      Ea----Eb


  Ea=    E1;
  Eb=min(E2,br0);
  if(Eb>Ea)      result0=A0*(pow(Ea,-g0+1)-pow(Eb,-g0+1))/(g0-1);

  if(verbose==1) //AWS20100817
  *txt_stream<<"SourcePopulation::integratepowbr: Ea="<<Ea<<" Eb="<<Eb<<" result0="<<result0<<endl;

  //      E1-------------E2
  //        br0------------br1
  //        Ea-----------Eb

  Ea=max(E1,br0);
  Eb=min(E2,br1);
  if(Eb>Ea)      result1=A1*(pow(Ea,-g1+1)-pow(Eb,-g1+1))/(g1-1);

  if(verbose==1) //AWS20100817
  *txt_stream<<"SourcePopulation::integratepowbr: Ea="<<Ea<<" Eb="<<Eb<<" result1="<<result1<<endl;


  //      E1-------------E2
  //          br1------------>>inf
  //          Ea---------Eb

  Ea=max(E1,br1);
  Eb=    E2     ;
  if(Eb>Ea)      result2=A2*(pow(Ea,-g2+1)-pow(Eb,-g2+1))/(g2-1);

  if(verbose==1) //AWS20100817
  *txt_stream<<"SourcePopulation::integratepowbr: Ea="<<Ea<<" Eb="<<Eb<<" result2="<<result2<<endl;

  result=result0+result1+result2;

  if(verbose==1) //AWS20100817
  *txt_stream<<"SourcePopulation::integratepowbr: E1="<<E1<<" E2="<<E2<<" result ="<<result <<endl;

  return result;
 
}
*/
///////////////////////////////////////////////////////////////////
 double SourcePopulation::shapepowbr(double A0,double g0,double g1,double g2, double br0,double br1, double E)
{

  // shape of power-law with two breaks, from E1 to E2
  // A0 E^-g0      A0 input
  // A1 E^-g1      defined by continuity at br0
  // A2 E^-g2      defined by continuity at br1

  double A1,A2;
  double result;
  



  A1=A0*pow(br0,-g0+g1); //continuity at br0
  A2=A1*pow(br1,-g1+g2); //continuity at br1



  // sketch of ranges
 
  //0<<---------br0-----br1------------>>inf
  //      g0        g1         g2 
 
  //      


  
  if(E <br0)      result=A0*pow(E,-g0);
  if(E>=br0)      result=A1*pow(E,-g1);
  if(E> br1)      result=A2*pow(E,-g2);

  if(verbose==1) //AWS20100817
  *txt_stream<<"SourcePopulation::shapepowbr:  E="<<E<<" result="<<result<<endl;



  return result;
 
}
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
double   SourcePopulation::flux_cutoff_powerlaw(double E                , double g, double Ecutoff,double Emin_ref, double Emax_ref, double flux_ref)
{
  // computes flux at an energy       given the flux over a reference range.
  // for a cutoff powerlaw E^-g exp(-E/Ecutoff)
  // using: integral x_to_inf [t^-g e^-ct] = c^(g-1) integral cx_to_inf [t^-g e^-t] = c^(g-1) \Gamma(a,cx) with -g=a-1 -> a=1-g
  // From gsl manual:
  // double gsl_sf_gamma_inc (double a, double x)
  // is the unnormalized incomplete Gamma Function \Gamma(a,x) = \int_x^\infty dt t^{a-1} \exp(-t) for a real and x >= 0. 

  // tests in test_gsl_gamma_function.cc

  double a=1.0-g;
  double c=1.0/Ecutoff;
  
  double x3=Emin_ref;
  double x4=Emax_ref;

  if(E>10*Ecutoff){double result=0; return result;} // to avoid gsl underflow
 

  // normalization to reference range
  double s_ref =  pow(c,(g-1)) *gsl_sf_gamma_inc(a,c*x3)
                - pow(c,(g-1)) *gsl_sf_gamma_inc(a,c*x4);

  
  double result = flux_ref / s_ref * pow(E,-g) * exp(-E/Ecutoff);

  return result;
 
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

double SourcePopulation::flux_cutoff_powerlaw(double Emin, double Emax, double g, double Ecutoff,double Emin_ref, double Emax_ref, double flux_ref)
{
  // computes flux in an energy range given the flux over a reference range.
  // for a cutoff powerlaw E^-g exp(-E/Ecutoff)
  // using: integral x_to_inf [t^-g e^-ct] = c^(g-1) integral cx_to_inf [t^-g e^-t] = c^(g-1) \Gamma(a,cx) with -g=a-1 -> a=1-g
  // From gsl manual:
  // double gsl_sf_gamma_inc (double a, double x)
  // is the unnormalized incomplete Gamma Function \Gamma(a,x) = \int_x^\infty dt t^{a-1} \exp(-t) for a real and x >= 0. 

  // tests in test_gsl_gamma_function.cc

  if(Emax>10*Ecutoff){double result=0; return result;} // to avoid gsl underflow

  double a=1.0-g;
  double c=1.0/Ecutoff;
  double x1=Emin;
  double x2=Emax;
  double x3=Emin_ref;
  double x4=Emax_ref;

  // integral over required range
  double s     =  pow(c,(g-1)) *gsl_sf_gamma_inc(a,c*x1)
                - pow(c,(g-1)) *gsl_sf_gamma_inc(a,c*x2);

  // normalization to reference range
  double s_ref =  pow(c,(g-1)) *gsl_sf_gamma_inc(a,c*x3)
                - pow(c,(g-1)) *gsl_sf_gamma_inc(a,c*x4);

  
  double result = flux_ref * s/s_ref;

  return result;
 
}