/************************************************/
/*												*/
/*	Model.c - Model handler routines.			*/
/*												*/
/*	Revised 7/6/94	by RMR						*/
/*												*/
/************************************************/

#include <stdio.h>
#include "Snakes.h"

/************************************************/
/*		Model Placement Routines				*/
/************************************************/

void ClearModel(ModelPtr model)
{
	model->sideAParticles = 0;
	model->sideBParticles = 0;
	model->numberOfSprings = 0;
	model->sideADone = false;
	model->sideBDone = false;
	model->totalKineticEnergy = 0.0;
	model->totalForceMagnitude = 0.0;
	model->averageWidth = 0.0;
	SetVector(&(model->sumInternalForces), 0.0, 0.0);
	SetVector(&(model->sumExternalForces), 0.0, 0.0);
}

static void AddSpring(ModelPtr model, Particle *from, Particle *to, SpringType type)
{
	short n;
	
	n = model->numberOfSprings;
	if (n >= MAXSPRINGS)
		DoError("\pExceeded max # of springs!");
	model->springList[n].from = from;
	model->springList[n].to = to;
	model->springList[n].type = type;
	model->springList[n].lengthAtRest = Distance(from->position, to->position);
	model->numberOfSprings += 1;
	SetVector(&(from->velocity), 0.0, 0.0);
	SetVector(&(from->force), 0.0, 0.0);
	SetVector(&(to->velocity), 0.0, 0.0);
	SetVector(&(to->force), 0.0, 0.0);
}

void BuildSpringsFromParticles(ModelPtr model)
{
	short i, maxA, maxB;
	
	model->numberOfSprings = 0;
	maxA = model->sideAParticles;
	maxB = model->sideBParticles;
	for (i = 0; i < MAXRUNGS; i++)
	{
		if (i < maxA && i < maxB)
		{
			AddSpring(model, &model->sideAList[i], &model->sideBList[i], tRung);
			if (i > 0)
			{
				AddSpring(model, &model->sideAList[i], &model->sideAList[i-1], tEdge);
				AddSpring(model, &model->sideBList[i], &model->sideBList[i-1], tEdge);
				AddSpring(model, &model->sideAList[i], &model->sideBList[i-1], tBrace);
				AddSpring(model, &model->sideBList[i], &model->sideAList[i-1], tBrace);
			}
		}
		else if (i < maxA && i >= maxB)
		{
			if (i > 0)
				AddSpring(model, &model->sideAList[i], &model->sideAList[i-1], tEdge);
		}
		else if (i >= maxA && i < maxB)
		{
			if (i > 0)
				AddSpring(model, &model->sideBList[i], &model->sideBList[i-1], tEdge);
		}
		else
			break;
	}
}

static Particle *PtInParticle(Point p)
{
	short		i, j, k, d;
	Rect		r;
	Particle	*particle;
	
	d = PARTICLE_RADIUS;
	for (i = 0; i < gRefModel->sideAParticles; i++)
	{
		particle = &gRefModel->sideAList[i];
		j = (short)particle->position.x;
		k = (short)particle->position.y;
		SetRect(&r, j - d, k - d, j + d, k + d);
		if (PtInRect(p, &r))
			return particle;
	}
	for (i = 0; i < gRefModel->sideBParticles; i++)
	{
		particle = &gRefModel->sideBList[i];
		j = (short)particle->position.x;
		k = (short)particle->position.y;
		SetRect(&r, j - d, k - d, j + d, k + d);
		if (PtInRect(p, &r))
			return particle;
	}
	for (i = 0; i < gTestModel->sideAParticles; i++)
	{
		particle = &gTestModel->sideAList[i];
		j = (short)particle->position.x;
		k = (short)particle->position.y;
		SetRect(&r, j - d, k - d, j + d, k + d);
		if (PtInRect(p, &r))
			return particle;
	}
	for (i = 0; i < gTestModel->sideBParticles; i++)
	{
		particle = &gTestModel->sideBList[i];
		j = (short)particle->position.x;
		k = (short)particle->position.y;
		SetRect(&r, j - d, k - d, j + d, k + d);
		if (PtInRect(p, &r))
			return particle;
	}
	particle = nil;
	return nil;
}

