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

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

hcu

こんな感じですかね.

もうちょっと頑張ったら, これの断面図をグロテスクに描いたりもできますが,
まぁアクセス数次第ですかね.

レッツ!!!等値面!!!!!!!


コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です