Plot for body launched at alpha to horizon with circular velocity

kepler2alpha-1.py

import os
import numpy as np
from matplotlib.patches import Ellipse
import matplotlib as mpl
from matplotlib import pyplot as plt
from matplotlib import rc, rcParams
from matplotlib.ticker import AutoMinorLocator, MultipleLocator
from matplotlib.patches import Arc

# Trigonometry {{{
def cos(x):
    return np.cos( np.radians(x) )
def sin(x):
    return np.sin( np.radians(x) )
def tg(x):
    return np.tan( np.radians(x) )
#}}}

# Const. {{{
lw=1.3
blue = '#4477EE'
red = '#EE4466'
gray = '#333333'
lightgray = '#888888'
#}}}

# Plot setup {{{
from matplotlib import pyplot as plt
from matplotlib import rc, rcParams
rcParams['font.size'] = 16.
rcParams['text.usetex'] = True
#rc('text.latex', preamble=r'\usepackage{esvect}')
rcParams['font.sans-serif'] = ['Computer Modern']
rc('text.latex', preamble=r'\usepackage[T2A]{fontenc} \usepackage[utf8]{inputenc} \usepackage{bm}')
#rc('text.latex', preamble=r'\usepackage[utf8]{inputenc}')

rc('axes', linewidth=1.2)
rc('xtick.major', size=4, width=1.1)
rc('ytick.major', size=4, width=1.1)
fig, ax = plt.subplots()
fig.subplots_adjust(left=.15, bottom=.15, right=.95, top=.95)
#}}}

# Set ellipse parameters {{{
w = 10
a = w/2

h = 6
b = h/2

c = (a**2-b**2)**0.5
e = c/a
angle = 90
x0 = 0
y0 = 0
#}}}

# Plot ellipse {{{
ell = mpl.patches.Ellipse(
        xy=(x0,y0),
        width=w, height=h,
        fill=False,
        angle = angle,
        lw=lw,
        color='k',
        zorder=11
        )

ax.add_patch(ell)
#}}}

# Plot semimajor and semiminor axes {{{
# SemiMajor axis
xl = x0 - a*cos(angle)
xr = x0 + a*cos(angle)
x = [xl, xr]
dy = a * (2 * (1-cos(angle)))**0.5 * cos(angle/2)
yl = y0 - dy
yr = y0 + dy
y = [yl, yr]
line,=ax.plot(x,y, color=gray, lw=lw)

# Minor
dx = b * (2*(1-cos(angle)))**0.5 * cos(angle/2)
xl = x0 + dx
xr = x0 - dx
x = [xl,xr]
y = [y0-b*cos(angle), y0 + b*cos(angle)]
line,=ax.plot(x,y, color=gray, lw=lw)
#}}}

# Plot Focuses Foci {{{
dxf = c*cos(angle)
dyf = dy*c/a

#F₁
ax.plot(x0 - dxf, y0 - dyf, 'o', ms=5, color=gray, mec=gray)
#F₂
#ax.plot(x0 + dxf, y0 + dyf, 'o', ms=3, color=gray, mec=gray)
#}}}

# Plot circle {{{
r = 5
x0c = x0
y0c = y0-dyf
circ = plt.Circle((x0c, y0c), r, color='k', 
        fill=True, 
        fc='#CCCCCC',
        ec='#AAAAAA',
        lw=1.3)
ax.add_artist(circ)
#}}}

# Plot tangent to circle and text K {{{
# Tangent point
xt = x0c - 3 # -(r+x0c) < xt < (r+x0c)

D = 4*y0c**2 - 4*(y0c**2 + (xt-x0c)**2 - r**2)
if y0c>0: yt = ( 2*y0c - D**0.5 ) / 2
else:    yt = ( 2*y0c + D**0.5 ) / 2

# dy/dx = tga
if xt==x0c-r or xt==x0c+r:
    tga=1e99
else:
    tga = -( (xt-x0c) / (r**2-(xt-x0c)**2)**0.5 )
    #if yt>0: tga = +( (xt-x0c) / (r**2-(xt-x0c)**2)**0.5 )
    #else:    tga = -( (xt-x0c) / (r**2-(xt-x0c)**2)**0.5 )

h=2.5
x = np.linspace( xt-1.6*h, xt+h, 2)
y = tga*(x - xt) + yt

