Last updated: Aug. 31, 2010

You're th visitor.

マイスナー効果

マイスナー効果による磁束線排除の様子を図示したものを作製しました。この例では超伝導体は球形の完全反磁性体（比透磁率が-1）とし、「平行磁場中に球形の磁性体をおいた場合の磁場分布は球の中心に適当な大きさの磁気ダイポールをおいた場合の磁場分布に等しい」ことを用いて磁場分布を計算しています。磁束線は適当なスタート地点から磁場の方向に少しずつ折れ線を書いていく方法で近似的に描画しました。

図をかっこよくするために、超伝導体の描画には3次元グラフの表面の色づけの項で説明した「光の当たっている物体に見える色のつけ方」を利用しています。

下のスクリプトを実行すると、こちらのようなPostScriptファイルが出来ます。上の図はそのPostScriptファイルを切り取ったものです。

なお、スクリプトの実行にはgnuplot4.3以降が必要です。また、このスクリプトを初めて実行する際は必ずMODE変数を"MAKETABLE"にした上で実行してください。2回目以降は図の大きさを変えたりするだけであればMODE変数を"PLOT"にして（つまりMODE = "MAKETABLE"の行をコメントアウトして）実行すれば、所要時間が短くてすみます。

#version 4.3

############################################################
####   以下のどちらかをコメントアウトしてください。    #####
############################################################

MODE = "PLOT"        # プロットのみ
                     #（MAKEDATAで一回は実行していないとエラーがでる）
MODE = "MAKEDATA"    # データファイルを作ってプロットする
                     #（時間かかるが一度はこのモードで実行しないと
                     #  データファイルが作れないのでエラーとなる）

############################################################


set angles radians

mu_ = pi*4e-7 #真空の透磁率
a=1.0         #超伝導体の半径
H0=1.0        #外部磁場 （y方向）
mus = 0.0     #比透磁率
m = 4.0*pi*a**3*mu_*(mus-1)/(mus+2)*H0           #仮想的なダイポールの大きさ

Hr(r,theta)=   m/(2.0*pi*mu_)/r**3*cos(theta)    #ダイポールの磁場のr成分
Ht(r,theta)=   m/(4.0*pi*mu_)/r**3*sin(theta)    #ダイポールの磁場のtheta成分

Hxr(r,theta) = Hr(r,theta)*sin(theta)+Ht(r,theta)*cos(theta)  #ダイポールの磁場のx成分
Hyr(r,theta) = Hr(r,theta)*cos(theta)-Ht(r,theta)*sin(theta)  #ダイポールの磁場のy成分

Hx(x,y)=   Hxr(sqrt(x**2+y**2), atan2(x,y))                   # 全磁場のx成分
Hy(x,y)=H0+Hyr(sqrt(x**2+y**2), atan2(x,y))                   # 全磁場のy成分

atany(x,y)=(y==0.0 ? pi*0.5 : (atan(x/y) >= 0.0 ? atan(x/y) : atan(x/y)+pi))

vlength = 0.0001  # 磁束線を描く際のステップ幅（短いほど正確だが計算時間がかかる）

paramfile      = "meissner_param.dat"  # 磁束線を描く際のパラメーターtを保存しておくファイル
outputfile     = "meissner_output.dat" # 磁束線を保存するファイル
theta_phi_file = "sphere.dat"          # 球を描くためのtheta,phi座標を保存しておくファイル

#
# データファイルの作成 （MODE が "MAKEDATA" の場合のみ実行される）
#
if (MODE eq "MAKEDATA") set samples 100000
if (MODE eq "MAKEDATA") set table paramfile
if (MODE eq "MAKEDATA") plot 0
if (MODE eq "MAKEDATA") unset table

if (MODE eq "MAKEDATA") set table outputfile
if (MODE eq "MAKEDATA") set multiplot

