Vec3D new_dir = new Vec3D(velocity);
    new_dir.normalize();
    Vec3D new_up = new Vec3D(0.0, 1.0, 0.0);
    new_up.normalize();
    Vec3D crossP = new_dir.cross(new_up);
    crossP.normalize();
    
    float dotP  = new_dir.dot(new_up);
    float angle = new_dir.angleBetween(new_up, true); 
    
    
    pushMatrix();
    // update location
    translate(location.x, location.y, location.z);
    // orientation to velocity
    rotate(-angle, crossP.x, crossP.y, crossP.z);
    
    // draw your stuff here
    
    popMatrix(); 

    // Rotation vectors
    // use to perform orientation to velocity vector
    Vec3D new_dir = new Vec3D(velocity);
    float r = sqrt(new_dir.x * new_dir.x + new_dir.y * new_dir.y + new_dir.z * new_dir.z);
    float theta = atan2(new_dir.y, new_dir.x);
    float phi = acos(new_dir.z / r);
    
    pushMatrix();
    // update location
    translate(location.x, location.y, location.z);
    // orientation to velocity
    rotateZ(theta);
    rotateY(phi);
    rotateX(HALF_PI);
    
    // draw your stuff here
    
    popMatrix(); 

that doesn't account for is "roll" (to use the airplane analogy).  with only two points (well, a position and a velocity vector, if I've correctly followed your several threads) you can't nail down all three implied axes, all you can get is two. (which is enough to "point at", but not enough to determine where "up/roll" should be)

typically if you need all three axes you'll have to carry around your own pitch/roll/yaw or a frenet frame or orthonormal TNB vectors or etcetera, potentially converting back and forth to quats if you need to slerp them, and WON'T be attempting to "derive" them from velocity-position alone.  

you can try faking that third local axes ("up") by hard-coding it (perhaps as the world y-axis) sometimes.  then you can cross your "at" vector (vel-pos) with "up" to get "right" and you now have local TNB.  keep in mind that crossing only guarantees perpendicular, not WHICH perpendicular, so you typically cross back, maybe iterate, etc, yuck -- it's still only a partial solution.  so you can expect to experience problems as your velocity approaches being colinear with whatever hard-coded value you choose.