Last updated: Aug. 31, 2010

You're th visitor.

マイスナー効果（日本語説明つき）

磁場を超伝導転移後に印加しても超伝導転移前に印加してもマイスナー効果が現れることを説明する図です。「最先端科学の体験型学習講座」Experienced-based Learning Course for Advanced Science（ELCAS）のブックレット用に作成しました。描画の基本的な部分は「マイスナー効果」をご覧ください。

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

なお、このスクリプトの実行結果は使用するバージョンによって異なる場合があるようです。上の図は角藤さんのところからダウンロードしたgnuplot4.5を使用して描画しました。手元のgnuplot4.2.6では常伝導状態の球が現れないということが起こりました。

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

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

MODE = "PLOT"        # プロットのみ（MAKEDATAで一回は実行していないとエラーがでる）
MODE = "PRINT"        # プロットのみ（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 ( -1 "pink", -0.000001 "dark-red", 0 "dark-blue" , 1 "cyan")
set cbrange [-1:1]

set style arrow 1 size graph 0.02,20,50   filled linewidth 8

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


set multiplot

# 超伝導状態＋磁場なしの描画

set tmargin screen 0.8
set bmargin screen 0.2
set rmargin screen 0.4
set lmargin screen -0.2

set arrow 1 from screen 0.3,0.6 rto screen 0.1,0.1 as 1 lc rgb "blue"
set label 1 at screen 0.345,0.67 right font "GothicBBB-Medium-90pv-RKSJ-H,36" "冷却" tc rgb "blue"

splot theta_phi_file using (a*sin($1)*cos($2)):(a*sin($1)*sin($2)):(a*cos($1)):(1+a*sin($1)*cos($2))*(-0.5) with pm3d notitle

if (MODE eq "PLOT") pause -1


# 超伝導状態＋磁場ありの描画

set tmargin screen 0.8
set bmargin screen 0.2
set rmargin screen 1.2
set lmargin screen 0.6

set arrow 1 from screen 0.3,0.4 rto screen 0.1,-0.1 as 1 lc rgb "dark-green"
set label 1 at screen 0.34,0.33 right font "GothicBBB-Medium-90pv-RKSJ-H,36" "磁場印加" tc rgb "dark-green"

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




if (MODE eq "PLOT") pause -1


# 常伝導状態＋磁場ありの描画


set tmargin screen 0.4
set bmargin screen -0.2
set rmargin screen 0.8
set lmargin screen 0.2

set arrow 1 from screen 0.6,0.3 rto screen 0.1,0.1 as 1 lc rgb "blue"
set label 1 at screen 0.66,0.38 right font "GothicBBB-Medium-90pv-RKSJ-H,36" "冷却" tc rgb "blue"

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

if (MODE eq "PLOT") pause -1



# 常伝導状態＋磁場なしの描画

set tmargin screen 1.2
set bmargin screen 0.6
set rmargin screen 0.8
set lmargin screen 0.2

set arrow 1 from screen 0.6,0.7 rto screen 0.1,-0.1 as 1 lc rgb "dark-green"
set label 1 at screen 0.65,0.62 right font "GothicBBB-Medium-90pv-RKSJ-H,36" "磁場印加" tc rgb "dark-green"

splot theta_phi_file using (a*sin($1)*cos($2)):(a*sin($1)*sin($2)):(a*cos($1)):(1+a*sin($1)*cos($2))*(+0.5) with pm3d notitle


unset multiplot

# 画像ファイルの作成
if (MODE eq "PRINT" && ( exist("TERM")==0 || TERM eq "WIN") ) pause -1;\
     set out "./Meissner_Japanese.eps";\
     set term post eps color solid "GothicBBB-Medium-90pv-RKSJ-H" 20 size 7in,7in;\
     TERM = "PS";\
     reread

undefine TERM
set out
set term pop

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