The following text lists the FORTRAN code, written by Allen W. Shafter (1997), following Keith Horne (1991), that utilizes Fourier-filtered back-projection techniques to create a Doppler tomogram from spectroscopic data. This author
revised the code to include a Gaussian cut-off term to the high-pass ramp filter, and FITS header information.
program backproj
dimension data(751,54), d(751,54), bp(65,65), ph(54)
dimension dd(751)
character*40 fname,inlist,ofn,pfname
character*1 dummy
write (*,'(''ENTER FILENAME OF INPUT LIST: ''$)')
read (*,550) inlist
open (4,file=inlist,status='old')
550 format (a40)
write (*,'(''ENTER INITIAL WAVELENGTH (A): ''$)')
read (5,*) w0
write (*,'(''ENTER DISPERSION (A/ch): ''$)')
read (5,*) dis
write (*,'(''ENTER CUT-OFF FREQUENCY: ''$)')
read (5,*) wc
write (*,'(''ENTER EMISSION LINE CENTRAL WAVELENGTH (A): ''$)')
read (5,*) xlc
write (*,'(''ENTER DIMENSION (even number): ''$)')
read (5,*) numm
write (*,'(''ENTER VELOCITY EXTENT OF MAP (km/s): ''$)')
read (5,*) extent
write (*,'(''ENTER PHASE FILENAME: ''$)')
read (*,550) pfname
write (*,'(''ENTER OUTPUT FILENAME: ''$)')
read (*,550) ofn
num=numm+1
numo2=numm/2
xmult=extent/float(numo2)
xnyq=1./(2.*dis)
cf=100.
open (7,file=pfname)
open (2,file=ofn)
pi=3.1415926
c=2.9979e5
j=0
99 continue
j=j+1
read (4,550,end=11) fname
read (7,*) ph(j)
open (3,file=fname,status='old')
7 read (3,105) dummy
105 format(a1)
if (dummy.ne.' ') go to 7
i=0
1 i=i+1
read (3,*,end=99) data(i,j)
go to 1
11 continue
jnum=j-1
close (3)
close (4)
ii=int((xlc-w0)/dis)-32
if (ii.lt.0) then
print *, 'ERROR: spectral subset too big. try inum=128.'
go to 999
endif
inum=64
nn=inum/2
isign=1
jsign=-1
do j=1,jnum
do i=1,inum
dd(i)=data(i+ii,j)
end do
call realft(dd,nn,isign)
do i=1,inum
xi=1./(float(i)*dis)
dd(i)=dd(i)*(xi/xnyq)*exp((-1.*((xi/wc)**2))/2.0)
end do
call realft(dd,nn,jsign)
do i=1,inum
d(i,j)=dd(i)/float(nn)
end do
end do
C HEADER INFORMATION
C **************************************************************
write (2,22) 'SIMPLE ','=','T'
22 format (A8,A1,T11,A20)
23 format (A8,A1,T11,I20)
24 format (A8,A1,T11,A1,A18,T30,A1)
25 format (A8,A1,T11,G20.5)
26 format (A3)
write (2,23) 'BITPIX ','=', 8
write (2,23) 'NAXIS ','=', 2
write (2,23) 'NAXIS1 ','=', num
write (2,23) 'NAXIS2 ','=', num
write (2,24) 'OBJECT ','=', '''', object, ''''
write (2,24) 'TELESCOP','=', '''', 'MLO 1m' ,''''
write (2,24) 'OBSERVER','=', '''', 'Etzel & Dewey' ,''''
C HEADER FOR COORDINATES OF IMAGE
write (2,25) 'CRVAL1 ','=', (-1.*extent)
write (2,23) 'CRPIX1 ','=', 1
write (2,25) 'CDELT1 ','=', xmult
write (2,24) 'CRTYPE1 ','=', '''','VELO',''''
write (2,25) 'CRVAL2 ','=', (-1.*extent)
write (2,23) 'CRPIX2 ','=', 1
write (2,25) 'CDELT2 ','=', xmult
write (2,24) 'CRTYPE2 ','=', '''','VELO',''''
write (2,24) 'FILENAME','=', '''', 'DOPPLER TOMOGRAM', ''''
write (2,26) 'END'
C **************************************************************
do l=1,num
do k=1,num
vx=xmult*float(k-1-numo2)
vy=xmult*float(l-1-numo2)
sum=0.
do j=1,jnum
phi=ph(j)
v=vy*sin(2.*pi*phi)-vx*cos(2.*pi*phi)
ich=int((xlc/dis)*v/c)+32
sum=sum+d(ich,j)/((d(1,j)+d(inum,j))/2.)
end do
bp(k,l)=sum
end do
end do
write (2,100) ((bp(k,l), k=1,num), l=1,num)
100 format (6e12.5)
close (2)
999 stop
end
subroutine realft(data,n,isign)
real*8 wr,wi,wpr,wpi,wtemp,theta
dimension data(*)
theta=6.28318530717959d0/2.0d0/dble(n)
c1=0.5
if (isign.eq.1) then
c2=-0.5
call four1(data,n,+1)
else
c2=0.5
theta=-theta
endif
wpr=-2.0d0*dsin(0.5d0*theta)**2
wpi=dsin(theta)
wr=1.0d0+wpr
wi=wpi
n2p3=2*n+3
do 11 i=2,n/2+1
i1=2*i-1
i2=i1+1
i3=n2p3-i2
i4=i3+1
wrs=sngl(wr)
wis=sngl(wi)
h1r=c1*(data(i1)+data(i3))
h1i=c1*(data(i2)-data(i4))
h2r=-c2*(data(i2)+data(i4))
h2i=c2*(data(i1)-data(i3))
data(i1)=h1r+wrs*h2r-wis*h2i
data(i2)=h1i+wrs*h2i+wis*h2r
data(i3)=h1r-wrs*h2r+wis*h2i
data(i4)=-h1i+wrs*h2i+wis*h2r
wtemp=wr
wr=wr*wpr-wi*wpi+wr
wi=wi*wpr+wtemp*wpi+wi
11 continue
if (isign.eq.1) then
h1r=data(1)
data(1)=h1r+data(2)
data(2)=h1r-data(2)
else
h1r=data(1)
data(1)=c1*(h1r+data(2))
data(2)=c1*(h1r-data(2))
call four1(data,n,-1)
endif
return
end
subroutine four1(data,nn,isign)
real*8 wr,wi,wpr,wpi,wtemp,theta
dimension data(*)
n=2*nn
j=1
do 11 i=1,n,2
if(j.gt.i)then
tempr=data(j)
tempi=data(j+1)
data(j)=data(i)
data(j+1)=data(i+1)
data(i)=tempr
data(i+1)=tempi
endif
m=n/2
1 if ((m.ge.2).and.(j.gt.m)) then
j=j-m
m=m/2
go to 1
endif
j=j+m
11 continue
mmax=2
2 if (n.gt.mmax) then
istep=2*mmax
theta=6.28318530717959d0/(isign*mmax)
wpr=-2.d0*dsin(0.5d0*theta)**2
wpi=dsin(theta)
wr=1.d0
wi=0.d0
do 13 m=1,mmax,2
do 12 i=m,n,istep
j=i+mmax
tempr=sngl(wr)*data(j)-sngl(wi)*data(j+1)
tempi=sngl(wr)*data(j+1)+sngl(wi)*data(j)
data(j)=data(i)-tempr
data(j+1)=data(i+1)-tempi
data(i)=data(i)+tempr
data(i+1)=data(i+1)+tempi
12 continue
wtemp=wr
wr=wr*wpr-wi*wpi+wr
wi=wi*wpr+wtemp*wpi+wi
13 continue
mmax=istep
go to 2
endif
return
end