쌍극자형 지구 자기장 그리기(Python) (Dipole model of the Earth's magnetic field coding)

안녕하십니까, 현재 천문학을 전공하고 있는 필명 AstroPenguin이라고 합니다.

천문학이라 하면 망원경으로 별을 보는 학문이라고 생각하시는 분들이 많겠지만

현대의 천문학은 유도된 공식을 python을 통해 프로그래밍하여 그래프로 나타내기도 하고

문자형식의 파일을 python을 통해 불러와 천체의 이미지를 확인하기도 합니다.

 

예전의 천문학과 달리 현재는 프로그래밍과 시뮬레이션을 통한 연구가 많이 이루어지고 있습니다.

제가 실습을 통해 이러한 과정이 어떻게 이루어지고 있는지 보여드리도록 하겠습니다.

이번에는 쌍극자형 지구 자기장에 대해 알아보고 python을 통해 그려보는 시간을 가져보도록 하겠습니다.

 

지구의 내부에는 핵이 존재합니다. 철과 니켈 등 금속 물질로 이루어진 액체 상태의 외핵이 온도차에 의해 대류를 합니다. 자유전자를 가진 금속물질이 이동을 하면 유도 전류가 발생하는데, 이 때문에 외핵 내부에 자기장이 형성됩니다. 그리고 지구는 항상 자전을 하는데, 이로 인해 외핵이 원운동하게 되며 거듭 유도 전류가 발생하여 지구 자기장이 형성되는 원리입니다.

 

이러한 지구 자기장은 쌍극자의 형태를 가집니다. 쌍극자란 말 그대로 +와 -의 양쪽의 전하를 모두 가진 물질을 말합니다. 그리고 +와 -극 사이에서 전류나 자기장 유도되면서 다음과 같은 자기력선의 형태가 발생합니다.

이를 우주에서 본다면 다음과 같은 형태입니다. 원래대로라면 양쪽이 대칭인 형태여야하지만, 강력한 에너지를 지닌 태양풍에 의해 태양을 바라보는 지구 자기장은 태양풍의 힘을 받아 압축이 됩니다.

관측적인 지구의 자기장을 표현하려면 태양풍의 영향도 고려해야하지만,

현재 제 대학생의 수준으로는 한계가 있어 태양풍이 없을 때의, 양쪽이 대칭적인

쌍극자형 지구 자기장을 포현해보도록 하겠습니다.

지금부터 시작하겠습니다!!!

 

우선 사용할 기능인 numpy와 matplotlib.pyplot를 불러옵니다.

import numpy as np
import matplotlib.pyplot as plt

지구와 지구 자기력선 사이의 거리인 r은 다음과 같은 관계가 존재합니다.

r  = L X RE

r을 지구 반지름인 RE의 실수배로 표현할 수 있습니다.

그러므로 기준이 되는 임의의 req를 RE에 6.6(정지궤도위성의 고도)(많이 사용되는 L값)을 곱하여 사용합니다.

r과 req, 그리고 세타와의 관계를 이용하여 r을 구하는 함수를 정의합니다.

def r(theta):
    RE = 6378000
    req = 6.6 * RE
    return req*(np.cos(theta)**2.0)

반지름을 뜻하는 r에 sin과 cos을 곱해준다면 무슨 일이 일어날까요?

바로 r과 세타를 x와 y로 표현할 수 있게 됩니다.

r(theta)에 각각 sin과 cos을 곱해주어 x와 y를 선언합니다.

그리고 좀 전에 req를 6.6 X 지구반지름(RE)으로 임의로 지정해주었는데, 사실은 req의 값은 매우 많습니다.

               ( 자기력선이 존재하는 거리 = L X 지구반지름(RE) , L은 임의의 실수) (즉 L은 여러 값을 가진다)

이를 표현하기 위해 저는 for문을 사용하여 r(theta)를 임의의 실수로 나누어주는 반복문을 돌렸습니다.

                ( r(theta)는 req에 비례하기 때문 )

그 후 plot을 하고 색깔을 지정해주어 자기력선을 표현합니다.

n = 8
for i in range(n):
    theta = np.linspace(0.0, 2.0*np.pi,100)
    ni = np.linspace(0,2,8)
    x = (r(theta)/(ni[i]+1))*np.cos(theta)
    y = (r(theta)/(ni[i]+1))*np.sin(theta)
    plt.plot(x,y,c='blue')
    x2 = (r(theta)*(ni[i]+1))*np.cos(theta)
    y2 = (r(theta)*(ni[i]+1))*np.sin(theta)
    plt.plot(x2,y2,c='blue') 

 

마지막 과정은 그래프의 정가운데에 파란색 점을 찍어 지구를 표현하는 과정입니다.

dot_x와 dot_y를 통해 그래프의 중점을 지정해주고

matplotlib의 scatter기능을 통해 점을 찍어줍니다.

이 때 s = 1800으로 놓는 이유는 자기력선간의 간격과 지구의 크기를 실제와 같게 맞춰주기 위해서입니다.

그리고 마지막으로 c = 'blue' 또는 'b'를 통해 색깔을 파란색으로 지정해줍니다.

순서를 바꾸어서 헷갈리는 분들도 계실텐데 plt.title와 plt.xlim등을 이용해 그래프의 제목과 출력 범위를 설정하는것 잊지 말아주세요!!

dot_x = [0.0]
dot_y = [0.0]
plt.title('Earth dipole magnetic field')
plt.xlim(-0.5e8, 0.5e8)
plt.ylim(-2e7, 2e7)
plt.scatter(dot_x,dot_y,s=1800,c='blue')
plt.show()

두근두근!! 이제 실행을 하시면 쌍극자 모양 지구 자기장이 나옵니다~~~(박수)

출력된 쌍극자형 지구자기장

저는 r에 관한 공식과 L에 관하여 지구의 쌍극자를 표현하였지만,

r벡터나 각도 세타 등 여러 가지 표현법이 존재합니다. 사람마다 생각하는 바가 다르기 때문에

코드를 짜는 방식도 각각 다를 것입니다.

여러분도 한 번 각자만의 생각으로 지구 자기장을 표현하심이 어떨까요?

참고할 수 있는 코딩이 잘 된 예시를 아래 남겨드리겠습니다!

scipython.com/blog/visualizing-the-earths-magnetic-field/

Visualizing the Earth's dipolar magnetic field