MATLABを用いた3次元等値面の描き方

さて, 数値計算で得られた計算データを3次元で美しく描きたくなると思います.
なると思います.
MATLABで等値面を描くと, とても綺麗になるのでおすすめです.
ということで, 描き方について紹介したいと思います.
まぁ, 手始めに楕円体ですかね.

program main
implicit none

integer(8)::na,nb,nc,intangle
real(8)::x,y,z,V,theta
integer(8),parameter::no=50
real(8)::a,b,c
real(8)::ss(1:3)

Open(50,FILE='data.dat',STATUS='unknown')

a=0.6d0
b=1.2d0
c=0.6d0

do nc=-no,no
write(*,*)nc+no,"/",no*2
do na=-no,no
do nb=-no,no

ss(1)=dble(na)/dble(no)
ss(2)=dble(nb)/dble(no)
ss(3)=dble(nc)/dble(no)

V=(b**2d0 *c**2d0 *SS(1)**2d0)+(c**2d0 *a**2d0 *SS(2)**2d0)+(a**2d0 *b**2d0 *SS(3)**2d0) !楕円体の式

write(50,'(3f10.5,2f11.7)')ss(1),ss(2),ss(3),V
end do
end do
end do

close(50)

end do
end program main

今回は, 101x101x101で計算しますが, これらの範囲は自由に変更できます.
実際に使う際は, vの式を任意の式に変更するだけです.

x, y, zの回す順番だけは, このプログラムコードに合わせて下さい.

さて, data.datができました. 次に, MATLABで読み込めるようにファイルを変換します.

program main
implicit none

integer(8)::nx,ny,nz,xx,yy,zz
real(8),allocatable::x(:,:,:)
real(8),allocatable::y(:,:,:)
real(8),allocatable::z(:,:,:)
real(8),allocatable::v(:,:,:)

nx=50
ny=50
nz=50

open(62,file="data.dat",status="old")!読み込みファイル
open(64,file="x.txt",status="unknown")!出力ファイル
open(65,file="y.txt",status="unknown")!出力ファイル
open(66,file="z.txt",status="unknown")!出力ファイル
open(67,file="v.txt",status="unknown")!出力ファイル

allocate(x(-nx:nx,-ny:ny,-nz:nz))
allocate(y(-nx:nx,-ny:ny,-nz:nz))
allocate(z(-nx:nx,-ny:ny,-nz:nz))
allocate(v(-nx:nx,-ny:ny,-nz:nz))

do zz=-ny,ny
do yy=-nx,nx
do xx=-nz,nz
read(62,*)x(xx,yy,zz),y(xx,yy,zz),z(xx,yy,zz),v(xx,yy,zz)
end do
end do
end do

do xx=-nx,nx
write(64,*) x(xx,:,:)
write(65,*) y(xx,:,:)
write(66,*) z(xx,:,:)
write(67,*) v(xx,:,:)
end do

deallocate(x)
deallocate(y)
deallocate(z)
deallocate(v)

close(62)
close(64)
close(65)
close(66)
close(67)

end program

これで, x.txt, y.txt, z.txt, v.txtが出力されます.
nx, ny, nzの値は, 範囲設定を変えた場合は適宜それに合わせて書き換えて下さい.
さて, 準備ができたので, MATLABを起動しましょう.

x.txt, y.txt, z.txt, v.txtの階層に移動して,
次のコマンドを入力してこれらのファイルを読み込みます.
次のコマンドはスクリプトファイルにしちゃうと楽ですね.

xn=50*2+1 %xのデータ数
yn=50*2+1 %yのデータ数
zn=50*2+1 %zのデータ数
x=load('x.txt')
x2=reshape(x,xn,yn,zn)
y=load('y.txt')
y2=reshape(y,xn,yn,zn)
z=load('z.txt')
z2=reshape(z,xn,yn,zn)
v=load('v.txt')
v2=reshape(v,xn,yn,zn)

コードの内容は,
まず, 先ほどのデータを読み込みます.
次に, 読み込んだデータをちょっと整理してます.

さて, あとは描くだけです！！
次のコマンドを入力していきましょう.

E=0.5 %Eが描きたい等値面の値
p = patch(isosurface(x2,y2,z2,v2,E));
isonormals(x2,y2,z2,v2,p)
set(p,'FaceColor','green','EdgeColor','none'); %ここで色変更可
daspect([1,1,1])
view(30,30); axis([-1.5 1.5 -1.5 1.5 -1.5 1.5]) %axisで表示範囲を設定
camlight ('headlight') %光の当てる方向
camlight ('right') %光の当てる方向
lighting gouraud

これで, 綺麗な楕円体が描けると思います.