yini=-4
xini=-1.75
ev=10
if (MODE eq "MAKEDATA") \
     plot paramfile every ev using (xx = ($0==0 ? xini : xx+dx), yy = ($0==0 ? yini : yy+dy), \
     dx=(Hx(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     dy=(Hy(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     xx):(yy) w l

yini=-4
xini=-1.5
if (MODE eq "MAKEDATA") \
     plot paramfile every ev using (xx = ($0==0 ? xini : xx+dx), yy = ($0==0 ? yini : yy+dy), \
     dx=(Hx(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     dy=(Hy(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     xx):(yy) w l

yini=-4
xini=-1.25
if (MODE eq "MAKEDATA") \
     plot paramfile every ev using (xx = ($0==0 ? xini : xx+dx), yy = ($0==0 ? yini : yy+dy), \
     dx=(Hx(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     dy=(Hy(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     xx):(yy) w l

yini=-4
xini=-1.0
if (MODE eq "MAKEDATA") \
     plot paramfile every ev using (xx = ($0==0 ? xini : xx+dx), yy = ($0==0 ? yini : yy+dy), \
     dx=(Hx(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     dy=(Hy(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     xx):(yy) w l

yini=-4
xini=-0.75
if (MODE eq "MAKEDATA") \
     plot paramfile every ev using (xx = ($0==0 ? xini : xx+dx), yy = ($0==0 ? yini : yy+dy), \
     dx=(Hx(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     dy=(Hy(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     xx):(yy) w l

yini=-4
xini=-0.5
ev = 3
if (MODE eq "MAKEDATA") \
     plot paramfile every ev using (xx = ($0==0 ? xini : xx+dx), yy = ($0==0 ? yini : yy+dy), \
     dx=(Hx(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     dy=(Hy(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     xx):(yy) w l

yini=-4
xini=-0.25
if (MODE eq "MAKEDATA") \
     plot paramfile every ev using (xx = ($0==0 ? xini : xx+dx), yy = ($0==0 ? yini : yy+dy), \
     dx=(Hx(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     dy=(Hy(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     xx):(yy) w l

yini=-4
xini=-0.05
ev=1
if (MODE eq "MAKEDATA") \
     plot paramfile every ev using (xx = ($0==0 ? xini : xx+dx), yy = ($0==0 ? yini : yy+dy), \
     dx=(Hx(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     dy=(Hy(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     xx):(yy) w l

yini=-4
xini=0.05
if (MODE eq "MAKEDATA") \
     plot paramfile every ev using (xx = ($0==0 ? xini : xx+dx), yy = ($0==0 ? yini : yy+dy), \
     dx=(Hx(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     dy=(Hy(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     xx):(yy) w l

yini=-4
xini=0.25
ev=3
if (MODE eq "MAKEDATA") \
     plot paramfile every ev using (xx = ($0==0 ? xini : xx+dx), yy = ($0==0 ? yini : yy+dy), \
     dx=(Hx(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     dy=(Hy(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     xx):(yy) w l

yini=-4
xini=0.5
if (MODE eq "MAKEDATA") \
     plot paramfile every ev using (xx = ($0==0 ? xini : xx+dx), yy = ($0==0 ? yini : yy+dy), \
     dx=(Hx(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     dy=(Hy(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     xx):(yy) w l

ev=10
yini=-4
xini=0.75
if (MODE eq "MAKEDATA") \
     plot paramfile every ev using (xx = ($0==0 ? xini : xx+dx), yy = ($0==0 ? yini : yy+dy), \
     dx=(Hx(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     dy=(Hy(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     xx):(yy) w l

yini=-4
xini=1.0
if (MODE eq "MAKEDATA") \
     plot paramfile every ev using (xx = ($0==0 ? xini : xx+dx), yy = ($0==0 ? yini : yy+dy), \
     dx=(Hx(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     dy=(Hy(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     xx):(yy) w l

yini=-4
xini=1.25
if (MODE eq "MAKEDATA") \
     plot paramfile every ev using (xx = ($0==0 ? xini : xx+dx), yy = ($0==0 ? yini : yy+dy), \
     dx=(Hx(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     dy=(Hy(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     xx):(yy) w l

yini=-4
xini=1.5
if (MODE eq "MAKEDATA") \
     plot paramfile every ev using (xx = ($0==0 ? xini : xx+dx), yy = ($0==0 ? yini : yy+dy), \
     dx=(Hx(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     dy=(Hy(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     xx):(yy) w l

yini=-4
xini=1.75
if (MODE eq "MAKEDATA") \
     plot paramfile every ev using (xx = ($0==0 ? xini : xx+dx), yy = ($0==0 ? yini : yy+dy), \
     dx=(Hx(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     dy=(Hy(xx,yy)/sqrt(Hx(xx,yy)**2+Hy(xx,yy)**2)*vlength*ev), \
     xx):(yy) w l
unset multiplot
unset table


#
# theta-phi座標ファイルの作成
#
set xrange [0:pi]
set yrange [-1*pi:1*pi]
set samples 50
set isosamples 50
set table theta_phi_file
splot 0.0
unset table


#
# 磁束線と超伝導体の描画 
#
set view map
set pm3d depthorder
set size ratio 1
set palette defined ( 0 "dark-blue" , 1 "cyan")

set xrange [-4:4]
set yrange [-4:4]
unset xtics
unset ytics
unset colorbox

# 背景を黒にする
set object 1 rect from screen 0,0 to screen 1,1 \
           fc rgb "black" fillstyle solid 1.0

splot outputfile using 1:(stringcolumn(3)eq "i" ? $2 : 1/0):(0) every 10 w l lw 2 lc rgb "green" notitle,\
      theta_phi_file using (a*sin($1)*cos($2)):(a*sin($1)*sin($2)):(a*cos($1)):(a*sin($1)*cos($2)) with pm3d notitle


pause -1



# 画像ファイルの作成
set out "Meissner.eps"
set term post color solid size 7in,7in
replot
set out
set term pop

ついでに、常伝導状態の図は以下のスクリプトで作製できます。こちらのようなPostScriptファイルが出来ます。

a=1.0      # 超伝導体の半径
theta_phi_file = "sphere.dat"  # 球を描くためのtheta,phi座標を保存しておくファイル

set angles radians

#
# theta-phi座標ファイルの作成
#
set xrange [0:pi]
set yrange [-1*pi:1*pi]
set samples 50
set isosamples 50
set table theta_phi_file
splot 0.0
unset table


#
# 磁束線と超伝導体の描画 
#
set view map
set pm3d depthorder
set size ratio 1
set palette defined ( 0 "dark-blue" , 1 "cyan")

set xrange [-4:4]
set yrange [-4:4]
set urange [-4:4]
unset xtics
unset ytics
unset colorbox


# 背景を黒くする
set object 1 rect from screen 0,0 to screen 1,1 fc rgb "black" fillstyle solid 1.0

splot theta_phi_file using (a*sin($1)*cos($2)):(a*sin($1)*sin($2)):(a*cos($1)):(a*sin($1)*cos($2)) with pm3d notitle,\
     "+" using (-1.75):1:(0) w l lw 2 lc rgb "green" notitle,\
     "+" using (-1.50):1:(0) w l lw 2 lc rgb "green" notitle,\
     "+" using (-1.25):1:(0) w l lw 2 lc rgb "green" notitle,\
     "+" using (-1.00):1:(0) w l lw 2 lc rgb "green" notitle,\
     "+" using (-0.75):1:(0) w l lw 2 lc rgb "green" notitle,\
     "+" using (-0.50):1:(0) w l lw 2 lc rgb "green" notitle,\
     "+" using (-0.25):1:(0) w l lw 2 lc rgb "green" notitle,\
     "+" using ( 1.75):1:(0) w l lw 2 lc rgb "green" notitle,\
     "+" using ( 1.50):1:(0) w l lw 2 lc rgb "green" notitle,\
     "+" using ( 1.25):1:(0) w l lw 2 lc rgb "green" notitle,\
     "+" using ( 1.00):1:(0) w l lw 2 lc rgb "green" notitle,\
     "+" using ( 0.75):1:(0) w l lw 2 lc rgb "green" notitle,\
     "+" using ( 0.50):1:(0) w l lw 2 lc rgb "green" notitle,\
     "+" using ( 0.25):1:(0) w l lw 2 lc rgb "green" notitle,\
     "+" using ( 0.00):1:(0) w l lw 2 lc rgb "green" notitle

pause -1

# 画像ファイルの作成
set out "Meissner_normal.eps"
set term post color solid size 7in,7in
replot
set out
set term win

• トップページ
• 自己紹介・経歴
• 研究紹介
  • 最近の成果
• 論文・学会発表
  • 論文リスト
  • 発表リスト
• 教育
  • 教育
  • 課題研究Q4
  • 課題演習B4
• アウトリーチ
• プログラミング
  • xyzzy
    • wgnuplot-mode
    • kyotex-mode
    • eng-comp
  • gnuplot
  • JavaScript
    • プレゼンタイマー
• その他
  • リンク集
  • ニュース一覧
  • 更新履歴
• English Top Page