#ifndef __USER_ULIB_FLT_H
#define __USER_ULIB_FLT_H


class Ulib_Flt {
	//
  // eq. of Flt.fis (FISH Library)
	//
  // PURPOSE:  These functions used when attempting to eliminate
  //           "floaters" (particles with less than [flt_def] contacts)
  //           from the initial packed assembly.
  private:
    double _pc_r_mult;  // Initial radius multiplier
    double _pc_f_tol;   // Fraction of mean contact force to aim for
    double _pc_ncut;    // Number contacts to qualify as a floater
    double _pc_relax;   // ... account for Kn in series
    double _pc_hyst;    // Hysteresis factor on force reduction

    void _pc_rif(void);
      //
    void _pc_flag_floaters(void);
		  // Flag floaters & expand
    void _pc_shrink(void);
		  // shrink active balls
  public:
    Ulib_Flt(void) {};
    ~Ulib_Flt(void) {};
    // static function
    static void flt_eliminate(void);
      // ----- Entry-level call to invoke floater-reduction algorithm.
      // INPUT: flt_def    - balls with < [flt_def] contacts are "floaters"
      //        flt_remain - proportion of total balls that can be floaters
    static void pc_zap_floaters(void);
      // ----- Eliminate the floaters in the assembly.
      // CAUTION:  Friction will be set to zero! ... re-specify at end
      //           Ball color indices will be changed!
      //           All balls will be freed
      //           Set bonds immediately after exit.
      // Can monitor via internal variables:
      //    _pc_n_float  ; Number of floaters
      //    _pc_n_active ; Number of floaters being adjusted
      //    _pc_max_act_force    ; Max floater force
      //    _pc_av_contact_force ; Average contact force
      //;        flt_def    - balls with < [flt_def] contacts are "floaters"
    static void flt_num(void);
      // ----- Returns number of floaters (balls with < [flt_def] contacts)
      //       in model.
      // INPUT:  flt_def - balls with < [flt_def] contacts are "floaters"
      // OUTPUT: flt_num - number of floaters
    static void flt_item(void);
      // ----- Generates a plotitem for visualizing the floaters.
      // Visualization Mapping:
      //    Floaters are outlined by filled circles that use color(0).
};
#endif

void IcgCppLibary::flt_eliminate(void)
{
  _debug = 1      // always dump run-tracking info. to screen
  pc_zap_floaters();
  command
    prop color 0  // reset all ball color indices to 0
  end_command
}

void IcgCppLibary::pc_zap_floaters(void)
{
  // --- parameters ---
  _pc_r_mult = 1.3      // Initial radius multiplier
  _pc_f_tol  = 0.1      // Fraction of mean contact force to aim for
  _pc_ncut   = flt_def  // Number contacts to qualify as a floater
  _pc_relax  = 1.5      // ... account for Kn in series
  _pc_hyst   = 0.9      // Hysteresis factor on force reduction
  // ------------------
  _pc_remain = int(flt_remain * md_numballs)
  //
  command
    hist _pc_n_float          // Number of floaters (internal to this alg.)
    hist _pc_n_active         // Number of floaters being adjusted
    hist _pc_max_act_force    // Max floater force
    hist _pc_av_contact_force // Average contact force
  end_command
  _pc_old_count = 0

	for (int _pc_outer=0;i<10;i++) {
    if (flt_num <= _pc_remain) {
      _pc_rif();
			return;
    }
    _pc_flag_floaters(); // flag and expand floaters
    if (_pc_outer > 1) { // skip this check on first pass
      if (_pc_n_float>=_pc_old_count) { // (got itself into a corner)
        _pc_ncut = max(2,_pc_ncut-1);
        flt_def = _pc_ncut;
      }
    }

    command
      cyc 200 // get clusters to equilibrium
    endCommand

    section
    for (int i=0;i<100;i++) {
      _pc_shrink()
      if (_debug) {
        _s_flt = string(_pc_n_float);
        _s_alt = string(_pc_n_active);
        oo = out(' floaters: '+_s_flt+' active: '+_s_alt);
      }
      if (_pc_n_active==0) {
        _pc_old_count = _pc_n_float;
        break;
      }
      command
        cycle 100
      endCommand
    }
    _pc_rif();
  }
  _pc_rif();
}

void IcgCppLibary::flt_num(void)
{
  flt_num_ = 0
  bp = Udata::ball_head();
  while (bp) {
    if (md_numcon_bp < flt_def)
      flt_num_ = flt_num_ + 1;
		bp = Udata::b_next(bp);
  }
  flt_num = flt_num_
}

void IcgCppLibary::flt_item(void)
{
  plot_item
  Array _flt_pt1(dim)

  stat = set_line_width(0)
  stat = set_color( 0 )

  bp = Udata::ball_head();
  while (bp) {
    if (md_numcon_bp<flt_def) {
      _flt_pt1(1) = b_x(bp);
      _flt_pt1(2) = b_y(bp);
#if DIM == 3
//''      if dim=3 then
        _flt_pt1(3) = b_z(bp)
//      end_if
#endif
      stat = fill_circle( _flt_pt1, b_rad(bp) )
    }
    bp = Udata::b_next(bp);
  }
}

