P2 - P1 will give you a direction vector but you also need an up vector or else you have infinite number of solutions. Imagine you airplane rolling around its tail to nose axis. Most libraries give you a lookat function which returns a rotation matrix which can be converted to a quaternion. In DirectX you can chain XMQuaternionRotationMatrix(XMMatrixLookAtLH(..)); If you want faster code search for better implementations.

3 hours ago, mmakrzem said:

How do I calculate this quaternion given P1(x,y,z) and P2(x,y,z) ?

That's not possible - you can't specify orientation from a single direction alone - you need at least two directions.

You could use an up vector for the second, which causes flipping issues if the other direction becomes close to that. If the up vector is editable by the user, e.g. by dragging a point around, flipping issues can be avoided.

You could use euler angles to represent orientation like anyone else does.

You could use your approach to get a rotation of the object like this: P2-P1 represents target front direction of object. Calculate axis and angle rotation between current and target front dir and apply it to the object. So you can modify an existing orientation with a single direction, and there are no flipping issues, but there is no direct control of the twist around the front dir.

in your example you give the distance between the tail and nose is 10.  Does your p1 and p2 define the axis of alignment or does the air-plane stretch?  Otherwise it seems you're just readjusting the axis of orientation.  Like @JoeJ mentioned you'll need to have an 'up' Axis, otherwise you'll achieve the result your looking for but each time the plane will seem to be rotated about it's relative z axis differently.  And like Ivan says, does the library you're using support lookAt() functionality?  Or is it a variation of the lookAt function you're attempting to calculate?

Would something like this work??

p3 = p2.sub( p1 ).normalize()  // gives you a normalised axis of direction from p2 to p1.

// assuming the tail of the plane is middle of the air-plane creation model i.e At( 0 , 0 , 0 ) in the modelling program used to make the air-plane.
airplane.position = p2;
airplane.lookAt( p3 ); // that will adjust the lootAt axis accordingly.  

The syntax will be different of course.

If I've understood this correctly I'm going to assume the vector that goes from P2 to P1 is the new direction the object is supposed to face.
Let's call this vector: V1 = P2 - P1.

Then we assume, as others have already said, that the model has it's default orientation along the "up" vector in world coordinates or at least along one cardinal axis, X, Y or Z.
Let's assume Z so that our original direction is: V0 = Z = [0, 0, 1]
 

Now to calculate a quaternion Q that will align the model currently facing V0 with V1.

I had this same problem a while ago and found this:
http://stackoverflow.com/questions/1171849/finding-quaternion-representing-the-rotation-from-one-vector-to-another
That I adapted into my own code which written in pseudo C++ would be like:

Quaternion alignDirection(Vector3 V0, Vector3 V1){
    Vector3 A = V0.cross(V1);
    
    float length_0 = V0.norm();
    float length_1 = V1.norm();
    float dp = V0.dot(V1);
    
    float w = sqrt(length_0*length_0 * length_1*length_1) + dp;

    Quaternion Q(w, A.x, A.y, A.z);
    Q.normalize();
    return Q;
}

Take note on the quaternion initialization in your framework as they sometimes differ.
A is the "vector" of the quaternion values.

I don't get how it works. But I tested it before I realized this function already existed in the math library I use.
This will only align the model though. Then you would also need to scale and translate it so that the nose touches P2 and the tail P1.

(First post.)

14 minutes ago, UfU__X said:

Quaternion alignDirection(Vector3 V0, Vector3 V1){
    Vector3 A = V0.cross(V1);
    
    float length_0 = V0.norm();
    float length_1 = V1.norm();
    float dp = V0.dot(V1);
    
    float w = sqrt(length_0*length_0 * length_1*length_1) + dp;

    Quaternion Q(w, A.x, A.y, A.z);
    Q.normalize();
    return Q;
}

Take note on the quaternion initialization in your framework as they sometimes differ.
A is the "vector" of the quaternion values.

I don't get how it works. But I tested it before I realized this function already existed in the math library I use.
This will only align the model though. Then you would also need to scale and translate it so that the nose touches P2 and the tail P1.

(First post.)

This function calculates the rotation from one direction to another. You could look at how to calculate an axis and angle between two directions, how to convert axis and angle to quat and then you get the code you posted.

There are issues, however:

You calculate lengths (length_0 and length_1) just to square the result later. You should calculate only the squared length to avoid 2 unnecessary square roots.

If V0 and V1 are opposite directions, everything becomes zero. You either need to handle this in this function or in normalize(), but the latter is unlikely - infinites may creep in an cause nasty bugs

17 minutes ago, JoeJ said:

You calculate lengths (length_0 and length_1) just to square the result later. You should calculate only the squared length to avoid 2 unnecessary square roots.

Oh that's true. There is probably a squared norm function in the library (eigen) I should use instead.

18 minutes ago, JoeJ said:

If V0 and V1 are opposite directions, everything becomes zero. You either need to handle this in this function or in normalize(), but the latter is unlikely - infinites may creep in an cause nasty bugs

That's why I use the function in eigen instead. I'm just assuming it would be able to treat such cases better if not completely.
But I actually left a comment to myself in the code that normalize "should" treat such cases but I'm not sure if I tested if it did.

I forgot about another special case: Zero input vectors. Here's how i handle them both:

	void FromVecToVecNonUniformLength (const sVec3 &dir0, const sVec3 &dir1)
    {
        sScalar dot = dir0.Dot(dir1);
        (*this)[3] = sqrt(dir1.SqL() * dir0.SqL()) + dot; // set w
        if ((*this)[3] >= FP_EPSILON)
        { 
            *((sVec3*)this) = dir0.Cross(dir1); // set xyz
            Normalize();
        }
        else
        {
            if (dot<0) *this = sQuat (1.0f, 0, 0, 0); // input 180 degrees apart -> 180 degrees rotation around x
            else *this = sQuat (0, 0, 0, 1.0f); // zero vector input -> zero rotation output
        }
    }
	

Code is ugly because i cast quat to vector, which only works with quats in xyzw format.

Doing such checks in Normalize would be a waste of performance for most usecases, and Eigen maybe puts the responsibility to the user because it is necessary to pick an arbitary axis for a 180 degrees rotation.