Matplotlib을 이용한 탄성파 자료처리 그림 그리기

이전 글에서는 SU 명령어들을 이용해 탄성파 자료처리 결과 확인용 그림을 그리는 방법을 살펴보았습니다. 이번에는 Python의 Matplotlib을 이용하여 그린 그림 예제들을 보겠습니다. 그림은 IPython Processing 모듈을 이용해 그렸으며, 그릴 때 사용한 코드는 github에서 볼 수 있습니다.

속도모델, 구조보정 영상

우선, 다음과 같이 이진파일로부터 2차원 속도모델과 구조보정 결과를 그릴 수 있습니다. 기본적으로 속도모델은 컬러, 구조보정 영상은 흑백으로 그리도록 했지만, 필요에 따라 코드를 수정해서 색상을 바꿀 수 있습니다. 색상을 바꾸고 싶을 경우 imshow 함수의 cmap 인자를 이용하면 됩니다.

%matplotlib inline
from pkprocess import *
import numpy as np

vel = np.fromfile("marm16km.drt", dtype=np.float32)

nx = 576
nz = 188
h = 0.016
vel.shape = (nx, nz)

plot_vel(vel, h)

png


plot_mig(vel,h)

png

공통송신원 모음, 스펙트럼

그리고 SU 파일로부터 공통송신원 모음이나 F-X, F-K 스펙트럼을 그릴 수 있습니다. 공통송신원 모음은 Wiggle trace 또는 이미지로 그릴 수 있고, 이미지 색상은 cmap으로 조절 가능합니다.


su = read_su("marm3000.su")

plot_wiggle(su, perc=97)

min=-616.05078125 max=613.4453125

png


plot_image(su, perc=97)

min=-616.05078125 max=613.4453125

png


plot_image(su, perc=97, cmap='bwr')

min=-616.05078125 max=613.4453125

png


specfx(su)

dt=0.004, fmax=125.0

png


specfk(su)

dt=0.004, fmax=125.0

dx=0.025, kmax=20.0

png

위의 그림들 모두 Matplotlib으로 그렸으므로, 수정이 필요할 경우 Matplotlib 문서를 참고하여 수정해서 사용하시면 되겠습니다.

탄성파 자료처리 그림 그리기

탄성파 자료처리 결과 그림을 쉽게 그리는 방법을 살펴보겠습니다.

탄성파 자료처리를 하다 보면 결과물을 그림으로 확인해야 하는 경우가 많습니다. 특별히 노력해서 그려야 하는 그림도 있지만 속도모델, 공통송신원모음 등 대부분의 그림은 거의 비슷한 명령으로 그릴 수 있습니다. 개인적으로 논문이나 발표자료에 넣을 그림을 그릴 때 Seismic Un*x(SU)를 많이 이용하는데, 몇 가지 자주 그리는 그림들을 쉽게 그릴 수 있도록 파이썬 모듈을 만들었습니다. 모듈은 gpl라이브러리에 포함되어 있습니다. 최근 python 3 용으로 수정하였습니다.

먼저 속도모델을 예로 들어보겠습니다. 그림을 그리기 위한 코드는 다음과 같습니다.


from gpl.psplot import plot

vel="marm16km.drt"
opt="n1=188 d1=0.016 d2=0.016 d1num=1 d2num=2"

plot.velocity_color("vel_color.png",vel,opt)

위 코드는 gpl.psplot 모듈에서 plot을 가져오고, marm16km.drt 파일로부터 opt 문자열의 옵션을 이용하여 vel_color.png 파일을 생성하는데, 컬러로 된 속도모델 그림으로 만들라는 코드입니다.

velocity_color는 그림 종류를 지정하는 명령인데, 현재 다음과 같은 명령들을 지원합니다.

  • velocity(target, source, option, unit=”km/s”)
  • velocity_color(target, source, option, unit=”km/s”)
  • gradient(target, source, option)
  • gradient_color(target, source, option)
  • migration(target, source, option)
  • contour(target, source, option)
  • seismogram(target, source, option)
  • spectrum(target, source, option)

