1 // Compute geodetic longitude, latitude and altitude

2 // from Cartesian orbit vector in ECR coordinates

3 

4 // Uses geodetic approximation from Patt and Gregg, IJRS, 1993.

5 

6 // Arguments

7 //

8 // Name Type I/O Description

9 // ---- ---- --- -----------

10 // nRecords int I number of records

11 // orb double I 3 x nRecords array of orbit vectors

12 // lon double O nRecords array of longitude

13 // lat double O nRecords array of latitude

14 // alt double O nRecords array of altitude

15 

16 // Liang Hong, Jan 14, 2016, Ported from IDL

17 // V0.1, Feb 1, 2016

18 

19 #include <math.h>

20 

21 int orb2lla(int nRecords, double orb[], double lon[], double lat[], double alt[]) {

22  // Define constants

23  double re = 6378.137; // Earth equatorial radius (km)

24  double rem = 6371; // Earth mean radius (km)

25  double f = 1.0 / 298.257; // Earth flattening factor

26 

27  int i;

28  double omf2, omf2p, rad, pxy, temp, rl;

29  omf2 = (1.0 - f)*(1.0 - f);

30 

31  for (i = 0; i < nRecords; i++) {

32  // Compute longitude

33  lon[i] = (180.0 / M_PI) * atan2(orb[1 + i * 6], orb[0 + i * 6]);

34 

35  // Compute geodetic latitude

36  rad = sqrt(orb[0 + i * 6] * orb[0 + i * 6] + orb[1 + i * 6] * orb[1 + i * 6] + orb[2 + i * 6] * orb[2 + i * 6]);

37  omf2p = (omf2 * rem + rad - rem) / rad;

38  pxy = orb[0 + i * 6] * orb[0 + i * 6] + orb[1 + i * 6] * orb[1 + i * 6];

39  temp = sqrt(orb[2 + i * 6] * orb[2 + i * 6] + omf2p * omf2p * pxy);

40  lat[i] = (180.0 / M_PI) * asin(orb[2 + i * 6] / temp);

41 

42  // Compute altitude

43  rl = re * (1.0 - f) / sqrt(1.0 - (2.0 - f) * f * pxy / (rad * rad));

44  alt[i] = rad - rl;

45  }

46  return 0;

47 }