static void PlaceParticle(Point p, ModelPtr model, Boolean sideB)
{
	short 			*numberOfParticles;
	Particle 		*particleBase;

	numberOfParticles = sideB ? &model->sideBParticles : &model->sideAParticles;
	particleBase = sideB ? model->sideBList : model->sideAList;
	if (*numberOfParticles >= MAXRUNGS)
		DoError("\pExceeded number of particles allowed!");
	particleBase[*numberOfParticles].position.x = (Scalar)p.h;
	particleBase[*numberOfParticles].position.y = (Scalar)p.v;
	particleBase[*numberOfParticles].velocity.x = (Scalar)0.0;
	particleBase[*numberOfParticles].velocity.y = (Scalar)0.0;
	particleBase[*numberOfParticles].force.x = (Scalar)0.0;
	particleBase[*numberOfParticles].force.y = (Scalar)0.0;
	particleBase[*numberOfParticles].velocity.x = (Scalar)0.0;
	*numberOfParticles += 1;
	if (*numberOfParticles == MAXRUNGS)
	{
		if (sideB)
			model->sideBDone = true;
		else
			model->sideADone = true;
	}
	BuildSpringsFromParticles(gRefModel);
	BuildSpringsFromParticles(gTestModel);
	DrawPlaceDialog();
}

static void DragParticle(Particle *particle, Point p, WindowPtr wind)
{
	GrafPtr	temp;
	Point	pt;
	Rect	r;
	short	d;
	
	GetPort(&temp);
	SetPort(wind);
	d = PARTICLE_RADIUS;
	RGBBackColor(&kBlack);
	RGBForeColor(&kBlack);
	pt = p;
	SetRect(&r, pt.h - d, pt.v - d, pt.h + d, pt.v + d);
	FrameOval(&r);
	while (StillDown())
	{
		GetMouse(&pt);
		SetRect(&r, pt.h - d, pt.v - d, pt.h + d, pt.v + d);
		RGBForeColor(&kWhite);
		FrameOval(&r);
		SkelPause(2);
		SetRect(&r, pt.h - d, pt.v - d, pt.h + d, pt.v + d);
		RGBForeColor(&kBlack);
		FrameOval(&r);
	}
	particle->position.x = (Scalar)pt.h;
	particle->position.y = (Scalar)pt.v;
	SetPort(temp);
}

void PictPoint(Point p, WindowPtr wind)
{
	Particle *particle;
	
	if (!gPictArray.valid)
		DoError("\pPict Array not valid!");
	particle = PtInParticle(p);
	switch (gModelStatus)
	{
		case sPlaceRefA:
			if 	(particle != nil 		|| 
				gRefModel->sideADone	||
				gRefModel->sideAParticles >= MAXRUNGS)
				SysBeep(10);
			else
				PlaceParticle(p, gRefModel, false);
			break;
			
		case sPlaceRefB:
			if 	(particle != nil		|| 
				gRefModel->sideBDone	||
				gRefModel->sideBParticles >= MAXRUNGS)
				SysBeep(10);
			else
				PlaceParticle(p, gRefModel, true);
			break;
			
		case sPlaceTestA:
			if 	(particle != nil 		|| 
				gTestModel->sideADone	||
				gTestModel->sideAParticles >= MAXRUNGS)
				SysBeep(10);
			else
				PlaceParticle(p, gTestModel, false);
			break;
			
		case sPlaceTestB:
			if 	(particle != nil		|| 
				gTestModel->sideBDone	||
				gTestModel->sideBParticles >= MAXRUNGS)
				SysBeep(10);
			else
				PlaceParticle(p, gTestModel, true);
			break;
			
		default:
			if (particle != nil)
				DragParticle(particle, p, wind);
	}
	KickPictWindow();
}


/************************************************/
/*		Model Placement Dialog Routines			*/
/************************************************/

