3
$\begingroup$

I'm calculating the acceleration due to the harmonics in ECEF frame. Gravity potential in spherical harmonics is shown here (I just removed "$1+$", as considering only harmonics effect).

$U_{har}=\frac{\mu}{r}[\sum_{i=2}^d\sum_{j=0}^o (\frac{R_{eq}}{r})^iP_{ij}(\sin\phi)(S_{ij}\sin{j\lambda}+C_{ij}\cos{j\lambda})]$,

where $d$ and $o$- degree and order, $\phi$ and $\lambda$- latitude and longitude respectively.

I compare the results with GMAT. For degree 2 and order 0 (J2) the error in propagation was 5m. But for degree/order=8 the error is 350km!

The steps:

  1. In a loop $i\in[2,2]$, $j\in[0,0]$

    • Calculate the $P_i=\frac{d^{i+j}}{d\mu^{i+j}}(\mu^2-1)^i$
    • Calculate the Legendre polynom $P_{ij}=\frac{(1-\mu^2)^{\frac{j}{2}}}{i!*2^i}P_i$
    • Calculate the $sum+=P_{ij}(\frac{z}{r})*(C\cos(j*atan(\frac{y}{x}))+S\sin(j*atan(\frac{y}{x})))(\frac{R_{eq}}{r})^i$
  2. Calculate the potential $U_{har}=\mu\frac{sum}{r}$
  3. Calculate the $a_x$: $f=\frac{dU_{har}}{dx}$
  4. Calculate the value in the point $f(r_x,r_y,r_z)$

As it seen, the latitude is $asin(\frac{z}{r})$ and longitude is $atan(\frac{y}{x})$

The coefficients (JGM-3):

2 0 -0.10826360229840e-02 0.0
2 1 -0.24140000522221e-09 0.15430999737844e-08
2 2 0.15745360427672e-05 -0.90386807301869e-06

I have written a code on Julia language, which builds the expression (depending on the degree and order).

using SatelliteToolbox
using SymEngine

path="C:/xampp/htdocs/sat_prop/";
JGM_coeff_file=string(path,"coeff.txt");

const date  = DatetoJD(2100,01,01,0,0,0)
const degree = 8

y = [-4617E+03, 1709E+03, -5040E+03]

const Req = 6378136.3
const GMe = 398600.4415E+9

function harmonics(dy,y,dU,date)
  dy= [
        dU[1](y[1],y[2],y[3]),
        dU[2](y[1],y[2],y[3]),
        dU[3](y[1],y[2],y[3])
      ]
end

function potential()

  @vars x y z myu

  CS=open(readdlm,JGM_coeff_file)

  longitude=atan( y/x );
  r=(x^2+y^2+z^2)^(1/2)

  U_sum=0
  for i=2:degree
      for j=0:degree

          index=1+j; for ll=2:i-1 index+=ll+1; end

          P_i=(myu^2-1)^i
          for k=1:i+j P_i=diff(P_i,myu) end

          P_ij=(((1-myu^2)^(j/2))/(factorial(i)*2^i))*P_i

          if(P_ij!=0)
            U_sum+= P_ij(z/r)*(CS[index,3]*cos(j*longitude)+CS[index,4]*sin(j*longitude))*(Req/r)^i
          end
       end
  end

  U=GMe*(U_sum)/r

  return lambdify(expand(diff(U,x)),[x,y,z]),lambdify(expand(diff(U,y)),[x,y,z]),lambdify(expand(diff(U,z)),[x,y,z])
end


dU=potential();
dy=zero(y)
@time harmonics(dy,y,dU,date)
@time harmonics(dy,y,dU,date)
$\endgroup$
  • $\begingroup$ In what frame is the X, Y, and Z coordinates? They need to be in ECEF. $\endgroup$ – ChrisR Jul 29 '18 at 18:59
  • 1
    $\begingroup$ +1 Great editing, this looks much better, very nice! I'll try to take a careful look today, thanks for taking the time to dig in with MathJax! $\endgroup$ – uhoh Jul 29 '18 at 23:17
  • 1
    $\begingroup$ @Leeloo, okay, thanks for the edits. To be honest, I'm not sure how to answer this question sadly, as I've only used the Pines equations to compute the harmonics. $\endgroup$ – ChrisR Jul 30 '18 at 0:03
  • $\begingroup$ @Leeloo the equations are several pages long. A good reference is Chapter 2 of the PhD dissertation of Brandon A. Jones 2010. It's available freely on the University of Colorado website. The filename should be called "bajones2010" I think. He also mentions two other formulations which should speed up calculations (Fantino 2008), but they seemed significantly more complex to me. I've implemented Pines here, in Rust, as part of a toolkit I'm writing to automate most of the analysis I do, & with the fidelity of GMAT. $\endgroup$ – ChrisR Jul 30 '18 at 18:05
  • $\begingroup$ @Leeloo I'll take a look, that looks very promising! Possibly helpful here as well: How can I verify my reconstructed gravity field of Ceres from spherical harmonics? $\endgroup$ – uhoh Aug 1 '18 at 8:20
2
+50
$\begingroup$

Be careful with arc tangent

I am unfamiliar with Julia, but in most languages $\arctan(y/x)$ returns a value between $[-\pi/2, \pi/2]$. It does this because it doesn't know if $y$ or $x$ is positive or negative. So when you calculate longitude as:

$$ \lambda = \arctan(\frac{y}{x}) = \arctan(\frac{4617\times10^3}{-1709\times 10^3}) = -69.7^\circ$$

That is in the 4th quadrant, which is $180^\circ$ off from where the longitude actually is (the 2nd quadrant when $y$ is positive and $x$ is negative). This doesn't matter for the $j=0$ terms like J2, but it inverts the tesseral terms, since $\sin( \theta + 180 ) = - \sin \theta$ and $\cos( \theta + 180) = -\cos ( \theta )$. That would explain the discrepancy in agreement when adding more terms.

Use a arctan2

Most languages have a special function (usually called "atan2(x,y)" or some variant) that takes two parameters: an x and y parameter. If you don't have one, you will have to use some if statements to consider which quadrant the longitude is actually in.

$\endgroup$
  • $\begingroup$ You're absolutely right. However, there're some troubles with atan2 in Julia, so I used longitude=atan(y/x)+pi*sign(y)*(1-sign(x))/2. Is it OK? Assumed, that x and y are never 0. $\endgroup$ – Leeloo Aug 3 '18 at 8:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.