なかのにっき

sourceforge の Octave Repository にある。debian だと octave-forge を入れれば使える。

measured date    :  5: 9 PM 11/21/ 2
data title       : %T   T ALL RANGE  
operating mode   :     spectrum scan
system responce  :            medium
bandpass(nm)     :           2.0(nm)
UV-VIS change WL :         340.0(nm)
baseline         :            system
Lamp chg control :              auto
Higher bound WL  :         900.0(nm)
Lower  bound WL  :         185.0(nm)
scan speed       :         120.0(nm/s)

wavelength(nm), Intensity ratio
       900.000,       0.8973543
       899.600,       0.9003937
        :
       185.200,       0.0121768
[EOF]

というようなデータファイル (DOS 形式) を

#!/usr/bin/perl

while(<>){/^wavelength/ && last; }

while(<>){
	chomp;
	s/\r//;
	(my $wltmp, $trtmp) = split /,/;
	unshift @wl, $wltmp;
	unshift @tr, $trtmp;
}

printf "%d\n", scalar(@wl);

while(@wl){
	printf "%20f%20f\n", pop(@wl), pop(@tr);
}

という perl script (chop.pl) を通し、さらに

#!/usr/bin/octave -qf

LOADPATH=":/usr/share/octave//";

n = scanf("%d", "C");
for i = 1:n
	[wl(i) tr(i)] = scanf("%f%f", "C");
endfor

sr = sgolayfilt(tr, 5, 7);

nout = 0;
for i = 1:n
	# [400nm, 800nm], 2nm step
	if (ceil( wl(i) / 2) * 2 == wl(i) && wl(i) <= 800 && wl(i) >= 400)
		nout += 1;
	endif
endfor
printf("%5d\n", nout);

for i = 1:n
	if (ceil( wl(i) / 2) * 2 == wl(i) && wl(i) <= 800 && wl(i) >= 400)
		printf("%20f%20f\n", wl(i), sr(i));
	endif
endfor

という octave script (prefit.m) を通して出来上がり。 出力行カウントに 2 回ループをまわしてるのがマヌケだが、まあいいや。 こいつを fortran プログラムに食わせれば良いだろう。

なお octave-forge のスクリプトを使うときは LOADPATH 指定が必要。

• 最初の : (DEFAULT_LOADPATH に置き換わる)
• 最後の // (kpathsea 的に下層も再帰的に探す)

の 2 点がミソ ( このへん 参照)。