# Plot tangent line
ax.plot(x,y, color=gray, lw=1.1)
ax.text(xt-h*1.25, -4.75, '$\\textrm{Горизонт}$' + '\n' + u'$\\textrm{в}$ $\\textrm{точке}$ $B$', ha='right')

## Add K letter
#xk = -1
#yk = tga*(xk - xt) + yt
#ax.plot(xk, yk, 'o', ms=4, color=gray, mec=gray)
#ax.text(xk-0.45, yk+0.55, '$K$', color='k')
#}}}

# Annotate ellipse ABCD O F₁ F₂ {{{
#ax.plot(xf2, yf2, 'o', ms=5, color=gray, mec=gray)

ax.plot(x0,      y0, 'o', ms=4, color=gray, mec=gray)
ax.text(x0+0.15, y0-0.77, '$O$', color='k', ha='left')

ax.plot(x0,      y0+a-c, 'o', ms=4, color=gray, mec=gray)
ax.text(x0+0.15, y0+a-c+0.2, '$L$', color='k', ha='left')

ax.plot(xt, yt, 'o', ms=4, color=gray, mec=gray, zorder=102)
ax.text(xt-1.05, yt-0.05, '$B$', color='k', va='bottom', ha='left')

ax.plot(x0,      y0+a, 'o', ms=4, color=gray, mec=gray)
ax.text(x0+0.15, y0+a+0.25, '$C$', color='k', ha='left')

ax.plot(xt+2*b, yt, 'o', ms=4, color=gray, mec=gray)
ax.text(xt+2*b+0.20, yt-0.05, '$D$', color='k', va='bottom', ha='left')

ax.plot(x0,      y0-a, 'o', ms=4, color=gray, mec=gray)
ax.text(x0+0.15, y0-a-0.95, '$A$', color='k', ha='left')

ax.text(x0+0.20, y0-dyf+0.20, '$F_1$', color='k', ha='left')
#}}}

# Plot V, annotate V, α {{{
hl=0.3
hw=0.2
z=2.5
theta1=np.degrees(np.arctan(tga))
print(90-theta1)
ax.add_patch(Arc((xt, yt), z, z, theta2=90, theta1=np.degrees(np.arctan(tga)), lw=1.1))
ax.arrow( xt, yt, 0, 4, head_length=hl, head_width=hw, fc=red, ec=red, lw=1.4, zorder=99)
ax.text(xt-0.85, yt+4.1, r'$\bm{v}$', color=red, fontsize=20)
ax.text(xt+0.6, yt+z/2, '$\\alpha$', color='k')
#}}}

## Plot parabola BCD for comparison {{{
## B coordinates
#c1=5
#xb=xt
#yb=yt+c1
#print(xb,yb)
#a1 = - yb/xb**2
#b1=0
#x_vertex = -b1/(2*a1)
#y_vertex = (4*a1*c1 - b1**2) / (4*a1)
#xp = np.linspace(-3,3,30)
#yp = a1*xp**2 + b1*xp + c1
#line,=ax.plot(xp,yp, lw=1.2, color=lightgray)
#line.set_dashes([1.3, 2.3])
##}}}

# Plot r1 {{{
print(a,b,c)
dx = b * (2*(1-cos(angle)))**0.5 * cos(angle/2)
xl = x0 - c * cos(angle)
xr = x0 - dx
x = [xl,xr]

dy = c * (2 * (1-cos(angle)))**0.5 * cos(angle/2)
yl = y0 - dy
yr = y0 + b*cos(angle)
y = [yl,yr]
ax.plot(x,y, color=gray, lw=1.2)
#}}}

## Plot r2{{{
#dx = b * (2*(1-cos(angle)))**0.5 * cos(angle/2)
#xl = x0 + c * cos(angle)
#xr = x0 - dx
#x = [xl,xr]
#
#dy = c * (2 * (1-cos(angle)))**0.5 * cos(angle/2)
#yl = y0 + dy
#yr = y0 + b*cos(angle)
#y = [yl,yr]
#ax.plot(x,y, color=gray, lw=1.2)
##}}}

# Set limits etc {{{
xmax = x0 + w/2
ymax = y0 + h/2
kx = 2.5
ky = kx
#ax.set_xlim([-xmax*kx,xmax*kx])
#ax.set_ylim([-ymax*ky,ymax*ky])
#ax.set_axisbelow(False)
plt.axis('equal')
ax.set_xlim([-12, 9])
ax.set_ylim([-12, 9])