void DrawPlaceDialog(void)
{
	GrafPtr	temp;
	
	GetPort(&temp);
	SetPort(gPlaceDialog);
	RGBForeColor(&kBlack);
	RGBBackColor(&kWhite);
	switch (gModelStatus)
	{
		case sPlaceRefA:
			SkelSetDlogRadioButtonSet(gPlaceDialog, iPlaceRefA, iPlaceTestB, iPlaceRefA);
			break;
			
		case sPlaceRefB:
			SkelSetDlogRadioButtonSet(gPlaceDialog, iPlaceRefA, iPlaceTestB, iPlaceRefB);
			break;
			
		case sPlaceTestA:
			SkelSetDlogRadioButtonSet(gPlaceDialog, iPlaceRefA, iPlaceTestB, iPlaceTestA);
			break;
			
		case sPlaceTestB:
			SkelSetDlogRadioButtonSet(gPlaceDialog, iPlaceRefA, iPlaceTestB, iPlaceTestB);
	}
	DrawParticleStatus(gPlaceDialog, iTextRefA, gRefModel, false);
	DrawParticleStatus(gPlaceDialog, iTextRefB, gRefModel, true);
	DrawParticleStatus(gPlaceDialog, iTextTestA, gTestModel, false);
	DrawParticleStatus(gPlaceDialog, iTextTestB, gTestModel, true);
	SkelDrawButtonOutline(SkelGetDlogCtl(gPlaceDialog, iPlaceClose));
	SetPort(temp);
}

void SetPlaceDone(void)
{
	switch (gModelStatus)
	{
		case sPlaceRefA: 
			gRefModel->sideADone = true;
			if (!gRefModel->sideBDone)
				gModelStatus = sPlaceRefB;
			else if (!gTestModel->sideADone)
				gModelStatus = sPlaceTestA;
			else if (!gTestModel->sideBDone)
				gModelStatus = sPlaceTestB;
			break;

		case sPlaceRefB: 
			gRefModel->sideBDone = true; 
			if (!gRefModel->sideADone)
				gModelStatus = sPlaceRefA;
			else if (!gTestModel->sideADone)
				gModelStatus = sPlaceTestA;
			else if (!gTestModel->sideBDone)
				gModelStatus = sPlaceTestB;
			break;

		case sPlaceTestA: 
			gTestModel->sideADone = true; 
			if (!gTestModel->sideBDone)
				gModelStatus = sPlaceTestB;
			else if (!gRefModel->sideADone)
				gModelStatus = sPlaceRefA;
			else if (!gRefModel->sideBDone)
				gModelStatus = sPlaceRefB;
			break;

		case sPlaceTestB: 
			gTestModel->sideBDone = true;
			if (!gTestModel->sideADone)
				gModelStatus = sPlaceTestA;
			else if (!gRefModel->sideADone)
				gModelStatus = sPlaceRefA;
			else if (!gRefModel->sideBDone)
				gModelStatus = sPlaceRefB;
	}
}

void PlaceDelOne(void)
{
	switch (gModelStatus)
	{
		case sPlaceRefA:
			gRefModel->sideADone = false;
			if (gRefModel->sideAParticles > 0)
				gRefModel->sideAParticles -= 1;
			break;
			
		case sPlaceRefB:
			gRefModel->sideBDone = false;
			if (gRefModel->sideBParticles > 0)
				gRefModel->sideBParticles -= 1;
			break;
			
		case sPlaceTestA:
			gTestModel->sideADone = false;
			if (gTestModel->sideAParticles > 0)
				gTestModel->sideAParticles -= 1;
			break;
			
		case sPlaceTestB:
			gTestModel->sideBDone = false;
			if (gTestModel->sideBParticles > 0)
				gTestModel->sideBParticles -= 1;
	}
}

static void PlaceClear(void)
{
	switch (gModelStatus)
	{
		case sPlaceRefA:
			gRefModel->sideADone = false;
			gRefModel->sideAParticles = 0;
			break;
			
		case sPlaceRefB:
			gRefModel->sideBDone = false;
			gRefModel->sideBParticles = 0;
			break;
			
		case sPlaceTestA:
			gTestModel->sideADone = false;
			gTestModel->sideAParticles = 0;
			break;
			
		case sPlaceTestB:
			gTestModel->sideBDone = false;
			gTestModel->sideBParticles = 0;
	}
}

static pascal Boolean PlaceDlogFilter(DialogPtr dlog, EventRecord *evt, short *item)
{
	if (dlog == gPlaceDialog && evt->what == updateEvt)
	{
		DrawPlaceDialog();
		return false;
	}
	return SkelDlogMapKeyToButton(dlog, evt, item, iPlaceClose, iPlaceClose);
}

