class CellInfo {

public:

  void operator() (Cell &c,std::ostream &os) const {

    os << "Cell " << c.index << " says: " << endl;

    os << "c.nodes.size() = " << c.nodes.size() << endl;

      cerr << (*i)->Index() << " ";

    }

    cerr << endl;

  }

};

      Vector startpoint = Vector(sp.x(), sp.y()) / Cell::Factor() - Cell::Offset();

      Vector endpoint = Vector(ep.x(), ep.y()) / Cell::Factor() - Cell::Offset();

      NodeSet *node_set = new NodeSet;

	(*it)->setBrush(QBrush("purple"));

	Cell &c=(*it)->getCell();

	// sometimes the cell hasn't properly divided yet before the

	// next division is called?  so check for it?  let's find a way

	// to do this later. Note that this functionality currently

	// might result in a segmentation fault for users who are

      mesh.TestIllegalWalls();

      // Do the actual cutting and removing

      if (intersected_cells.size()) {

	mesh.CutAwayBelowLine( startpoint, endpoint ); 

	// Correct flags of nodeset

	  (*i)->SetSAM();

	  (*i)->SetBoundary();

	}

	// Make cells attached to nodeset part of the boundary

	// This is a temporary solution for the following:

	//   If a cell attached to a NodeSet divides (with a division plane intersecting the node set), 

	//   we must insert a new node into the node set.

	// For now, we add a layer of "virtual cells" inbetween. 

	list<Cell *> cells_attached_to_nodeset = node_set->getCells();

	  (*c)->SetBoundary(Cell::SAM);

	}

#ifdef QDEBUG

  Vector centroid=Centroid();

  double prev_cross_z=(axis * (centroid - *(nodes.back()) ) ).z ;

  ItList new_node_locations;

    // cross product to detect position of division

    Vector cross = axis * (centroid - *(*i));

    if (cross.z * prev_cross_z < 0 ) {

  for (list<Wall *>::iterator w=walls.begin(); w!=walls.end(); w++) {

#ifdef QDEBUG

    qDebug() << "Before apoptosis, wall " << (*w)->Index() << " says: c1 = "

	     << (*w)->c1->Index() << ", c2 = " << (*w)->c2->Index() << endl;

#endif

  }

    bool illegal_flag = false;

    if ((*w)->c1 == (*w)->c2 ) 

      illegal_flag=true;

    if ((*w)->c1 == this) {

    }

  }

  walls.remove(0);

  // Unregister me from my nodes, and delete the node if it no longer belongs to any cells

  list<Node *> superfluous_nodes;

    Node &no(*(*n));

    // locate myself in the node's owner list

    list<Neighbor>::iterator cellpos;

    bool cell_found=false;

      if (nb->cell == this) {

	cellpos = nb;

	cell_found = true;

	break;

      }

    }

void Cell::ConstructConnections(void)

