Python

쌍극좌표계 작도

coding art 2025. 3. 18. 16:23
728x90

쌍극 좌표계는 서로 부호가 반대이며 일정 거리만큼 떨어져 있는 크가가 같은 전하 배치에 의해서 생성되는 전기장 계산 즉 dipole 문제로 잘 알려져 있다.

 

 

 

아래는 쌍극좌표계 코드 bipolarplot.py 의 아나콘다에서의 실행 결과이다.

 

 

 

 

# bipolarplot.py

 

import numpy as np
import matplotlib.pyplot as plt

def bpplot(a, n):

# 1. 상면에 2개의 원 작도
    alpha1 = 2.0 ; alpha2 = 0.8
    m = 2*n+1
    x = [] ; y = [] ; my =[] ; xs = [] ; ys = [] ; mys = []
    angle = []
    
    for i in range(m):
    # 파이 N 등분 후 -pi 에서 0까지 완쪽 원 0에서 +pi까지 우측원 연산
        theta = i*np.pi/n - np.pi
    # (xx, yy)는 alpha1 즉 작은 원의 좌표, (xxx, yyy)는 큰 원의 좌표
        xx = a*np.sin(theta)/(np.cosh(alpha1)-np.cos(theta))
        xxx = a*np.sin(theta)/(np.cosh(alpha2)-np.cos(theta))
        yy = a*np.sinh(alpha1)/(np.cosh(alpha1)-np.cos(theta))
        yyy = a*np.sinh(alpha2)/(np.cosh(alpha2)-np.cos(theta))
   
    # angle: 쌍극좌표계 각 
        angle.append(theta)
        x.append(xx)
        xs.append(xxx)
        y.append(yy)
        my.append(-yy)
        ys.append(yyy)
        mys.append(-yyy)
        
    print('쌍극좌표계 각 = ', angle)
    print('쌍극좌표계 x좌표 = ', x)
    print('쌍극좌표계 y좌표 = ', y)
    

# 2. 상면의 2개의 원과 직교하는 원들의 부분 작도  
    alphamin = 5 ; alphamax = 10
    etha1 = (1/6)*np.pi ; etha2 = -(1/6)*np.pi
    etha3 = (6/7)*np.pi ; etha4 = -(6/7)*np.pi
    
    etha5 = (1/4)*np.pi ; etha6 = -(1/4)*np.pi
    etha7 = (3/4)*np.pi ; etha8 = -(3/4)*np.pi 
    
    delalpha = alphamax - alphamin
    
    mm = 2*n+1
    xa = [] ; ya = [] ; mya = [] ; xsa = [] ; ysa = [] ; mysa = []
    za = [] ; wa = [] ; mwa = [] ; zsa = [] ; wsa = [] ; mwsa = []
    ga = [] ; ha = [] ; mha = [] ; gsa = [] ; hsa = [] ; mhsa = []
    
    plt.figure()

    for i in range(mm):
        
        xsi = i*delalpha/mm
        xxa = a*np.sin(etha1)/(np.cosh(xsi)-np.cos(etha1))
        xxxa = a*np.sin(etha2)/(np.cosh(xsi)-np.cos(etha2))
        yya = a*np.sinh(xsi)/(np.cosh(xsi)-np.cos(etha1))
        yyya = a*np.sinh(xsi)/(np.cosh(xsi)-np.cos(etha2))
        
        zza = a*np.sin(etha3)/(np.cosh(xsi)-np.cos(etha3))
        zzza = a*np.sin(etha4)/(np.cosh(xsi)-np.cos(etha4))
        wwa = a*np.sinh(xsi)/(np.cosh(xsi)-np.cos(etha3))
        wwwa = a*np.sinh(xsi)/(np.cosh(xsi)-np.cos(etha4))
        
        gga = a*np.sin(etha5)/(np.cosh(xsi)-np.cos(etha5))
        ggga = a*np.sin(etha6)/(np.cosh(xsi)-np.cos(etha6))
        hha = a*np.sinh(xsi)/(np.cosh(xsi)-np.cos(etha5))
        hhha = a*np.sinh(xsi)/(np.cosh(xsi)-np.cos(etha6))
        
        # 아래 위 대칭 원과 오뚜기 형상 작도
        xa.append(xxa)  ; xsa.append(xxxa)
        ya.append(yya)  ; mya.append(-yya)
        ysa.append(yyya); mysa.append(-yyya)
        
        # 마주보는 호 작도 계산
        za.append(zza)  ; zsa.append(zzza)
        wa.append(wwa)  ; mwa.append(-wwa)
        wsa.append(wwwa); mwsa.append(-wwwa)
        
        ga.append(gga)  ; gsa.append(ggga)
        ha.append(hha)  ; mha.append(-hha)
        hsa.append(hhha); mhsa.append(-hhha)
        
        # 그래프 박스 작도       
        xlow = -4.5 ; xhigh = +4.5 ; ylow = -3.0 ; yhigh = +3.0
        plt.xlim(xlow, xhigh) ; plt.ylim(ylow, yhigh)
        
        # x, y 축 작도
        plt.plot((xlow, xhigh), (0.0, 0.0)) ; plt.plot((0.0, 0.0), (ylow, yhigh))
        
        # xsi 상수인 원 아래 위 대칭 작도
        plt.plot(x, y)  ; plt.plot(x, my) ; plt.plot(xs, ys); plt.plot(xs, mys)
        
        # ehta 상수인 원 반만 작도 후 mya, mysa 사용 대칭 작도
        plt.plot(xa,ya)   ; plt.plot(xa, mya) ; plt.plot(xsa, ysa); plt.plot(xsa, mysa)
        
        # mwa, mwsa 사용 작은 호 마주 작도
        plt.plot(za,wa)   ; plt.plot(za, mwa)
        plt.plot(zsa, wsa); plt.plot(zsa, mwsa)
        
        # ehta 상수(pi/4)인 원 반만 작도 후 mha, mhsa 사용 아래 대칭 작도
        plt.plot(ga, ha)   ; plt.plot(ga, mha)
        plt.plot(gsa, hsa) ; plt.plot(gsa, mhsa)

    #d_R = 1.543
    #R = a * 1.0/ np.arcsin(alpha1)
    #d = a * 1/np.tanh(alpha1)

    #print('(d/R) = ', d_R)    #print(' Radius = ', R)
    #print('Distance = ', d)    #plt.plot(x,y)    #plt.plot(xs , ys)

a = 1.0
n = 50
bp = bpplot(a, n)