#include "sopnamsp.h"
#include "machdefs.h"
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <unistd.h>

#include "constcosmo.h"
#include "cosmocalc.h"
#include "geneutils.h"
#include "integfunc.h"
#include "schechter.h"
#include "planckspectra.h"

inline double rad2deg(double trad) {return trad/M_PI*180.;}
inline double rad2min(double trad) {return trad/M_PI*180.*60.;}
inline double rad2sec(double trad) {return trad/M_PI*180.*3600.;}
inline double deg2rad(double tdeg) {return tdeg*M_PI/180.;}
inline double min2rad(double tmin) {return tmin*M_PI/(180.*60.);}
inline double sec2rad(double tsec) {return tsec*M_PI/(180.*3600.);}

void usage(void);
void usage(void) {
  cout<<"cmvdefsurv [-r] -x adtx,atxlarg [-y adty,atylarg] -z dred,redlarg redshift"<<endl
      <<" -x adtx,atxlarg : resolution en Theta_x (arcmin), largeur (degre)"<<endl
      <<" -y adty,atylarg : idem selon y, si <=0 meme que x"<<endl
      <<" -z dred,redlarg : resolution en redshift, largeur en redshift"<<endl
      <<" -P : on donne -x -y -z en Mpc aulieu d\'angles et de redshift"<<endl
      <<" -O surf,tobs : surface effective (m^2) et temps d\'observation (s)"<<endl
      <<" -A Tbruit : temperature de bruit (K)"<<endl
      <<" -S Tsynch,indnu : temperature (K) synch a 408 Mhz, index d\'evolution"<<endl
      <<"                   (indnu==0 no evolution with freq.)"<<endl
      <<" -M  : masse de HI de reference (MSol), si <=0 mean schechter"<<endl
      <<" -F  : HI flux factor to be applied for our redshift"<<endl
      <<" -1 : on ne mesure qu\'une polarisation"<<endl
      <<" redshift : redshift moyen du survey"<<endl
      <<endl;
}

