/*
  Standalone program to estimate the sensitivity of a pulsar timing
  array to a gravitational wave background.

  Based on Verbiest et al., MNRAS 2009 and arXiv:0906.4246 (Ph.D. Thesis)

  Joris P.W. Verbiest, MPIfR, 12 Feb 2010
 */
/*
  26 April 2010: Added in self-noise as explained by KJ Lee.
 */
// g++ -o PTAS PTASens.C -Wall
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <string.h>
#include <vector>
#include <math.h>

using namespace std;

#define MAX_NPSR 60

int verb = 0;

void GetCoordinates(char *Name,long double *dec,long double *RA);
long double GetHD(char *Name1,char *Name2);
void help();

int main(int argc,char *argv[]){
  int ct; // generic counter
  char infile[100]="Standard.dat";
  char outfile[100]="Sensitivity.dat";
  long double alpha = -2.0/3.0;
  long double minAmp=1e-16,maxAmp=1e-11;
  int Namp=300;

  // Get Essential Data from the Command Line
  for(ct = 1;ct<argc;ct++){
    if(strcmp(argv[ct],"-file")==0){
      strcpy(infile,argv[++ct]);
    }else if(strcmp(argv[ct],"-alpha")==0){
      sscanf(argv[++ct],"%Lg",&alpha);
    }else if(strcmp(argv[ct],"-minAmp")==0){
      sscanf(argv[++ct],"%Lg",&minAmp);
    }else if(strcmp(argv[ct],"-maxAmp")==0){
      sscanf(argv[++ct],"%Lg",&maxAmp);
    }else if(strcmp(argv[ct],"-namp")==0){
      sscanf(argv[++ct],"%d",&Namp);
    }else if(strcmp(argv[ct],"-outfile")==0){
      strcpy(outfile,argv[++ct]);
    }else if(strcmp(argv[ct],"-h")==0){
      help();
      exit(0);
    }
  }

  char PulsarNames[MAX_NPSR][20];
  long double RMS[MAX_NPSR];
  long double Tspan[MAX_NPSR];
  int Nobs[MAX_NPSR];
  int Npsr = 0;
  FILE *fin;
  fin = fopen(infile,"r");
  if(fin==NULL){
    cout << "Problem opening file. Will output help information and "
         << "quit.\n";
    help();
    exit(1);
  }
  // Read Essential Data from Input File
  while (!feof(fin)){
    fscanf(fin,"%s %Lg %Lg %d\n",PulsarNames[Npsr],&RMS[Npsr],&Tspan[Npsr],
           &Nobs[Npsr]);
    RMS[Npsr]/=(365.25*24.0*3600.0);
    Npsr++;
  }
  fclose(fin);

  long double FNyq[Npsr];
  long double FCorn[Npsr];

  // Determine the Nyquist Frequency for each pulsar
  for(ct = 0;ct<Npsr;ct++){
    FNyq[ct] = (long double)Nobs[ct]/(2.0*Tspan[ct]);
  }

  int pct1,pct2;
  long double Amplitude = 1e-20;
  long double LogAmp = log10(Amplitude);
  int Nfree;
  long double minLength;
  long double Freq;
  long double SigNSq[Npsr];
  long double SigGSq[Npsr];
  long double SigFrac[Npsr];
  long double Xi;
  vector <long double> Sens;
  int sensCt=0;
  long double AmpIncrease = (log10(maxAmp)-log10(minAmp))/(long double)(Namp-1);
  // Loop over the GWB Amplitudes
  for(LogAmp=log10(minAmp);LogAmp<log10(maxAmp);LogAmp += AmpIncrease){
    Amplitude = pow(10.0,LogAmp);
    // Loop over the pulsars
    for(ct=0;ct<Npsr;ct++){
      // Determine the Corner frequency
      FCorn[ct] = pow(24.0*Tspan[ct]/(long double)(Nobs[ct])
                      *pow(M_PI*RMS[ct]/Amplitude,2.0),1.0/(2.0*alpha-3.0));
      if(FCorn[ct]>FNyq[ct]) 
        FCorn[ct] = FNyq[ct];

      SigNSq[ct] = 0.0;
      SigGSq[ct] = 0.0;
      // Determine the integrated Noise and GW power spectra:
      for(Freq = 1.0/Tspan[ct];Freq<=FNyq[ct];Freq+=1.0/Tspan[ct]){
        //SigNSq[ct] += pow(Freq/FCorn[ct],2.0*alpha-3.0)
        //  /pow(1+pow(Freq/FCorn[ct],2.0*alpha-3.0),2.0);
        SigNSq[ct] += pow(Freq/FCorn[ct],2.0*alpha-3.0)
          /pow(1.0+pow(Freq/FCorn[ct],2.0*alpha-3.0),2.0)+
          pow(Freq/FCorn[ct],4.0*alpha-6.0)
          /pow(1.0+pow(Freq/FCorn[ct],2.0*alpha-3.0),2.0);
        //SigGSq[ct] += pow(Freq/FCorn[ct],4.0*alpha-6.0)
        //  /pow(1+pow(Freq/FCorn[ct],2.0*alpha-3.0),2.0);
        SigGSq[ct] += pow(Freq/FCorn[ct],4.0*alpha-6.0)
          /pow(1.0+pow(Freq/FCorn[ct],2.0*alpha-3.0),2.0);
      }
      // Only the fraction of these comes up in the Sensitivity Equation:
      SigFrac[ct] = SigNSq[ct]/SigGSq[ct];
    }

    // Now loop over the pulsar pairs to add up the sensitivity:
    Sens.push_back(0.0);
    for(pct1=0;pct1<Npsr-1;pct1++){
      for(pct2=pct1+1;pct2<Npsr;pct2++){
        Xi = GetHD(PulsarNames[pct1],PulsarNames[pct2]);
        minLength = Tspan[pct1];
        if(Tspan[pct2]<Tspan[pct1])
          minLength = Tspan[pct2];
        Nfree = (int)(2.0*minLength*sqrt(FCorn[pct1]*FCorn[pct2]));
        if(Nfree<1)
          Nfree = 1;
        Sens[sensCt] += pow(Xi,2.0)*Nfree
          /(1.0+pow(Xi,2.0)+SigFrac[pct1]+SigFrac[pct2]
            +SigFrac[pct1]*SigFrac[pct2]);
      }
    }
    Sens[sensCt] = sqrt(Sens[sensCt]);
    sensCt++;
  }

  // And output the result:
  FILE *fout;
  fout = fopen(outfile,"w");
  ct = 0;
  for(LogAmp=log10(minAmp);LogAmp<log10(maxAmp);LogAmp+=AmpIncrease){
    Amplitude = pow(10.0,LogAmp);
    fprintf(fout,"%Lg\t%Lg\n",Amplitude,Sens[ct]);
    ct++;
  }
  fclose(fout);
  return(0);
}

