Skip to main content
Space Exploration

Return to Question

Notice removed Draw attention by CommunityBot
Bounty Ended with cms's answer chosen by CommunityBot
added 172 characters in body
Source Link
Leeloo
Leeloo
  • 743
  • 2
  • 12
  • 30

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

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
deleted 52 characters in body
Source Link
Leeloo
Leeloo
  • 743
  • 2
  • 12
  • 30
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

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

          legendre=P_ij=(((1-myu^2)^(j/2))/(factorial(i)*2^i))*legendre_part*P_i
          
          if(legendreP_ij!=0)
            U_sum+= legendreP_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)

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

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

          legendre=(((1-myu^2)^(j/2))/(factorial(i)*2^i))*legendre_part
          
          if(legendre!=0)
            U_sum+= legendre(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)

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)
added 102 characters in body
Source Link
Leeloo
Leeloo
  • 743
  • 2
  • 12
  • 30

To test the algorithm, I'm trying to calculate the acceleration (derivative of $U$), considering only J2, however the result is wrong (II compare with the equation, described here, p.9). The result of $a_x$ isresults with -3e-3GMAT, but should be about. For degree 2 and order 0 -1e-2(J2) the error in propagation was 5m. But for degree/order=8 the error is 350km!

  • 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(\sin\phi)$$P_{ij}=\frac{(1-\mu^2)^{\frac{j}{2}}}{i!*2^i}P_i$
  • Calculate the $sum+=P_{ij}*(C\cos(j\lambda)+S\sin(j\lambda))(\frac{R_{eq}}{r})^i$$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$

I suppose that latitude and longitude should be expressed through x,yAs it seen,z and differentiated then. Tried to replace the latitude withis $\text{asin}(\frac{z}{r})$,$asin(\frac{z}{r})$ and longitude withis $\text{atan}(\frac{y}{x})$, but it didn't solve the problem.$atan(\frac{y}{x})$

#Initialize theusing librariesSatelliteToolbox
using SymEngine
using Geodesy
path="C:/xampp/htdocs/sat_prop/";
#Initialize theJGM_coeff_file=string(path,"coeff.txt");

const variablesdate  myu-= sinDatetoJD(latitude2100,01,01,0,0,0)
const degree = 8

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

const longitudeReq = 6378136.3
@varsconst xGMe = 398600.4415E+9

function harmonics(dy,y,dU,date)
 z myudy= lambda
[
# Degree 2, order 0- only J2
d=2; o=0;
 dU[1](y[1],y[2],y[3]),
# C       dU[2](y[1],y[2],y[3]), 
 S       dU[3](y[1],y[2],y[3])
 gravitational parameter and equatorial radius ]
C=-0.10826360229840e-02end
S=0
GMefunction =potential()

 398600.4415E+9 @vars x y z myu
Req 
 = 6378136.3CS=open(readdlm,JGM_coeff_file)

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

  U_sum=0
  for i=2:ddegree
      for j=0:odegree

        #calculate the Legendreindex=1+j; polynomialfor ll=2:i-1 index+=ll+1; end

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

          legendre=(((1-myu^2)^(j/2))/(factorial(i)*2^i))*legendre_part(myu)
          
        #Calculate the sumif(legendre!=0)
        U_sum+=legendre*    U_sum+= legendre(C*cosz/r)*(j*lambdaCS[index,3]*cos(j*longitude)+S*sin+CS[index,4]*sin(j*lambdaj*longitude))*(Req/r)^i
          end
       end
  end

#Calculate the potential
U=GMe*(U_sum)/r

