#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_integration.h>

#define H0 70
#define omegav 0.7
#define omegam 0.3
#define omegak 0.0
#define zmax 100
#define zmin 0
#define conv_factor 1
#define c 1     

/************************** List of FUNCTIONS *********************************************/

double Hfunc(double z)
          {  
           double H;
	   H=H0*pow((omegam*pow((1+z),3)+omegak*pow((1+z),2)+omegav),0.5);
           /*printf("hubble test passed succesfully: %g\n",H);*/
           return H;
          }
  
double SFRfunc(double z, int SFR_index)
      {  
      double SFR;
      if (SFR_index == 1)
              {
             SFR=0.32*(exp(3.4*z)/(exp(3.8*z)+45))*(pow((1+z),1.5)*Hfunc(z));
              }
      else if (SFR_index==2)
              {
             SFR=0.16*(exp(3.4*z)/(exp(3.4*z)+22))*(pow((1+z),1.5)*Hfunc(z));
              }
      else if (SFR_index==3)
              {
             SFR=0.22*(exp(3.05*z-0.4)/(exp(2.93*z)+15))*(pow((1+z),1.5)*Hfunc(z));
              }
      else if (SFR_index==4) 
              { 
                if ((z>=0) && (z<=1.04))
                   {
		   SFR=0.015*pow((1+z),3.28)*pow((0.3*pow((1+z),3)+0.7),0.5)*Hfunc(z);
                   }
                else if ((z>=1.04) && (z<=4.48))
                   {
                   SFR=0.188*pow((1+z),-0.26)*pow((0.3*pow((1+z),3)+0.7),0.5)*Hfunc(z);
                   }
                else /*if ((z>=4.48) && (z<=6))*/
                   {
                   SFR=97.72*pow((1+z),-8.0)*pow((0.3*pow((1+z),3)+0.7),0.5)*Hfunc(z);
                   }                          
              }

      else if (SFR_index==5)
              {
             SFR=(0.2331*(pow((1+z),3.7))/pow(1+0.075*pow(1+z,3.7),1.84))*Hfunc(z);
              } 
      else if (SFR_index==6)
              {
            SFR=(0.0103+0.088*z)/((1+pow(z/2.4,2.8))*pow(0.3*pow((1+z),3)+0.7,0.5))*Hfunc(z);
              } 
      else  if (SFR_index==7)
              {
             SFR=(0.014+0.11*z)/((1+pow(z/1.4,2.2))*pow(0.3*pow((1+z),3)+0.7,0.5))*Hfunc(z);
              } 
      else 
              { 
                if ((z>=0) && (z<=0.993))
                   {
                   SFR=0.02*pow(1+z,3.30)*pow(0.3*pow((1+z),3)+0.7,0.5)*Hfunc(z);
                   }
                else if ((z>=0.993) && (z<=3.80))
                   {
                   SFR=0.19*pow(1+z,-0.0549)*pow(0.3*pow((1+z),3)+0.7,0.5)*Hfunc(z);
                   }
                else /*if ((z>=3.80) && (z<=10))*/
                   {
                   SFR=223.9*pow(1+z,-4.46)*pow(0.3*pow((1+z),3)+0.7,0.5)*Hfunc(z);
                   }
                /*else break*/                        
              }
              printf("SFR test passed succesfully: %g\n",SFR);       
   return SFR;
      } 


double inv_H(double z, void *params_ptr) 
      {
      double params = *(double*)params_ptr;
      /*printf("H_inv test passed succesfully: %g\n", params*pow(Hfunc(z), -1.0)); */
      double inv_H = params*pow(Hfunc(z), -1.0);  
      return inv_H;
      }


double Dlfunc(double z)
       { 
        double Dl;
        double int_invH, err; /*declare the ouptput of the integration and its error*/
        double params = 1.0;
        gsl_integration_workspace* w=gsl_integration_workspace_alloc(10000);
	/*gsl_function F = {&H_inv,0};*/ 
        gsl_function F;
	F.function = &inv_H;
	F.params = &params;
        /*printf("integration_before_test passed succesfully")*/
        gsl_integration_qag (&F, zmin, z, 0, 1e-7, 50, 6, w, &int_invH, &err);
        Dl=c*(1+z)*int_invH ;
        printf("Dl test passed : %g\n", Dl);
        gsl_integration_workspace_free (w);/*clean workspace for the second call of Dlfunc*/
        return Dl;
       }


double dvdzfunc(double z)
       { 
        double dvdz;
        dvdz=4*M_PI*c*pow(Dlfunc(z),2)/(pow(1+z,2)*Hfunc(z));
        return dvdz;
       }


double p_unor( double z, void *SFR_index_ptr) 
{
  int SFR_index = *(int*)SFR_index_ptr;
  return conv_factor*SFRfunc(z, SFR_index)*dvdzfunc(z)*pow(1+z,-1);
}


double p_simfunc(double z, int SFR_index)
       {
        double p_sim ;
        double int_p_unor, err;
        /*normalize p_unor*/
        gsl_integration_workspace*w=gsl_integration_workspace_alloc(10000);
        gsl_integration_workspace_free(w);
	gsl_function F = { &p_unor, &SFR_index };
        gsl_integration_qag (&F, zmin, zmax, 0, 1e-7, 50, 6,w, &int_p_unor, &err);
        p_sim=p_unor(z,&SFR_index)/int_p_unor;
        gsl_integration_workspace_free (w);
        return p_sim;
       }
       
/************************************* MAIN ***************************************************/
int main (void)
{
 /********************************* simulated REDSHIFT *************************************/
   int SFR_index=1;       
   int i,i_sim=1; 
   int i_max=10000;
   double z_rand, z_test, z_accept;
   printf("program started");
   const gsl_rng_type * T;  
   gsl_rng * r;             
   gsl_rng_env_setup();     
   T = gsl_rng_ranlxd2;     
   r = gsl_rng_alloc (T);  
   gsl_vector *u=gsl_vector_alloc (10000);  
   gsl_vector *g=gsl_vector_alloc (10000);                                                     
   printf("generator declared");
   
   for(i=0;i<i_max;i++)
     {     /*use Rejection method*/
       z_test=gsl_rng_uniform(r);  
       z_rand=gsl_rng_uniform(r)*zmax;  
       if (z_test > p_simfunc(z_rand, SFR_index))  
       	 break;
       printf("rejection 1 test passed succesfully");
       z_accept = z_rand;  
       gsl_vector_set(g,i,z_rand); 
       gsl_vector_set(u,i_sim,z_accept); 
       i_sim = i_sim +1;
     }
   for(i=0;i<i_max;i++) /*display simulated redshift distribution before  rejection*/
     {
       printf("%f \n", gsl_vector_get(g,i));
     }
   for(i=0;i<i_sim;i++) /*display simulated redshift distribution after rejection*/
     {
       printf("%f \n", gsl_vector_get(u,i));
     }
     
   gsl_rng_free(r);  /*cleanup the generator workspace, and the input vectors*/
   gsl_vector_free(u);
   gsl_vector_free(g);
   return 0 ;
}