squirmgrid.cpp

本程序模拟细胞的自我繁殖

    SquirmCellList reaction_candidates;

    POSITION pos1,pos2;

    SquirmCell *cell,*other;
    C2dVector force,to;
    float to_len2;
	float wall_range;

    const float XY_DIST = 2.0F*(float)sqrt(MAX_RANGE*MAX_RANGE/2.0F);

    for(cx=0;cx<this->slots_X;cx++)
    {
        for(cy=0;cy<this->slots_Y;cy++)
        {
            InitSearchArea(cx,cy,search_slots,central_cells,nearby_cells);
            pos1 = central_cells->GetHeadPosition();
            while(pos1)
            {
                cell = central_cells->GetNext(pos1);
                reaction_candidates.RemoveAll();

		        // -- sum the forces acting on this cell --

		        force.x = 0.0F;
                force.y = 0.0F;

                // a special case: unbonded atoms don't have any physical forces on them
                if(cell->GetBondedCells().IsEmpty() && cell->GetState()==0)
                {
                    // no forces
                }
                else {
		            // the repulsive forces between nearby cells
		            pos2 = nearby_cells.GetHeadPosition();
		            while(pos2)
		            {
			            other = nearby_cells.GetNext(pos2);
                        if(other==cell) continue; // no interaction with self!

				        to.x = cell->location.x;
				        to.y = cell->location.y;
                        to.Subtract(other->location);

                        // for speed (minor saving), we first check if this cell is within 
                        // an octagon before going for the full circle test
                        if(fabs(to.x)+fabs(to.y) > XY_DIST) continue;

                        to_len2 = Length2(to);

                        if(to_len2<REACTION_RANGE2)
                        {
                            // we don't react if the atom is too close. this was introduced
                            // when we turned off the volume exclusion force for atoms in 
                            // state 0

                            // check whether anything is too near to this atom (slow)
                            SquirmCellList too_close;
                            this->GetAllWithinRadius(other->location,RADIUS*1.5F,too_close);
                            if(too_close.GetCount()<=1)
                                reaction_candidates.AddTail(other);
                        }

                        if((cell->GetBondedCells().IsEmpty() && cell->GetState()==0) ||
                            (other->GetBondedCells().IsEmpty() && other->GetState()==0))
                        {
                            // no forces
                        }
                        else {
				            // atoms that are too close get repelled
                            to.Mult(compute_force(to_len2,PHYS_RANGE,0.4F));
				            force.Add(to);
                        }
		            }
                }

		        // pull forces from bonded cells
		        {
			        pos2 = cell->GetBondedCells().GetHeadPosition();
			        while(pos2)
			        {
				        other = cell->GetBondedCells().GetNext(pos2);
				        if(other!=cell)
				        {
                            to.x = other->location.x;
                            to.y = other->location.y;
                            to.Subtract(cell->location);
                            to.Mult(spring_force(Length2(to)));
                            force.Add(to);
				        }
                        else
                            SquirmError("Cell bonded with self!?");
			        }
		        }

		        // forces from walls
			    wall_range = 1.0F*RADIUS;
		        if(cell->location.x<wall_range) // left wall
		        {
			        to.x = cell->location.x;
                    to.y = 0.0F;
                    // length of to vector is easy to compute (sign must be positive)
                    to.Mult(wall_range/cell->location.x-1.0F);
                    force.Add(to);
		        }
		        if(cell->location.x>size.x-wall_range) // right wall
		        {
                    to.x = cell->location.x-size.x;
                    to.y = 0.0F;
                    // length of to vector is easy to compute (sign must be positive)
                    to.Mult(wall_range/(size.x-cell->location.x)-1.0F);
                    force.Add(to);
		        }
		        if(cell->location.y<wall_range) // top wall
		        {
                    to.x = 0.0F;
                    to.y = cell->location.y;
                    // length of to vector is easy to compute (sign must be positive)
                    to.Mult(wall_range/cell->location.y-1.0F);
                    force.Add(to);
		        }
		        if(cell->location.y>size.y-wall_range) // bottom wall
		        {
                    to.x = 0.0F;
                    to.y = cell->location.y-size.y;
                    // length of to vector is easy to compute (sign must be positive)
                    to.Mult(wall_range/(size.y-cell->location.y)-1.0F);
                    force.Add(to);
		        }

		        cell->velocity.Add(force);

		        // we must limit the velocity else atoms can pass through each other
                float speed2 = Length2(cell->velocity);
		        if(speed2>MAX_VELOCITY*MAX_VELOCITY)
			        cell->velocity.Mult(MAX_VELOCITY/(float)sqrt(speed2));

                // apply any reactions that are possible
                int ret = this->chemistry.react(cell,reaction_candidates,REACTION_RANGE2);

                // observer: make a note if a certain reaction occurred (highly specific)
                if(this->do_output && ret==0 && cell->GetState()==10 && cell->GetType()=='e')
                {
                    CString molecule;
                    molecule += cell->GetType();
                    SquirmCell *a,*found;
                    SquirmCellList so_far;
                    so_far.AddTail(cell);
                    a=cell;
                    bool ok=true,found_next;
                    while(ok)
                    {
                        if(a->GetBondedCells().GetCount()!=2) ok=false;
                        else { 
                            POSITION pos = a->GetBondedCells().GetHeadPosition();
                            found_next=false;
                            while(pos)
                            {
                                found = a->GetBondedCells().GetNext(pos);
                                if(found->GetState()==1 && so_far.Find(found)==NULL)
                                {
                                    a=found;
                                    molecule += a->GetType();
                                    found_next=true;
                                    so_far.AddTail(found);
                                    break;
                                }
                            }
                            if(!found_next || a->GetType()=='f') ok=false;
                        }
                    }
                    fprintf(this->output,"%d,%d,%s\n",this->iterations,ret,molecule);
                }
            }
        }
    }
}

