さて, 数値計算で得られた計算データを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
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
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)
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
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
これで, 綺麗な楕円体が描けると思います.
こんな感じですかね.
もうちょっと頑張ったら, これの断面図をグロテスクに描いたりもできますが,
まぁアクセス数次第ですかね.
レッツ!!!等値面！！！！！！！