plt.grid()
#ax.xaxis.set_major_locator( MultipleLocator(3))
#ax.yaxis.set_major_locator( MultipleLocator(2))
#ax.tick_params(top='off', right='off', which='both')

#plt.axis('off')

plt.tick_params(
    axis='both',
    which='both',
    bottom=False,
    top=False,
    labelbottom=False)
plt.tick_params(
    axis='y',
    which='both',
    left=False,
    right=False,
    labelleft=False)


#}}}

## Save the figure {{{
#plt.savefig('k1.png')
##}}}

plt.show()

kepler2alpha-2.py

import math, os
import numpy as np
from matplotlib.patches import Ellipse
import matplotlib as mpl
from matplotlib import pyplot as plt
from matplotlib import rc, rcParams
from matplotlib.ticker import AutoMinorLocator, MultipleLocator


def rotatePoint(centerPoint,point,angle): #{{{
    angle = math.radians(angle)
    temp_point = point[0]-centerPoint[0] , point[1]-centerPoint[1]
    temp_point = ( temp_point[0]*math.cos(angle)-temp_point[1]*math.sin(angle) , temp_point[0]*math.sin(angle)+temp_point[1]*math.cos(angle))
    temp_point = temp_point[0]+centerPoint[0] , temp_point[1]+centerPoint[1]
    return temp_point
#}}}

# Trigonometry {{{
def cos(x):
    return math.cos( math.radians(x) )
def sin(x):
    return math.sin( math.radians(x) )
#}}}

# Const. {{{
res=50
lw=1.3
blue = '#4477EE'
lightblue = '#78B7D8'
red = '#EE4466'
gray = '#333333'
#}}}

# Plot setup {{{
from matplotlib import pyplot as plt
from matplotlib import rc, rcParams
rcParams['font.size'] = 16.
rcParams['text.usetex'] = True
#rc('text.latex', preamble=r'\usepackage{esvect}')
rcParams['font.sans-serif'] = ['Computer Modern']
rc('text.latex', preamble=r'\usepackage[T2A]{fontenc} \usepackage[utf8]{inputenc} \usepackage{bm}')
#rc('text.latex', preamble=r'\usepackage[utf8]{inputenc}')

rc('axes', linewidth=1.2)
rc('xtick.major', size=4, width=1.1)
rc('ytick.major', size=4, width=1.1)
fig, ax = plt.subplots()
fig.subplots_adjust(left=.15, bottom=.15, right=.95, top=.95)
#}}}

# Set ellipse parameters {{{
w = 10
a = w/2

h = 6
b = h/2

c = (a**2-b**2)**0.5
e = c/a
angle = 90
x0 = 0
y0 = 0
#}}}

# Plot ellipse {{{
ell = mpl.patches.Ellipse(
        xy=(x0,y0), 
        width=w, height=h, 
        fill=False,
        angle = angle,
        lw=lw,
        color='k',
        zorder=11
        )

ax.add_patch(ell)
#}}}

# Plot semimajor and semiminor axes {{{
# SemiMajor axis
xl = x0 - a*cos(angle)
xr = x0 + a*cos(angle)
x = [xl, xr]
dy = a * (2 * (1-cos(angle)))**0.5 * cos(angle/2)
yl = y0 - dy
yr = y0 + dy
y = [yl, yr]
line,=ax.plot(x,y, color=gray, lw=lw)

# Minor
dx = b * (2*(1-cos(angle)))**0.5 * cos(angle/2)
xl = x0 + dx
xr = x0 - dx
x = [xl,xr]
y = [y0-b*cos(angle), y0 + b*cos(angle)]
line,=ax.plot(x,y, color=gray, lw=lw)
#}}}

## Plot Focuses (Foci) {{{
dxf = c*cos(angle)
dyf = dy*c/a
xf1 = x0+dxf
yf1 = y0+dyf
xf2 = x0-dxf
yf2 = y0-dyf
#ax.plot(xf1, yf1, 'o', ms=5, color=gray, mec=gray)
#ax.plot(xf2, yf2, 'o', ms=5, color=gray, mec=gray)
#}}}