static pascal void PlaceDlogSelect(DialogPtr dlog, short item)
{
	if (dlog == gPlaceDialog)
	{
		switch (item)
		{
			case iPlaceClose:
				ClosePlaceDialog();
				return;
				
			case iPlaceDone: SetPlaceDone(); break;
			case iPlaceDelOne: PlaceDelOne(); break;
			case iPlaceClear: PlaceClear(); break;
			case iPlaceRefA: gModelStatus = sPlaceRefA; break;
			case iPlaceRefB: gModelStatus = sPlaceRefB; break;
			case iPlaceTestA: gModelStatus = sPlaceTestA; break;
			case iPlaceTestB: gModelStatus = sPlaceTestB; break;
			default: return;
		}
		BuildSpringsFromParticles(gRefModel);
		BuildSpringsFromParticles(gTestModel);
		KickPictWindow();
		DrawPlaceDialog();
	}
}

static pascal void PlaceDlogClobber(void)
{
	DisposeDialog(gPlaceDialog);
}

Boolean InitPlaceDialog(void)
{
	return SkelDialog(gPlaceDialog, PlaceDlogFilter, PlaceDlogSelect, nil, PlaceDlogClobber);
}

void OpenPlaceDialog(void)
{
	if (!gPictArray.valid)
		DoError("\pPict Array not valid!");
	if (((WindowPeek)gExecDialog)->visible)
		CloseExecDialog();
	gModelStatus = sPlaceRefA;
	ShowWindow(gPlaceDialog);
	SelectWindow(gPlaceDialog);
	DrawPlaceDialog();
	BuildSpringsFromParticles(gRefModel);
	BuildSpringsFromParticles(gTestModel);
	KickPictWindow();
}

void ClosePlaceDialog(void)
{
	if (((WindowPeek)gPlaceDialog)->visible)
	{
		HideWindow(gPlaceDialog);
		gModelStatus = sIdle;
	}
	AdjustMenus(true);
}


/************************************************/
/*		Execution Dialog Routines				*/
/************************************************/

static Boolean TerminationCondition(void)
{
	Boolean b1, b2, b3, b4;
	
	b1 = gRefModel->totalKineticEnergy < gRefModel->kineticLimit;
	b2 = gRefModel->totalForceMagnitude < gRefModel->forceLimit;
	b3 = gTestModel->totalKineticEnergy < gTestModel->kineticLimit;
	b4 = gTestModel->totalForceMagnitude < gTestModel->forceLimit;
	return (b1 && b2 && b3 && b4);
}

static void ComputeAverageWidth(ModelPtr model)
{
	short 	i, rungs;
	Vector	v1, v2, v3, v4, diff;
	Scalar 	area2x, sumArea2x, sumSideA, sumSideB;
	
	rungs = MinShort(model->sideAParticles, model->sideBParticles);	
	if (rungs == 0)
		model->averageWidth = 0.0;
	else if (rungs == 1)
	{
		v1 = model->sideAList[0].position;
		v2 = model->sideBList[0].position;
		VectorDifference(v1, v2, &diff);
		model->averageWidth = Magnitude(diff);
	}
	else
	{
		sumArea2x = 0.0;
		sumSideA = 0.0;
		sumSideB = 0.0;
		for (i = 1; i < rungs; i++)
		{
			v1 = model->sideAList[i-1].position;
			v2 = model->sideAList[i].position;
			v3 = model->sideBList[i-1].position;
			v4 = model->sideBList[i].position;
			area2x = ((v4.x - v1.x) * (v2.y - v3.y)) +
					 ((v3.x - v2.x) * (v4.y - v1.y));
			sumArea2x += (area2x < 0.0 ? -area2x : area2x);
			VectorDifference(v1, v2, &diff);
			sumSideA += Magnitude(diff);
			VectorDifference(v3, v4, &diff);
			sumSideB += Magnitude(diff);
		}
		model->averageWidth = sumArea2x / (sumSideA + sumSideB);
	}
}