long double GetHD(char *Name1,char *Name2){
  long double dec1,dec2,RA1,RA2;
 
  GetCoordinates(Name1,&dec1,&RA1);
  GetCoordinates(Name2,&dec2,&RA2);

  long double angle = acos(cos(dec1)*cos(dec2)*cos(RA1-RA2)+
                      sin(dec1)*sin(dec2));
  long double xx = 0.5*(1.0-cos(angle));
  long double Xi = 1.5*xx*log(xx)-xx/4.0+0.5;

  return(Xi);
}

void GetCoordinates(char *Name,long double *dec,long double *RA){
  int length=strlen(Name);
  int xtra=0;
  if(Name[0]=='J' || Name[0]=='B')
    xtra = 1;

  long double hours = 0.0;
  long double temp = 0.0;
  *RA = 0.0L;

  // Determine the right ascension
  temp = (long double)(Name[xtra]-'0');
  hours += 10.0*temp;
  temp = (long double)(Name[++xtra]-'0');
  hours += temp;

  temp = (long double)(Name[++xtra]-'0');
  hours += temp/6.0;
  temp = (long double)(Name[++xtra]-'0');
  hours += temp/60.0;
  *RA = hours/12.0*M_PI;
  if(verb>20)
    cout << "Right ascension: " << *RA << " hours: " << hours << endl;

  // Declination:
  int PosNeg = 1;
  if(Name[++xtra]=='-')
    PosNeg = -1;
  else if(Name[xtra]!='+')
    cout << "Don't understand declination sign in Name: " 
         << Name << endl;
  
  *dec = 0.0;
  long double degrees = 0.0;
  temp = (long double)(Name[++xtra]-'0');
  degrees += 10.0*temp;
  temp = (long double)(Name[++xtra]-'0');
  degrees += temp;
  
  if(length>8){
    temp = (long double)(Name[++xtra]-'0');
    degrees += temp/6.0;
    temp = (long double)(Name[++xtra]-'0');
    degrees += temp/60.0;
  }
  degrees *= (long double)PosNeg;
  *dec = degrees/180.0*M_PI;
  if(verb>20)
    cout << "Declination: " << *dec << " degrees: " << degrees << endl;
}  

void help(){
  cout << " _______________________________________________________________\n";
  cout << "/PTAS - a standalone program to estimate PTA sensitivity       \\\n"
       << "|            curves.                                            |\n";
  cout << "|                                                               |\n";
  cout << "| USAGE and defaults:                                           |\n"
       << "|      PTAS -h -file Standard.dat -alpha -0.66666 -minAmp 1e-16 |\n"
       << "|                -maxAmp 1e-11 -namp 300                        |\n"
       << "|                -outfile Sensitivity.dat                       |\n"
       << "|                                                               |\n";
  cout << "| OPTIONS:                                                      |\n"
       << "|      -h      : this help information                          |\n"
       << "|      -file   : input filename (see notes below)               |\n"
       << "|      -outfile: output filename (see notes below)              |\n"
       << "|      -alpha  : GWB spectral index (default -2/3)              |\n"
       << "|      -minAmp : minimum GWB amplitude                          |\n"
       << "|      -maxAmp : maximum GWB amplitude                          |\n"
       << "|      -namp   : number of GWB amplitudes                       |\n";
  cout << "| NOTES:                                                        |\n"
       << "|       -file: The input file must contain a four-column line   |\n"
       << "|              for each pulsar of the PTA. The four columns must|\n"
       << "|              contain respectively:                            |\n"
       << "|              1: the pulsar name                               |\n"
       << "|              2: the timing rms in seconds                     |\n"
       << "|              3: the data span in years                        |\n"
       << "|              4: the number of TOAs.                           |\n"
       << "|                                                               |\n"
       << "|       -outfile: The outfile will contain two columns:         |\n"
       << "|                 1: the GWB amplitude                          |\n"
       << "|                 2: the PTA sensitivity.                       |\n";
  cout << "| Have a nice day!                                              |\n"
       << "\\_______________________________________________"
       <<"________________/\n";
  }