위의 명령들은 SU를 이용해 해당 그림을 그리라는 명령으로, contourpscontour를 사용하고 나머지는 psimage를 사용합니다. 입력 파일이 SU 파일이라면 supscontour 또는 supsimage를 사용합니다.

인자들 중 target은 출력 파일, source는 입력 파일, option은 그림 그릴 때 사용할 옵션입니다. 그림 종류에 따라 기본적으로 몇 가지 옵션이 들어가있는데, n1, d1, d2와 같이 입력 파일에 따라 달라지는 옵션을 option에 넣어주면 됩니다. 그리고 기본 옵션을 덮어쓰고 싶은 경우에도 option에 추가해줍니다.

속도모델의 단위는 기본적으로 km/s로 지정해 놓았는데, 필요에 따라 바꿔서 사용할 수 있습니다. g/cc로 바꾸면 밀도 모델을 그릴 수도 있겠죠. migration은 snapshot을 그릴 때 사용할 수도 있습니다.

SU 명령은 기본적으로 eps 파일을 생성합니다. target을 eps 외의 다른 파일(png, tiff, jpg 등)로 지정하면ImageMagickconvert 명령을 이용해 eps 파일을 변환합니다.

따라서 본 모듈의 모든 기능을 이용하려면 Python, SU, ImageMagick이 필요합니다.

터미널에서 위의 코드를 실행했을 때 나오는 메시지는 다음과 같습니다.

psimage height=1.0 width=2.65 d2s=0.5 lwidth=0.1 lstyle="vertright" lheight=1.0
label2="Distance (km)" ghls="0.33,0.5,1" bps=24 whls="0,0.5,1" legend=1
bhls="0.67,0.5,1" labelsize=8 label1="Depth (km)" d1s=0.5  n1=188 d1=0.016
d2=0.016 d1num=1 d2num=2 < marm16km.drt > vel_color.eps

psimage: bclip=5.5 wclip=1.5

// adding velocity unit (km/s)

// fixing bounding box
Original:  %%BoundingBox: 66 41 353 207
Updated:   %%BoundingBox: 85 104 324 202

// converting .eps to .png ..

내용을 살펴보면 다음 순서로 실행됩니다.

  1. SU의 psimage 명령을 이용해 속도모델 eps 파일을 생성합니다. 옵션은 컬러 속도모델에 맞춰서 들어갑니다. 참고로, 그림 크기는 Geophysics 논문 기준에 맞춘 것입니다.
  2. km/s라는 단위를 넣어줍니다(postscript 수정).
  3. 그림 여백을 조절합니다(bounding box 수정).
  4. eps 파일을 png 파일로 수정합니다.

그리고, 결과물인 vel_color.png은 다음과 같습니다. output(vel_color.png)

아래 코드와 다른 그림 예시를 올리니 필요한 그림에 해당하는 명령을 사용하시면 되겠습니다.


from gpl.psplot import plot

vel="marm16km.drt"
opt="n1=188 d1=0.016 d2=0.016 d1num=1 d2num=2"

plot.velocity("vel.png",vel,opt+"lbeg=1.5 lend=5.5 lfnum=1.5")
plot.velocity_color("vel_color.png",vel,opt)
plot.velocity_color("density_color.png",vel,opt,unit="g/cc")
plot.gradient("grad.png",vel,opt)
plot.gradient_color("grad_color.png",vel,opt)
plot.migration("mig.png",vel,opt)
plot.contour("contour.png",vel,opt)

seismo="marm3000.su"
opt2="f2=0 d2=0.025 d1s=0.5 d2s=0.5"
plot.seismogram("seismo.png",seismo,opt2)

spec="marm3000fx.su"
plot.spectrum("spec.png",spec,opt2)

output(vel.png)