void DrawExecDialog(void)
{
	GrafPtr	temp;
	
	GetPort(&temp);
	SetPort(gExecDialog);
	RGBForeColor(&kBlack);
	RGBBackColor(&kWhite);
	if (TerminationCondition())
	{
		SkelSetDlogStr(gExecDialog, iExecStatus, "\pConverged.");
		SysBeep(10);
	}
	else 
		switch (gModelStatus)
		{
			case sPaused: 
				SkelSetDlogStr(gExecDialog, iExecStatus, "\pPaused..."); 
				break;
				
			case sRunning: 
				SkelSetDlogStr(gExecDialog, iExecStatus, "\pRunning..."); 
				break;
	
			default: 
				SkelSetDlogStr(gExecDialog, iExecStatus, "\pHalted...");
		} 
	if (gModelStatus == sRunning)
		SkelSetDlogRadioButtonSet(gExecDialog, iExecStep, iExecCont, iExecCont);
	else if (gModelStatus == sPaused)
		SkelSetDlogRadioButtonSet(gExecDialog, iExecStep, iExecCont, iExecStep);
	else
		SkelSetDlogRadioButtonSet(gExecDialog, iExecStep, iExecCont, 0);
	ComputeAverageWidth(gRefModel);
	ComputeAverageWidth(gTestModel);
	ScalarToPString(gRefModel->averageWidth, gTempStr, "%6.2f");
	SkelSetDlogStr(gExecDialog, iRefWidth, gTempStr);
	ScalarToPString(gTestModel->averageWidth, gTempStr, "%6.2f");
	SkelSetDlogStr(gExecDialog, iTestWidth, gTempStr);
	ScalarToPString(gRefModel->totalKineticEnergy, gTempStr, "%6.5f");
	SkelSetDlogStr(gExecDialog, iRefKECur, gTempStr);
	ScalarToPString(gRefModel->totalForceMagnitude, gTempStr, "%6.5f");
	SkelSetDlogStr(gExecDialog, iRefSFCur, gTempStr);
	ScalarToPString(gTestModel->totalKineticEnergy, gTempStr, "%6.5f");
	SkelSetDlogStr(gExecDialog, iTestKECur, gTempStr);
	ScalarToPString(gTestModel->totalForceMagnitude, gTempStr, "%6.5f");
	SkelSetDlogStr(gExecDialog, iTestSFCur, gTempStr);
	ScalarToPString(gRefModel->kineticLimit, gTempStr, "%6.5f");
	SkelSetDlogStr(gExecDialog, iRefKELim, gTempStr);
	ScalarToPString(gRefModel->forceLimit, gTempStr, "%6.5f");
	SkelSetDlogStr(gExecDialog, iRefSFLim, gTempStr);
	ScalarToPString(gTestModel->kineticLimit, gTempStr, "%6.5f");
	SkelSetDlogStr(gExecDialog, iTestKELim, gTempStr);
	ScalarToPString(gTestModel->forceLimit, gTempStr, "%6.5f");
	SkelSetDlogStr(gExecDialog, iTestSFLim, gTempStr);
	SkelDrawButtonOutline(SkelGetDlogCtl(gExecDialog, iExecStop));
	SetPort(temp);
}

static pascal Boolean ExecDlogFilter(DialogPtr dlog, EventRecord *evt, short *item)
{
	if (dlog == gExecDialog && evt->what == updateEvt)
	{
		DrawExecDialog();
		return false;
	}
	return SkelDlogMapKeyToButton(dlog, evt, item, iExecStop, iExecStop);
}

static pascal void ExecDlogSelect(DialogPtr dlog, short item)
{
	if (dlog == gExecDialog)
		switch (item)
		{
			case iExecStop:
				HideWindow(dlog);
				gModelStatus = sIdle;
				AdjustMenus(true);
				break;

			case iExecStep:
				gModelStatus = sPaused;
				SkelSetDlogRadioButtonSet(gExecDialog, iExecStep, iExecCont, iExecStep);
				KickPictWindow();
				AdjustMenus(true);
				break;

			case iExecCont:
				gModelStatus = sRunning;
				SkelSetDlogRadioButtonSet(gExecDialog, iExecStep, iExecCont, iExecCont);
				KickPictWindow();
				AdjustMenus(true);
		}
}

static pascal void ExecDlogClobber(void)
{
	DisposeDialog(gExecDialog);
}

Boolean InitExecDialog(void)
{
	return SkelDialog(gExecDialog, ExecDlogFilter, ExecDlogSelect, nil, ExecDlogClobber);
}