////

void IcgCppLibary::_pc_rif(void)
{
  cp = Udata::contact_head(); // remove inhibit flags
  while (cp) {
    c_inhibit(cp) = 0;
		cp = Udata::c_next(cp);
  }
#if DIM == 2
  command
    free x y spin
  end_command
#elif DIM == 3
  command
    free x y z xspin yspin zspin
  end_command
#endif
}

void IcgCppLibary::_pc_flag_floaters(void)
{
  cp = Udata::contact_head(); // ... get mean contact force
  _pc_mcf  = 0.0
  _pc_n_cf = 0
  while (cp) {
    _pc_mcf  = _pc_mcf + abs(Udata::c_nforce(cp))
    _pc_n_cf = _pc_n_cf + 1
    cp       = Udata::c_next(cp);
  }
  if (_pc_n_cf > 0)
    _pc_mcf = _pc_mcf / float(_pc_n_cf);
  _pc_av_contact_force = _pc_mcf;
  _pc_target           = _pc_f_tol * _pc_mcf;
  _pc_n_float          = 0;
  bp = Udata::ball_head(); //... identify floaters
  while (bp) {
    b_color(bp) = 0;
#if DIM == 2
    b_xvel(bp)  = 0.0;
    b_yvel(bp)  = 0.0;
    b_rvel(bp)  = 0.0;
    b_xfix(bp)  = 0;
    b_yfix(bp)  = 0;
    b_rfix(bp)  = 0;
#elif DIM == 3
    b_xvel(bp)  = 0.0;
    b_yvel(bp)  = 0.0;
    b_zvel(bp)  = 0.0;
    b_rxvel(bp)  = 0.0;
    b_ryvel(bp)  = 0.0;
    b_rzvel(bp)  = 0.0;
    b_xfix(bp)  = 0;
    b_yfix(bp)  = 0;
    b_zfix(bp)  = 0;
    b_rxfix(bp)  = 0;
    b_ryfix(bp)  = 0;
    b_rzfix(bp)  = 0;
#endif
    _pc_count   = 0;
		cp = Udata::b_clist(bp);
    while (cp) {
      if (Udata::c_nforce(cp)) 
        _pc_count ++;
			b1    = Udata::c_ball1(cp);
      b2    = Udata::c_ball2(cp);
      if (b1==bp)
        cp = Udata::c_b1clist(cp);
      else
        cp = Udata::c_b2clist(cp);
    }
    if (_pc_count<_pc_ncut) {
			Udata::b_color(bp) = 1;
			Udata::b_fric(bp)  = 0.0;
			Udata::b_rad(bp)   = Udata::b_rad(bp) * _pc_r_mult;
      _pc_n_float ++;
		} else {
#if DIM == 2
		  Udata::b_xfix(bp) = 1;
			Udata::b_yfix(bp) = 1;
			Udata::b_rfix(bp) = 1;
#elif DIM == 3
			Udata::b_xfix(bp)  = 1;
			Udata::b_yfix(bp)  = 1;
			Udata::b_zfix(bp)  = 1;
			Udata::b_rxfix(bp)  = 1;
			Udata::b_ryfix(bp)  = 1;
			Udata::b_rzfix(bp)  = 1;
#endif
    }
    bp = Udata::b_next(bp);
  }
  // --- set inhibit flags ---
  cp = Udata::contact_head();
	while (cp) {
    b1 = Udata::c_ball1(cp);
    b2 = Udata::c_ball2(cp);
#if DIM == 2
    if (Udata::b_rfix(b1)==1)
      if (Udata::b_rfix(b2)==1)
				Udata::c_inhibit(cp) = 1;
#elif DIM == 3
    if (Udata::b_rxfix(b1)==1)
      if (Udata::b_rxfix(b2)==1)
				Udata::c_inhibit(cp) = 1;
			cp = Udata::c_next(cp);
  }
}

void IcgCppLibary::_pc_shrink(void)
{ 
  _pc_max_act_force = 0.0;
  _pc_n_active      = 0;
  bp = Udata::ball_head();
  while (bp) {
    if (Udata::b_color(bp)==1) {
      _pc_count = 0;
      _pc_sum   = 0.0;
      cp = Udata::b_clist(bp);
      while (cp) {
        if (c_nforce(cp)>0.0) {
          _pc_sum   = _pc_sum + Udata::c_nforce(cp);
          _pc_count = _pc_count + 1;
        }
        b1    = Udata::c_ball1(cp);
        b2    = Udata::c_ball2(cp);
        if (b1==bp) 
          cp = Udata::c_b1clist(cp);
        else
          cp = Udata::c_b2clist(cp);
      }
      if (_pc_count>1) {
        _pc_sum /= float(_pc_count);
        if (_pc_sum > _pc_target) {
          _pc_delta_r  = _pc_relax*(_pc_sum-_pc_hyst*_pc_target)/b_kn(bp);
          b_rad(bp)    = b_rad(bp) - _pc_delta_r;
          _pc_n_active = _pc_n_active + 1;
        }
      }
    }
    _pc_max_act_force = max(_pc_max_act_force,_pc_sum);
    bp = Udata::b_next(bp);
  }
}
