Leeway calibration (Part 2)

In the previous post, I have presented a method to calibrate the leeway of a sailboat. Here, I present the details of the required calculations, in the form of a C-language routine. Beware that if you want to program this in Excel, the Excel ATAN2 function needs the arguments in reverse order : (x, y) instead of (y,x).

#define PI 3.14159265
#define DEG_TO_RAD ((double)(PI/180.0))
#define RAD_TO_DEG ((double) (180.0/PI))

// inputs from downwind leg
double sog2;  // GPS SOG of downwind leg
double cogm2;  // GPS COG (corrected to magnetic) of downwind leg
double stw2;  // measured boat speed (also real STW) of downwind leg
double heading2; // heading of downwing leg

// inputs from upwind leg
double sog1;  // GPS SOG of upwind leg
double cogm1;  // GPS COG (corrected to magnetic) of upwind leg
double heading1; // heading of upwind leg

// outputs
double leeway;  // leeway angle of upwind leg;
double soc;   // speed of current (common to both legs)
double doc;   // direction of current (common to both legs)
double stw1;  // real speed through water of upwind leg;
double dir_stw1;    // direction of stw1 in compass coordinates

// intermediate values
double cogm2_rad; // cogm2 in radians in cartesian coordinates
double sog2x;  // horizontal component of sog2 in cartesian coordinates
double sog2y;  // vertical component of sog2 in cartesian coordinates
double heading2_rad; // heading2 in radians in cartesian coordinates
double stw2x;  // horizontal component of stw2 in cartesian coordinates
double stw2y;  // vertical component of stw2 in cartesian coordinates
double socx;  // horizontal component of soc in cartesian coordinates
double socy;  // vertical component of soc in cartesian coordinates
double doc_rad;  // doc in radians in cartesian coordinates
double cogm1_rad; // cogm1 in radians in cartesian coordinates
double sog1x;  // horizontal component of sog1 in cartesian coordinates
double sog1y;  // vertical component of sog1 in cartesian coordinates
double stw1x;  // horizontal component of stw1 in cartesian coordinates
double stw1y;  // vertical component of stw1 in cartesian coordinates
double dir_stw1_rad // direction of stw1 in radians in cartesian coordinates


cogm2_rad = (90.0 - cogm2) * DEG_TO_RAD;
sog2x = sog2 * cos(cogm2_rad);
sog2y = sog2 * sin(cogm2_rad);
heading2_rad = (90.0 - heading2) * DEG_TO_RAD;
stw2x = stw2 * cos(heading2_rad);
stw2y = stw2 * sin(heading2_rad);
socx = sog2x - stw2x;
socy = sog2y - stw2y;

/******* this section is optional *************
soc = sqrt(socx * socx + socy * socy); 
doc_rad = atan2(socy, socx);
if(isnan(doc_rad))   // singularity
{
 if(socy < 0.0)
  doc = 180.0;
 else doc = 0.0;
}
else
{
 doc = 90.0 - doc_rad * RAD_TO_DEG;
 if(doc < 0.0)
  doc += 360.0;
 else if(doc >= 360.0)
  doc -= 360.0;
}
************************************************/

cogm1_rad = (90.0 - cogm1) * DEG_TO_RAD;
sog1x = sog1 * cos(cogm1_rad);
sog1y = sog1 * sin(cogm1_rad);
stw1x = sog1x - socx;
stw1y = sog1y - socy;

/******* this section is optional *************
stw1 = sqrt(stw1x * stw1x + stw1y * stw1y);
***********************************************/

dir_stw1_rad = atan2(stw1y, stw1x);
if(isnan(dir_stw1_rad))   // singularity
{
 if(stw1y < 0.0)
  dir_stw1 = 180.0;
 else dir_stw1 = 0.0;
}
else
{
 dir_stw1 = 90.0 - dir_stw1_rad * RAD_TO_DEG;
 if(dir_stw1 < 0.0)
  dir_stw1 += 360.0;
 else if(dir_stw1 >= 360.0)
  dir_stw1 -= 360.0;
}
leeway = dir_stw1 - heading1;