void OpenExecDialog(void)
{
	if (!gPictArray.valid)
		DoError("\pPict Array not valid!");
	if (!((WindowPeek)gGradWindow)->visible)
		OpenGradWindow();
	if (!((WindowPeek)gGradWindow)->visible)
		return;
	if (((WindowPeek)gPlaceDialog)->visible)
		ClosePlaceDialog();
	gModelStatus = sRunning;
	ShowWindow(gExecDialog);
	SelectWindow(gExecDialog);
	gRefModel->totalKineticEnergy = 999.0;
	gRefModel->totalForceMagnitude = 999.0 ;
	gTestModel->totalKineticEnergy = 999.0 ;
	gTestModel->totalForceMagnitude = 999.0 ;
	DrawExecDialog();
}

void CloseExecDialog(void)
{
	if (((WindowPeek)gExecDialog)->visible)
	{
		HideWindow(gExecDialog);
		gModelStatus = sIdle;
	}
	AdjustMenus(true);
}


/************************************************/
/*		Model Drawing Routines					*/
/************************************************/

static void DrawParticle(Particle particle, short wsize)
{
	short j, k, d;
	Rect r;

	j = (short)particle.position.x;
	k = (short)particle.position.y;
	if (gOptions->showParticleWindows)
	{
		d = wsize / 2;
		SetRect(&r, j - d, k - d, j + d, k + d);
		FrameRect(&r);
	}
	else
	{
		d = PARTICLE_RADIUS;
		SetRect(&r, j - d, k - d, j + d, k + d);
		FrameOval(&r);
	}
}

void DrawModel(ModelPtr model, WindowPtr window)
{
	GrafPtr 		temp;
	Vector 			v;
	short 			i;
	
	GetPort(&temp);
	SetPort(window);
	RGBBackColor(&kBlack);
	RGBForeColor(&model->particleColor);
	for (i = 0; i < model->sideAParticles; i++)
		DrawParticle(model->sideAList[i], model->windowSize);
	for (i = 0; i < model->sideBParticles; i++)
		DrawParticle(model->sideBList[i], model->windowSize);
	for (i = 0; i < model->numberOfSprings; i++)
	{
		v = model->springList[i].from->position;
		MoveTo(v.x, v.y);
		switch(model->springList[i].type)
		{
			case tRung: RGBForeColor(&model->rungColor); break;
			case tEdge: RGBForeColor(&model->edgeColor); break;
			case tBrace: RGBForeColor(&model->braceColor); break;
			default: RGBForeColor(&kBlack);
		}
		v = model->springList[i].to->position;
		LineTo(v.x, v.y);
	}
	SetPort(temp);
}


/************************************************/
/*		Model Execution Routines				*/
/************************************************/

static Vector InternalForces(ModelPtr model, Particle *p, Vector *sumForces)
{
	short		i, n;
	Spring		*spring;
	Particle	*other;
	Scalar		stiffness, friction;
	Scalar		relDist, relSpeed, springForce, frictForce;
	Vector 		relPos, relVel, relDir, force;
	
	sumForces->x = 0.0;
	sumForces->y = 0.0;
	n = model->numberOfSprings;
	spring = model->springList;
	for (i = 0; i < n; i++)
	{
		if (spring[i].from == p)
			other = spring[i].to;
		else if (spring[i].to == p)
			other = spring[i].from;
		else continue;
		switch (spring[i].type)
		{
			case tRung:  
				stiffness = model->rungStiffness;  
				friction = model->rungFriction; 
				break;
			
			case tEdge:  
				stiffness = model->edgeStiffness;  
				friction = model->edgeFriction; 
				break;
			
			case tBrace: 
				stiffness = model->braceStiffness; 
				friction = model->braceFriction; 
				break;
			
			default:     
				stiffness = 0.0;  			 
				friction = 0.0;
		}
		VectorDifference(p->position, other->position, &relPos);
		Direction(relPos, &relDir);
		relDist = Magnitude(relPos);
		if  (relDist == 0.0) relDist = 1.0;
		springForce = (spring[i].lengthAtRest - relDist) * stiffness;
		VectorDifference(p->velocity, other->velocity, &relVel);
		relSpeed = DotProduct(relVel, relDir);
		frictForce = relSpeed * friction;
		VectorByScalar(relDir, (springForce - frictForce), &force);
		VectorSum(*sumForces, force, sumForces);
	}
}