{

  // Tie up the nodes of this cell, assuming they are correctly ordered

  //cerr << "Constructing connections of cell " << index << endl;

    //cerr << "Connecting node " << *i << endl;

    //cerr << "Node " << *i << endl << " = " << *(*i) << endl;

    // 1. Tidy up existing connections (which are part of this cell)

    if ((*i)->owners.size()>0) {

      list<Neighbor>::iterator neighb_with_this_cell=

  // check each edge for intersection with the line

  ItList new_node_locations;

#ifdef QDEBUG

  qDebug() << "Cell " << Index() << " is doing DivideOverGivenLine" << endl;

#endif

    Vector v3 = *(*i);

    list<Node *>::iterator nb=i;

    nb++;

    if (nb == nodes.end()) {

      nb = nodes.begin();

    // this construction with "TINY" should simulate open/closed interval <0,1]

    if ( ( TINY < ua && ua < 1.+TINY ) && ( TINY < ub && ub < 1.+TINY ) ) {

      // yes, intersection detected. Push the location to the list of iterators

      new_node_locations.push_back(nb);

    } 

#ifdef QDEBUG

  if (new_node_locations.size()<2) { 

    qDebug() << "Line does not intersect with two edges of Cell " << Index() << endl;

    qDebug() << "new_node_locations.size() = " << new_node_locations.size() << endl;

    return false;

  }

#endif

}

// Core division procedure

void Cell::DivideWalls(ItList new_node_locations, const Vector from, const Vector to, bool fix_cellwall, NodeSet *node_set)

{

  Wall *div_wall[4];

  double orig_length[2];

  for (int i=0;i<4;i++) { div_wall[i]=0; orig_length[i/2] = 0.; }

  // construct new Nodes at the intersection points

  // unless they coincide with existing points

    // intersection between division axis

    // and line from this node to its predecessor

    // method used: http://astronomy.swin.edu.au/~pbourke/geometry/lineline2d/

    Vector v1 = from;

      // find first non-self duplicate in the owners: 

      // cells owning the same two nodes

      // share an edge with me

      owners.sort( mem_fun_ref( &Neighbor::Cmp ) );

      list<Neighbor>::iterator c;

	c=adjacent_find(c,owners.end(),neighbor_cell_eq);

	if (c->cell->Index() != this->Index() || c==owners.end()) break;

      }

      if (c!=owners.end())

	neighbor_cell = c->cell;

  daughter->AddWall( wall );

  //cerr << "Correct walls of cell " << Index() << " and daughter " << daughter->Index() << endl;

  // Move Walls to daughter cell

  list <Wall *> copy_walls = walls;

    //cerr << "Doing wall, before:  " << **w << endl;

    //  checks the nodes of the wall and gives it away if appropriate

    (*w)->CorrectWall ( );

  }

  //cerr << "Cell " << index << " has been dividing, and gave birth to Cell " << daughter->index << endl;

  // now reconstruct neighbor list for all "broken" neighbors

    ((Cell *)(*i))->ConstructNeighborList();

  }

  ConstructNeighborList();

  daughter->ConstructNeighborList();

  daughter->div_counter=(++div_counter);

}

// Move the whole cell

void Cell::Move(const Vector T) {

    *(*i)+=T;

  }

}

double Cell::Displace(double dx, double dy, double dh)

{

  // for the length constraint, collect all edges to this cell's nodes,

  // which are not part of the cell

  // the length of these edges will change

  double old_length=0.;

      if (n->getCell()!=this) {

	length_edges.push_back( pair <Node *,Node *> (*i, n->nb1) );

	length_edges.push_back( pair <Node *,Node *> (*i, n->nb2) );

	old_length += 

	  DSQR(Node::target_length-(*(*i)-*(n->nb1)).Norm())+

	  DSQR(Node::target_length-(*(*i)-*(n->nb2)).Norm());

    }

  }

  // calculate area energy difference of neighboring cells

  // (this cells' shape remains unchanged)

  double old_area_energy=0., old_length_energy=0.;

    old_area_energy += DSQR((*i)->Area()-(*i)->TargetArea());

    old_length_energy += DSQR((*i)->Length()-(*i)->TargetLength());

  }

  Move(movement);

  double new_area_energy=0., new_length_energy=0.;

    cellareas.push_back((*i)->CalcArea());

    new_area_energy += DSQR(cellareas.back()-(*i)->TargetArea());

    new_length_energy += DSQR((*i)->CalcLength()-(*i)->TargetLength());

  }

  double new_length=0;

  }

  dh += (new_area_energy - old_area_energy) + (new_length_energy - old_length_energy) * lambda_celllength +

    par.lambda_length * (new_length - old_length);

  if (dh<0 || RANDOM()<exp(-dh/par.T)) {

    // update areas of cells

    //cerr << "neighbors: ";

    list<CellBase *>::const_iterator nb_it = neighbors.begin();

      ((Cell *)(*nb_it))->area = *ar_it;

      (*nb_it)->SetIntegrals(); 

    }

    //cerr << endl;

// Get energy level of whole cell (excluding length constraint?)

double Cell::Energy(void) const

