# Untitled

unknown
plain_text
a year ago
2.1 kB
3
Indexable
Never
```import numpy as np
import matplotlib.pyplot as plt

def cart2pol(x, y):
rho = np.sqrt(x**2 + y**2)
phi = np.arctan2(y, x)
return(rho, phi)

def pol2cart(rho, phi):
x = rho * np.cos(phi)
y = rho * np.sin(phi)
return(x, y)

def derivative(i, x, y):

DXP=x[i+1]-x[i]
DXM=x[i]-x[i-1]
DYP=y[i+1]-y[i]+1.E-30
DYM=y[i]-y[i-1]+1.E-30
ELP=1.
ELM=1.
ELP=np.sqrt(DXP**2+DYP**2)
ELM=np.sqrt(DXM**2+DYM**2)
DXX=(DXP*ELM+DXM*ELP)/(ELP+ELM)
DYY=(DYP*ELM+DYM*ELP)/(ELP+ELM)
DXDY=DXX/(DYY+1.E-30)

return DXDY

x = []; y = []
with open("P1.dat", 'r') as f:
for line in f:
split = line.split("  ")
x.append(float(split[1]))
y.append(float(split[2]))

x2 = []; y2 = []
with open("P2.dat", 'r') as f:
for line in f:
split = line.split("  ")
x2.append(float(split[1]))
y2.append(float(split[2]))

plt.plot(x,y)
plt.plot(x2,y2)
plt.show()

x = np.array(x)
y = np.array(y)
r1, t1 = cart2pol(x, y)
angle1 = 2*np.pi/4

t1 = [t1 + angle1*i for i in range(4)]
x1 = np.array([]); y1 = np.array([])

plt.figure()
for i in range(4):
xm, ym = pol2cart(r1, t1[i])
x1 = np.concatenate((x1,xm)); y1 = np.concatenate((y1,ym))

r1w = np.max(r1) - 0.01
theta = np.linspace(0, 2*np.pi, 100)
xc1 = r1w * np.cos(theta)
yc1 = r1w * np.sin(theta)

r2w = r1w * 5/4
xc2 = r2w * np.cos(theta)
yc2 = r2w * np.sin(theta) + (r1w+r2w)
E = r1w + r2w

plt.plot(xc1, yc1)
plt.plot(xc2, yc2)
plt.plot(x1, y1)
plt.show()

phi1 = np.linspace(0, 4*np.pi, 100)

plt.figure()
for PHI_1 in phi1:
PHI_2 = 4/5 * PHI_1

M21 = [[np.cos(PHI_2-PHI_1), np.sin(PHI_2-PHI_1), 0, -E*np.sin(PHI_2)],
[-np.sin(PHI_2-PHI_1), np.cos(PHI_2-PHI_1), 0, -E*np.cos(PHI_2)],
[0, 0, 1, 0],
[0, 0, 0, 1]]
M21 = np.matrix(M21)

tmp = np.array(M21 * np.matrix([x1, y1, np.zeros(np.size(x1)), np.ones(np.size(x1))]))

plt.plot(tmp[0], tmp[1])

plt.show()
print("S")

```