Plotting satellite orbits with Gnuplot

orbits.gnuplot

r(p,e,f) = p/(1.0 + e * cos(f))
 
orbit_x(p,e,i,an,argp,f) = \
    r(p,e,f) * (\
    cos(f) * ((cos(argp) * cos(an)) - (sin(argp) * sin(an) * cos(i))) + \
    sin(f) * (-(cos(argp) * sin(an) * cos(i)) - (sin(argp) * cos(an))))
orbit_y(p,e,i,an,argp,f) = \
    r(p,e,f) * (\
    cos(f) * ((sin(argp) * cos(an) * cos(i)) + (cos(argp) * sin(an))) + \
    sin(f) * ((cos(argp) * cos(an) * cos(i)) - (sin(argp) * sin(an))))
orbit_z(p,e,i,an,argp,f) = \
    r(p,e,f) * (\
    cos(f) * (sin(argp) * sin(i)) + \
    sin(f) * (cos(argp) * sin(i)))
 
set view equal xyz
set parametric
 
set key off
set grid
 
p1 = 1.0
e1 = 0.0
i1 = 0.0
an1 = 0.0
argp1 = 0.0
 
p2 = 1.0
e2 = 0.7
i2 = pi/4
an2 = 0.0
argp2 = 0.0
 
splot \
    orbit_x(p1,e1,i1,an1,argp1,u), \
    orbit_y(p1,e1,i1,an1,argp1,u), \
    orbit_z(p1,e1,i1,an1,argp1,u) \
    title "Orbit 1", \
    orbit_x(p2,e2,i2,an2,argp2,u), \
    orbit_y(p2,e2,i2,an2,argp2,u), \
    orbit_z(p2,e2,i2,an2,argp2,u) \
    title "Orbit 2"