Untitled

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")