static void ExternalForces(ModelPtr model, Particle *p, Vector *sumForces)
{
	long	x, xmin, xmax, y, ymin, ymax, numPixels;
	ValPtr	grad;
	Scalar	wsize, relDist2, pixMass, massSum, aveMass;
	Vector 	pixPos, relPos, forceSum, force, damping, gravity;
	
	grad = (ValPtr)LockArray(&gGradArray);
	wsize = (Scalar)(model->windowSize);
	xmin = (long)(p->position.x - (wsize / 2.0) + 0.5);
	if (xmin < 0) xmin = 0;
	xmax = (long)(p->position.x + (wsize / 2.0) + 0.5);
	if (xmax >= gPictWidth) xmax = gPictWidth - 1;
	ymin = (long)(p->position.y - (wsize / 2.0) + 0.5);
	if (ymin < 0) ymin = 0;
	ymax = (long)(p->position.y + (wsize / 2.0) + 0.5);
	if (ymax >= gPictHeight) ymax = gPictHeight - 1;
	numPixels = 0;
	massSum = 0.0;
	SetVector(&forceSum, 0.0, 0.0);
	for (y = ymin; y <= ymax; y++)
		for (x = xmin; x <= xmax; x++)
		{
			SetVector(&pixPos, (Scalar)x, (Scalar)y);			
			VectorDifference(pixPos, p->position, &relPos);
			relDist2 = relPos.x * relPos.x + relPos.y * relPos.y;
			pixMass = (Scalar)grad[y * gPictWidth + x];
			massSum += pixMass;
			VectorByScalar(relPos, pixMass / (relDist2 + 1.0), &force);
			VectorSum(forceSum, force, &forceSum);
			numPixels++;
		}
	UnlockArray(&gGradArray);
	if (numPixels == 0L) numPixels = 1L;
	aveMass = massSum / (Scalar)numPixels;
	if (aveMass == 0.0) aveMass = 1.0;
	VectorByScalar(p->velocity, gGlobalParams->viscousFriction, &damping);
	VectorByScalar(forceSum, model->particleMass, &forceSum);
	VectorByScalar(forceSum, gGlobalParams->gravityConstant / aveMass, &gravity);
	VectorByScalar(p->velocity, gGlobalParams->viscousFriction, &damping);
	VectorDifference(gravity, damping, sumForces);
}

void UpdateModel(ModelPtr model)
{
	Scalar			totalEnergy = 0.0;
	Scalar			totalForces = 0.0;
	Scalar			timeIncr, speed;
	Vector			iforce, eforce, vec;
	Particle		*p;
	short			i;

	timeIncr = gGlobalParams->timeIncrement;
	for (i = 0; i < model->sideAParticles; i++)
	{
		p = &(model->sideAList[i]);
		if (gOptions->firstRungDoesntMove && i == 0)
			SetVector(&(p->force), 0.0, 0.0);
		else
		{
			InternalForces(model, p, &iforce);
			ExternalForces(model, p, &eforce);
			VectorSum(iforce, eforce, &(p->force));
		}
		totalForces += Magnitude(p->force);
	}
	for (i = 0; i < model->sideBParticles; i++)
	{
		p = &(model->sideBList[i]);
		if (gOptions->firstRungDoesntMove && i == 0)
			SetVector(&(p->force), 0.0, 0.0);
		else
		{
			InternalForces(model, p, &iforce);
			ExternalForces(model, p, &eforce);
			VectorSum(iforce, eforce, &(p->force));
		}
		totalForces += Magnitude(p->force);
	}
	for (i = 0; i < model->sideAParticles; i++)
	{
		p = &(model->sideAList[i]);
		if (gOptions->firstRungDoesntMove && i == 0)
			SetVector(&(p->velocity), 0.0, 0.0);
		else
		{
			VectorByScalar(p->force, timeIncr / model->particleMass, &vec);
			VectorSum(p->velocity, vec, &(p->velocity));
		}
		speed = Magnitude(p->velocity);
		totalEnergy += (Scalar)(0.5 * model->particleMass * speed * speed);
	}
	for (i = 0; i < model->sideBParticles; i++)
	{
		p = &(model->sideBList[i]);
		if (gOptions->firstRungDoesntMove && i == 0)
			SetVector(&(p->velocity), 0.0, 0.0);
		else
		{
			VectorByScalar(p->force, timeIncr / model->particleMass, &vec);
			VectorSum(p->velocity, vec, &(p->velocity));
		}
		speed = Magnitude(p->velocity);
		totalEnergy += (Scalar)(0.5 * model->particleMass * speed * speed);
	}
	for (i = 0; i < model->sideAParticles; i++)
	{
		p = &(model->sideAList[i]);
		if (!(gOptions->firstRungDoesntMove && i == 0))
		{
			VectorByScalar(p->velocity, timeIncr, &vec);
			VectorSum(p->position, vec, &(p->position));
			if (Constrain(&(p->position), gPictRect))
				SetVector(&(p->velocity), 0.0, 0.0);
		}
	}
	for (i = 0; i < model->sideBParticles; i++)
	{
		p = &model->sideBList[i];
		if (!(gOptions->firstRungDoesntMove && i == 0))
		{
			VectorByScalar(p->velocity, timeIncr, &vec);
			VectorSum(p->position, vec, &(p->position));
			if (Constrain(&(p->position), gPictRect))
				SetVector(&(p->velocity), 0.0, 0.0);
		}
	}
	model->totalKineticEnergy = totalEnergy;
	model->totalForceMagnitude = totalForces;
}

