This is the function where the error occurs. As you can see there are a
lot of calculations in it, they all work fine except the one mentioned
(which seems to be the easiest one).
Fx = un.x * dEdp;
Fy = un.y * dEdp;
This always gives the wrong results... Everything before it works and
everything after it works too (except that everything is calculated with
the wrong values for Fx and Fy).
But I guess it is too complex for a quick fix...
Thanks anyway
Jens
int GBE_MoveNode( int index, Coords * movedir )
{
int i, l, k, moveflag = 0;
double E1, E2, E3, E4, vangle[3],lenL[3] , lenK[3];
double truetimestep = 3.1536e10;
float Fx,Fy;
double lenP, lenV, lenF, mobility1, mobility2, mobility3 , switchd;
Coords p1, p2, gvector, un, sigma1, sigma2, sigma3, L[3], F;
Coords newxy, xynb, xy, prev, V;
float dEdp;
int nb[3];
char cc[255];
mobility1 = 1e-12;
mobility2 = 1e-12;
mobility3 = 1e-12;
switchd = ElleSwitchdistance() / 10;
//Get position of first node
ElleNodePosition( index, & p1 );
ElleNodePrevPosition( index, & prev );
//this is the first energy we need: E((x+dx),y)
p2.x = p1.x + switchd;
p2.y = p1.y;
ElleSetPosition( index, & p2 );
GetNodeEnergy( index, & p2, & E1 );
//this is the second energy we need:E((x-dx),y)
p2.x = p1.x - switchd;
p2.y = p1.y;
ElleSetPosition( index, & p2 );
GetNodeEnergy( index, & p2, & E2 );
//this is the third energy we need: E(x,(y+dy))
p2.x = p1.x;
p2.y = p1.y + switchd;
ElleSetPosition( index, & p2 );
GetNodeEnergy( index, & p2, & E3 );
//this is the fourth energy we need: E(x,(y-dy))
p2.x = p1.x;
p2.y = p1.y - switchd;
ElleSetPosition( index, & p2 );
GetNodeEnergy( index, & p2, & E4 );
//Reset node position to starting values
ElleSetPosition( index, & p1 );
ElleSetPrevPosition( index, & prev );
//So now we can calculate the gradient vector (P)
gvector.x = ( E1 - E2 ) / ( 2 * switchd );
gvector.y = ( E3 - E4 ) / ( 2 * switchd );
if ( gvector.x == 0.0 && gvector.y == 0.0 )
{
movedir->x = 0.0;
movedir->y = 0.0;
return ( 0 );
}
else
{
//reverse gradient
gvector.x *= -1;
gvector.y *= -1;
//length of gradient
lenP = sqrt( ( gvector.x * gvector.x ) + ( gvector.y * gvector.y ) );
//unit vector
un.x = gvector.x / lenP;
un.y = gvector.y / lenP;
dEdp = ( ( gvector.x / gvector.y ) + ( gvector.y / gvector.x ) ) /
lenP;
////////
//////// wrong results in this calculation
////////
Fx = un.x * dEdp;
Fy = un.y * dEdp;
////////
////////
////////
lenF = sqrt( ( Fx * Fx ) + ( Fy * Fy ) );
ElleNeighbourNodes( index, nb );
for ( i = 0, k = 0; i < 3; i++ )
{
if ( nb != NO_NB && ElleNodeIsActive( index ) )
{
ElleNodePosition( nb, & xynb );
L[k].x = p1.x - xynb.x;
L[k].y = p1.y - xynb.y;
lenL[k] = sqrt( ( L[k].x * L[k].x ) + ( L[k].y * L[k].y ) );
vangle[k] = fabs(90- ( acos( ( Fx * L[k].x + Fy * L[k].y ) / (
lenL[k] * lenF ) ) ));
lenK[k] = lenL[k] * fabs( sin( vangle[k] * ( 180 / 3.1415926535
) ) );
k++;
}
}
sigma1.x = Fx / ( lenK[0] + ( mobility1 / mobility2 ) * lenK[1] );
sigma1.y = Fy / ( lenK[0] + ( mobility1 / mobility2 ) * lenK[1] );
V.x = sigma1.x * (mobility1/switchd/switchd);
V.y = sigma1.y * (mobility1/switchd/switchd);
lenV = sqrt( ( V.x * V.x ) + ( V.y * V.y ) );
movedir->x = p1.x - V.x;
movedir->y = p1.y - V.y;
moveflag = 1;
return ( moveflag );
}
}
)