Re: Problems with multiplications of doubles and/or floats

From: J.K. Becker (jkbecker_at_becker.de)
Date: 04/15/04 

Date: Thu, 15 Apr 2004 11:11:15 +0200

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[i] != NO_NB && ElleNodeIsActive( index ) )
       {
         ElleNodePosition( nb[i], & 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 );
   }
}

Relevant Pages

  • Re: (N+1)^N/N^N
    ... >and yet it seems to be FOUNDATIONAL in many calculations. ... David C. Ullrich ... Prev by Date: ...
    (sci.math)
  • Re: Totals Queries
    ... or you can uncheck the Show checkbox on the rate if ... you only need to use it in calculations. ... Prev by Date: ...
    (microsoft.public.access.queries)
  • Format a field for Time, but as duration?
    ... What is the best way to format the field in the tabel so that the data will ... be used to summarise without adding calculations to it? ... Prev by Date: ...
    (microsoft.public.access.tablesdbdesign)
  • Re: Varying a Latitude/Longitude
    ... longitude. ... I'll bet one of them can do these calculations, ... jumped out at me when I very quickly skimmed over the one-liner ... Prev by Date: ...
    (comp.lang.perl.misc)
  • add two matrix
    ... or in another exemple: ... but this way of making the calculations, with For loops is very slow, ... Prev by Date: ...
    (comp.soft-sys.matlab)