functions - FunctionInterpolation Errors / Question re Evaluation Order and Options

I have using Mathematica functions that takes a Cartesian coordinate relative to the Earth (xyz) and converts it to a latitude, longitude, and altitude (lla). And here it is:

xyz2lla = First@GeoPosition@GeoPositionXYZ[#, "WGS84"] &

I'm using it to convert a satellite viewing line (line of sight or los) from xyz to lla. For instance, if I have a satellite at the following observation point (obs, meters) looking in the direction of look:

obs = {2.560453600382259, 5.245110323032143, -3.819772142191310} 1*^6;
look = {-0.233218833096895, -0.814561997858160, -0.531128729720249};

It has the line of sight:

obs+d look

The transformation xyz2lla is smooth, so I was hoping to use function interpolation:

f = FunctionInterpolation[{xyz2lla[obs + # look]} &[d], 
                          {d, 2000000, 3700000}
    ]

And while this works, and I get the following lat, lon, alt functions from it:

:

I also get the following errors that I'm wondering why I get:

GeoPositionXYZ::invcoord: "\!\(\"{2.560453600382259*^6 - 0.233218833096895*d, 5.245110323032143*^6 - 0.81456199785816*d, -3.81977214219131*^6 - 0.531128729720249*d}\"\) is not a valid coordinate specification."

Thread::tdlen: Objects of unequal length in {-212500.,-70833.3,70833.3,212500.}^{} cannot be combined. >>

Thread::tdlen: "Objects of unequal length in {12.2667 +3.14159\ I,11.1681 +3.14159\ I,11.1681,12.2667}\ {}\\n cannot be combined."

FunctionInterpolation::nreal: Near d = 2.2125`*^6, the function did not evaluate to a real number.

Does anyone have any insight into these errors?

Also, how can I adjust the quality of the interpolation / how many points are investigated?

Answer

I don't think InterpolatingFunction is intended to work on vector functions. The doc page doesn't say anything about it. Try for instance

f = FunctionInterpolation[{x, x}, {x, 0, 6}]

And

Table[f[x], {x, 0, 6, 1}]

returns

(*
==> {0., 0.9999652778, 1.998576389, 3., 4.004131944, \
4.994409722, 6.}
*)

So, no two dimensional output.

It's probably better to come up with three separate interpolating functions for each of the three components of the output.

---

Having said that, it doesn't look like this is the end of the problems.

xyz2llaPhi = GeoPosition[GeoPositionXYZ[obs + # look, "WGS84"]][[1, 1]] &;
fPhi = FunctionInterpolation[xyz2llaPhi[d], {d, 2000000, 3700000}]

yields

Whereas

xyz2llaPhi[200000]

yields a nice numerical result:

(*

==> -34.86629487
*)

It looks like the combination of FunctionInterpolation and GeoPosition isn't a healthy one. Somehow, the index d is held unevaluated.

A workaround would be to generate a table of values and then use ListInterpolation.