# Plot r1 {{{
dx = b * (2*(1-cos(angle)))**0.5 * cos(angle/2)
xl = x0 - c * cos(angle)
xr = x0 - dx
x = [xl,xr]

dy = c * (2 * (1-cos(angle)))**0.5 * cos(angle/2)
yl = y0 - dy
yr = y0 + b*cos(angle)
y = [yl,yr]
ax.plot(x,y, color=gray, lw=1.2)

xper=yper=[]
xper = np.append(xper, x)
yper = np.append(yper, y)
#}}}

# Plot r2{{{
dx = b * (2*(1-cos(angle)))**0.5 * cos(angle/2)
xl = x0 - c * cos(angle)
xr = x0 + dx
x = [xl,xr]

dy = c * (2 * (1-cos(angle)))**0.5 * cos(angle/2)
yl = y0 - dy
yr = y0 - b*cos(angle)
y = [yl,yr]
ax.plot(x,y, color=gray, lw=1.2)
xper = np.append(xper, x)
yper = np.append(yper, y)
#}}}

# Get perimeter of area to be filled {{{
xs = np.linspace(0,5,res)
y1 =  b * (1 - xs**2/a**2)**0.5
y2 = -b * (1 - xs**2/a**2)**0.5

xn=yn=xn1=yn1=xn2=yn2=[]
for x,y in zip(xs,y1):
    (xr, yr) = rotatePoint( (x0,y0), (x,y), angle)
    xn = np.append (xn, xr)
    yn = np.append (yn, yr)
    xn1 = np.append (xn1, xr)

xper = np.append(xn[::-1], xper)
yper = np.append(yn[::-1], yper)

for x,y in zip(xs,y2):
    (xr, yr) = rotatePoint( (x0,y0), (x,y), angle)
    xn = np.append (xn, xr)
    yn = np.append (yn, yr)
    xn2 = np.append (xn2, xr)
    yn2 = np.append (yn2, yr)

xper = np.append(xn2, xper)
yper = np.append(yn2, yper)
#}}}
plt.fill(xper, yper, facecolor=lightblue, edgecolor='None')

# Annotate ABCD  {{{
dx = 0.82
dy = 0.16
ad = {
'A':   [x0-a, y0,   'left',  'bottom' ],
'B':   [x0,   y0+b, 'right', 'top' ],
'C':   [x0+a, y0,   'right',  'top' ],
'D':   [x0,   y0-b, 'left',  'top' ],
'O':   [x0,   y0,   'right', 'top' ],
'F_1': [x0-c, y0,   'right', 'bottom' ],
#'F_2': [x0+c, y0,   'right', 'top' ],
}

for p in ad:
    (x,y) = rotatePoint((x0,y0),(ad[p][0], ad[p][1]),angle)
    ax.plot(x,y, 'o', ms=4, color=gray)
    if ad[p][2] == 'left':  x-=dx
    if ad[p][2] == 'right': x+=dx
    if ad[p][3] == 'top':   y+=dy
    if ad[p][3] == 'bottom':y-=3.9*dy

    ax.text(x ,y, '${}$'.format(p), ha=ad[p][2])#, va=ad[p][3])

#}}}

# Set limits etc {{{
xmax = x0 + w/2
ymax = y0 + h/2
kx = 2.5
ky = kx
#ax.set_xlim([-xmax*kx,xmax*kx])
#ax.set_ylim([-ymax*ky,ymax*ky])
k = 3.3
ax.set_xlim([-k,    k])
ax.set_ylim([-2*k,2*k])

#ax.autoscale()
#plt.axis('off')
#plt.axis('equal')
ax.set_aspect('equal')
ax.grid()

plt.tick_params(
    axis='both',
    which='both',
    bottom=False,
    top=False,
    labelbottom=False)
plt.tick_params(
    axis='y',
    which='both',
    left=False,
    right=False,
    labelleft=False)

ax.xaxis.set_major_locator( MultipleLocator(2))
ax.yaxis.set_major_locator( MultipleLocator(2))

plt.tick_params(
    axis='both',
    which='both',
    bottom=False,
    top=False,
    labelbottom=False)
plt.tick_params(
    axis='y',
    which='both',
    left=False,
    right=False,
    labelleft=False)
#}}}

## Save the figure {{{
##plt.savefig('k2.png')
##}}}


plt.show()