17 December 2004
atomtab.f
In case you want this without the hassle of getting the whole VPFIT setup,
the source for the Fortran77 program for converting a file of Jason Prochaska's
wavelengths and oscillator strengths is given below. It does what it says
as a series of write statements at the beginning of the listing.
*
* program to compare atom.dat file with Prochaska
data file, and replace
* any wavelengths or oscillator strengths with his
values
*
character*132 inchstr,jnchstr
double precision ddiff,dtol,difx,difp
double precision dv(24),djv(8)
integer iv(24),ijv(8)
character*32 cv(24),cjv(8)
* Prochaska values
character*8 cjion(1000)
double precision djwv(1000),djfik(1000)
integer lnj(1000)
* atom.dat values
double precision drv(24,1000)
double precision rwv(1000),rfk(1000)
character*132 iachstr(1000)
integer lchstr(1000)
character*32 crv(24,1000)
integer lrv(24,1000),lastrv(1000)
logical lcont
*
write(6,*) '__________________________________________________
'
write(6,*) ' '
write(6,*) 'This program takes a VPFIT-based
atom.dat file,'
write(6,*) 'and a Prochaska file of wavelengths
and oscillator'
write(6,*) 'strengths, and compares the two.
Where differences'
write(6,*) 'are found in the nearest wavelengths
for a given ion'
write(6,*) 'the wavelength is set to the
Prochaska value in the'
write(6,*) 'output file, and similarly if
the oscillator'
write(6,*) 'strengths for an ion at similar
wavelengths differ'
write(6,*) 'then the Prochaska table value
is written to the'
write(6,*) 'output file. Upper level transition
rates and atomic'
write(6,*) 'masses are left as thay were,
and the resulting'
write(6,*) 'output file is readable by VPFIT'
write(6,*) ' '
write(6,*) 'Note that only lines which appear
in the input'
write(6,*) 'atom.dat will appear in the output
file. Lines'
write(6,*) 'from the Prochaska file which
have no counterpart'
write(6,*) 'in atom.dat are omitted from
the output.'
write(6,*) ' '
write(6,*) 'Values modified from the original
atom.dat are'
write(6,*) 'labelled Prochaska at the end
of the line. These'
write(6,*) 'should be checked to make sure
they are right. In'
write(6,*) 'particular please check that
lines are not'
write(6,*) 'duplicated in the output file
replacing those with'
write(6,*) 'close wavelengths in the original
atom.dat '
write(6,*) ' '
write(6,*) '__________________________________________________
'
write(6,*) ' '
dtol=2.0d-1
do k=1,1000
rwv(k)=0.0d0
rfk(k)=0.0d0
end do
*
write(6,*) 'VPFIT atomic data file?'
read(5,'(a)') inchstr
ln=lastchpos(inchstr)
open(unit=1,name=inchstr(1:ln),status='old',err=999)
* read everything in
lcont=.true.
nat=1
do while (lcont)
read(1,'(a)',end=991) iachstr(nat)
lchstr(nat)=lastchpos(iachstr(nat))
call dsepvar(iachstr(nat),24,dv,iv,cv,lastrv(nat))
if(cv(1)(2:2).eq.' ') then
cv(1)=cv(1)(1:1)//cv(2)
do j=2,23
dv(j)=dv(j+1)
cv(j)=cv(j+1)
iv(j)=iv(j+1)
end do
lastrv(nat)=lastrv(nat)-1
end if
do j=1,lastrv(nat)
crv(j,nat)=cv(j)
drv(j,nat)=dv(j)
lrv(j,nat)=lastchpos(cv(j))
end do
if(cv(1)(1:2).eq.'EN') then
lcont=.false.
end if
goto 992
991 lcont=.false.
992 nat=nat+1
end do
nat=nat-1
close(unit=1)
*
write(6,*) 'Prochaska atomic data file?'
read(5,'(a)') inchstr
ln=lastchpos(inchstr)
open(unit=2,name=inchstr(1:ln),status='old',err=999)
* read it all in
read(2,'(a)') jnchstr
call dsepvar(jnchstr,8,djv,ijv,cjv,njv)
nrec=ijv(1)
if(nrec.gt.1000) then
write(6,*) 'Using only first
1000 lines from Prochaska list'
nrec=1000
end if
do j=1,nrec
read(2,'(a)',end=997) jnchstr
call dsepvar(jnchstr,8,djv,ijv,cjv,njv)
cjion(j)=cjv(2)
lnj(j)=lastchpos(cjion(j))
djwv(j)=djv(1)
djfik(j)=djv(4)
end do
997 close(unit=2)
*
write(6,*) 'Destination atomic data file?'
read(5,'(a)') inchstr
ln=lastchpos(inchstr)
open(unit=3,name=inchstr(1:ln),status='unknown',err=999)
*
* cycle through atom.dat material (k=atom.dat, j=Jason
P.)
do kat=1,nat
* if special character, do nothing
if(crv(1,kat)(1:2).ne.'??'.and.crv(1,kat)(1:2).ne.'<>'.and.
: crv(1,kat)(1:2).ne.'__'.and.crv(1,kat)(1:2).ne.'>>'.and.
: crv(1,kat)(1:2).ne.'<<')
then
ln=lrv(1,kat)
jv=1
jvbest=0
ddiff=1.0d30
do while(jv.le.nrec.and.lnj(jv).gt.0)
if(ln.eq.lnj(jv))
then
if(crv(1,kat)(1:ln).eq.cjion(jv)(1:ln)) then
if(abs(drv(2,kat)-djwv(jv)).lt.ddiff) then
jvbest=jv
ddiff=abs(drv(2,kat)-djwv(jv))
end if
end if
end if
jv=jv+1
end do
* jvbest is closest relevant
wavelength. Is it close enough?
jrep=0
if(ddiff.lt.dtol)
then
* possibly replace
value, unless they are identical
if(ddiff.gt.1.0d-6)
then
rwv(kat)=djwv(jvbest)
end if
if(abs(drv(3,kat)-djfik(jvbest)).gt.1.0d-4*drv(3,kat))
then
*
possibly replace oscillator strength
rfk(kat)=djfik(jvbest)
rwv(kat)=djwv(jvbest)
end if
end if
end if
end do
*
* now have a set of possible replacements, but should
only replace closest
do k=1,nat
* if(crv(1,k)(1:4).eq.'CII*')
then
* write(6,*) crv(1,k)(1:6),drv(2,k),drv(3,k),drv(4,k),rwv(k),
* : rfk(k)
* end if
krep=0
if(rwv(k).ne.0.0d0.or.rfk(k).ne.0.0d0)
then
* might want to replace
something
krep=k
difx=0.0d0
if(rwv(k).gt.0.0d0)
then
difx=abs(rwv(k)-drv(2,k))*1000.0d0
end if
if(rfk(k).gt.0.0d0)
then
difx=difx
+ abs(1.0d0-rfk(k)/drv(3,k))*10.0d0
end if
do k2=1,nat
if(crv(1,k).eq.crv(1,k2).and.rwv(k).eq.rwv(k2).and.
:
rfk(k).eq.rfk(k2)) then
*
possible duplicate
difp=0.0d0
if(rwv(k).gt.0.0d0) then
difp=abs(rwv(k)-drv(2,k2))*1000.0d0
end if
if(rfk(k).gt.0.0d0) then
difp=difp + abs(1.0d0-rfk(k)/drv(3,k2))*10.0d0
end if
if(difp.lt.difx) then
krep=k2
difx=difp
end if
end if
end do
* now have closest, so replace
that one and eliminate the others
if(rwv(krep).gt.0.0d0)
then
drv(2,krep)=rwv(krep)
end if
if(rfk(krep).gt.0.0d0)
then
drv(3,krep)=rfk(krep)
end if
do k2=1,nat
if(k2.ne.krep.and.crv(1,k2).eq.crv(1,krep).and.
: rwv(k2).eq.rwv(krep).and.rfk(k2).eq.rfk(krep))
then
rwv(k2)=0.0d0
rfk(k2)=0.0d0
end if
end do
* it may visit krep again,
but no longer matters
end if
end do
*
* appropriate replacements now done, so write out
the results
do k=1,nat
if(rwv(k).gt.0.0d0.or.rfk(k).gt.0.0d0)
then
if(lastrv(k).ge.5)
then
if(crv(5,k)(1:1).ne.'
') then
jnchstr=' '//crv(5,k)(1:lrv(5,k))
lxj=lastchpos(jnchstr)
if(lastrv(k).ge.6) then
do j=6,lastrv(k)
lx=lastchpos(cv(j))
jnchstr=jnchstr(1:lxj)//' '//crv(j,k)(1:lrv(j,k))
lxj=lastchpos(jnchstr)
end do
jnchstr=jnchstr(1:lxj)//' Prochaska'
lxj=lxj+10
end if
else
jnchstr=' Prochaska'
lxj=1
end if
end if
if(crv(1,k)(2:2).eq.'I'.or.crv(1,k)(2:2).eq.'V')
then
crv(1,k)=crv(1,k)(1:1)//'
'//crv(1,k)(2:6)
end if
write(3,'(a,f10.4,1pe12.3,e12.3,a)')
crv(1,k)(1:7),drv(2,k),
: drv(3,k),drv(4,k),jnchstr(1:lxj)
else
* just copy the input line
write(3,'(a)') iachstr(k)(1:lchstr(k))
end if
end do
998 write(6,*) 'Finished OK'
996 stop
999 write(6,*) 'disaster'
goto 996
end
SUBROUTINE DSEPVAR(CHST,NVAR,RV,IV,CV,NV)
c
c Separate a string of variables, separated by blanks,
commas or &,
c into variable arrays. ' or " are used to delineate character
c variables which contain any special characters (space,
comma,
c ' or "), so ' and " cannot appear in the same string.
Successive
c commas are treated as zero or blank, as appropriate.
Variables
c which do not translate to real or integer are set to
zero.
c Note - spaces after a character before commas are treated
c as another separator, so
c 1 ,,2 becomes
1 0 0 2
c and
c 1,, 2
" 1 0 2
c
c input: chst character
variable containing data
c nvar number of
variables expected
c
c output: cv character
variables
c rv double precision
variables
c iv integer variables
c nv actual number
of variables found
c
c
character*132 cx
character*(*) cv(*),chst
character*1 chdg,ctab
character*3 cdum
dimension iv(*)
double precision rv(*)
double precision xxx
* tab character
ctab=' '
c single quote fix
cdum='x''x'
chdg=cdum(2:2)
c
c
c string lengths
lstr=len(chst)
if(lstr.le.0) then
c no data at all
cv(1)=' '
rv(1)=0.0d0
iv(1)=0
nv=0
return
end if
c lv=len(cv(1))
c
kk=1
c
c find first non-blank character
k1=1
1 continue
do while (chst(k1:k1).eq.' ')
if(k1.ge.lstr) goto 3
k1=k1+1
end do
c
c check that it is not a special character
if(chst(k1:k1).eq.chdg.or.chst(k1:k1).eq.'"') then
c
c character string in quotes - search for matching
quote
k2=k1+1
do while (chst(k2:k2).ne.chst(k1:k1).and.
1 k2.le.lstr)
k2=k2+1
end do
if(k2.le.k1+1) then
cv(kk)=' '
else
cv(kk)=chst(k1+1:k2-1)
end if
rv(kk)=0.0d0
iv(kk)=0
k1=k2+1
kk=kk+1
c if separator added, go one more step
if(chst(k1:k1).eq.','.or.chst(k1:k1).eq.'&'.or.
: chst(k1:k1).eq.ctab)
k1=k1+1
else
if(chst(k1:k1).eq.','.or.chst(k1:k1).eq.'&'.or.
: chst(k1:k1).eq.ctab)
then
c comma separator, straight after previous
separator
cv(kk)=' '
rv(kk)=0.0d0
iv(kk)=0
kk=kk+1
k1=k1+1
else
c
c some other character - search for
the next
c space or comma or '&'
k2=k1+1
do while (chst(k2:k2).ne.' '.and.chst(k2:k2).ne.','
1 .and.chst(k2:k2).ne.'&'.and.chst(k2:k2).ne.ctab
: .and.k2.le.lstr)
k2=k2+1
end do
k2m=k2-1
cv(kk)=chst(k1:k2m)
kvar=k2m-k1
c check if a reasonable integer
if(kvar.gt.9) goto 91
cx=' '
if(kvar.lt.20) then
cx(20-kvar:20)=chst(k1:k2m)
end if
c read as an integer
read(cx,'(i20)',err=91) iv(kk)
rv(kk)=iv(kk)
goto 92
c
c integer failed - try as a real
91 cx=' '
if(kvar.lt.30) then
cx(30-kvar:30)=chst(k1:k2m)
end if
read(cx,'(f30.0)',err=93) rv(kk)
xxx=rv(kk)
if(xxx.lt.0.0) xxx=-xxx
if(xxx.le.2.0d9) then
iv(kk)=int(rv(kk))
else
iv(kk)=0
end if
goto 92
c
c character string only thing that
works
93 rv(kk)=0.0d0
iv(kk)=0
92 k1=k2+1
kk=kk+1
end if
end if
if(k1.le.lstr.and.kk.le.nvar) goto 1
c
c end of data
3 nv=kk-1
c
c zero or blank off nv+1 to nvar
if(nv.lt.nvar) then
do i=kk,nvar
cv(i)=' '
rv(i)=0.0d0
iv(i)=0
end do
end if
return
end
integer function lastchpos(chstr)
* finds the position of the last non-blank character in
* the character string chstr
character*(*) chstr
lastchpos=len(chstr)
do while(lastchpos.gt.1.and.
: chstr(lastchpos:lastchpos).eq.'
')
lastchpos=lastchpos-1
end do
if(lastchpos.eq.1.and.chstr(1:1).eq.' ') then
lastchpos=0
end if
return
end