スカラーインデックス(n_s)について

計算したいn_sの数は「nn」。これは、driversub.f 中で、nn=num_slopeとして決められる。num_slopeはsubroutine run_cmbfastの引数。nnはcommon /initialps/　に属していて、これでcmbflat等に渡される。

cmbflatが計算するclts等は、２次元配列で、clts(l0max,nnmax)のように定義されている。これらはいずれもcmbfast.incで

c l0max is used to dimension arrays. It must be bigger than the
c maximum l to be calculated by at least 300.
c
parameter (l0max=5300)


c nnmax sets the maximum number of spectral indeces that
c can be calculated.

parameter(nnmax=51)

のように決められている。１回Boltzmann方程式を解くごとに、最大nnmax個のスカラーインデックスに対応するC_lが計算することができるようになっている。しかし、そのスカラーインデックスは前もって指定しておかなくてはならず、brentを使うときのようにchi2を計算した結果を見て次に計算するスカラーインデックスの値を決めるというアルゴリズムにはそのまま使えない。

cmbflat中には、以下のような長いループがあって、kで指定される各モードについてBoltzmann方程式が積分されたり、line of sight integralがされたりしている（中略と書いてあるところで）。スカラーとテンサーそれぞれあるが、いずれも積分し終わったあとに、C_lを計算する準備として、初期パワースペクトル（powersflat, powertflatを呼んで用意する）をかけたcl, cpl, ccl, ckkl, ctkl （スカラーの場合）を計算している。

c Begin k-loop
do 200 k=kkmin+myid,kkmax,numprocs
akt=ak1(k)


　　　　　　＜中略＞

c Adding to calculate the integral over k.
c Scalar case
c
do in=1,nn
do j=1,l0
call powersflat(akt,in,apowers)
apowers=apowers/akt
ckj=apowers*dl(j)*dl(j)*dak1(k)
cl(j,in)=cl(j,in)+ckj
cpkj=apowers*dpl(j)*dpl(j)*dak1(k)
cpl(j,in)=cpl(j,in)+cpkj
cckj=apowers*dl(j)*dpl(j)*dak1(k)
ccl(j,in)=ccl(j,in)+cckj
ckkkj=apowers*dkl(j)*dkl(j)*dak1(k)
ckkl(j,in)=ckkl(j,in)+ckkkj
ctkkj=apowers*dl(j)*dkl(j)*dak1(k)
ctkl(j,in)=ctkl(j,in)+ctkkj
end do
end do



　　　　　　　　＜中略＞

c Tensor case.
c
c Adding to calculate the integral over k.
c Tensor case

do in=1,nn
do j=1,l0
call powertflat(akt,in,apowert)
apowert=apowert/akt


ctkj=apowert*dtl(j)*dtl(j)*dak1(k)
ctl(j,in)=ctl(j,in)+ctkj

ctekj=apowert*dtel(j)*dtel(j)*dak1(k)
ctel(j,in)=ctel(j,in)+ctekj

ctbkj=apowert*dtbl(j)*dtbl(j)*dak1(k)
ctbl(j,in)=ctbl(j,in)+ctbkj

ctckj=apowert*dtl(j)*dtel(j)*dak1(k)
ctcl(j,in)=ctcl(j,in)+ctckj

end do
end do
c
end if

200 continue


スカラーの場合で言うと、ckj, cpkj, cckj,ckkkj,ctkkj (kとjの足を持っているようなもん）を２次元配列に保存して、cl等をk-loopの外で計算するようにすれば、 n_s だけ変えるときのC_lを、積分の結果をそのままつかってapowersの計算だけし直すことで計算できる。