Calculating in-between vector angles

Question

I have an elliptical particles system and I am trying to calculate the angle between the major axis of each particle and the vector of the distance from the center of the next particle. This picture should illustrate nicely, I hope.

In this example I initialize two particles at (0,0) and (10,-10). The angle of each particle relative to the Cartesian system is 45 degrees. So these 2 are parallel to each other. I want to calculate the angles that I have marked as 90 degree angles in the picture.

For this, I have made this function, using the dot product:

    double temp = 0;
    double result = 0;
    double cosPhi = 0;
    double dotProd = 0;
    double cosTheta = cos(current->getTheta());
    //  If Theta = PI/2, there is a tiny difference between this rounded value
    // and the exact value of PI/2 due to double precision representation. This
    // gives cos(PI/2) = 6*10^-17 instead of 0. So, manually set cos(PI/2) = 0.
    if(almostEqual(cosTheta,0))
        cosTheta = 0;
    double cellLengthX = current->getLength()*cosTheta;
    double cellLengthY = current->getLength()*sin(current->getTheta());
    double dist = Distance(current, next);
    if(dist == 0)
        return 0.0;
    //cout << "Dist " << dist << endl;
    // calculate the vector of distance 
    // next - current, because the end of the vector is the next cell,
    // while the start of the vector is the current cell.
    double distX = next->getCurrX() - current->getCurrX();
    double distY = next->getCurrY() - current->getCurrY();
    //cout <<"DistX " << distX <<" DistY " << distY << endl;
    // Dot product
    dotProd = cellLengthX*distX + cellLengthY*distY;
    //cout << "DotP " << dotProd << endl; 
    // angle from formula : vector_a*vector_b = |a||b|cos(phi)
    cosPhi = dotProd/(current->getLength()*dist);
    cout << "aCos " << acos(cosPhi) << endl;
    if(cosPhi == 1) // 0 degrees
    {
        return 0;
    }
    else if(cosPhi == -1) // 180 degrees 
    {
        return M_PI;
    }
    else if(cosPhi == 0) // 90 or 270 degrees
    {
        if(distX>0 || distY>0)
            return M_PI_2; // 90
        else
            return 3*M_PI_2; // 270
    }
    else if(-distX<=0 && -distY <0) // 1st quadrant -distX<=0
    {
        //cout << "Angle 1st :" << acos(cosPhi) << endl;
        return acos(cosPhi);
    }
    else if(-distX>=0 && -distY >0) // 3rd quadrant    -distX>=0
    {
        temp = acos(cosPhi);
        result = M_PI_2 + temp;
        //cout << "Angle 3rd :" << result << endl;
        return result;
    }
    else if(-distX >0 && -distY <=0) // 2nd quadrant  -distY <=0
    {
        //cout << "Angle 2nd :" << acos(cosPhi) << endl;
        return acos(cosPhi);
    }
    else if(-distX <0 && -distY >=0) // 4th quadrant    -distY >=0
    {
        temp = acos(cosPhi);
        result = 3*M_PI_2 + temp;
        //cout << "Angle 4th :" << result << endl;
        return result;
    } 

This method worked in principle, but there is something wrong when I judge in which quadrant to assign the return of acos(cosPhi). As I have read from others on SO (and found out myself) there is a problem with acos because you don't know in which quadrant to assign the final angle.

So I made a method using atan2(y,x):

    double x1, y1, x2, y2, phi1, phi2, theta, omega;
    omega = 0;
    theta = current->getTheta();
    phi1 = 0;
    phi2 = 0;
    x1 = current->getLength()*cos(theta);
    y1 = current->getLength()*sin(theta);
    x2 = next->getCurrX() - current->getCurrX();
    y2 = next->getCurrY() - current->getCurrY();
    //phi1 = atan2(y1,x1);
    //phi2 = atan2(y2,x2);
    omega = atan2(y2-y1,x2-x1);//phi2 - phi1;
    cout << phi1 << "   " << phi2 <<  "   " << omega << endl;
    return omega;

But what I get back is -90 degrees for the first particle (I guess it's OK) and -180 degrees for the second particle which is not OK. Any help would be really appreciated.

Answer 1

It's hard to be sure, based on your description, but... Try this:

double x2, y2, theta, omega;

theta = current->getTheta();

x2 = next->getCurrX() - current->getCurrX();
y2 = next->getCurrY() - current->getCurrY();

omega = atan2(y2,x2) - theta;
cout << omega << endl;
return omega;