output(vel_color.png)

output(density_color.png)

output(grad.png)

output(grad_color.png)

output(mig.png)

output(contour.png)

output(seismo.png)

output(spec.png)

Muting scripts

SU 파일 뮤팅 스크립트를 공개하였습니다. SU 파일을 손으로 뮤팅할 때 모든 shot gather를 뮤팅하는 것이 아니라 건너 뛰면서 뮤팅한 후 건너뛴 부분은 interpolation하도록 만들었습니다. 해상 스트리머 탐사 자료에 사용할 수 있습니다. Github에서 받을 수 있습니다.

 

참고 문헌

Kim, A., D. Ryu, W. Ha, and C. Shin, 2015, An efficient first arrival picking procedure for marine streamer data, Geosystem Engineering,  18(5), 251-258.

IPython Processing

학부 탄성파 자료 처리 수업에서 사용하기 위한 파이썬 패키지입니다. 기존의 공개된 Matlab 패키지를 파이썬으로 변환하였습니다. 간단한 2차원 자료 처리(time processing)가 가능합니다.

육상, 해상 자료 처리 예제와 함께 공개하였고, Github에서 받으실 수 있습니다.

참고 문헌

하완수, 2015, 대화식 탄성파 자료 처리 수업을 위한 파이썬 패키지 개발, 한국자원공학회지, 52(4), 414-421.

SU 파일 입출력을 위한 포트란 라이브러리

포트란 사용자들이 Seismic Un*x의 SU 파일을 읽고, 수정하고, 쓸 수 있도록 하는 입출력 라이브러리를 작성하여 공개하였습니다.

Github에서 받으실 수 있습니다.

매뉴얼은 ReadTheDocs를 참고하세요.

참고문헌

하완수, 2015, SU 파일 입출력을 위한 포트란 라이브러리 개발, 한국자원공학회지, 52(1), 81-90.

Postscript language editing

Postscript로 만들어진 .ps 또는 .eps 파일은 앞의 글에서 보셨던 것처럼, 일반적인 text 편집기로 편집할 수 있는 ascii 파일입니다. 파일의 내용은 출력물을 만들어내는 postscript 언어죠. 따라서 postscript 언어를 알면 eps 그림 파일도 마음대로 편집할 수 있습니다. Postscript language를 배우고 싶으신 분은 Adobe site에 가셔서 매뉴얼을 받아보시면 됩니다. 여기서는 앞의 글에서 만들었던 파일에서 legend unit의 위치를 바꾸는 법만 살펴보도록 하겠습니다. 앞에서 보았던 그림은 다음과 같습니다.

before editing

before editing

위 그림에서 오른쪽 끝에 있는 “m/s”를 legend(scale bar) 위로 옮겨봅시다. 결과는 다음과 같습니다.(옮긴 후 bounding box도 바꿔줬습니다.)

after editing unit

after editing unit

위의 결과를 얻기 위해서는 .eps 파일을 열어서 ‘m/s’라는 문자열을 찾아 지워줍니다. 파일 끝에서 약간 앞에 두 개가 있을겁니다. 그런 후 아래의  코드를 .eps 파일 끝부분의 ‘showpage’ 명령 앞에 넣어줍니다.

%%%%% changed the position of unit
GS
270 190 TR
NP
/Helvetica findfont 8 scalefont setfont
0 0 0 setrgbcolor
21.96 -6.462 M
(m/s) SW exch -0.5 mul
exch -0.5 mul RM (m/s) SH
S
GR
%%%%%

%%%%%는 comment이고, ‘m/s’라는 문자열의 위치는 ‘GS’ 아래에 있는 두 개의 숫자(x좌표, y좌표)로 조정합니다. 그림 크기에 따라 위치는 달라집니다.

Postscript bounding box

SU(Seismic Un*x)에 있는 psimage로 Marmousi 속도모델을 그리면 다음과 같습니다.(그림 겉부분의 회색은 그림에 포함되어 있지 않은 부분으로, 경계를 표시하기 위해 넣었습니다.)

