cerr << "Cell " << c.Index() << " is _not_ a source cell.\n";

    }

  }

};

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;

  }

};

double PINSum(Cell &c) {

  return c.Chemical(1) + c.SumTransporters(1);// + c.ReduceCellAndWalls<double>( complex_PijAj );

}

class DrawCell {

#ifdef QDEBUG

	qDebug() << "No cells detected :-(" << endl;

#endif

	return;

      }

      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

	// quickly dragging and releasing division lines...

	scene()->update();

#ifdef QDEBUG

	qDebug() << "Dividing Cell " << c.Index() << endl;

#endif

      mesh.AddNodeSet(node_set);

#ifdef QDEBUG

      qDebug() << "Done DivideOverGivenLine" << endl;

#endif

      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

	qDebug() << "Done CutAwayBelowLine" << endl;

#endif

	mesh.TestIllegalWalls();

	mesh.RepairBoundaryPolygon();

#ifdef QDEBUG

	qDebug() << "Done RepairBoundaryPolygon" << endl;

  // better look for intersection with a simple line intersection algorithm as below?

  // this leads to some exceptions: e.g. dividing a horizontal rectangle.

  // leaving it like this for the time being

  if (dead) return;

  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 ) {

      new_node_locations.push_back(i);

    }		

    prev_cross_z=cross.z;

  }

{

  // First kill walls

#ifdef QDEBUG

  qDebug() << "This is cell " << Index() << endl;

  qDebug() << "Number of walls: " << walls.size() << endl;

#endif

  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) {

      // invert wall?

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

      (*w)->c2 = m->boundary_polygon;

      Node *n1 = (*w)->n1;

      (*w)->n1 = (*w)->n2;

    } else {

#ifdef QDEBUG

      qDebug() << "Not killing wall." << endl;

      qDebug() << "Wall " << (*w)->Index() << " says: c1 = " 

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

#endif

    }

  }

  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;

      }

    }

    if (!cell_found) {

      // I think this cannot happen; but if I am wrong, unpredictable things would happen. So throw an exception.

      throw ("Cell not found in CellBase::Apoptose()\n\rPlease correct the code to handle this situation.");

    }

  }

  // mark cell as dead

  MarkDead();

}

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=

	// remove myself from the neighbor list of the node

	find_if((*i)->owners.begin(),

		(*i)->owners.end(),

		bind2nd(mem_fun_ref( &Neighbor::CellEquals ),this->Index() )  );

      if (neighb_with_this_cell!=(*i)->owners.end()) 

	(*i)->owners.erase(neighb_with_this_cell);

}

bool Cell::DivideOverGivenLine(const Vector v1, const Vector v2, bool fix_cellwall, NodeSet *node_set )

{

  if (dead) return false;

  // 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();

    }

    Vector v4 = *(*nb);

    double denominator = 

      (v4.y - v3.y)*(v2.x - v1.x) - (v4.x - v3.x)*(v2.y - v1.y);

      ((v4.x - v3.x)*(v1.y - v3.y) - (v4.y - v3.y)*(v1.x -v3.x))/denominator;

    double ub = 

      ((v2.x - v1.x)*(v1.y-v3.y) - (v2.y- v1.y)*(v1.x - v3.x))/denominator;

    //cerr << "Edge " << *i << " to " << *nb << ": ua = " << ua << ", ub = " << ub << ":  ";

    // 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)

{

  if (dead) return;

  bool boundary_touched_flag=false;

  // Step 0: keep some data about the parent before dividing

  int new_node_flag[2];

  Edge div_edges[2];

  int nnc=0;

  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;

    Vector v2 = to;

    Vector v3 = *(**i);

    // get previous node

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

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

	   back_inserter(owners));

      copy((div_edges[i].second)->owners.begin(),

	   (div_edges[i].second)->owners.end(),

	   back_inserter(owners));

      // 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;

      else 

	neighbor_cell = 0;

    }

    if (neighbor_cell /* && !neighbor_cell->BoundaryPolP() */) {

  // Add a "Cell Wall" for diffusion algorithms

  Wall *wall = new Wall( new_node_ind[0], new_node_ind[1], this, daughter );

  AddWall( wall );

  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 << "and after: " << **w << endl;

  }

  // Correct tranporterconcentrations of divided walls

  for (int i=0;i<4;i++) {

    if (div_wall[i]) {

      div_wall[i]->SetLength();

      div_wall[i]->CorrectTransporters(orig_length[i/2]);

    }

  }

  //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();

  m->plugin->OnDivide(&parent_info, daughter, this);

  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)

{

  // Displace whole cell, add resulting energy to dh,

  // and accept displacement if energetically favorable

  // 

  // Method is called if a "fixed" node is displaced

  Vector movement(dx,dy,0);

  vector< pair<Node *, Node *> > length_edges;

  vector<double> cellareas;

  cellareas.reserve(neighbors.size());

  // 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;

  } else {

    Move ( -1*movement);

  }

void Cell::Displace (void)

{

  Displace(par.mc_cell_stepsize*(RANDOM()-0.5),par.mc_cell_stepsize*(RANDOM()-0.5),0);

}

// 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;

  // area constraint

  energy += DSQR(CalcArea() - target_area);

  // cell length constraint

  energy += lambda_celllength * DSQR(Length() - target_length);

  // The (obvious) O(N*N) algorithm

  // Compare each edge against each other edge

  // An O(N log(N)) algorithm by Shamos & Hoey (1976) supposedly exists;

  // it was mentioned on comp.graphics.algorithms

  // 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 = nodes.begin();

	} 

	Vector v2 = *(*nb);

	Vector v3 = *(*j);

	nb=j;

	nb++;

  neighbor_of_moving_node[0]=*(*nb); 

  nb=moving_node_ind_pos;

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

    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();

      } 

      if (*i == moving_node_ind || *nb == moving_node_ind) {

	// do not compare to self

	continue;

      }

      Vector v3 = *(*i);

*/

void Cell::ConstructWalls(void)

{

  return;

  if (dead) return;

  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;

      }

    }

    // One Wall for each neighbor, so we should be able to correctly construct neighbor lists here.

    if (duplicates[0]==this) {

void BoundaryPolygon::Draw(QGraphicsScene *c, QString tooltip)

{

  // Draw the BoundaryPolygon on a QCanvas object

  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);

  if (!tooltip.isEmpty())

    p->setToolTip(tooltip);

  p->show();

}

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

    for (int c=0;c<NChem();c++) {

      double phi = (*i)->length * ( D[c] ) * ( ((Cell *)(*i)->c2)->chem[c] - chem[c] );

#ifdef QDEBUG

      if ((*i)->c1!=this) {

	qDebug() << "Warning, bad cells boundary: " << (*i)->c1->Index() << ", " << index << endl;

  if (DeadP()) { 

#ifdef QDEBUG

    qDebug() << "Cell " << index << " not drawn, because dead." << endl;

#endif

    return;

  }

  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) );

  }

  QColor cell_color;

  m->plugin->SetCellColor(this,&cell_color);

  p->setPolygon(pa);

  // construct an ellipse

  QGraphicsEllipseItem *disk = new QGraphicsEllipseItem ( -1*mag, -1*mag, 2*mag, 2*mag, 0, c );

  disk->setBrush( QColor("forest green") );

  disk->setZValue(5);

  disk->show();

  Vector centroid=Centroid();

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

}