/*  Bulirsch-Stoer integrator  

    (Unfortunately this is nearly a direct translation of a 
     FORTRAN code.  The original FORTRAN was ably written by
     Michael Henon. Jack Wisdom)

    must define EQUATIONS_OF_MOTION, MAX_NUMBER_EQUATIONS,
    INTEGRATION_EPSILON, number_equations,
    and possibly NAME in calling routine

*/

#define SUCCESS 1
#define FAILURE 0

#ifndef EQUATIONS_OF_MOTION
#define EQUATIONS_OF_MOTION(state, deriv) \
  equations_of_motion( &(state), &(deriv)); 
#endif

int integrate(state,official_delta_t)
     State *state;
     double *official_delta_t;

{
  State local_state; //creates local versions of state to manipulate
  State state_deriv;
  static int lt_BS[10] = {1,2,3,4,6,8,12,16,24,32};
  static double temp_BS[MAX_NUMBER_EQUATIONS][12];
  double d_BS[6], delta_t_BS;

  int i_BS, m_BS, m1_BS, k_BS, i1max_BS, i1_BS;
  double xa_BS, xb_BS, varm_BS, h_BS, hd_BS, eta2_BS, dta_BS;
  double yb_BS, c_BS, b1_BS, den_BS, dtn_BS, b_BS, var_BS, varma_BS;
  double foo_BS, flt_BS;

/* reset state variables--local_state set equal to state */
  local_state.t = state->t; 
  for(i_BS=0;i_BS<number_equations;i_BS++)
    local_state.x[i_BS] = state->x[i_BS];
  delta_t_BS = *official_delta_t;

  xa_BS = local_state.t;

  EQUATIONS_OF_MOTION(local_state,state_deriv)
 
  for(i_BS=0; i_BS<number_equations; i_BS++)
    {
      temp_BS[i_BS][1] = ABS(local_state.x[i_BS]);
      if(temp_BS[i_BS][1] < INTEGRATION_EPSILON)
	temp_BS[i_BS][1] = INTEGRATION_EPSILON;
      temp_BS[i_BS][4] = state_deriv.x[i_BS];
      temp_BS[i_BS][0] = local_state.x[i_BS];
    }

 outside:
  xb_BS = delta_t_BS + xa_BS;
  for(m_BS=0; m_BS<10; m_BS++)
    {
      flt_BS = lt_BS[m_BS];
      varm_BS = 0.0;
      m1_BS = MIN(m_BS, 6);
      if(m1_BS != 0)
	  for(k_BS=0; k_BS<m1_BS; k_BS++)
	    {
	      d_BS[k_BS] = flt_BS/lt_BS[m_BS-k_BS-1];
	      d_BS[k_BS] *= d_BS[k_BS];
	    }
      h_BS = delta_t_BS / flt_BS;
      hd_BS = 0.5 * h_BS;
      for(i_BS=0; i_BS<number_equations; i_BS++)
	{
	  temp_BS[i_BS][3] = temp_BS[i_BS][0];
	  local_state.x[i_BS] = temp_BS[i_BS][0] + hd_BS*temp_BS[i_BS][4];
	}
      i1max_BS = 2*flt_BS - 1;
      local_state.t = xa_BS;
      for(i1_BS=0; i1_BS<i1max_BS; i1_BS++)
	{
	  local_state.t += hd_BS;

  	  EQUATIONS_OF_MOTION(local_state,state_deriv)

	  for(i_BS=0; i_BS<number_equations; i_BS++)
	    {
	      foo_BS = local_state.x[i_BS];
	      temp_BS[i_BS][1] = MAX(temp_BS[i_BS][1], ABS(foo_BS));
	      if((temp_BS[i_BS][1]) == (0.0))
		{
		  printf("difsys: temp_BS[i_BS][1] = 0\n");
		  return(FAILURE);
		}
	      eta2_BS = temp_BS[i_BS][3] + h_BS*state_deriv.x[i_BS];
	      temp_BS[i_BS][3] = local_state.x[i_BS];
	      local_state.x[i_BS] = eta2_BS;
	    }
	}

      local_state.t = xb_BS;
      EQUATIONS_OF_MOTION(local_state,state_deriv)

      for(i_BS=0; i_BS<number_equations; i_BS++)
	{
	  dta_BS = temp_BS[i_BS][11];
	  yb_BS = (temp_BS[i_BS][3] + local_state.x[i_BS] 
		   + hd_BS*state_deriv.x[i_BS]) / 2.0;
	  c_BS = yb_BS;
	  temp_BS[i_BS][11] = yb_BS;
	  if(m1_BS != 0)
	    {
	      for(k_BS=0; k_BS<m1_BS; k_BS++)
		{
		  b1_BS = d_BS[k_BS] * dta_BS;
		  den_BS = b1_BS - c_BS;
		  dtn_BS = dta_BS;
		  if(den_BS != 0)
		    {
		      b_BS = (c_BS - dta_BS) / den_BS;
		      dtn_BS = c_BS * b_BS;
		      c_BS = b1_BS * b_BS;
		    }
		  dta_BS = temp_BS[i_BS][11-k_BS-1];
		  temp_BS[i_BS][11-k_BS-1] = dtn_BS;
		  yb_BS += dtn_BS;
		}
	      var_BS = ABS(temp_BS[i_BS][2] - yb_BS) / temp_BS[i_BS][1];
	      varm_BS = MAX(varm_BS, var_BS);
	    }
	  temp_BS[i_BS][2] = yb_BS;
	}
      
      if(m_BS < 3)
	varma_BS = varm_BS; 
      else if(varm_BS <= INTEGRATION_EPSILON) /* check for convergence or divergence */
	goto end;
      else if(varm_BS >= varma_BS)
	break;
    }
  
  delta_t_BS /= 2.0;
  goto outside;
  
 end: local_state.t = xb_BS;
  for(i_BS=0; i_BS<number_equations; i_BS++)
    local_state.x[i_BS] = temp_BS[i_BS][2];

/* reset state variables */
  state->t = local_state.t; 
  for(i_BS=0;i_BS<number_equations;i_BS++)
    state->x[i_BS] = local_state.x[i_BS];
  *official_delta_t = delta_t_BS * 1.5 * pow(0.6, (double)(m_BS-m1_BS));

  return(SUCCESS);
}