void RelaxModel(void)
{
	BuildSpringsFromParticles(gRefModel);
	BuildSpringsFromParticles(gTestModel);
	KickPictWindow();
}

void ExecuteOneCycle(void)
{
	UpdateModel(gRefModel);
	UpdateModel(gTestModel);
	ShowPictWindow();
	ShowGradWindow();
	DrawExecDialog();
}


void LoadModel(void)
{
	StandardFileReply 	reply;
	SFTypeList			typeList;
	short				numTypes, srcFile, oldRefNum;
	FInfo				myInfo;
	
	typeList[0] = kModelType;
	numTypes = 1;
	StandardGetFile(nil, numTypes, typeList, &reply);
	if (!reply.sfGood)
		return;	
	if (FSpGetFInfo(&(reply.sfFile), &myInfo) != noErr)
	{
		DoWarning("\pCannot get finder info!");
		return;
	}
	if (myInfo.fdCreator == kCreatorID)
	{
		if ((srcFile = FSpOpenResFile(&(reply.sfFile), fsRdWrPerm)) == -1)
			DoWarning("\pCannot open model file!");
		else
		{
			oldRefNum = CurResFile();
			UseResFile(srcFile);
			LoadPreference((Ptr)gGlobalParams, sizeof(GlobalParams), kGlobalParamType, kGlobalParams);
			LoadPreference((Ptr)gRefModel, sizeof(Model), kModelType, kRefModel);
			LoadPreference((Ptr)gTestModel, sizeof(Model), kModelType, kTestModel);
			CloseResFile(srcFile);
			UseResFile(oldRefNum);
		}
	}
	BuildSpringsFromParticles(gRefModel);
	BuildSpringsFromParticles(gTestModel);
	KickPictWindow();
}

void SaveModelAs(void)
{
	StandardFileReply 	reply;
	short				myRefNum, oldRefNum;
	OSErr				err;
	
	StandardPutFile("\pSave Model...", "\p", &reply);
	if (!reply.sfGood)
		return;
	SetCursor(*gWatchCursor);
	FSpCreateResFile(&(reply.sfFile), kCreatorID, kModelType, 0);
	if ((myRefNum = FSpOpenResFile(&(reply.sfFile), fsRdWrPerm)) == -1)
		DoError("\pCannot open resource file!");
	UseResFile(myRefNum);
	if ((err = ResError()) != noErr)
		DoError("\pCannot use resource fork!");
	SavePtrRsrc((Ptr)gGlobalParams, sizeof(GlobalParams), kGlobalParamType, kGlobalParams, "\pGlobal Parameters");
	SavePtrRsrc((Ptr)gRefModel, sizeof(Model), kModelType, kRefModel, "\pReference Model");
	SavePtrRsrc((Ptr)gTestModel, sizeof(Model), kModelType, kTestModel, "\pTest Model");
	UpdateResFile(myRefNum);
	CloseResFile(myRefNum);
	SetCursor(&arrow);
}