#ifndef NeutronMultiFunction_H
#define NeutronMultiFunction_H

#include<cmath>

class NeutronMulti{


 public:
  NeutronMulti(double the_a1, double the_a2,double the_a3,
	       double the_t2, double the_t3){
    a1=the_a1;
    a2=the_a2;
    a3=the_a3;
    t2=the_t2;
    t3=the_t3;
  }

  // Define functions to get memeber values.
  // (Getter Functions)
     double Ga1(){return a1;}
     double Ga2(){return a2;}
     double Ga3(){return a3;}
     double Gt2(){return t2;}
     double Gt3(){return t3;}
     char nan[5];
     
     double operator()(double x) const{
       double b1 =0;
       double b2 = -1.27580114;
       double b3 =  1.363936055;
       double c1 = 0.481385;
       double c2 = 0.5219742;
       double c3 = 1.0387373;
       double t1 =  1.777;
       double x1 = a1*std::sqrt(x)/(2.*std::pow((3.14*t1),(3./2.)));
       double x2 = a2*std::sqrt(x)/(2.*std::pow((3.14*t2),(3./2.)));
       double x3 = a3*std::sqrt(x)/(2.*std::pow((3.14*t3),(3./2.)));
       double x4 = exp(-(x - b1*std::sqrt(x) + c1)/t1);
       double x5 = exp(-(x - b2*std::sqrt(x) + c2)/t2);
       double x6 = exp(-(x - b3*std::sqrt(x) + c3)/t3);
       double y = x1*x4 + x2*x5 + x3*x6;
       
       if(y ==!nan) 
       return y;
     }
  private:
     // Indics 1 -CN , 2-First Fragment , 3- Second Fragment
     // t = temp of nuclei
     double a1,a2,a3;
     double t1,t2,t3;
};

#endif 



 int main(){
   NeutronMulti myexample(0.00100954,0.00100954,1,0.0000000,2.01972);
  cout<<myexample(0.0)<<endl;


}