{

  double energy = 0.;

  double length_contribution = 0.;

      if (n->getCell()==this) {

	length_contribution += 

	  DSQR(Node::target_length-(*(*i)-*(n->nb2)).Norm());

      }

    }

  }

  // wall elasticity constraint

  energy += par.lambda_length * length_contribution;

  // But I haven't been able to lay my hand on the paper.

  // Let's try whether we need it....

  // method used: http://astronomy.swin.edu.au/~pbourke/geometry/lineline2d/

    list<Node *>::const_iterator j=i; 

    ++j;

      {

	Vector v1 = *(*i);

	list<Node *>::const_iterator nb=i;

	nb++;

	if (nb == nodes.end()) {

  }

  nb--;

  neighbor_of_moving_node[1]=*( *nb ); 

    for (int j=0;j<2;j++) { // loop over the two neighbors of moving node

      list<Node *>::const_iterator nb=i;

      nb++;

      if (nb == nodes.end()) {

	nb = nodes.begin();

      } 

  walls.clear();

  neighbors.clear();

  // Get "corner points; i.e. nodes where more than 2 cells are connected

  list<Node *> corner_points;

    // look for nodes belonging to >2 cells

    if ((*i)->owners.size()>2) {

      // push onto list

      corner_points.push_back(*i);

  // Construct Walls between corner points

  // previous one in list

  list<Node *>::const_iterator nb = (--corner_points.end());

  // loop over list, 

    if (nb==corner_points.end()) nb=corner_points.begin();

    // add owning cells to a list

    list<Cell *> owning_cells;

    Node &n(*(*i));

      owning_cells.push_back(j->cell);

    }

    Node &n2(*(*nb));

      owning_cells.push_back(j->cell);

    }

    // sort cell owners

    owning_cells.sort( mem_fun( &Cell::Cmp ));

    // find duplicates

    vector<Cell *> duplicates;

    list<Cell *>::const_iterator prevj = (--owning_cells.end());

    }

    if (duplicates.size()==3) { // ignore cell boundary (this occurs only after the first division, I think)

      vector<Cell *>::iterator dup_it = find_if(duplicates.begin(),duplicates.end(),mem_fun(&Cell::BoundaryPolP) );

      if (dup_it!=duplicates.end()) 

	duplicates.erase(dup_it);

      else {

	return;

  CellItem* p = new CellItem(this, c);

  QPolygonF pa(nodes.size());

  int cc=0;

    Node *i=*n;

  }

  p->setPolygon(pa);

  p->setPen(par.outlinewidth>=0?QPen( QColor(par.cell_outline_color),par.outlinewidth):QPen(Qt::NoPen));

  p->setBrush( Qt::NoBrush );

  p->setZValue(1);

}

void Cell::Flux(double *flux, double *D)

{

  // loop over cell edges

    // leaf cannot take up chemicals from environment ("no flux boundary")

    if ((*i)->c2->BoundaryPolP()) continue;

    // flux depends on edge length and concentration difference

  CellItem* p = new CellItem(this, c);

  QPolygonF pa(nodes.size());

  int cc=0;

    Node *i=*n;

    pa[cc++] = QPoint((int)((offset[0]+i->x)*factor),

		      (int)((offset[1]+i->y)*factor) );

  }

  Vector centroid=Centroid();

  disk -> setPos((offset[0]+centroid.x)*factor,(offset[1]+centroid.y)*factor);

}

void Cell::DrawNodes(QGraphicsScene *c) const {

    Node *i=*n;

    NodeItem *item = new NodeItem ( &(*i), c );

    item->setColor();

    item->setZValue(5);

    item->show();

  }

}