original eps file

original eps file

여기서 사용한 명령은 다음과 같습니다.

psimage par='../marm8m.txt' label1="Depth (km)" label2="Offset (km)" labelsize=8 height=1.0 width=2.4 legend=1 lstyle=vertright lwidth=0.1 lheight=1 units="m/s" < ../marm8m.drt > marm.eps

이 때, psimage는 그림 주위로 지나치게 넓은 공간을 만들어 줍니다. Bounding box 정보가 정확하지 않기 때문이죠. 이 상태로는 eps 파일을 다른 그림파일로 변환하여 paper에 넣거나 power point 발표자료에 넣기에 좋지 않습니다(물론 자르기 crop 기능을 이용할 수도 있기는 하죠).

이 공간을 없애기 위해서는 아래 명령을 이용합니다.

gs -sDEVICE=bbox -dNOPAUSE -dBATCH marm.eps

그럼 다음과 같은 결과를 보여줍니다.

GPL Ghostscript 8.63 (2008-08-01)
Copyright (C) 2008 Artifex Software, Inc.  All rights reserved.
This software comes with NO WARRANTY: see the file PUBLIC for details.
Loading NimbusSanL-Regu font from /usr/share/fonts/default/Type1/n019003l.pfb… 2656772 1085343 2641408 1357198 2 done.
Loading NimbusSanL-Bold font from /usr/share/fonts/default/Type1/n019004l.pfb… 2673436 1178370 2661504 1363393 2 done.
%%BoundingBox: 87 107 327 200
%%HiResBoundingBox: 87.695997 107.509005 326.645990 199.601994

위의 결과에서 마지막 두 줄에 나온 것이 흰 공간을 없앤 bounding box의 크기입니다. 둘 중 하나를 쓰시면 됩니다. 네 개의 숫자는 각각 왼쪽 아래 x좌표, 왼쪽 아래 y좌표, 오른쪽 위 x좌표, 오른쪽 위 y좌표를 의미합니다. Eps 파일을 텍스트 편집기로 열어서 %%BoundingBox 라고 써진 줄을 찾아 bounding box 크기를 위의 정보로 고쳐주면 아래와 같은 결과를 얻을 수 있습니다.

after fixing bounding box

after fixing bounding box

또는 SU에 있는 psbbox 라는 프로그램을 이용할 수도 있습니다.

psbbox llx=87 lly=107 urx=327 ury=200 < marm.eps >marmfx.eps

Gpl에 있는 fixbbox 프로그램은 위의 과정을 자동으로 실행하는 Python 프로그램으로,

fixbbox <input eps file> <output eps file>

과 같이 실행할 수 있습니다.

Fortran option parser

gpl에서 사용하기 위해 개발한 포트란 option parser를 소개합니다. Option parser는 command line option들을 분석해서 프로그램에서 사용할 수 있도록 변수로 저장해주는 역할을 합니다. 리눅스나 맥에서 명령을 실행할 때

$ gcc -c -o main.e main.c
$ tar -zxvf file.tgz

와 같은 유닉스 표준 형태의 옵션이 있고,

$ suximage n1=101 d1=0.5 perc=99 < file.su

와 같이 Seismic Un*xMadagascar 등에서 사용하는 형태의 옵션도 있는데, 지금 소개해드리는 option parser는 두 번째 형태의 옵션을 지원합니다. Option parser는 한 번 작성한 프로그램의 재사용을 위해 매우 유용한 기능입니다. 우선 테스트 프로그램을 통해 사용 예를 살펴본 후 자세한 사용법을 보겠습니다.

테스트 프로그램