void SquirmGrid::MoveCells()
{
	// then, update the position of each cell using its velocity
	C2dVector new_loc;
    SquirmCell *cell;

    POSITION pos = this->cell_list.GetHeadPosition();
	while(pos)
	{
		cell = this->cell_list.GetNext(pos);

	    int old_slot_x = (int)floor(cell->location.x/slot_size);
	    int old_slot_y = (int)floor(cell->location.y/slot_size);

		new_loc = cell->location + cell->velocity;
		if(Valid(new_loc))
			cell->location = new_loc;
		else
        {
            if(cell->location.x<0.0F)
                cell->location.x = 0.0F;
            else if(cell->location.x>this->size.x)
                cell->location.x = this->size.x;
            if(cell->location.y<0.0F)
                cell->location.y = 0.0F;
            else if(cell->location.y>this->size.y)
                cell->location.y = this->size.y;
			cell->velocity = C2dVector(0,0);
        }

	    int new_slot_x = (int)floor(cell->location.x/slot_size);
	    int new_slot_y = (int)floor(cell->location.y/slot_size);

        if(old_slot_x!=new_slot_x || old_slot_y!=new_slot_y)
        {
            this->cell_grid[old_slot_x][old_slot_y].removeOccupant(cell);
            this->cell_grid[new_slot_x][new_slot_y].addOccupant(cell);
        }
	}
}

CString SquirmGrid::ProbeAt(float x,float y,float scale)
{ 
    /*if(this->probe) 
    {
        RemoveCell(probe); 
        probe->location = C2dVector(x,y); 
        PutCell(probe);
    }*/

    SquirmCellList dest;
    this->GetAllWithinRadius(C2dVector((x-OFFSET)/scale,(y-OFFSET)/scale),RADIUS,dest);
    CString text;
    if(!dest.IsEmpty())
    {
        SquirmCell *cell = dest.GetHead();
        text.Format("%c%d",cell->GetType(),cell->GetState());
        return text;
    }
    else
        return "";
}

float SquirmGrid::compute_force(float r2,float range,float magnitude)
{
    if(r2 < range*range)
        return magnitude * (range/(float)sqrt(r2) - 1.0F); // our force/distance curve
    else
        return 0.0F;
}

/* code leftover from when we implemented Ono's continuous lipids model

float SquirmGrid::force1(float r2)
{
    static const float RANGE = 10.0F*RADIUS;
    return compute_force(r2,RANGE,0.02F);
}

float SquirmGrid::force2(float r2)
{
    static const float RANGE = 3.0F*RADIUS;
    return compute_force(r2,RANGE,0.1F);
}

float SquirmGrid::force3(float r2)
{
    static const float RANGE = 2.0F*RADIUS;
    return compute_force(r2,RANGE,0.2F);
}*/

float SquirmGrid::spring_force(float r2)
{
    static const float RANGE = 2.0F*RADIUS;
    static const float RANGE2 = RANGE*RANGE;
    if(r2 > RANGE2) 
        return 0.4F * ( (float)sqrt(r2)/RANGE - 1.0F );
    else
        return 0.0F;
}

void SquirmGrid::EdgeEffect()
{
	// any unbonded cell near the edge becomes a random type in state 0 (occasionally)
    // (simulates bordering on a larger world)

    SquirmCell *cell;

    POSITION pos = this->cell_list.GetHeadPosition();
	while(pos)
	{
		cell = this->cell_list.GetNext(pos);

        if( cell->GetBondedCells().IsEmpty() && (cell->location.x<RADIUS || cell->location.y<RADIUS ||
          cell->location.x>this->GetSize().x-RADIUS || cell->location.y>this->GetSize().y-RADIUS) )
        {
            if(RandomN(100)==0)
            {
                cell->SetType(SquirmCellProperties::GetRandomType());

                // to introduce caustic particles, we occasionally set the state to -1
                // to remove this effect, comment out the next three lines
                if(RandomN(100)==0)
                    cell->SetState(-1);
                else
                    cell->SetState(0);
            }
        }
    }
}