int main(int narg,char *arg[])
{
 // --- Valeurs fixes
 // WMAP
 unsigned short flat = 0;
 double h100=0.71, om0=0.267804, or0=7.9e-05, ol0=0.73,w0=-1.;
 cout<<"\nh100="<<h100<<" Om0="<<om0<<" Or0="<<or0<<" Or0="<<or0<<" Ol0="<<ol0<<" w0="<<w0<<endl;
 // Schechter
 double h75 = h100 / 0.75;
 double nstar = 0.006*pow(h75,3.);  //
 double mstar = pow(10.,9.8/(h75*h75));  // MSol
 double alpha = -1.37;
 cout<<"nstar= "<<nstar<<"  mstar="<<mstar<<"  alpha="<<alpha<<endl;


 // --- Arguments
 bool inmpc = false;
 double adtx=1., atxlarg=90., adty=-1., atylarg=-1.;
 double dx=1.,txlarg=1000., dy=-1.,tylarg=1000., dz=1.,tzlarg=100.;
 int nx,ny,nz;
 double redshift = 1., dred=0.01, redlarg=0.3;
 double tobs = 3600., surfeff = 400000.;
 double Tbruit=75.;
 // a 408 MHz (Haslam) + evol index a -2.6
 double Tsynch408=60., nuhaslam=0.408, indnu = -2.6;
 double mhiref = -1.; // reference Mass en HI (def integ schechter)
 double hifactor = 1.;
 double facpolar = 1.; // si on ne mesure qu'un polarisation 0.5

 // --- Decodage arguments
 char c;
  while((c = getopt(narg,arg,"hP1x:y:z:A:S:O:M:F:")) != -1) {
  switch (c) {
  case 'P' :
    inmpc = true;
    break;
  case 'x' :
    sscanf(optarg,"%lf,%lf",&adtx,&atxlarg);
    break;
  case 'y' :
    sscanf(optarg,"%lf,%lf",&adty,&atylarg);
    break;
  case 'z' :
    sscanf(optarg,"%lf,%lf",&dred,&redlarg);
    break;
  case 'O' :
    sscanf(optarg,"%lf,%lf",&surfeff,&tobs);
    break;
  case 'A' :
    sscanf(optarg,"%lf",&Tbruit);
    break;
  case 'S' :
    sscanf(optarg,"%lf,%lf",&Tsynch408,&indnu);
    break;
  case 'M' :
    sscanf(optarg,"%lf",&mhiref);
    break;
  case 'F' :
    sscanf(optarg,"%lf",&hifactor);
    break;
  case '1' :
    facpolar = 0.5;
    break;
  case 'h' :
  default :
    usage(); return -1;
  }
 }
 if(optind>=narg) {usage(); return-1;}
 sscanf(arg[optind],"%lf",&redshift);
 if(redshift<=0.) {cout<<"Redshift "<<redshift<<" should be >0"<<endl; return -2;}

 // --- Initialisation de la Cosmologie
 cout<<"\n--- Cosmology for z = "<<redshift<<endl;
 CosmoCalc univ(flat,true,2.*redshift);
 double perc=0.01,dzinc=redshift/100.,dzmax=dzinc*10.; unsigned short glorder=4;
 univ.SetInteg(perc,dzinc,dzmax,glorder);
 univ.SetDynParam(h100,om0,or0,ol0,w0);
 univ.Print(redshift);

 double dang    = univ.Dang(redshift);
 double dtrcom  = univ.Dtrcom(redshift);
 double dlum    = univ.Dlum(redshift);
 double dloscom = univ.Dloscom(redshift);
 double dlosdz  = univ.Dhubble()/univ.E(redshift);
 cout<<"dang="<<dang<<" dlum="<<dlum<<" dtrcom="<<dtrcom
     <<" dloscom="<<dloscom<<" dlosdz="<<dlosdz<<" Mpc"<<endl;

 cout<<"\n1\" -> "<<dang*sec2rad(1.)<<" Mpc = "<<dtrcom*sec2rad(1.)<<" Mpc com"<<endl;
 cout<<"1\' -> "<<dang*min2rad(1.)<<" Mpc = "<<dtrcom*min2rad(1.)<<" Mpc com"<<endl;
 cout<<"1d -> "<<dang*deg2rad(1.)<<" Mpc = "<<dtrcom*deg2rad(1.)<<" Mpc com"<<endl;

 cout<<"dz=1 -> "<<dlosdz<<" Mpc com"<<endl;

 cout<<"1 Mpc los com -> dz = "<<1./dlosdz<<endl;
 cout<<"1 Mpc transv com -> "<<rad2sec(1./dtrcom)<<"\" = "
     <<rad2min(1./dtrcom)<<" \' = "<<rad2deg(1./dtrcom)<<" d"<<endl;

 // --- Mise en forme dans les unites appropriees
 if(adty<=0.) adty=adtx;
 if(atylarg<=0.) atylarg=atxlarg;
 if(inmpc) {
   dx = adtx; txlarg = atxlarg; nx = int(txlarg/dx+0.5);
   dy = adty; txlarg = atxlarg; ny = int(tylarg/dy+0.5);
   dz = dred; tzlarg = redlarg; nz = int(tzlarg/dz+0.5);
   adtx = dx/dtrcom; atxlarg = adtx*nx;
   adty = dy/dtrcom; atylarg = adty*ny;
   dred = dz/dlosdz; redlarg = dred*nz;
 } else {
   adtx = min2rad(adtx); atxlarg = deg2rad(atxlarg); nx = int(atxlarg/adtx+0.5);
   adty = min2rad(adty); atylarg = deg2rad(atylarg); ny = int(atylarg/adty+0.5);
   nz = int(redlarg/dred+0.5);
   dx = adtx*dtrcom; txlarg = dx*nx;
   dy = adty*dtrcom; tylarg = dy*ny;
   dz = dred*dlosdz; tzlarg = dz*nz;
 }
 double Npix = (double)nx*(double)ny*(double)nz;

 double redlim[2] = {redshift-redlarg/2.,redshift+redlarg/2.};
 if(redlim[0]<=0.)
   {cout<<"Lower redshift limit "<<redlim[0]<<" should be >0"<<endl; return -3;}
 double dtrlim[2] = {univ.Dtrcom(redlim[0]),univ.Dtrcom(redlim[1])};
 double loslim[2] = {univ.Dloscom(redlim[0]), univ.Dloscom(redlim[1])};
 double dlumlim[2] = {univ.Dlum(redlim[0]),univ.Dlum(redlim[1])};

 cout<<"\n---- Type de donnees: inmpc = "<<inmpc<<endl;
 cout<<"---- Line of Sight: Redshift = "<<redshift<<endl
     <<"dred = "<<dred<<"  redlarg = "<<redlarg<<endl
     <<"  dz = "<<dz<<" Mpc  redlarg = "<<tzlarg<<" Mpc, nz = "<<nz<<" pix"<<endl;
 cout<<"---- Transverse X:"<<endl
     <<"adtx = "<<rad2min(adtx)<<"\',  atxlarg = "<<rad2deg(atxlarg)<<" d"<<endl
     <<"  dx = "<<dx<<" Mpc,  txlarg = "<<txlarg<<" Mpc, nx = "<<nx<<" pix"<<endl;
 cout<<"---- Transverse Y:"<<endl
     <<"adty = "<<rad2min(adty)<<"\',  atylarg = "<<rad2deg(atylarg)<<" d"<<endl
     <<"  dy = "<<dy<<" Mpc,  tylarg = "<<tylarg<<" Mpc, ny = "<<ny<<" pix"<<endl;
 cout<<"---- Npix total = "<<Npix<<" -> "<<Npix*sizeof(double)/1.e6<<" Mo"<<endl;

 // --- Cosmolographie Transverse
 cout<<"\n--- Transverse"<<endl;
 cout<<"dang comoving = "<<dtrcom<<" Mpc (com) var_in_z ["
                         <<dtrlim[0]<<","<<dtrlim[1]<<"]"<<endl;

 cout<<"... dx = "<<dx<<" Mpc (com), with angle "<<adtx*dtrcom<<endl
     <<"       with angle var_in_z ["<<adtx*dtrlim[0]<<","<<adtx*dtrlim[1]<<"]"<<endl;
 cout<<"... largx = "<<txlarg<<" Mpc (com), with angle "<<atxlarg*dtrcom<<endl
     <<"          with angle var_in_z ["<<atxlarg*dtrlim[0]<<","<<atxlarg*dtrlim[1]<<"]"<<endl;

 cout<<"... dy = "<<dy<<" Mpc (com), with angle "<<adty*dtrcom<<endl
     <<"       with angle var_in_z ["<<adty*dtrlim[0]<<","<<adty*dtrlim[1]<<"]"<<endl;
 cout<<"... largy = "<<tylarg<<" Mpc (com), with angle "<<atylarg*dtrcom<<endl
     <<"          with angle var_in_z ["<<atylarg*dtrlim[0]<<","<<atylarg*dtrlim[1]<<"]"<<endl;

 // --- Cosmolographie Line of sight
 cout<<"\n--- Line of Sight"<<endl;
 cout<<"los comoving distance = "<<dloscom<<" Mpc (com) in ["
     <<loslim[0]<<","<<loslim[1]<<"]"<<endl
     <<"                                    diff = "
     <<loslim[1]-loslim[0]<<" Mpc"<<endl;

 cout<<"...dz = "<<dz<<" Mpc (com), with redshift approx "<<dred*dlosdz<<endl;
 cout<<"...tzlarg = "<<tzlarg<<" Mpc (com), with redshift approx "<<redlarg*dlosdz<<endl;

 // --- Solid Angle & Volume
 cout<<"\n--- Solid angle"<<endl;
 double angsol = AngSol(adtx/2.,adty/2.,M_PI/2.);
 cout<<"Elementary solid angle = "<<angsol<<" sr = "<<angsol/(4.*M_PI)<<" *4Pi sr"<<endl;
 double angsoltot = AngSol(atxlarg/2.,atylarg/2.,M_PI/2.);
 cout<<"Total solid angle = "<<angsoltot<<" sr = "<<angsoltot/(4.*M_PI)<<" *4Pi sr"<<endl;

 cout<<"\n--- Volume"<<endl;
 double dvol = dx*dy*dz;
 cout<<"Pixel volume comoving = "<<dvol<<" Mpc^3"<<endl;
 double vol = univ.Vol4Pi(redlim[0],redlim[1])/(4.*M_PI)*angsoltot;
 cout<<"Volume comoving = "<<vol<<" Mpc^3 = "<<vol/1.e9<<" Gpc^3"<<endl
     <<"Pixel volume comoving = vol/Npix = "<<vol/Npix<<" Mpc^3"<<endl;

 // --- Fourier space: k = omega = 2*Pi*Nu
 cout<<"\n--- Fourier space"<<endl;
 cout<<"Array size: nx = "<<nx<<",  ny = "<<ny<<",  nz = "<<nz<<endl;
 double dk_x = 2.*M_PI/(nx*dx), knyq_x = M_PI/dx;
 double dk_y = 2.*M_PI/(nx*dy), knyq_y = M_PI/dy;
 double dk_z = 2.*M_PI/(nz*dz), knyq_z = M_PI/dz;
 cout<<"Resolution: dk_x = "<<dk_x<<" Mpc^-1  (2Pi/dk_x="<<2.*M_PI/dk_x<<" Mpc)"<<endl
     <<"            dk_y = "<<dk_y<<" Mpc^-1  (2Pi/dk_y="<<2.*M_PI/dk_y<<" Mpc)"<<endl;
 cout<<"Nyquist:    kx = "<<knyq_x<<" Mpc^-1  (2Pi/knyq_x="<<2.*M_PI/knyq_x<<" Mpc)"<<endl
     <<"            ky = "<<knyq_y<<" Mpc^-1  (2Pi/knyq_y="<<2.*M_PI/knyq_y<<" Mpc)"<<endl;
 cout<<"Resolution: dk_z = "<<dk_z<<" Mpc^-1  (2Pi/dk_z="<<2.*M_PI/dk_z<<" Mpc)"<<endl;
 cout<<"Nyquist:    kz = "<<knyq_z<<" Mpc^-1  (2Pi/knyq_z="<<2.*M_PI/knyq_z<<" Mpc)"<<endl;

 // --- Masse de HI
 cout<<"\n--- Mass HI"<<endl;
 Schechter sch(nstar,mstar,alpha);
 sch.SetOutValue(1);
 cout<<"nstar= "<<nstar<<"  mstar="<<mstar<<"  alpha="<<alpha<<endl;
 cout<<"mstar*sch(mstar) = "<<sch(mstar)<<" Msol/Mpc^3/Msol"<<endl;
 int npt = 10000;
 double lnx1=log10(1e-6), lnx2=log10(1e+14), dlnx=(lnx2-lnx1)/npt;
 double masshimpc3 = IntegrateFuncLog(sch,lnx1,lnx2,0.001,dlnx,10.*dlnx,6);
 cout<<"Mass density: "<<masshimpc3<<" Msol/Mpc^3"<<endl;

 double masshipix = masshimpc3*dvol;
 double masshitot = masshimpc3*vol;
 cout<<"Pixel mass = "<<masshipix<<" Msol"<<endl
     <<"Total mass in survey = "<<masshitot<<" Msol"<<endl;
 if(mhiref<=0.)  mhiref = masshipix;

 // --- Survey values
 double unplusz = 1.+redshift;
 double nuhiz = Fr_HyperFin_Par / unplusz;  // GHz
 // dnu = NuHi/(1.+z0-dz/2) - NuHi/(1.+z0+dz/2)
 //     = NuHi*dz/(1.+z0)^2 * 1/[1-(dz/(2*(1+z0)))^2]
 double dnuhiz = Fr_HyperFin_Par *dred/(unplusz*unplusz)
               / (1.- (dred/.2/unplusz)*(dred/.2/unplusz));
 cout<<"\n--- Observation:"<<endl
     <<"      surf_eff="<<surfeff<<" m^2, tobs="<<tobs<<" s"<<endl
     <<"      nu="<<nuhiz<<" GHz, dnu="<<dnuhiz*1.e3<<" Mhz"<<endl;
 cout<<"dang lumi = "<<dlum<<" in ["<<dlumlim[0]<<","<<dlumlim[1]<<"] Mpc"<<endl;

 // --- Power emitted by HI
 cout<<"\n--- Power from HI for M = "<<mhiref<<" Msol at "<<nuhiz<<" GHz"<<endl;
 cout<<"flux factor = "<<hifactor<<" at redshift = "<<redshift<<endl;

 double fhi = hifactor*FluxHI(mhiref,dlum);
 cout<<"FluxHI("<<dlum<<" Mpc) all polar:"<<endl
     <<"  Flux= "<<fhi<<" W/m^2 = "<<fhi/Jansky2Watt_cst<<" Jy.Hz"<<endl
     <<"      in ["<<hifactor*FluxHI(mhiref,dlumlim[0])
     <<","<<hifactor*FluxHI(mhiref,dlumlim[1])<<"] W/m^2"<<endl;
 double sfhi = fhi / (dnuhiz*1.e+9) / Jansky2Watt_cst;
 cout<<"If spread over "<<dnuhiz<<" GHz, flux density = "<<sfhi<<" Jy"<<endl;

 // --- Signal analysis
 cout<<"\n--- Signal analysis"<<endl;
 PlanckSpectra planck(T_CMB_Par);
 planck.SetApprox(1);  // Rayleigh spectra
 planck.SetVar(0); // frequency
 planck.SetUnitOut(0); // output en W/....
 planck.SetTypSpectra(0); // radiance W/m^2/Sr/Hz

 double hnu = h_Planck_Cst*(nuhiz*1.e+9);
 cout<<"Facteur polar = "<<facpolar<<",   h.nu="<<hnu<<" J"<<endl;

 double psig = facpolar * fhi * surfeff;
 double esig = psig * tobs;
 double nsig = esig / hnu;
 double tsig = psig / k_Boltzman_Cst / (dnuhiz*1.e+9);
 double ssig = psig / surfeff / (dnuhiz*1.e+9) / Jansky2Watt_cst;
 cout<<"Signal("<<mhiref<<" Msol): P="<<psig<<" W -> E="<<esig<<" J -> "<<nsig<<" ph"<<endl;
 cout<<"     Antenna temperature: tsig="<<tsig<<" K"<<endl;
 cout<<"     flux density = "<<ssig<<" Jy"<<endl;

 double facnoise = facpolar * tobs * surfeff * angsol * (dnuhiz*1.e+9);

 double tsynch = Tsynch408;
        if(fabs(indnu)>1.e-50) tsynch *= pow(nuhiz/nuhaslam,indnu);
 planck.SetTemperature(tsynch);
 double psynch = facpolar * planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1.e+9);
 double esynch = psynch * tobs;
 double nsynch = esynch / hnu;
 double ssynch = psynch / surfeff / (dnuhiz*1.e+9) / Jansky2Watt_cst;
 cout<<"Synchrotron: T="<<Tsynch408<<" K ("<<nuhaslam<<" GHz), "
     <<tsynch<<" K ("<<nuhiz<<" GHz)"<<endl
     <<"             P="<<psynch<<" W -> E="<<esynch<<" J -> "<<nsynch<<" ph"<<endl;
 cout<<"     flux density = "<<ssynch<<" Jy"<<endl;

 double tcmb = T_CMB_Par;
 planck.SetTemperature(tcmb);
 double pcmb = facpolar * planck(nuhiz*1.e+9) * surfeff * angsol * (dnuhiz*1.e+9);
 double ecmb = pcmb * tobs;
 double ncmb = ecmb / hnu;
 double scmb = pcmb / surfeff / (dnuhiz*1.e+9) / Jansky2Watt_cst;
 cout<<"CMB: T="<<tcmb<<" K -> P="<<pcmb<<" W -> E="<<ecmb<<" J -> "<<ncmb<<" ph"<<endl;
 cout<<"     flux density = "<<scmb<<" Jy"<<endl;

 // --- Noise analysis
 cout<<"\n--- Noise analysis"<<endl;
 double pbruit = k_Boltzman_Cst * Tbruit * (dnuhiz*1.e+9);
 double ebruit = pbruit * tobs;
 cout<<"Bruit: T="<<Tbruit<<" K, P="<<pbruit<<" W -> E="<<ebruit<<" J"<<endl;

 double seq = (1./facpolar) * k_Boltzman_Cst * Tbruit / surfeff /Jansky2Watt_cst;
 cout<<"\nSystem equivalent flux density: "<<seq<<" Jy"<<endl;

 double slim = seq / sqrt(2.*(dnuhiz*1.e+9)*tobs);
 cout<<"Observation flux density limit: "<<slim<<" Jy"<<endl;

 double SsN = ssig/slim;
 cout<<"\nSignal to noise ratio = "<<SsN<<endl;

 double smass = mhiref/ssig*slim;
 cout<<"\nSigma noise equivalent = "<<smass<<" Msol"<<endl;

 return 0;
}