A hint... Using trigonometry to attach can work well, but make sure you use a Vincenty algorithm corrected for altitude to calculate the Lat,Lon, and don't assume a spherical earth. The code I've been using which works well is this, ported from javascript I found somewhere online:
void destVincenty(double lat1, double lon1, double bearing, double dist, double altitude,
double *lat2out, double *lon2out)
{
/* local variable definitions */
// WGS-84 ellipsiod
double a = RAD_EARTH_MAJOR + altitude, b = RAD_EARTH_MINOR + altitude, f = WGS_84_FLATTEN;
// double a = 6378137.0 + altitude, b = 6356752.3142 + altitude, f = 1 / 298.257223563;
double alpha1, sinAlpha, sinAlpha1, cosAlpha1, cosSqAlpha;
double sigma, sigma1, cos2SigmaM, sinSigma, cosSigma, deltaSigma, sigmaP;
double tanU1, cosU1, sinU1, uSq;
double A, B, C, L, lambda;
double tmp, lat2;
//double revAz; /* unused but retained for alg completeness */
/* code body */
alpha1 = bearing;
sinAlpha1 = sin(alpha1);
cosAlpha1 = cos(alpha1);
tanU1 = (1.0 - f) * tan(lat1);
cosU1 = 1.0 / sqrt((1.0 + tanU1*tanU1));
sinU1 = tanU1*cosU1;
sigma1 = atan2(tanU1, cosAlpha1);
sinAlpha = cosU1 * sinAlpha1;
cosSqAlpha = 1.0 -sinAlpha*sinAlpha;
uSq = cosSqAlpha * (a*a - b*b) / (b*b);
A = 1.0 + uSq / 16384.0 * (4096.0 + uSq*(-768.0 + uSq*(320.0 - 175.0 * uSq)));
B = uSq / 1024.0 * (256.0 + uSq*(-128.0 + uSq*(74.0 - 47.0 * uSq)));
sigma = dist / (b*A);
sigmaP = 2.0 * M_PI;
while (fabs(sigma - sigmaP) > 1e-12) {
cos2SigmaM = cos(2.0 * sigma1 + sigma);
sinSigma = sin(sigma);
cosSigma = cos(sigma);
deltaSigma = B*sinSigma*(cos2SigmaM + B / 4.0 * (cosSigma*(-1.0 + 2.0 * cos2SigmaM*cos2SigmaM) - B / 6.0 * cos2SigmaM*(-3.0 + 4.0 * sinSigma*sinSigma)*(-3.0 + 4.0 * cos2SigmaM*cos2SigmaM)));
sigmaP = sigma;
sigma = dist / (b*A) + deltaSigma;
}
tmp = sinU1*sinSigma - cosU1*cosSigma*cosAlpha1;
lat2 = atan2(sinU1*cosSigma + cosU1*sinSigma*cosAlpha1,
(1 - f)*sqrt(sinAlpha*sinAlpha + tmp*tmp));
lambda = atan2(sinSigma*sinAlpha1,
cosU1*cosSigma - sinU1*sinSigma*cosAlpha1);
C = f / 16.0 * cosSqAlpha*(4.0 + f*(4.0 - 3.0 * cosSqAlpha));
L = lambda - (1 - C)*f*sinAlpha*(sigma + C*sinSigma*(cos2SigmaM + C*cosSigma*(-1.0 + 2.0 * cos2SigmaM*cos2SigmaM)));
// final bearing
// revAz = atan2(sinAlpha, -tmp);
*lat2out = lat2;
*lon2out = lon1 + (L);
}
Hope it's helpful,
Farley