Untitled
plain_text
a month ago
2.1 kB
2
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")