void Cell::DrawIndex(QGraphicsScene *c) const {

  DrawText( c, QString("%1").arg(index));

  Call this function after the Monte Carlo updates, and before doing the reaction-diffusion iterations.

*/

void Cell::SetWallLengths(void)

{

    // Step 1: find the path of nodes leading along the Wall.

    // A Wall often represents a curved cell wall: we want the total

    // length _along_ the wall here...

    list<Node *>::const_iterator first_node_edge = find(nodes.begin(), nodes.end(), (*de)->n1);

    list<Node *>::const_iterator second_node_edge_plus_1 = ++find(nodes.begin(), nodes.end(), (*de)->n2);

    double sum_length = 0.;

    // Now, walk to the second node of the edge in the list of nodes

      list<Node *>::const_iterator prev_n = n; 

      --prev_n;

      // Note that Node derives from a Vector, so we can do vector calculus as defined in vector.h 

      sum_length += (*(*prev_n) - *(*n)).Norm(); 

      //cerr << "Node " << *prev_n << " to " << *n << ", cumulative length = " << sum_length << endl;

    }

#endif

  // Add Wall to Cell's list

  walls.push_back( w );

  // Add wall to Mesh's list if it isn't there yet

  }

}

//! Remove Wall w from the list of Walls

list<Wall *>::iterator Cell::RemoveWall( Wall *w )

{

  // remove wall from Mesh's list

  // remove wall from Cell's list

  return walls.erase (find( walls.begin(), walls.end(), w ));

}

void Cell::EmitValues(double t)

{

  //  cerr << "Attempting to emit " << t << ", " << chem[0] << ", " << chem[1] << endl;

  emit ChemMonValue(t, chem);

}

/* finis */

{

  if (ch>=NChem()) {

    stringstream error;

    error << "SetChemical: value ch = " << ch << " is out of range\n";

    throw error.str().c_str();

  }

    (*w)->setTransporter(this, ch, conc);

  }

}

ostream &CellBase::print(ostream &os) const

{

  } 

  os << endl;

  os << " [ area = " << area << " ]";

  os << " [ walls = ";

    os << (*i)->n1->Index() << " -> " << (*i)->n2->Index() << ", " <<  (*i)->c1->Index() << " | " << (*i)->c2->Index() << ", ";

  }

  os << " ] ";

  os << "div_counter = " << div_counter << endl;

  os << "cell_type = " << cell_type << endl;

  os << endl;

double CellBase::CalcArea(void) const

{

  double loc_area=0.;

    list<Node *>::const_iterator i_plus_1=i; i_plus_1++;

    if (i_plus_1==nodes.end())

      i_plus_1=nodes.begin();

    loc_area+= (*i)->x * (*i_plus_1)->y;

Vector CellBase::Centroid(void) const

{

  double area=0.;

  double integral_x_dxdy=0.,integral_y_dxdy=0.;

    list<Node *>::const_iterator i_plus_1=i; i_plus_1++;

    if (i_plus_1==nodes.end())

      i_plus_1=nodes.begin();

    area+= (*i)->x * (*i_plus_1)->y;

  intgrl_xx=0.; intgrl_xy=0.; intgrl_yy=0.;

  intgrl_x=0.; intgrl_y=0.;

  area=0.;

  list<Node *>::const_iterator nb;

  list<Node *>::const_iterator i=nodes.begin();

    nb = i; nb++; if (nb==nodes.end()) nb=nodes.begin();

    area+=(*i)->x*(*nb)->y;

    area-=(*nb)->x*(*i)->y;

    intgrl_xx+= 

  double my_intgrl_xx=0., my_intgrl_xy=0., my_intgrl_yy=0.;

  double my_intgrl_x=0., my_intgrl_y=0., my_area=0.;

  my_area=0.;

  list<Node *>::const_iterator nb;

  list<Node *>::const_iterator i=nodes.begin();

    nb = i; nb++; if (nb==nodes.end()) nb=nodes.begin();

    my_area+=(*i)->x*(*nb)->y;

    my_area-=(*nb)->x*(*i)->y;

    my_intgrl_xx+= 

     << " " << intgrl_x << " " << intgrl_y << " " << source << " ";

  cellvec.Dump(os);

  os << " " << source_conc << " " << source_chem;

  os << endl;

}

void CellBase::UnfixNodes(void)

{

    (*i)->Unfix();

  }

}

void CellBase::FixNodes(void)

{

    (*i)->Fix();

  }

}

// returns true if cell is at border

bool CellBase::AtBoundaryP(void) const

{

  return at_boundary;