Trilateration (2D) algorithm implementation

I am trying to implement the Trilateration process in 2D. The wikipedia article relating to this: Tilateration

I have found a nice question here at this site, where the algorithm is well explained:artifical intelligence

After all, I tried to implement the algorithm in c++. Unfortunately I faced some problems... Let's see my implementation. It is only a function: The first inputs are three vector, each representing a 2D point with X,Y coordinates. The other (r1,r2,r3) input variables stand for the distance/radius of each point.

#include <iostream>
#include <fstream>
#include <sstream>
#include <math.h> 
#include <vector>
using namespace std;

std::vector<double> trilateration(double point1[], double point2[], double point3[], double r1, double r2, double r3) {
    std::vector<double> resultPose;
    //unit vector in a direction from point1 to point 2
    double p2p1Distance = pow(pow(point2[0]-point1[0],2) + pow(point2[1]-point1[1],2),0.5);
    double exx = (point2[0]-point1[0])/p2p1Distance;
    double exy = (point2[1]-point1[1])/p2p1Distance;
    //signed magnitude of the x component
    double ix = exx*(point3[0]-point1[0]);
    double iy = exy*(point3[1]-point1[1]);
    //the unit vector in the y direction. 
    double eyx = (point3[0]-point1[0]-ix*exx)/pow(pow(point3[0]-point1[0]-ix*exx,2) + pow(point3[1]-point1[1]-iy*exy,2),0.5);
    double eyy = (point3[1]-point1[1]-iy*exy)/pow(pow(point3[0]-point1[0]-ix*exx,2) + pow(point3[1]-point1[1]-iy*exy,2),0.5);
    //the signed magnitude of the y component
    double jx = eyx*(point3[0]-point1[0]);
    double jy = eyy*(point3[1]-point1[1]);
    //coordinates
    double x = (pow(r1,2) - pow(r2,2) + pow(p2p1Distance,2))/ (2 * p2p1Distance);
    double y = (pow(r1,2) - pow(r3,2) + pow(iy,2) + pow(jy,2))/2*jy - ix*x/jx;
    //result coordinates
    double finalX = point1[0]+ x*exx + y*eyx;
    double finalY = point1[1]+ x*exy + y*eyy;
    resultPose.push_back(finalX);
    resultPose.push_back(finalY);
    return resultPose;
}

As I mentioned I followed this article. I am of the opinion that the problem lies at the part where the y coordinate is calculated. I am also not sure about last part, where I calculate finalX, finalY...

My main function is the following:

int main(int argc, char* argv[]){
    std::vector<double> finalPose;
    double p1[] = {4.0,4.0};
    double p2[] = {9.0,7.0};
    double p3[] = {9.0,1.0};
    double r1,r2,r3;
    r1 = 4;
    r2 = 3;
    r3 = 3.25;
    finalPose = trilateration(p1,p2,p3,r1,r2,r3);
    cout<<"X:::  "<<finalPose[0]<<endl;
    cout<<"Y:::  "<<finalPose[1]<<endl; 
    //x = 8, y = 4.1

}

The result should be around X~8 and Y~4.1, but I got X = 13.5542 and Y=-5.09038

So my problem is and question is: I have problem with dividing the calculations for x and y. I think I could solve the algorithm till x, after that I have problems with calculating y.

The calculation is the following for y: y = (r12 - r32 + i2 + j2) / 2j - ix / j

I do not know which i and j I should use here since I have two i (ix,iy) and two j(jx,jy). As you can see I used iy and jy but at the end of the line I used ix due to multiplication with x. Thanks in advance!

It's a little unclear, and perhaps incorrect, in the linked SO answer that the values of i and j are scalar values and computed a bit differently than the other vector quantities. More explicitly you should have:

i = e_x · (P3 - P1) = e_xx (P3_x - P1_x) + e_xy (P3_y - P1_y) = ix + iy

j = e_y · (P3 - P1) = e_yx (P3_x - P1_x) + e_yy (P3_y - P1_y) = jx + jy

Note that · is the dot product of two vectors here. Thus, in your code there should be no ix, iy, jx or jy.

Also, in your calculation of y you should change the denominator of /2*j to:

 / (2*j)

otherwise you are multiplying by j instead of dividing. Making these changes gives me the result of [7.05, 5.74] which is closer to your expected values.