表面波の伝搬計算とクラドニ図形の作成

楕円形もしくは長方形の板における表面波の反射を計算してクラドニ図形を作成するProcedureです。各時間の表面波の伝搬も保存しています。

Function SurfaceWave(RorC,radius,ratio)

Variable RorC	//境界条件が長方形か円かの選択 0なら楕円、それ以外なら長方形
Variable radius,ratio	//長方形及び楕円の長辺の長さと縦横比
Variable freq,m,k,att
	freq	= 0.3
	m		= 0.5
	k		= 1
	att		= 2e-4
Variable amp,wid	 //初期条件
	amp		= 1
	wid		= 0.3
Variable Num_l,l_xy,l_start,d_l, Num_t,d_t
	Num_l		= 250
	l_xy		= 10
	l_start	= -l_xy/2
	Num_t		= 500//1000
	d_t			= 0.5
	d_l		= l_xy/Num_l
Variable Start_x,Start_y,cell_x,cell_y
	Start_x = 0
	Start_y = 0
	cell_x	= (Start_x-l_start)*d_l
	cell_y	= (Start_y-l_start)*d_l
Variable f_spr,force,accer
Variable i,j,t

String FilmName = "film_Mat"
	Make/o/d/n = (Num_l,Num_l,Num_t) $(FilmName)
	Wave MatWave = $(FilmName)
		setscale/p x,l_start,d_l,MatWave
		setscale/p y,l_start,d_l,MatWave
			MatWave = 0
	Make/o/d/n = (Num_l,Num_l) film_start,film_z,film_v,Chladni
		setscale/p x,l_start,d_l,film_start,film_z,film_v,Chladni
		setscale/p y,l_start,d_l,film_start,film_z,film_v,Chladni
			film_v			= 0
			film_start	= amp*exp(-((sqrt((x-start_x)^2+(y-start_y)^2))/wid)^2)		//Gaussian分布を持つ初期条件を仮定
			film_z			= film_start
			Chladni		=0

print "Start: "+time()
Variable posi_x,posi_y,posi_e
		for (t = 0 ; t < Num_t ; t += 1)	
			for(i = 1; i < Num_l-2; i += 1)
			for(j = 1; j < Num_l-2; j += 1)	
				f_spr	= -k*(4*film_z[i][j]-film_z[i-1][j]-film_z[i+1][j]-film_z[i][j+1]-film_z[i][j-1])
					force	= f_spr-film_v[i][j]*att
					accer	= force/m
					film_v[i][j] = film_v[i][j]+accer*d_t
			endfor
			endfor

			for(i = 0 ; i < Num_l ; i += 1)
			for(j = 0 ; j < Num_l ; j += 1)
				film_z[i][j] += film_v[i][j] * d_t
				posi_x	= 	ABS(l_start+i*d_l)
				posi_y	= 	ABS(l_start+j*d_l)
				posi_e	= sqrt((posi_x/radius)^2+(posi_y/(radius/ratio))^2)
				if(RorC == 0)
					if(posi_e>1)				
						film_z[i][j] = 0
					endif
				else
					if((posi_x>radius)||(posi_y>radius*ratio))
						film_z[i][j] = 0
					endif
			endif
				MatWave[i][j][t] = film_v[i][j]
			endfor
			endfor
			Chladni += ABS(film_z)

			film_z += film_start*cos(freq*t*d_t)	//CW波源を仮定。初期条件からの表面波の広がりのみを考える場合はsかうじょ

		endfor
		Chladni	= 1/Chladni
		killwaves film_v,film_z
print "Finish: "+time()
end

保存された表面波の伝搬は三次元行列形式で保存されていて表示しにくいため以下の表示用/動画出力用procedureで確認しています。

Function Disp(t)

Variable t
Variable i,j

Variable range =1
Variable Num_l,Num_t,d_t,d_l,l_xy
	l_xy	= 10
	Num_l	= 250
	d_l	= l_xy/Num_l

String DName = "Film_Mat"
	Wave DWave = root:$(DName)

Make/o/d/n = (Num_l,Num_l) XYplane
	setscale/p x,-l_xy/2,d_l,XYPlane
	setscale/p y,-l_xy/2,d_l,XYPlane

	for(i = 0 ; i < Num_l; i += 1)
	for(j = 0 ; j < Num_l; j += 1)
		XYplane[i][j] = DWave[i][j][t]
	endfor
	endfor
	
Display/k=1/W=(0,0,300,300)
AppendImage XYPlane
ModifyImage XYplane ctab= {-range,range,RedWhiteBlue,0}
setaxis bottom -l_xy/2,l_xy
setaxis left -l_xy/2,l_xy
ModifyGraph nticks=0,btLen=1
set()
end

Function Makemovie()
    String MoviePath = "C:Users:PCUser:Desktop:Chladni.mov"
    //String WindowName = "Pict0000"
   Disp(0)
    NewMovie/F=20/CF=10/O as MoviePath
   Variable i
   for (i=0 ; i<700 ; i+=2)
		Disp(i)
      Doupdate
		AddMovieFrame
		killwindow #
   endfor
    CloseMovie
    //PlayMovie as MoviePath
End