#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <float.h>
#include "nrutil.h"
#include "nr.h"


#define C 3.00e8
#define NMAX 5000
#define DIMIN 1e-6

/* Calculate cosmological distance in Lambda-Universe */
/* Ver 0.91 */
/* 2003/5/27 K. Motohara */
/* kmotohara@ioa.s.u-tokyo.ac.jp */

/* Taken from Carroll et al. (1992) ARA&A 30, 499 */

/* 2003/5/27 Modified problem in if-clause of double */
/* 2003/5/27 Modified to use numerical-recipies-in-C */
/* 2003/5/27 Fixed Fatal Bug in function sinhh */

double omegak, omegam, omegal;

float sinn(float x, float a);
float function(float z);
float nintegral(float a, float b, int n);
float limitintegral(float a, float b);
int main(int argc,char *argv[])
{
  float h0,z,d;
  float factor;
  float integral;

  if(argc != 6){
    fprintf(stderr,"Usage: %s (P[roper]/L[uminosity]/D[iameter]) H0(km/s/Mpc) Omega_M Omega_lambda z\n",argv[0]);
    exit(1);
  }

  h0=atof(argv[2])*1e3/3.0856776e22;
  omegam=atof(argv[3]);
  omegal=atof(argv[4]);
  omegak=1.0-omegam-omegal;
  z=atof(argv[5]);

  fprintf(stderr,"H0 %e\n",h0);
  fprintf(stderr,"Omegam %e\n",omegam);
  fprintf(stderr,"Omegal %e\n",omegal);
  fprintf(stderr,"Omegak %e\n",omegak);
  if(strcmp(argv[1],"P")==0) factor=1.0;
  else if(strcmp(argv[1],"L")==0) factor=1.0+z;
  else if(strcmp(argv[1],"D")==0) factor=1.0/(1.0+z);
  else {
    fprintf(stderr,"Usage: %s (P[roper]/L[uminosity]/D[iameter]) H0(km/s/Mpc) Omega_M Omega_lambda z\n",argv[0]);
    exit(1);
  }
    
  integral=qromb(&(function), 0, z);

  if(fabs(omegak)<1e-9){
    /*    fprintf(stderr,"omegak=%f\n",omegak); */
    d=integral/h0*factor*C;
    fprintf(stdout,"%e\n",d); /* output in (m) */
  }
  else{
    integral=integral*sqrt(fabs(omegak));
    /*    d=sinn(integral,omegak)/h0/sqrt(fabs(omegak)); */
    d=sinn(integral,omegak)/sqrt(fabs(omegak))/h0*factor*C; 
    fprintf(stdout,"%e\n",d);
  }
}

float sinn(float x, float a)
{
  if(a>0.0) return sinh(x);
  else return sin(x);
}

float function(float z)
{
  return 1/sqrt((1.0+z)*(1.0+z)*(1.0+omegam*z)-z*(2.0+z)*omegal);
}