#Calculate the derivative
f=expandreturn lambdify(expand(diff(U,x))

#The values of x,[x,y,zz]), myu and lambda
vector = [-4617E+03lambdify(expand(diff(U, 1709E+03y)), -5040E+03]
latlon=LLA[x,y,z]),lambdify(ECEFexpand(vectordiff(U,z)), wgs84[x,y,z])
end
f

dU=potential(vector[1]);
dy=zero(y)
@time harmonics(dy,vector[2]y,vector[3]dU,sin(latlon.lat*pi/180date)
@time harmonics(dy,latlon.lon*pi/180y,dU,date)

To test the algorithm, I'm trying to calculate the acceleration (derivative of $U$), considering only J2, however the result is wrong (I compare with the equation, described here, p.9). The result of $a_x$ is -3e-3, but should be about -1e-2.

  • 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(\sin\phi)$
  • Calculate the $sum+=P_{ij}*(C\cos(j\lambda)+S\sin(j\lambda))(\frac{R_{eq}}{r})^i$

I suppose that latitude and longitude should be expressed through x,y,z and differentiated then. Tried to replace latitude with $\text{asin}(\frac{z}{r})$, longitude with $\text{atan}(\frac{y}{x})$, but it didn't solve the problem.

#Initialize the libraries
using SymEngine
using Geodesy

#Initialize the variables  myu- sin(latitude), lambda- longitude
@vars x y z myu lambda

# Degree 2, order 0- only J2
d=2; o=0;

# C, S, gravitational parameter and equatorial radius
C=-0.10826360229840e-02
S=0
GMe = 398600.4415E+9
Req = 6378136.3


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

U_sum=0
for i=2:d
    for j=0:o
        #calculate the Legendre polynomial
        legendre_part=(myu^2-1)^i
        for k=1:i+j
            legendre_part=diff(legendre_part,myu)
        end
        legendre=(((1-myu^2)^(j/2))/(factorial(i)*2^i))*legendre_part(myu)

        #Calculate the sum
        U_sum+=legendre*(C*cos(j*lambda)+S*sin(j*lambda))*(Req/r)^i
    end
end

#Calculate the potential
U=GMe*(U_sum)/r

#Calculate the derivative
f=expand(diff(U,x))

#The values of x,y,z, myu and lambda
vector = [-4617E+03, 1709E+03, -5040E+03]
latlon=LLA(ECEF(vector), wgs84)

f(vector[1],vector[2],vector[3],sin(latlon.lat*pi/180),latlon.lon*pi/180)

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!

  • 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$

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

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

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

          legendre=(((1-myu^2)^(j/2))/(factorial(i)*2^i))*legendre_part
          
          if(legendre!=0)
            U_sum+= legendre(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)
Math formatting
Source Link
edit approved Aug 1, 2018 at 9:39
Litho
edit approved Aug 1, 2018 at 9:39
Litho
  • 2.1k
  • 1
  • 14
  • 17
Tweeted twitter.com/StackSpaceExp/status/1024581057440612352
Notice added Draw attention by Leeloo
Bounty Started worth 50 reputation by Leeloo
[Edit removed during grace period]
Source Link
Leeloo
Leeloo
  • 743
  • 2
  • 12
  • 30
Post Undeleted by Leeloo
Post Deleted by Leeloo
added 1 character in body
Source Link
Leeloo
Leeloo
  • 743
  • 2
  • 12
  • 30
added 1 character in body
Source Link
Leeloo
Leeloo
  • 743
  • 2
  • 12
  • 30
added 940 characters in body
Source Link
Leeloo
Leeloo
  • 743
  • 2
  • 12
  • 30
added 940 characters in body
Source Link
Leeloo
Leeloo
  • 743
  • 2
  • 12
  • 30
added 940 characters in body
Source Link
Leeloo
Leeloo
  • 743
  • 2
  • 12
  • 30
edited tags
Link
uhoh
uhoh
  • 151k
  • 56
  • 505
  • 1.6k
added 206 characters in body; added 1 character in body
Source Link
Leeloo
Leeloo
  • 743
  • 2
  • 12
  • 30
added 55 characters in body
Source Link
Leeloo
Leeloo
  • 743
  • 2
  • 12
  • 30
Source Link
Leeloo
Leeloo
  • 743
  • 2
  • 12
  • 30