How to convert VSOP's xyz coordinates to heliocentric spherical lbr coordinates? I ran VSOP2013 in Fortran, and I now have the following position $x, y, z$ as well as velocity $v_x, v_y, v_z$ data.

How do I convert the position to heliocentric ecliptic spherical coordinates $l, b, r$? 

PLANETARY EPHEMERIS VSOP2013  EMB     
Ecliptic   Heliocentric Coordinates:  X,Y,Z (au)  X',Y',Z' (au/d)  - Dynamical Frame J2000

Julian Day 2415021.0   
-0.2054467990    0.9615310525    0.0002141233   -0.0171063024   -0.0036576954   -0.0000011815
Looking at the numbers and comparing to JPL Horizons, the best match is to the position of the Earth-Moon Barycenter in heliocentric coordinates. That's of course consistent with the data being labeled as:  EMB, Ecliptic Heliocentric Coordinates.

If you can live with some small uncertainties due to the details of how coordinate transforms may evolve over time, I believe you can just apply some simple trigonometry. You can confirm here.

$$r = \sqrt{x^2 + y^2 + z^2}$$

$$l = tan^{-1}(y/x)$$

$$b = sin^{-1}(z/r)$$





Data in Question: for JD 2415021.0
position (AU):     [-0.2054467990,  0.9615310525,  0.0002141233]
velocity (AU/day): [-0.0171063024, -0.0036576954, -0.0000011815]

JPL Horizions Earth-Moon barycenter J2000 Heliocentric for JD 2415021.0
position (AU):     [-0.2054465857,  0.9615310983,  0.00021393097] 
velocity (AU/day): [-0.0171063032, -0.0036576916, -0.00000118253]

JPL Horizions Earth J2000 Heliocentric for JD 2415021.0
position (AU):     [-0.2054522688,  0.9615603123,  0.00021303115] 
velocity (AU/day): [-0.0171136285, -0.0036594066, -0.00000182603]

JPL Horizions Earth J2000 Solar System barycenter for JD 2415021.0
position (AU):     [-0.2022720210,  0.9679285296,  0.00010939456]
velocity (AU/day): [-0.0171209826, -0.0036556277, -0.00000165111]