다음의 테스트 프로그램은 다양한 종류의 변수를 command line option으로부터 읽어서 출력하는 프로그램입니다.
지원하는 자료형은 integer, single precision real, double precision real, logical, character입니다(complex는 아직까지 필요가 없어서 안 넣었습니다). 단일 변수 뿐 아니라 각각의 배열도 지원합니다. 필수 입력 옵션과 기본값을 가진 옵션을 구분하여, 필수 입력 옵션들 중 하나라도 command line option에 없을 경우 도움말을 출력하고 프로그램을 종료합니다.

        program test_optparse
        use gpl_optparse
        implicit none
        integer,parameter:: mxlen=100,mxstr=100
        integer:: i,io,ia(mxlen)
        integer:: j,ni,nf,nb,nd,ns
        real:: f,fo,fa(mxlen)
        logical:: b,bo,ba(mxlen)
        real(kind=8):: d,d_o,da(mxlen)
        character(len=mxstr) :: s,so,sa(mxlen)

        ! required parameters
        call from_par('i',i,'integer number')
        call from_par('f',f,'float: single precision real number')
        call from_par('d',d,'double precision')
        call from_par('b',b,'boolean/logical')
        call from_par('s',s,'string')

        ! optional parameters
        call from_par('io',io,1,'1','optional integer')
        call from_par('fo',fo,1.0,'1.0','optional float')
        call from_par('do',d_o,1.0d0,'1.0d0','optional double')
        call from_par('bo',bo,.true.,'T','optional boolean')
        call from_par('so',so,'','empty','optional string')

        ! array, required parameters
        call from_par('ia',ia,ni,'integer array')
        call from_par('fa',fa,nf,'float array')
        call from_par('da',da,nd,'double array')
        call from_par('ba',ba,nb,'boolean array')
        call from_par('sa',sa,ns,'string array')

        call help_par()
        !call report_par()
        
        print*,'i=',i
        print*,'f=',f
        print*,'d=',d
        print*,'b=',b
        print*,'s=',trim(s)
        print*,''
        print*,'ia=',ia(1:ni)
        print*,'fa=',fa(1:nf)
        print*,'ba=',ba(1:nb)
        print*,'da=',da(1:nd)
        print*,'sa=',(trim(sa(j))//'_',j=1,ns)
        print*,''
        print*,'io=',io
        print*,'fo=',fo
        print*,'do=',d_o
        print*,'bo=',bo
        print*,'so=',so
        end program

실행 결과

위의 프로그램을 컴파일하여 아래와 같이 실행하면 변수들이 제대로 입력된 것을 확인할 수 있습니다. 문자열의 경우 문자열을 둘러싼 따옴표는 제거하고 변수에 저장합니다. 배열은 쉼표를 기준으로 배열의 원소를 나눕니다. 기본값을 가지고 있는 변수의 경우 command line option에 있으면 주어진 값을 저장하고 없으면 기본값을 저장합니다.

$ ./test_optparse s=test.dat i=2 f=2.5 d=3.0 b=T ia=1,2,3,4 fa=1.0,2.0,3.0,4.0 da=1.,2.,3. ba=T,F,T sa=st,'ri',"ng" io=5 so="gpl"
 i=           2
 f=   2.50000000    
 d=   3.0000000000000000     
 b= T
 s=test.dat
 
 ia=           1           2           3           4
 fa=   1.00000000       2.00000000       3.00000000       4.00000000    
 ba= T F T
 da=   1.0000000000000000        2.0000000000000000        3.0000000000000000     
 sa=st_ri_ng_
 
 io=           5
 fo=   1.00000000    
 do=   1.0000000000000000     
 bo= T
 so=gpl

만약 기본값이 없는 필수 옵션 중 하나라도 command line option에 빠져있다면 다음과 같이 도움말을 표시하고 프로그램을 종료합니다.

$ ./test_optparse 
 Required parameters:
     [i] i=             : integer number
     [f] f=             : float: single precision real number
     [d] d=             : double precision
     [b] b=             : boolean/logical
     [s] s=             : string
     [I] ia=            : integer array
     [F] fa=            : float array
     [D] da=            : double array
     [B] ba=            : boolean array
     [S] sa=            : string array
 Optional parameters:
     [i] io=1           : optional integer
     [f] fo=1.0         : optional float
     [d] do=1.0d0       : optional double
     [b] bo=T           : optional boolean
     [s] so=empty       : optional string

위의 도움말에서 i=integer, f=real, d=real(kind=8), b=logical, s=character(len=?)을 의미하고, 각각의 대문자는 배열을 의미합니다.

사용법

그럼 실제 모듈의 사용법을 알아보겠습니다. 모듈은 use gpl_optparse 또는 간편하게 use gpl로 불러올 수 있습니다. 가장 중요한 from_par 서브루틴은 다음과 같이 세 가지 방식으로 사용할 수 있습니다.

call from_par('parname',variable,'help message') !! 필수 옵션
call from_par('parname',variable_arr,len_arr,'help message') !! 필수 배열
call from_par('parname',variable,default,'default message','help message') !! 기본값이 있는 변수

Command line option은 parname=value 형태로 입력 받게 됩니다. 위의 서브루틴들에서 'parname'은 이 때 사용되는 이름이고, valuevariable에 저장됩니다. ’help message’는 도움말 출력시 보여주는 변수 설명입니다.

배열의 경우 parname=value1,value2,value3과 같은 형태로 입력받고, 입력받은 값은 variable_arr에 저장됩니다. 이 때 입력받은 원소는 len_arr 개입니다.

기본값이 있는 변수의 경우 default는 기본값이고, 'default message'는 도움말 출력시 기본값을 보여주기 위한 문자열입니다.

함수 오버로딩을 사용하였기 때문에 정수, 실수 등의 자료형에 상관없이 위의 서브루틴들을 이용할 수 있습니다.

기타 사용할 수 있는 서브루틴들과 함수는 다음과 같습니다.

call help_par()
call help_header('msg before the parameter help msg')
call help_footer('msg after the parameter help msg')
call force_help()

call report_par()

if(given_par('parname')) call do_something()
call from_parfile('parfile.txt')

에서 help_par는 도움말 출력을 위한 서브루틴입니다. 도움말과 관련된 서브루틴들로는 도움말을 보강하기 위한 call help_header('msg'), call help_footer('msg'), 필수 옵션이 주어졌는지와 무관하게 도움말을 출력하기 위한 call force_help()가 있습니다. report_par는 입력받은 변수들을 출력해서 확인하기 위한 서브루틴입니다. 이 외에도 옵션이 주어졌는지 확인하기 위한 logical function given_par('parname') 함수, 옵션들을 저장해놓은 텍스트파일로부터 옵션들을 읽어들이기 위한 call from_parfile('parfile.txt')와 같은 명령이 있습니다.

참고로, command line option에 par=parfile.txt와 같은 옵션이 있으면 parfile.txt에서 먼저 변수를 읽은 후 command line option을 읽습니다. parfile.txt에 주어진 변수가 command line option에 다시 나오면 command line에 주어진 값을 사용합니다.

gpl_optparse 모듈의 소스코드는 gplGitHub에서 받으실 수 있습니다.

Polymorphic Fortran & C

포트란 함수 오버로딩과 중복

앞서 포트란 함수 오버로딩에 관한 글을 올렸습니다. 포트란 모듈과 인터페이스를 이용하면 서로 다른 자료형을 인자로 받는 함수나 서브루틴이라도 같은 이름으로 사용할 수가 있었습니다. 그런데, 동적 자료형을 지원하는 언어와 달리, 포트란에서는 자료형별로 서브루틴들을 따로 만든 후에 같은 이름으로 호출하였습니다. 동적 자료형 언어에서는 함수 자체를 한 번만 작성하면 되지요. 여기에서 포트란 코드 작성에 중복이 발생하게 됩니다. 이러한 중복을 제거하기 위해 파이썬으로 만든 스크립트가 polyfc (Polymorphic Fortran & C)입니다. 예전에 Forpedo에서 아이디어를 얻었는데, Forpedo를 사용하려니 좀 복잡해서 Python 연습도 할 겸 gpl용으로 만들었습니다.

예제

이해하기 쉽게 예제를 살펴보겠습니다. 배열의 내용을 출력하는 서브루틴을 작성하려고 합니다. 배열은 integer, real, real(kind=8), complex, complex(kind=8) 다섯 종류의 자료형을 지원하려고 합니다. 그럼 서브루틴을 자료형에 따라 총 5개를 작성해야 하는데, 선언부만 다르고 나머지는 동일하거나 거의 비슷하게 됩니다. 그래서 일종의 템플릿 서브루틴을 만들고, 필요한 부분만 바꿔가며 서브루틴을 복제하려고 합니다. 그러면 거의 비슷한 코드를 중복해서 작성하는 수고를 덜 수 있겠죠. 아래는 polyfc에 입력으로 들어가는 템플릿 파일입니다.

module polyfc_example

!@interface print_array
contains
!@template print_array ifdcz
    subroutine print_array_<name>(arr)
    <type>,intent(in):: arr(:)
    integer i
    do i=1,size(arr)
        print*, i, arr(i)
    enddo
    end subroutine
!@end
end module

포트란 주석을 이용하여 인터페이스가 들어갈 부분과 템플릿 부분을 표시하고, 자료형에 따라 바뀌어야 하는 부분은 <name>, <type>으로 표시하였습니다. 이 외에 필요할 경우 <kind><esize>도 지원합니다. 바뀌며 들어가는 부분은 다음 표를 보면 알 수 있습니다.

Name Fortran type C type Fotran kind esize
i integer int (kind=4) 4
f real float (kind=4) 4
d real(kind=8) double (kind=8) 8
c complex float complex (kind=4) 8
z complex(kind=8) double complex (kind=8) 16
b logical (not supported) (kind=4) 4
s character(len=*) char (len=*) 1

그래서 polyfc input.f90 > output.f90과 같이 실행했을 때 얻게 되는 파일은 다음과 같습니다.

! This file was generated from pfc.example.f90 by Polyfc at Mon Jul 28 22:12:32 2014.
! Do not edit this file directly.

module polyfc_example

    interface print_array
        module procedure print_array_i
        module procedure print_array_f
        module procedure print_array_d
        module procedure print_array_c
        module procedure print_array_z
    end interface print_array

contains

    subroutine print_array_i(arr)
    integer,intent(in):: arr(:)
    integer i
    do i=1,size(arr)
        print*, i, arr(i)
    enddo
    end subroutine
 
    subroutine print_array_f(arr)
    real,intent(in):: arr(:)
    integer i
    do i=1,size(arr)
        print*, i, arr(i)
    enddo
    end subroutine
 
    subroutine print_array_d(arr)
    real(kind=8),intent(in):: arr(:)
    integer i
    do i=1,size(arr)
        print*, i, arr(i)
    enddo
    end subroutine
 
    subroutine print_array_c(arr)
    complex,intent(in):: arr(:)
    integer i
    do i=1,size(arr)
        print*, i, arr(i)
    enddo
    end subroutine
 
    subroutine print_array_z(arr)
    complex(kind=8),intent(in):: arr(:)
    integer i
    do i=1,size(arr)
        print*, i, arr(i)
    enddo
    end subroutine
 
end module

자동으로 작성하면 손으로 복사했을 경우 생길 수 있는 오류도 피할 수 있겠죠. 단, 코드에 버그가 있을 때에는 output 파일에서 버그가 있는 곳을 찾아 input 파일의 해당 위치를 고쳐줘야 합니다.

참고로, polyfc 이름이 의미하듯, C도 지원합니다. 단, C에서는 <name><type>만 지원합니다. 물론 인터페이스도 지원하지 않습니다. C에서는 아래와 같은 방식으로 사용 가능합니다.

//@template ifdczbs
void abc_<name>(<type> def)
{
    ...
}
//@end

GNU Quick Plot (gnuqp)

Gnuplot은 리눅스에서 텍스트파일에 저장된 값을 빠르게 그림으로 그려주는 프로그램입니다. 다양한 기능을 가지고 있지만, 제 경우에는 주로 수치해석 후 결과 확인용으로 씁니다.
gnuplot으로 그림을 그릴 때에는 command line 상에서 gnuplot이라고 치고 들어가서 gnuplot 명령어들을 이용하여 그림을 그리고 q를 입력하여 빠져나옵니다.
그런데 간단히 결과를 확인해보기 위해서 gnuplot에 들어가서

p 'file1' w l,'file1' u 1:3 w l,'file1' u 1:4 w p

또는

set grid
set xrange[:10]
set log y
p 'file1' w l,'file2' w l,'file3' w l

과 같이 매번 치려니 귀찮다는 생각이 들었습니다. 그래서 gnuqp (GNU Quick Plot)를 만들었습니다. 이 script를 사용하면 command line 상에서 바로 gnuplot 명령어를 사용하여 그림을 그릴 수 있습니다. 사용 방법은 아래와 같습니다.

Usage :
gnuqp [options] filename1 [u 1:2] [w l], filename2 [u 1:2] [w l], filename3 ...

실행파일 이름, 몇 가지 setting 관련 옵션들, 이후에는 gnuplot의 plot 명령어를 입력합니다.

Required parameters :
filename1
Empty filename[2,3,...] will be replaced by the filename1

두 번째 위치부터는 파일명을 생략하면 첫 번째 파일명으로 대체합니다. 하나의 파일에서 여러 column들을 그릴 때 편리합니다.

Optional parameters :
u 1:2   : columns you want to plot
w [lp..]: line style- line, point, dot or impulse ..etc (default: w l)

plot 명령어의 옵션들 중에는 using (columns)과 with (line style)만 지원합니다. 그 외의 명령은 제가 잘 안 써서요^^.
위의 옵션을 주지 않았을 때 기본적으로 with line 옵션으로 그립니다.

-p      : do not run gnuplot. just print the gnuplot command
-c      : no comma seperation - the arguments are filenames seperated with a blank- use with glob pattern
-l       : set logscale y
-g      : set grid
-x[:10] : set xrange [:10]
-y[1:5] : set yrange [1:5]

위의 옵션들은 gnuplot의 setting을 간편하게 하기 위해 만들었습니다.

-p 옵션을 붙이면 gnuplot의 명령어만 출력하고 그림은 안 그립니다.

-c 옵션을 붙이면 파일들을 기본 옵션(with line)으로 그립니다. 이 때 파일명들 사이의 “,”를 생략하고 파일명만 씁니다. command line상에서 glob pattern을 이용하여 여러 그림을 그릴 수 있도록 하기 위한 옵션입니다. 예를 들면, 다음과 같은 경우죠.

./gnuqp.py -p -c file.00*
-> p 'file.0010' w l,'file.0020' w l,'file.0030' w l,'file.0040' w l,'file.0050' w l

나머지 gnuplot setting들은 위의 설명으로 충분할 것이라 생각합니다.
앞에 예를 들었던 명령어들을 gnuqp를 이용하여 실행한다면 다음과 같습니다.

(gnuplot)
p 'file1' w l,'file1' u 1:3 w l,'file1' u 1:4 w p
(q)

gnuqp file1, u 1:3, u 1:4 wp,

(gnuplot)
set grid
set xrange[:10]
set log y
p 'file1' w l,'file2' w l,'file3' w l
(q)

gnuqp file1, file2, file3 -g -l -x[:10]와 같이 실행할 수 있습니다. gnuqp는 gpl에 포함되어 있습니다.