GLOBALS_SECTION
 #include <df1b2fun.h>
 
  dvariable lambda(const dvariable& t, const dvariable& a0, const dvariable& a1, 
                   const dvariable& a2, const dvariable& b1, const dvariable& b2){
	dvariable f = 0.5*a0 + a1*cos(M_PI*t/12.0) + a2*cos(M_PI*t/6.0) + b1*sin(M_PI*t/12.0) + b2*sin(M_PI*t/6.0);
    return f*f;
  }

  dvariable sigmoid(const dvariable& x, const dvariable& a, const dvariable& b, const dvariable& catchTime, const dvariable& rho){
    return (1.0/(1.0+exp(a*(-x+(catchTime-b/rho))))) / (1.0/(1.0+exp(a*(-catchTime+(catchTime-b/rho))))) ; 
  }

  dvariable romberg(const dvariable& a ,const dvariable& b, const int n, const int p, const dvariable& a0, const dvariable& a1, const dvariable& a2, const dvariable& b1, const dvariable& b2, const dvariable& alpha, const dvariable& beta, const dvariable& rho){
	dvariable  two   = 2.0;
	dvar_vector r(0,n);
	
	r[0]  = ( lambda(a,a0,a1,a2,b1,b2)*sigmoid(a,alpha,beta,b,rho) + lambda(b,a0,a1,a2,b1,b2)*sigmoid(b,alpha,beta,b,rho) ) * (b - a) / two; 
		
	cout << " r[0] = (" << lambda(a,a0,a1,a2,b1,b2) << " * " << sigmoid(a,alpha,beta,b,rho) << " + " << lambda(b,a0,a1,a2,b1,b2) << " * " << sigmoid(b,alpha,beta,b,rho) << ")* " << "((" << b << " - " << a << ") / " << two << ") = " <<  r[0] << endl;
	
 	cout << " r[0] = (" << lambda(a,a0,a1,a2,b1,b2) << " * " << sigmoid(a,alpha,beta,b,rho) << " + " << lambda(b,a0,a1,a2,b1,b2) << " * " << sigmoid(b,alpha,beta,b,rho) << ")* " << ((b - a) / two) << " = " <<  r[0] << endl;
	 

        return r[0];
  }

DATA_SECTION
  number rho;
  number catchTime;

PARAMETER_SECTION
  init_bounded_number a0(0.00001,1.0);
  init_number a1;
  init_number a2;
  init_number b1;
  init_number b2;
  init_bounded_number logalpha(-3.0,-2.0);
  init_bounded_number logbeta(1.0,3.0);
  sdreport_number alpha;
  sdreport_number beta;
  objective_function_value loglik
 
PRELIMINARY_CALCS_SECTION
  a0=0.0296148776575; 
  logbeta = 2.23636105363; 
  logalpha = -2.58907688300;
  a1 = -0.00519618644045;
  a2 = -0.0324360688600;
  b1 = -0.146826760910;
  b2 = -0.0520084376292;
  rho = 0.3048282;
  catchTime = 2.5;
 
PROCEDURE_SECTION
  loglik=0.0;
  alpha = exp(logalpha);
  beta = exp(logbeta);

  dvariable intFrom = (log(0.99/0.01) - alpha*( catchTime - beta/rho ) )/-alpha;
  dvariable integral = romberg(intFrom,catchTime,14,2,a0,a1,a2,b1,b2,alpha,beta,rho);
  exit(1);