#!/bin/env python from pylab import * from numpy import * from matplotlib.patches import Ellipse, Wedge rc('text', usetex=True) # For the bold face deg2rad = pi/180.0 rad2deg = 180.0/pi def Cos(x): return cos(pi*x/180.0) def Sin(x): return sin(pi*x/180.0) def findendpoint(m, b, Vlen, x1, y1): A = m*m+1 B = -2*x1 - 2*y1*m + 2*m*b C = -(Vlen*Vlen - x1*x1 - y1*y1 + 2*y1*b - b*b) x2 = (-B - sqrt(B*B-4*A*C))/(2*A) y2 = m*x2 + b x3 = (-B + sqrt(B*B-4*A*C))/(2*A) y3 = m*x3 + b return x2, y2, x3, y3 axis('equal') #axis([-6,6,-6,6]) axis('off') ax = gca() ex = 0.0; ey = 0.0 # Set properties ellipse and draw it major = 10.0; minor= 9.0 e1 = Ellipse((0,0), major, minor, 0, fill=False) ax.add_patch(e1) # Plot focus positions # Position of focus = c^2 = a^2 - b^2 a = 0.5*major; b = 0.5*minor c = sqrt(a*a-b*b) x = [c]; y= [0] plot(x,y, 'ro') text(c, -0.1*minor, '(c,0)', horizontalalignment='center') text(c, 0.02*minor, 'S', horizontalalignment='center') plot([-c],y, 'ro') text(-c, -0.1*minor, '(-c,0)', horizontalalignment='center') # Draw minor and major axis X1 = -a Y1 = 0 X2 = a Y2 = 0 plot([X1,X2], [Y1,Y2], 'g--') Y1 = -b X1 = 0 Y2 = b X2 = 0 plot([X1,X2], [Y1,Y2], 'g--') # Select position of object (e.g. earth) theta = 40.0 x1 = a * Cos(theta) y1 = b * Sin(theta) #aw = 0.04 #a1 = Arrow(c, 0.0, x1-c,y1, width=0.2) #ax.add_patch(a1) plot([x1],[y1], 'bo') text(x1, y1-0.03*minor, 'P', verticalalignment='top') angle = arctan2(y1, x1-c)*rad2deg radius = 0.07*major w1 = Wedge((c,0), radius, 0, angle, fill=False) ax.add_patch(w1) angle /= 2.0 xa = 1.2*radius * Cos(angle) ya = 1.2*radius * Sin(angle) text(c+xa, ya, r'$\alpha$') # Draw the line connecting the focal points with object on ellipse m1 = y1/(x1+c) m2 = y1/(x1-c) X1 = -c; Y1 = 0 lam = 1.0 X2 = -c + lam*(x1+c) Y2 = 0 + lam*y1 plot([X1,X2], [Y1,Y2], 'b--') X1 = c; Y1 = 0 lam = 1.0 X2 = c + lam*(x1-c) Y2 = 0 + lam*y1 plot([X1,X2], [Y1,Y2], 'b--') # Find the bisector that is tangent to the ellipse, not crossing a1 = y1/(x1+c); a2 = -1; k1 = y1*c/(x1+c) b1 = y1/(x1-c); b2 = -1; k2 = -y1*c/(x1-c) h=sqrt(b1*b1+b2*b2) k=sqrt(a1*a1+a2*a2) A1 = h*a1-k*b1; A2 = h*a2 - k*b2; A3 = h*k1 - k*k2 m = -A1/A2; off = -A3/A2 X1 = x1-3 Y1 = m*X1 + off X2 = x1+2 Y2 = m*X2 + off plot([X1,X2], [Y1,Y2], 'r--') # Draw velocity vector along tangent at x1,y1 Vlen = 3.0 b = off x2, y2, x3,y3 = findendpoint(m, b, Vlen, x1, y1) arr1 = Arrow( x1, y1, x2-x1,y2-y1, width=0.2, facecolor='red' ) ax.add_patch(arr1) text(x2, y2, r'\bf{V}', horizontalalignment='left', verticalalignment='bottom') xdotted1 = x2; ydotted1 = y2 # Draw vector for radial velocity i.e. perpendicular to radius vector # i.e. r.dtheta/dt phi1 = arctan2( y1,x1-c) phi2 = arctan2(y2-y1, x2-x1) phi3 = phi1 + pi/2 - phi2 m2 = tan(phi1+pi/2) b2 = y1 - m2*x1 PL = Vlen*cos(phi3) x2, y2, x3, y3 = findendpoint(m2, b2, PL, x1, y1) arr2 = Arrow( x1, y1, x2-x1,y2-y1, width=0.2, edgecolor='g' ) ax.add_patch(arr2) text(x2, y2-0.03*minor, 'L', verticalalignment='top') xdotted2 = x2; ydotted2 = y2 # Draw velocity component dr/dt as an extension to the radius vector m3 = tan(phi1) b3 = y1 - m3*x1 PR = Vlen*sin(phi3) x2, y2, x3, y3 = findendpoint(m3, b3, PR, x1, y1) arr3 = Arrow(x1, y1, x3-x1,y3-y1, width=0.2, edgecolor='g') ax.add_patch(arr3) text(x3, y3, 'R', horizontalalignment='left') xdotted3 = x3; ydotted3 = y3 plot([xdotted2, xdotted1,xdotted3], [ydotted2, ydotted1,ydotted3], ':') # Draw velocity component perpendicular to major axis PQ = PR / sin(phi1) arr4 = Arrow( x1, y1, 0, PQ, width=0.2, facecolor='y' ) ax.add_patch(arr4) text(x1, y1+PQ, r'$Q$', verticalalignment='bottom') # Draw velocity component perpendicular to radius vector PL2 = PL - PR/tan(phi1) x2, y2, x3, y3 = findendpoint(m2, b2, PL2, x1, y1) arr5 = Arrow(x1, y1, x2-x1,y2-y1, width=0.2, facecolor='y') ax.add_patch(arr5) text(x2, y2-0.03*minor, r'$L_1$', verticalalignment='top') # The other bisector """ A1 = h*a1+k*b1; A2 = h*a2 + k*b2; A3 = h*k1 + k*k2 m = -A1/A2; off = -A3/A2 X1 = x1 - 3 Y1 = m*X1 + off X2 = x1 + 3 Y2 = m*X2 + off plot( [X1,X2], [Y1,Y2]) """ savefig("etermfig.png") show()