mercredi 4 mars 2015

C fmod function: Floating point error and optimization


I'm trying to calculate the true course from one point to anoter on the surface of the earth in as few CPU cycles as possible. The result should be a double 0 <= tc < 360, however in a few special cases i get the result 360 (should be reported as 0). I realize that this is due to machine precision when working with fmod and floating point number, but what will be the most efficient workaround of the problem?



#include <stdio.h>
#include <math.h>

#define EPS 1e-15 // EPS a small number ~ machine precision
#define R2D 57.295779513082320876798154814105 //multiply radian with R2D to get degrees
#define D2R 0.01745329251994329576923690768489 //multiply degrees with D2R to get radians
#define TWO_PI 6.283185307179586476925286766559 //2*Pi

/*----------------------------------------------------------------------------
* Course between points

* We obtain the initial course, tc1, (at point 1) from point 1 to point 2
* by the following. The formula fails if the initial point is a pole. We can
* special case this with as IF statement
----------------------------------------------------------------------------
Implementation
Argument 1: INPUT - Pointer to double containing Latitude of point 1 in degrees
Argument 2: INPUT - Pointer to double containing Longitude of point 1 in degrees
Argument 3: INPUT - Pointer to double containing Latitude of point 2 in degrees
Argument 4: INPUT - Pointer to double containing Longitude of point 2 in degrees

RETURN: Double containing initial course in degrees from point1 to point 2
--------------------------------------------------------------------------*/
double _stdcall CourseInitial (double *lat1, double *lon1, double *lat2, double *lon2)
{
double radLat1 = D2R * *lat1;
double radLon1 = D2R * *lon1;
double radLat2 = D2R * *lat2;
double radLon2 = D2R * *lon2;
double tc = 0;

if (cos(radLat1) < EPS) { // EPS a small number ~ machine precision
if (radLat1 > 0) {
tc = 180; // starting from N pole
} else {
tc = 0; // starting from S pole
}
} else {
// Calculate true course [180, 540]
tc = R2D * atan2(sin(radLon2-radLon1),
cos(radLat1) * tan(radLat2) - sin(radLat1) * cos(radLon2-radLon1)
) + 360;
}

//Return 0 <= true course < 360
return fmod(tc, 360);
}

int main(void)
{
double lat1 = 89;
double lon1 = 17;
double lat2 = 68;
double lon2 = -163;
double tc = 0;
tc = CourseInitial(&lat1, &lon1, &lat2, &lon2);
printf("The course from point 1 to 2 is: %.5f", tc);

return 0;

}


Output:



The course from point 1 to 2 is: 360.00000



Aucun commentaire:

Enregistrer un commentaire