program fluo c---------------------------------------------------------------------- c version 002 2 May, 1999 c---------------------------------------------------------------------- c written by Daniel Haskel c Department of Physics c Box 351560 c University of Washington c Seattle, WA USA 98195 c c phone (206) 543-0435 c e-mail haskel@phys.washington.edu c---------------------------------------------------------------------- c description of the code: c c Fluo normalizes, aligns, and corrects xmu scans for self c absorption effects in fluorescence data. The normalization is done c by subtracting a pre-edge line then extrapolating a post-edge line to c e0. The e0 intercept is used for edge-step normalization and alignment. c The first derivative is also used in the alignment procedure. The c self absorption correction is evaluated using tabuted cross sections in c McMaster tables. c---------------------------------------------------------------------- c sample input file: c c e0 = 37442 c central = ba edge=k c atmwei = true c angin=60 angout=30 c formin=uwexafs formout=uwexafs c out =bafile c atoms c ba 1 c fe 12 c o 19 c files c bafeo_amor.xmu 1 * ba edge amorphous 36 K c bafeo_amor.xmu 3 * ba edge * amorphous 300 K c ! febao_crys300.dat * fe k-edge crystalline sample c---------------------------------------------------------------------- c version history c 001 pre-pre-pre version. No alignment. c 002 Mcfluo incorporated to calculate Xsections mutf, mub, muc c from atom list with relative amounts, central atom + edge. c Energy alignment of several files incorporated. c 5-21-99 If fluo-corrections calculated, added another edge-step c normalization after correction is done c---------------------------------------------------------------------- parameter(nptx=2**11, nfilx=50, ndocx=25, iou=10, $ nwdx=5) parameter(epsi=1.e-5, epsi2=2.e-3, zero=0.e0, pi= 3.141592654 ) character*5 vnum parameter(vnum='0.02') dimension energy(nptx), xmuraw(nptx), xmuout(nptx) dimension foo1(nptx), foo2(nptx), foo3(nptx) dimension nkey(nfilx) dimension der(nptx) logical eefind, stfind, vaxflg, lfluo, langle, le0fin character*3 cnum, test character*10 filtyp, frmin, frmout, skey(nfilx), skeyo, tmpe0 character*78 outfil, file(nfilx), line, of, string, $ commnt(nfilx), out2 character*100 doc(ndocx) real mutf, mub, muc real maxder, xmuint external xmuint common /idata/ nxmu save /idata/ common /rdata/ e0, energy, xmuout save /rdata/ character*78 comm 400 format(a) 410 format(1x,f13.6,3x,f13.6) 420 format(75('=')) 430 format(i3) 432 format(i3.3) 435 format(f9.2) 440 format('Fluo --',a5,'Correcting for self absorption.'16x, $ 'by Daniel Haskel') vaxflg = .false. c VAX USERS: change the line above to .true. c VAX VAX VAX VAX VAX VAX VAX VAX VAX VAX VAX VAX VAX VAX VAX VAX VAX test = 'abc' write(line,420) call messag(line) write(string,440)vnum call messag(string) call messag(line) call noinit(nkey, predg1, predg2, enor1, enor2, $ foo1, foo2, foo3, $ of, file, skey, commnt, angin, angout, $ mutf, mub, muc) call normin(nfile, nkey, file, e0, $ enor1, enor2, predg1, predg2, $ frmin, frmout, of, skey, commnt, $ mutf, mub, muc, angin, angout, $ le0fin, lfluo, langle) do 50 i=1,nfile c --- read the data filtyp = ' ' nxmu = nptx ndoc = ndocx call inpdat(filtyp, frmin, file(i), vaxflg, skey(i), nkey(i), $ ndoc, doc, nxmu, energy, xmuraw, foo1, foo2, foo3 ) call case(test,filtyp) call triml(frmin) if ((frmin(1:2).eq.'uw').and.(filtyp.ne.'xmu')) then call messag(' *** ERROR: uwexafs input file must be '// $ 'of type XMU.') call messag(' Skipping data set.') goto 50 endif if (frmout.eq.' ') frmout = frmin c --- run time message about input data file if (lfluo) then comm = ' Normalizing and correcting data from "' elseif (i.eq.1) then comm = ' Normalizing data from "' elseif (i.ne.1) then comm = ' Normalizing and aligning data from "' endif ii = istrln(file(i)) ic = istrln(comm) if (frmin(1:2).eq.'uw') then write(cnum,430)nkey(i) call messag(comm(:ic)//file(i)(:ii)//'", nkey = '//cnum) else call messag(comm(:ic)//file(i)(:ii)//'"') endif c --- remove pre-edge and normalize. use step as initial guess c --- in all other data step = 0 stfind = .true. slopre = 0 bpre = 0 if (i.eq.1) then if ((le0fin).and.(e0.gt.zero)) then eefind=.false. else eefind=.true. endif call prenew(eefind, e0, predg1, predg2, enor1, enor2, $ nxmu, energy, xmuraw, step, stfind, slopre, $ bpre, xmuout) if (eefind) then write(tmpe0, 435) e0 ccc write(tmpe0, '(f7.2)') e0 call messag(' e0 not provided; Found by fluo to be'// $ ' e0 = '//tmpe0) endif e0ref = e0 c print*,'e0ref=', e0ref if (step.lt.epsi) step = 1.e0 do 20 j=1,nxmu xmuout(j) = xmuout(j) / step 20 continue h=1.e0 do 21 j=1,nxmu der(j)= dfridr(xmuint,energy(j),h,err) 21 continue maxder=0.e0 do 22 j=1,nxmu maxder=max(der(j),maxder) 22 continue inmax=1 do 23 j=1,nxmu if (abs((der(j)-maxder)).lt.epsi) inmax=j 23 continue do 24 j=inmax,nxmu c if (abs((der(j)-zero)).lt.epsi2) then if (der(j).lt.zero) then inzero=j goto 25 endif 24 continue 25 continue c print*,'inmax=',inmax c print*,'inzero=',inzero derref=energy(inzero) emader=energy(inmax) c print*,'derref=',derref c print*,'emader=',emader else eefind = .true. call prenew(eefind, e0, predg1, predg2, enor1, enor2, $ nxmu, energy, xmuraw, step, stfind, slopre, $ bpre, xmuout) c print*, e0 e0shft = e0 - e0ref c print*,'e0shft=', e0shft if (step.lt.epsi) step = 1.e0 do 30 mm=1,nxmu xmuout(mm) = xmuout(mm) / step energy(mm) = energy (mm) - e0shft 30 continue h=1.e0 do 31 mm=1,nxmu der(mm)=dfridr(xmuint,energy(mm),h,err) 31 continue maxder=0.e0 do 32 mm=1,nxmu maxder=max(der(mm),maxder) 32 continue inmax=1 do 33 mm=1,nxmu if (abs((der(mm)-maxder)).lt.epsi) inmax=mm 33 continue do 34 mm=inmax,nxmu c if (abs((der(mm)-zero)).lt.epsi2) then if (der(mm).lt.zero) then inzero=mm goto 35 endif 34 continue 35 continue c print*,'inmax=',inmax c print*,'inzero=',inzero emader=energy(inmax) c print*,'emader=',emader dersh= energy(inzero) - derref c print*,'dersh=', dersh do 36 mm=1,nxmu energy(mm)=energy(mm) - dersh 36 continue endif if (lfluo) then if (langle) then angin = angin * (pi)/(180.0) angout = angout * (pi)/(180.0) else angin = 45.0 *(pi)/(180.0) angout = 45.0 *(pi)/(180.0) endif if (angout.le.epsi) angout = epsi mutf = mutf* (sin(angin))/(sin(angout)) alpha = ((mutf/muc) + (mub/muc)) betha = ( (mutf/muc) + (mub/muc) + 1 ) do 70 j=1,nxmu xmuout(j) = xmuout(j)*(alpha)/(betha - xmuout(j)) 70 continue do 71 j=1,nxmu xmuraw(j) = xmuout(j) 71 continue eefind = .true. call prenew(eefind, e0, predg1, predg2, enor1, enor2, $ nxmu, energy, xmuraw, step, stfind, slopre, $ bpre, xmuout) if (step.lt.epsi) step = 1.e0 do 72 mm=1,nxmu xmuout(mm) = xmuout(mm) / step 72 continue endif if (of.ne.' ') then outfil = of out2 = of if (frmout.eq.'ascii') then write(cnum,432)i jj=index(outfil,'.') if (jj .eq. 0) then im=istrln(outfil) outfil = outfil(1:im)//'.'//cnum out2 = outfil(1:im)//'_der.'//cnum else outfil = outfil(1:jj)//cnum out2 = outfil(1:jj-1)//'_der.'//cnum endif endif else if (frmout.eq.'ascii') then write(cnum,432)i ind = index(file(i), '.') ist = istrln(file(i)) if (ind.ne.0) then outfil = file(i)(:ind-1)//'.'//cnum out2 = file(i)(:ind-1)//'_der.'//cnum else outfil = file(i)(:ist)//'.'//cnum out2 = file(i)(:ist)//'_der.'//cnum endif else ind = index(file(i), '.') ist = istrln(file(i)) if (ind.ne.0) then outfil = file(i)(1:ind-1)//'.nor' out2 = file(i)(1:ind-1)//'_der'//'.nor' else outfil = file(i)(1:ist)//'.nor' out2 = file(i)(1:ist)//'_der'//'.nor' endif endif endif filtyp = 'xmu' nkeyo = 0 skeyo = ' ' iexist = 0 call addocn(ndoc, doc, commnt(i), step, slopre, bpre) call outdat(filtyp, frmout, outfil, vaxflg, skeyo, nkeyo, $ ndoc, doc, nxmu, energy, xmuout, foo1, foo2, foo3, iexist) call addocd(ndoc, doc, commnt(i), maxder, emader) call outdat(filtyp, frmout, out2, vaxflg, skeyo, nkeyo, $ ndoc, doc, nxmu, energy, der, foo1, foo2, foo3, iexist) ii = istrln(outfil) if (frmout(1:2).eq.'uw') then call messag(' >>> Wrote uwexafs output to "'// $ outfil(:ii)//'"') else call messag(' >>> Wrote ascii output to "'// $ outfil(:ii)//'"') endif 50 continue call messag(' fluo is done. Have a nice day ') call messag(line) stop c end main program fluo end subroutine noinit(nkey, predg1, predg2, enor1, enor2, $ foo1, foo2, foo3, $ of, file, skey, commnt, angin, angout, mutf, $ mub, muc) parameter(nptx=2**11, nfilx=50) character*10 skey(nfilx) character*78 file(nfilx), of, commnt(nfilx) integer nkey(nfilx) real mutf, mub, muc dimension energy(nptx), xmuout(nptx) dimension foo1(nptx), foo2(nptx), foo3(nptx) common /idata/ nxmu common /rdata/ e0, energy, xmuout parameter(zero=0.e0, nulval=-99 999) c --------- integers ---------------------------------------------- do 10 i=1,nfilx nkey(i) = 0 10 continue c --------- reals ------------------------------------------------- e0 = real(nulval-1) predg1 = -200.e0 predg2 = -50.e0 enor1 = 100.e0 enor2 = 500.e0 angin = 45.0 angout = 45.0 mutf = 10000.0 mub = 10000.0 muc = 10000.0 do 20 i=1,nptx energy(i) = zero xmuout(i) = zero foo1(i) = zero foo2(i) = zero foo3(i) = zero 20 continue c --------- logicals ---------------------------------------------- c eefind = .false. c lalign = .true. c lderiv = .false. c --------- characters -------------------------------------------- of = ' ' do 510 i=1,nfilx file(i) = ' ' skey(i) = ' ' commnt(i) = ' ' 510 continue return c end subroutine noinit end subroutine normin(nfile, nkey, file, e0, $ enor1, enor2, predg1, predg2, $ frmin, frmout, of, skey, commnt, $ mutf, mub, muc, angin, angout, $ le0fin, lfluo, langle) c-------------------------------------------------------------- c copyright 1996 university of washington Dani Haskel and Bruce Ravel c-------------------------------------------------------------- implicit real(a-h,o-z) implicit integer(i-n) c implicit double precision (a-h,o-z) c--------------------------------------------------------------------- c this parses the lines of the command file looking for keywords then c reads in the value for that keyword from the next word. c--------------------------------------------------------------------- parameter(nwdx=20, nfilx=50) parameter(zero=0.e0, nulval=-99999, imel=20) dimension perc(imel) dimension nkey(nfilx) character*2 test, cent, edge, elem(imel) character*10 frmin, frmout, skey(nfilx) character*20 fname, errmes*60 character*78 file(nfilx), string, messg, of, words(nwdx), $ commnt(nfilx), temp logical inpnul, there, le0fin, langle, lfluo, latoms, $ latwei real mutf, muc, mub 1410 format(bn,f10.0) 1440 format(a) c--------------------------------------------------------------------- c initialize some things used only in this routine iline = 0 mel = 0 inpnul = .true. test = 'ab' le0fin = .false. langle = .false. lfluo = .false. latoms = .false. latwei = .false. frmin = ' ' frmout = ' ' c the value of the variable test must be in the same case as the c keyword names in the long block of elseif's c--------------------------------------------------------------------- c open fluo.inp fname = 'fluo.inp' call upper(fname) inquire(file=fname,exist=there) if (.not.there) then call lower(fname) inquire(file=fname,exist=there) endif if (.not.there) then call messag(' ') call messag('* * * WARNING * * *') call messag('"fluo.inp" not found. fluo stopping.') call messag(' ') stop endif open(unit=1,file=fname,status='old') c--------------------------------------------------------------------- c begin reading the input file, words cleared each time 101 continue nwds = nwdx do 110 iw=1,nwds words(iw)=' ' 110 continue read (1,1440,end=191)string call untab(string) call triml(string) c - denotes end of job if (string(1:1).eq.'-') goto 191 c skip a comment line if ((string(1:1).eq.'!').or.(string(1:1).eq.'*') $ .or.(string(1:1).eq.'%')) goto 101 c skip a blank line if (string.eq.' ') goto 101 c begin reading line, i counts words i=1 call bwords(string,nwds,words) inpnul = .false. iline = iline+1 do 120 ib = 1,nwds call triml(words(ib)) 120 continue c************************ input file parsing *************************** c read a word, identify it, assign the value from the following word(s) c increment i and come back. i points to position in string, when i c exceeds nwds go read a new line. 130 continue call case(test,words(i)) c skip a blank line if (words(i).eq.' ') then goto 101 c ignore everything after !,*,% elseif ((words(i)(1:1).eq.'!').or.(words(i)(1:1).eq.'*').or. $ (words(i)(1:1).eq.'%')) then goto 101 c finish early elseif (words(i)(1:3).eq.'end') then goto 191 c central atom elseif (words(i)(1:4).eq.'cent') then cent=words(i+1)(1:2) i=i+2 c edge of interest elseif (words(i).eq.'edge') then edge=words(i+1)(1:2) i=i+2 c atomic weight flag elseif (words(i).eq.'atmwei') then if ((words(i+1)(1:1).eq.'t').or.(words(i+1)(1:1).eq.'T')) $ then latwei = .true. i=i+2 endif c e0 elseif ((words(i).eq.'e0').or.(words(i).eq.'ee')) then call getrea(words(i), words(i+1), e0) le0fin = .true. i=i+2 c pre-edge boundaries elseif (words(i).eq.'pre1') then call getrea(words(i), words(i+1), predg1) i=i+2 elseif (words(i).eq.'pre2') then call getrea(words(i), words(i+1), predg2) i=i+2 c normalization boundaries elseif (words(i).eq.'nor1') then call getrea(words(i), words(i+1), enor1) i=i+2 elseif (words(i).eq.'nor2') then call getrea(words(i), words(i+1), enor2) i=i+2 elseif (words(i).eq.'mutf') then call getrea(words(i), words(i+1), mutf) lfluo = .true. i=i+2 elseif (words(i).eq.'mub') then call getrea(words(i), words(i+1), mub) i=i+2 elseif (words(i).eq.'muc') then call getrea(words(i), words(i+1), muc) i=i+2 elseif (words(i).eq.'angin') then call getrea(words(i), words(i+1), angin) langle = .true. i=i+2 elseif (words(i).eq.'angout') then call getrea(words(i), words(i+1), angout) i=i+2 c format elseif (words(i).eq.'format') then frmin = words(i+1) frmout = words(i+1) call triml(frmin) call case(test,frmin) call triml(frmout) call case(test,frmout) i=i+2 c formin elseif (words(i).eq.'formin') then frmin = words(i+1) call triml(frmin) call case(test,frmin) i=i+2 c formout elseif (words(i).eq.'formout') then frmout = words(i+1) call triml(frmout) call case(test,frmout) i=i+2 c outfile elseif (words(i)(1:3).eq.'out') then of = words(i+1) i=i+2 c read in atoms list elseif (words(i)(1:4).eq.'atom') then mel = 0 do 528 ix=1,imel elem(ix) = ' ' perc(ix) = zero 528 continue 529 read (1,1440,end=191)string call untab(string) call triml(string) if (string(1:1).eq.'-') goto 191 c skip a comment line if ((string(1:1).eq.'!').or.(string(1:1).eq.'*') $ .or.(string(1:1).eq.'%')) goto 529 c skip a blank line if (string.eq.' ') goto 529 nwds = nwdx call bwords(string,nwds,words) do 720 ib = 1,nwds call triml(words(ib)) 720 continue if (words(1)(1:3).eq.'fil') then nfile=0 goto 140 else lfluo = .true. latoms = .true. mel=mel+1 elem(mel)=words(1)(1:2) call getrea('atomic concentration',words(2),perc(mel)) c print*, elem(mel), perc(mel) goto 529 endif c read in data filenames elseif ((words(i)(1:3).eq.'fil').or.(words(i)(1:3).eq.'dat')) $ then nfile = 0 140 continue read(1,1440,end=191)string call untab(string) call triml(string) nnw = nwdx do 150 ii=1,nwdx words(ii)=' ' 150 continue call bwords(string,nnw,words) do 170 n=1,nnw call triml(words(n)) 170 continue if ((words(1)(:1).eq.'!').or.(words(1)(:1).eq.' ').or. $ (words(1)(:1).eq.'%').or.(words(1)(:1).eq.'*')) then goto 140 elseif ((words(1)(1:3).eq.'---').or.(words(1).eq.'end')) then goto 191 elseif ((words(1).eq.'e0').or.(words(1)(1:4).eq.'mess').or. $ (words(1).eq.'nofit').or.(words(1)(1:4).eq.'noal').or. $ (words(1).eq.'pre1').or.(words(1)(1:4).eq.'norm').or. $ (words(1).eq.'pre2').or.(words(1)(1:3).eq.'out').or. $ (words(1).eq.'nor1').or.(words(1)(1:3).eq.'fil').or. $ (words(1).eq.'nor2').or.(words(1).eq.'data').or. $ (words(1).eq.'emax').or.(words(1).eq.'format').or. $ (words(1).eq.'ee').or. $ (words(1).eq.'formin').or.(words(1).eq.'formout')) $ then goto 130 else nfile = nfile+1 file(nfile) = words(1) call triml(file(nfile)) c print*, ' frmin = ', frmin(1:5) nf = 2 if (frmin(1:2).eq.'uw') then nf = 3 if (istrln(words(2)).le.3) then call getint('nkey',words(2),nkey(nfile)) else skey(nfile) = words(2) endif endif c will require user supplied file comment to be in same line after '*' if (frmin.eq.' ') then if ((words(2)(1:1).ne.'*').and. $ (words(2)(1:1).ne.' ')) then nf = 3 if (istrln(words(2)).le.3) then call getint('nkey',words(2),nkey(nfile)) else skey(nfile) = words(2) endif endif endif c print*, ' nfile, nkey, skey = ', nkey(nfile), skey(nfile) c --- individual file comments if ((nnw.ge.nf).and.(words(nf)(1:1).ne.'!').and. $ (words(nf)(1:1).ne.'%').and. $ (words(nf)(1:1).ne.'#').and. $ (words(nf)(1:1).eq.'*')) then commnt(nfile) = string do 160 k=1,nf-1 kk = istrln(words(k)) temp = commnt(nfile)(kk+1:) commnt(nfile) = temp call triml(commnt(nfile)) 160 continue endif goto 140 endif else iunk = istrln(words(i)) messg = words(i)(:iunk)//' is an unknown keyword.' call messag(messg) i=i+2 endif c if read entire line then read next line else read next word in line if (i.gt.nwds) goto 101 goto 130 c done reading lines 191 continue if (inpnul) then errmes = ' Your input file has no useful information ' goto 666 endif c calculate McMaster cross sections if (latoms) then print*, 'Atom list given ...' call mcfluo (cent, edge, mel, elem, perc, $ latwei, mutf, mub, muc) c print*, mutf, mub, muc endif if (nfile.eq.0) then errmes = 'You specified no input data.' goto 666 endif 300 continue return 666 continue ii = istrln(errmes) call messag('* * * ERROR! '//errmes(:ii)) call messag(' fluo is stopping.') stop c end subroutine normin end subroutine prenew(eefind, ee, predg1, predg2, enor1, enor2, $ nxmu, energy, xmuraw, step, stfind, slopre, bpre, xmuout) c c pre-edge subtraction and normalization of exafs data. c c inputs: c eefind logical flag for picking e0 at the maximum derivative c ee e0, energy origin of data c predg1 region for picking pre-edge line c predg2 region for picking pre-edge line c enor1 region for picking normalization c enor2 region for picking normalization c nxmu length of array energy, xmuraw, and xmuout c energy array of energy points c xmuraw array of raw absorption points c step edge step for normalization (found if zero on input) c stfind logical flag to tell whether to find edge step or not c outputs: c ee e0, energy origin of data c predg1 region used for picking pre-edge line c predg2 region used for picking pre-edge line c enor1 region used for picking normalization c enor2 region used for picking normalization c step edge step for normalization c slopre slope of pre-edge line c bpre intercept of pre-edge line (value at energy = 0) c xmuout array of unnormalized absorption points c after pre-edge line subtracted c c requires subroutine polyft (and so also function nofx) c c copyright 1992 university of washington : matt newville c Modified by Dani , march 96. Changes include checking pre1 and pre2 c and nor1, nor2 to be in data range. Also e0 has to be in data range. c Criteria changed to at least 6 consecutive increasing points. c----------------------------------------------------------------- integer ninc parameter (ninc = 6) logical eefind, stfind, inc(ninc) dimension energy(nxmu), xmuraw(nxmu), xmuout(nxmu), anorm(2) real zero, tiny, one, onepls parameter (zero = 0.0, tiny = 1.e-10, one = 1.0, onepls = 1.01) c c if ee was not specified, it is found in the data from c these conditions: c 1) ee is the third increasing point in a row. c 2) ee has the largest derivative up to that point. c these conditions will almost always pick a point on the c absorption step, and will fail less often than simply c taking the point with maximum derivative. c if (nxmu.le.5) return if ((.not. eefind).and. ((ee.lt.energy(1)).or. $ (ee.gt.energy(nxmu)))) then ccccc if (ee.lt.energy(1)) then eefind =.true. print*, ' The value of e0 given is outside data range' print*, ' fluo is finding its own value' endif if (eefind) then do 100 j = 1, ninc inc(j) = .false. 100 continue dxde = zero demx = zero ntry = max(2, int(nxmu/2)) + 3 do 150 i = 2, ntry deltax = xmuraw(i) - xmuraw(i-1) deltae = max(tiny, (energy(i) - energy(i-1))) dxde = deltax/deltae inc(1) = (dxde.gt.zero).and.(deltae.gt.tiny) if ((dxde.gt.demx).and.inc(6).and.inc(5).and.inc(4) $ .and.inc(3).and.inc(2).and.inc(1)) then ee = energy(i) demx = dxde * onepls end if if (inc(6).and.(.not.inc(5)).and.(.not.inc(4)) $ .and.(.not.inc(3)).and.(.not.inc(2)) $ .and.(.not.inc(1)) $ .and.(demx.gt.zero)) go to 200 do 130 j = ninc, 2, -1 inc(j) = inc(j - 1) 130 continue 150 continue end if 200 continue cc print*, ' pre-edge: loop stopped at = ', i c at this point, ee is known. so we go on c to do the pre-edge and normalization c c get pre-edge line (slope and intercept) by linear fit if ( (predg1.eq.zero).and.(predg2.eq.zero) ) then predg2 = -50 predg1 = -200 end if pre1 = ee + predg1 pre2 = ee + predg2 if (pre1.lt.energy(1)) pre1 = energy(2) if (pre2.lt.energy(1)) pre2 = ( ee + pre1 ) /2 call polyft(pre1, pre2, energy, xmuraw, nxmu, 2, anorm) bpre = anorm(1) slopre = anorm(2) c c apply e0 shift, do pre-edge subtraction, and calculate c the energy relative to edge for easier normalization do 300 i = 1, nxmu xmuout(i) = xmuraw(i) - bpre - slopre*energy(i) energy(i) = energy(i) - ee 300 continue c normalization : make pre-edge 0.0 and post-edge 1.0 c if step size wasn't given, get it by extracting to ee a c line that best fits the data on the range (ee+enor1,ee+enor2) if (stfind) then if ( (enor1.eq.zero).and.(enor2.eq.zero) ) then enor1 = 100 enor2 = 300 end if if (enor2.gt.energy(nxmu)) enor2 = energy(nxmu) if (enor1.gt.energy(nxmu)) enor1 = ( ee + enor2 ) /2 call polyft(enor1, enor2, energy, xmuout, nxmu, 2, anorm) step = anorm(1) end if if (abs(step).lt.tiny) step = one c put energy array back to original values, c recover from bad pre-edge fits do 500 i = 1, nxmu energy(i) = energy(i) + ee c for "post-edge" pre-edge fit: if (predg1.gt.zero) then xmuout(i) = one + xmuout(i) end if 500 continue return c end subroutine prenew end subroutine addocn(ndoc, doc, commnt, step, slopre, bpre) parameter(ndocx=25) character*78 commnt character*100 doc(ndocx) c If file comment is given, need to make room for 2 new lines; if not c just need one document line. Copy line i to line i+2 (or i+1). The c first line of doc in not altered. nnd = 2 if (commnt.ne.' ') then do 10 i=min(ndoc, ndocx-nnd), nnd, -1 doc(i+2) = doc(i) 10 continue write(doc(nnd),400)commnt(2:) 400 format('fluo: ',a) write(doc(1+nnd),410)step, slopre, bpre 410 format('fluo: STEP = ',f7.4', pre-edge line: mux = ', $ f9.7' *E +',f9.5) else do 20 i=min(ndoc, ndocx-nnd), nnd, -1 doc(i+1) = doc(i) 20 continue write(doc(nnd),410)step, slopre, bpre endif return c end subroutine addoc end c---------------------------------------------------------------------- c input/output routines for data files c for the uwxafs programs c c input/output routines for data files for the uwxafs programs c c copyright 1992 university of washington, seattle, washington c written by matthew newville c department of physics, fm-15 c university of washington c seattle, wa usa 98195 c phone (206) 543-0435 c e-mail newville@u.washington.edu c c these routines are the basic input/output routines for getting c numerical and document data from files into the uwxafs programs. c there are currently two data formats supported: c c 1. 'uw' : a binary file format known as the uwxafs file handling c routines. this is very efficient way to store data, and c can store several (191) data sets in a single file. the c drawback is that the files are not extremely portable. c c 2. 'asc': these are column files in a format that is fairly easy c for anything to deal with. the files have several lines c of documents. if the first character of the document is c '#' this character will be removed. after the documents c is a line with minus signs for characters(3:6), then an c ignored line (for column labels), and then the data. up c to five columns are used. the expected order is: c x, real(y), imag(y), ampl(y), phase(y). c if any column representing y is zero, the appropriate c value will be calculated and returned. the files in this c format hold only one data set, and use more memory than c the uwxafs files, but are portable and convenient. c c other file types can be added without too much difficulty. c the routines listed here are: c inpdat : retrieve data and documents from a file c inpcol : retrieve data and documents from an ascii file c inpuwx : retrieve data and documents from a uwxafs file c outdat : write data and documents to a file c outcol : write data and documents to an ascii file c outuwx : write data and documents to a uwxafs file c c note: the fortran input/output unit number 11 is used for all c unit numbers in these routines. conflicts between these c routines will not happen, but conflicts may arise if c unit = 11 indicates an open file in a calling subprogram. c---------------------------------------------------------------------- subroutine inpdat(filtyp, format, filnam, vax, skey, nkey, $ ndoc, doc, ndata, xdata, yreal, yimag, yampl, yphas ) c c copyright 1992 university of washington : matt newville c c retrieve data and documents from a file acording c to the format specified by 'format'. c inputs: c filtyp file type to open. if may be ' ' c format file format (uwxafs, ascii, column) c filnam file name c vax logical flag for being on a vax machine (binary file) c skey symbolic key for record in uwxafs file c ndata maximum number of elements in data arrays c nkey numeric key for record in uwxafs file c ndoc maximum number of document lines to get c note: ndoc cannot be less than or equal to zero! c outputs: c skey symbolic key of record in uwxafs file c ndoc number of document lines returned c doc array of document lines c ndata number of elements in data arrays c xdata array of x values of data c yreal array of real part of y data values c yimag array of imaginary part of y data values c yampl array of amplitude part of y data values c yphas array of phase part of y data values c--------------------------------------------------------------------- character*(*) filtyp, format, skey, filnam, doc(*) character*10 type, symkey, form, formin, errmsg*80 dimension xdata(*), yreal(*), yimag(*) dimension yampl(*), yphas(*) logical vax integer irecl, ndatmx, ndocmx data irecl, ndocmx, ndatmx / 512 , 19, 4096/ c--------------------------------------------------------------------- c some initializations if (vax) irecl = 128 type = filtyp call triml(type) call upper(type) symkey = skey call triml(symkey) call upper(symkey) c c determine format of the input file formin = format call triml(formin) call smcase(formin, 'a') call testrf(filnam, irecl, form, ier) call smcase(form, 'a') if (ier.eq.-1) then call messag(' inpdat error: file not found ') elseif (ier.eq.-2) then call messag(' inpdat error: unknown file format = '//formin) elseif (ier.eq.-3) then call messag(' inpdat error: poorly formatted ascii data? ') elseif (ier.eq.-4) then call messag(' inpdat error: no data in ascii format file? ') end if if (ier.ne.0) then errmsg = ' for file ' // filnam ilen = istrln(errmsg) call messag( errmsg(1:ilen) ) stop endif if ((formin.ne.' ').and.(formin(1:2).ne.form(1:2))) then call messag(' inpdat warning: the requested format was'// $ ' incorrect!') call messag(' form = '//form(1:5) ) call messag(' formin = '//formin(1:5) ) end if c now call the appropriate routine to get the data, c according to the format. ndata = max(1, min(ndata, ndatmx) ) ndoc = max(1, min(ndoc , ndocmx) ) if (form(1:2).eq.'uw') then ndoc = ndocmx call inpuwx(type, filnam, skey, nkey, irecl, ndoc, doc, $ ndata, xdata, yreal, yimag, yampl, yphas ) elseif ((form(1:2).eq.'co').or.(form(1:2).eq.'as')) then call inpcol(filnam, ndoc, doc, $ ndata, xdata, yreal, yimag, yampl, yphas ) skey = 'ascii' call upper(skey) else call messag(' inpdat error: unknown file format = '// form) ilen = min(54, max(1, istrln(filnam))) errmsg = ' for file ' // filnam(1:ilen) call messag( errmsg(1:ilen+26) ) stop end if filtyp = type format = form c return c end subroutine inpdat end subroutine inpcol(filnam, ndoc, doc, ndata, $ xdata, yreal, yimag, yampl, yphas) c c copyright 1992 university of washington : matt newville c c open and get all information from a column file. document c lines are read until a line of '----', then a label line is c skipped and the column data are read in. the data is read c and stored in the following order: c xdata yreal yimag yampl yphas c inputs: c filnam file name containing data c ndoc maximum number of document lines to get c ndata maximum number of elements in data arrays c outputs: c ndoc number of document lines returned c doc array of document lines c ndata number of elements in data arrays c xdata array of x values of data c yreal array of real part of y data values c yimag array of imaginary part of y data values c yampl array of amplitude part of y data values c yphas array of phase part of y data values c--------------------------------------------------------------------- integer ilen , istrln, j, i, mxword, ndoc, ndata, iounit integer iexist, ierr, nwords, idoc, id real zero , small parameter( zero = 0.e0 , small = 1.e-6, mxword = 5) real xdata(*), yreal(*), yimag(*), yampl(*), yphas(*) real xinp(mxword) logical isdat character*(*) filnam, doc(*) character*30 words(mxword), line*100, status*10, file*60 external istrln, isdat c--------------------------------------------------------------------- 10 format(a) file = filnam ilen = istrln(file) if (ilen.le.0) then call messag( ' inpcol: no file name given') stop end if c initialize buffers do 80 j = 1, ndoc doc(j) = ' ' 80 continue do 100 i = 1, mxword words(i) = '0. ' xinp(i) = zero 100 continue do 120 j = 1, ndata xdata(j) = zero yreal(j) = zero yimag(j) = zero yampl(j) = zero yphas(j) = zero 120 continue c open data file iounit = 7 status ='old' call openfl(iounit, filnam, status, iexist, ierr) if ((iexist.lt.0).or.(ierr.lt.0)) go to 900 c c get documents from header: up to ndoc c read file header, save as document lines, c remove leading '#' and '%' both of which are c known to be extraneous comment characters. nwords = 5 idoc = 0 id = 1 200 continue read(iounit, 10, end = 950, err = 960) line call triml (line) c if line is '----', read one more line, go read numerical data if (line(3:6) .eq. '----') then read(iounit, 10, end = 950, err = 960) line goto 400 end if c remove leading '#' or '%' from line if ( (line(1:1).eq.'#').or.(line(1:1).eq.'%') ) then line(1:1) = ' ' call triml(line) c if the line is all numbers, then this is data! elseif (isdat(line)) then goto 410 end if c save line in doc if there's room if ((idoc .lt. ndoc) .and. (istrln(line).gt.0) ) then idoc = idoc + 1 doc(idoc) = line endif goto 200 c c read numerical data 400 continue nwords = 5 read(iounit, 10, end = 600, err = 980) line 410 continue call untab(line) call bwords(line,nwords,words) if (nwords.le.1) goto 600 do 450 i = 1, nwords c >>> changed for phit call sgetvl(words(i), xinp(i), ierr) c <<< if (ierr.ne.0) goto 600 450 continue xdata(id) = xinp(1) yreal(id) = xinp(2) yimag(id) = xinp(3) yampl(id) = xinp(4) yphas(id) = xinp(5) if (id.ge.ndata) go to 610 if ( id.ge.2 ) then if (abs(xdata(id) - xdata(id-1)).ge.small) id = id + 1 else id = id + 1 end if goto 400 600 continue id = id - 1 if (id.lt.1) go to 950 610 continue ndata = id if (idoc.le.0) then ndoc = 1 doc(1) = 'inpdat: no document line found' end if c make sure that all columns are filled: c if yampl and yphas are both zero, compute them from yreal, yimag c if yreal and yimag are both zero, compute them from yampl, yphas do 800 i = 1, ndata if ( ( (yampl(i).eq.zero).and.(yphas(i).eq.zero) ) .and. $ ( (yreal(i).ne.zero).or. (yimag(i).ne.zero) ) ) then yampl(i) = sqrt( yreal(i)**2 + yimag(i)**2 ) yphas(i) = atan2( yimag(i), yreal(i) ) if (i.gt.1) call pijump( yphas(i), yphas(i-1) ) elseif ( (yreal(i).eq.zero).and.(yimag(i).eq.zero) $ .and.(yampl(i).ne.zero) ) then yreal(i) = yampl(i) * cos ( yphas(i) ) yimag(i) = yampl(i) * sin ( yphas(i) ) end if 800 continue c close data file and return close(iounit) return c error handling c open file - error 900 continue call messag(' inpcol: error opening file '//file(1:ilen) ) go to 990 c end or error at reading documents 950 continue 960 continue call messag( ' inpcol: error reading file '//file(1:ilen) ) call messag(' during reading of documents.') go to 990 c error at reading numerical data 980 continue call messag( ' inpcol: error reading file '//file(1:ilen) ) call messag(' during reading of numerical data.') 990 continue close(iounit) stop c end error handling c end subroutine inpcol end subroutine inpuwx(ftypin, filein, skey, nkey, irecl, ndoc, $ documt, ndata, xdata, yreal, yimag, yampl, yphas ) c c copyright 1992 university of washington : matt newville c c open and get all information from a uwxafs file c c inputs: c ftypin file type to open, checked for compatibility, may be ' ' c filein file name containing data c skey symbolic key for record in data file (only one of these) c nkey numeric key for record in data file (two is needed ) c ndoc maximum number of document lines to get c outputs: c skey symbolic key of record in uwxafs file c ndoc number of document lines returned c docu array of document lines c ndata number of elements in data arrays c xdata array of x values of data c yreal array of real part of y data values c yimag array of imaginary part of y data values c yampl array of amplitude part of y data values c yphas array of phase part of y data values c c notes: c 1 the full 'noabort' error checking is done for the calls to c the uwxafs routines, which means that marginally useful c error messages will be given when one of the uwxafs c filehandling routines dies. c c 2 currently, the following file types are supported: c xmu, chi, rsp, env, rspep, rip c c 3 uwxafs file handling routines only do single precision. c this routine can be made implicit double precision if the c array buffer is maintained as single precision: c implicit double precision(a-h,o-z) c real buffer(maxpts) c--------------------------------------------------------------------- parameter( maxpts = 2048, zero = 0. ) character*(*) ftypin, skey, filein, documt(*) character*10 type, ftype, safefl*8, abrtfl*8 character*70 filnam, messg dimension xdata(*), yreal(*), yimag(*) dimension yampl(*), yphas(*) real buffer(maxpts) c--------------------------------------------------------------------- c initialize 10 format(a) 20 format(2x,2a) 30 format(2x,a,i3) safefl = ' ' abrtfl = 'noabort' call upper(abrtfl) ftype = ftypin filnam= filein call upper(skey) call triml(skey) call triml(ftype) call triml(filnam) ilen = max(1, istrln(filnam)) c note: uwxafs requires ftype to be upper case. call upper (ftype) do 100 i = 1,ndata xdata(i) = zero yreal(i) = zero yimag(i) = zero yampl(i) = zero yphas(i) = zero 100 continue do 110 i = 1, maxpts buffer(i) = zero 110 continue c call uwxafs file handling routines: c : open data file iounit = 11 call openrf(iounit, filnam, abrtfl, safefl, ftype, irecl, ier) if (ier.ne.0) then messg = 'inpuwx: error opening file ' call messag(messg//filnam(:ilen)) write (messg, '(9x,a,i4)') 'openrf error code ',ier call messag(messg) stop end if c : check file type call gftype(iounit, type, ier) if (ier.ne.0) then messg = 'inpuwx: error getting file type for ' call messag(messg//filnam(:ilen)) write (messg, '(9x,a,i4)') 'gftype error code ',ier call messag(messg) stop end if call upper(type) if (ftype.eq.' ') then ftype = type elseif (ftype.ne.type) then messg = 'inpuwx: incorrect file type for ' call messag(messg//filnam(:ilen)) messg = ' file type for this file is ' call messag(messg//type) messg = ' file type requested was ' call messag(messg//ftype) stop endif ftypin = ftype c : find out how many records there are in the file call gnie (iounit, nie, ier) if (nie.le.0) then messg = 'inpuwx: no data records in ' call messag(messg//filnam(:ilen) ) stop end if c : get skey if it wasn't given as input if (skey.eq.' ') then call gskey(iounit, nkey, skey, ier) if (ier.ne.0) then messg = 'inpuwx: error getting skey for ' call messag(messg//filnam(:ilen)) write (messg, '(9x,a,i4)') 'gskey error code ',ier call messag(messg) stop end if if (skey.eq.' ') then write (messg, '(1x,2a,i4)') 'inpuwx: found no skey ', $ 'for nkey =',nkey call messag(messg) call messag(' in file = '//filnam(:ilen)) stop end if end if c : get nkey if it wasn't given as input if (nkey.eq.0) then call gnkey(iounit, skey, nkey, ier) if (ier.ne.0) then messg = 'inpuwx: error getting nkey for ' call messag(messg//filnam(:ilen)) write (messg, '(9x,a,i4)') 'gnkey error code ',ier call messag(messg) stop end if end if c c : get documents : up to ndoc c first check how many document lines there are call gdlen(iounit, nkey, ndocln, ier) if (ier.ne.0) then messg = 'inpuwx: error getting document length for ' call messag(messg//filnam(:ilen)) write (messg, '(9x,a,i4)') 'gdlen error code ',ier call messag(messg) stop end if if (ndoc.gt.ndocln) ndoc = ndocln c then get the documents call getdoc(iounit, documt, ndoc, skey, nkey, ndsent, ier) if (ier.eq.6) then messg = 'inpuwx error: reading file ' call messag(messg//filnam(:ilen) ) messg = ' no skey or nkey given to specify record, ' call messag(messg) messg = ' or an incorrect skey or nkey given ' call messag(messg) stop elseif (ier.ne.0) then messg = 'inpuwx: error getting documents for ' call messag(messg//filnam(:ilen)) write (messg, '(9x,a,i4)') 'getdoc error code ',ier call messag(messg) stop end if ndoc = ndsent c : get data call getrec(iounit, buffer, maxpts, skey, nkey, nbuff, ier) if (ier.ne.0) then messg = 'inpuwx: error getting data for ' call messag(messg//filnam(:ilen)) write (messg, '(9x,a,i4)') 'getrec error code ',ier call messag(messg) stop end if c : close file call closrf(iounit,ier) if (ier.ne.0) then messg = 'inpuwx: error closing data file ' call messag(messg//filnam(:ilen)) write (messg, '(9x,a,i4)') 'closrf error code ',ier call messag(messg) stop end if c----------------------------------------------------------------- c finished with uwxafs routines, so now sort the data into c xdata, re(y), imag(y), ampl(y), phase(y) according to file type c c convert ftype to the case of this routine. c 'case' controls the the case of this routine call smcase (ftype, 'case') c- xmu: nbuff energy, then nbuff y-values if (ftype.eq.'xmu') then ndata = nbuff/2 do 400 i = 1, ndata xdata(i) = buffer(i) yreal(i) = buffer(ndata + i) yampl(i) = yreal(i) 400 continue c- chi: xmin, deltax, chi(kmin + i*deltak) elseif (ftype.eq.'chi') then ndata = nbuff - 2 do 500 i = 1, ndata xdata(i) = buffer(1) + (i-1)*buffer(2) yreal(i) = buffer(2 + i) yampl(i) = yreal(i) 500 continue c- env,rspep: kmin, deltak, phase, amplitude pairs (kmin + i*deltak) elseif ( (ftype.eq.'env').or.(ftype.eq.'rspep') ) then ndata = (nbuff - 1) / 2 do 600 i = 1, ndata xdata(i) = buffer(1) +(i-1)*buffer(2) yphas(i) = buffer(2*i+1) yampl(i) = buffer(2*i+2) yreal(i) = yampl(i) * cos ( yphas(i) ) yimag(i) = yampl(i) * sin ( yphas(i) ) 600 continue c rsp, rip: kmin, deltak, real, imaginary pairs (kmin + i*deltak) elseif ( (ftype.eq.'rsp').or.(ftype.eq.'rip') ) then ndata = (nbuff - 1) / 2 do 700 i = 1, ndata xdata(i) = buffer(1) +(i-1)*buffer(2) yreal(i) = buffer(2*i+1) yimag(i) = buffer(2*i+2) yampl(i) = sqrt( yreal(i)**2 + yimag(i)**2 ) yphas(i) = atan2( yimag(i), yreal(i) ) if (i.gt.1) call pijump( yphas(i), yphas(i-1) ) 700 continue else messg = 'inpuwx: unrecognized file type for ' call messag(messg//filnam(:ilen)) messg = ' file type for this file is ' call messag(messg//ftype) stop end if return c end subroutine inpuwx end subroutine outdat(filtyp, format, filnam, vax, skey, nkey, $ ndoc, doc, ndata, xdata, yreal, yimag, yampl, yphas, iexist) c c copyright 1992 university of washington : matt newville c c write data and documents to a file acording to the c format specified c inputs: c filtyp file type to open, may be ' '. c format file format (uwxafs, ascii, column) c filnam file name c vax logical flag for being on a vax machine (binary file) c ndoc number of document lines returned c doc array of document lines c ndata number of elements in data arrays c xdata array of x values of data c yreal array of real part of y data values c yimag array of imaginary part of y data values c yampl array of amplitude part of y data values c yphas array of phase part of y data values c iexist flag for whether to write redundant data to uwxafs file c iexist = 1 : do not write redundant data c iexist = 0 : do write redundant data c outputs: c skey symbolic key for record in uwxafs file c nkey numeric key for record in uwxafs file c ndoc number of document lines written c c--------------------------------------------------------------------- character*(*) filtyp, format, filnam, skey, doc(*) character*30 type, form dimension xdata(*), yreal(*), yimag(*) dimension yampl(*), yphas(*) logical vax integer irecl c--------------------------------------------------------------------- irecl = 512 if (vax) irecl = 128 c idoc = ndoc skey = ' ' form = format type = filtyp call upper(type) call triml(type) call triml(form) c convert form to the case of this routine. c 'case' controls the the case of this routine call smcase (form, 'case') c if (form(1:2).eq.'uw') then if ( (idoc.le.0).or.(idoc.gt.19) ) idoc = 19 call outuwx(type, filnam, skey, nkey, irecl, idoc, doc, $ ndata, xdata, yreal, yimag, yampl, yphas, iexist) elseif ( (form(1:3).eq.'col').or.(form(1:3).eq.'asc') ) then call outcol(type, filnam, idoc, doc, $ ndata, xdata, yreal, yimag, yampl, yphas) skey = 'ascii' call upper(skey) else call messag('outdat: unknown file format = '//form) stop end if c return c end subroutine outdat end subroutine outcol(filtyp, filnam, ndoc, doc, ndata, $ xdata, yreal, yimag, yampl, yphas) c c copyright 1992 university of washington : matt newville c c open and write all information to a column file. document lines are c written, followed by a line of '----', then a label line, and then c the data are written. the file type tells what to use for the label c and how many columns to write. it may be left blank. c c inputs: c filtyp file type to write (may be ' ' : used for label only) c filnam file name to write (' ' and '*' mean write to unit 6) c ndoc maximum number of document lines to write c documt array of document lines c ndata number of elements in data arrays c xdata array of x values of data c yreal array of real part of y data values c yimag array of imaginary part of y data values c yampl array of amplitude part of y data values c yphas array of phase part of y data values c outputs: c ndoc number of document lines written c c--------------------------------------------------------------------- parameter(zero = 0.000000) character*(*) filtyp, filnam, doc(*) character*80 filout, errmsg, type*10, status*10 dimension xdata(*), yreal(*), yimag(*) dimension yampl(*), yphas(*) c--------------------------------------------------------------------- type = filtyp call triml(type) c convert type to the case of this routine. call smcase(type, 'a') filout = filnam call triml(filout) if (ndata.le.0) ndata = 2 25 format(2a) 26 format(a,5x,a) c open data file c if file name is ' ' or '*', write to standard output (unit 6) iounit = 6 if ((filout.ne.' ').and.(filout.ne.'*')) then iounit = 0 status ='unknown' call openfl(iounit, filout, status, iexist, ierr) if ((ierr.lt.0).or.(iexist.lt.0)) go to 990 endif c c write documents to header jdoc = 0 mxl = 76 mxlp1 = mxl + 1 do 200, idoc = 1, ndoc call triml(doc(idoc)) ilen = istrln(doc(idoc)) if (ilen.ge.1) then jdoc = jdoc + 1 if (ilen.gt.mxl) then write(iounit, 25) '# ', doc(idoc)(1:mxl) write(iounit, 26) '# ', doc(idoc)(mxlp1:ilen) else write(iounit, 25) '# ', doc(idoc)(1:ilen) end if end if 200 continue ndoc = jdoc write(iounit, 25) '#----------------------------------', $ '-----------------------------------' c write line of minus signs and the label for columns if (type.eq.'xmu') then write(iounit, 25) '# energy xmu ', ' ' c >>> changed for phit elseif (type.eq.'2col') then write(iounit, 25) '# x y(x) ', ' ' c <<< elseif (type.eq.'chi') then write(iounit, 25) '# k chi(k) ', ' ' elseif (type.eq.'env') then write(iounit, 25) '# k real(chi(k)) ', $ 'imag(chi(k)) ampl(chi(k)) phase(chi(k))' elseif ( (type.eq.'rsp').or.(type.eq.'rspep') ) then write(iounit, 25) '# r real(chi(r)) ', $ 'imag(chi(r)) ampl(chi(r)) phase(chi(r))' else write(iounit, 25) '# x real(y(x)) ', $ 'imag(y(x)) ampl(y(x)) phase(y(x)) ' end if c c make sure that all of re(y), im(y), amp(y), and phase(y) are known do 300 i = 1, ndata if ( ( (yampl(i).eq.zero).and.(yphas(i).eq.zero) ) .and. $ ( (yreal(i).ne.zero).or. (yimag(i).ne.zero) ) ) then yampl(i) = sqrt( yreal(i)**2 + yimag(i)**2 ) yphas(i) = atan2( yimag(i), yreal(i) ) if (i.gt.1) call pijump( yphas(i), yphas(i-1) ) elseif ( (yreal(i).eq.zero).and.(yimag(i).eq.zero) $ .and.(yampl(i).ne.zero) ) then yreal(i) = yampl(i) * cos ( yphas(i) ) yimag(i) = yampl(i) * sin ( yphas(i) ) end if 300 continue c c write out data : some file types only write out a few columns 380 format(2x,e13.7,3x,e13.7) 390 format(2x,e13.7,2x,e13.7,2x,e13.7,2x,e13.7,2x,e13.7) if ( (type.eq.'xmu').or.(type.eq.'chi').or. $ (type.eq.'2col') ) then do 400 i = 1, ndata write(iounit, 380) xdata(i), yreal(i) 400 continue else do 450 i = 1, ndata write(iounit, 390) xdata(i), yreal(i), yimag(i), $ yampl(i), yphas(i) 450 continue end if c c close data file and return close(iounit) return 990 continue ilen = max(1, istrln(filnam)) errmsg = 'outcol: error opening file '//filnam(:ilen) imsg = istrln(errmsg) call messag(errmsg(:imsg)) stop c end subroutine outcol end subroutine outuwx(ftypin, filein, skey, nkey, irecl, ndoc, doc, $ ndata, xdata, yreal, yimag, yampl, yphas, iexist) c c write out data and documents to a uwxafs file c c inputs: c ftypin file type to write to, may be ' ' if filnam exists. c filein file name to write to c skey symbolic key of record in uwxafs file c ndoc number of document lines returned c doc array of document lines c ndata number of elements in data arrays c xdata array of x values of data c yreal array of real part of y data values c yimag array of imaginary part of y data values c yampl array of amplitude part of y data values c yphas array of phase part of y data values c iexist flag for whether to write redundant data to file c iexist = 1 : do not write redundant data c iexist = 0 : do write redundant data c c copyright 1992 university of washington : matt newville c----------------------------------------------------------------------- parameter(maxpts = 2048, maxdoc = 19, zero = 0.) character*(*) filein, ftypin, doc(*), skey character*10 skyout, ftype, type, filnam*70, messg*70 character*100 docout(maxdoc), abrtfl*8, safefl*8 dimension xdata(*), yreal(*), yimag(*) dimension yampl(*), yphas(*) real buffer(maxpts) c----------------------------------------------------------------------- c initialize 10 format(a) safefl = ' ' abrtfl = 'noabort' call upper(abrtfl) skyout = ' ' type = ' ' filnam = filein call triml(filnam) ilen = max(1, istrln(filnam)) ftype = ftypin call upper(ftype) c >>> new for phit if (ftype.eq.'2COL') ftype = 'XMU' c <<< do 60 i = 1, maxdoc docout(i) = ' ' 60 continue c output documents idoc = 0 80 continue idoc = idoc + 1 if ((idoc.ge.maxdoc).or.(idoc.gt.ndoc)) then idoc = idoc - 1 go to 100 end if docout(idoc) = doc(idoc) call triml(docout(idoc)) go to 80 100 continue ccccc ndoc = idoc c open data file to check file type iounit = 11 call openrf(iounit, filnam, abrtfl, safefl, ftype, irecl, ier) if (ier.ne.0) then messg = 'outuwx: error opening file '//filnam(:ilen) imsg = max(1, istrln(messg)) call messag(messg(:imsg)) write(messg, '(9x,a,i3)' ) 'openrf error code ',ier call messag(messg) stop end if c check file type call gftype(iounit, type, ier) call upper(type) c if file type was not given, close and the re-open data file c with file type just found, so we can write to file if (ftype.eq.' ') then ftype = type call closrf(iounit,ier) call openrf(iounit, filnam, abrtfl, safefl, ftype,irecl,ier) c if file type was given but it was wrong, stop elseif (ftype.ne.type) then messg = 'outuwx: incorrect file type for ' call messag(messg//filnam(:ilen)) messg = ' file type for this file is ' call messag(messg//type) messg = ' file type requested was ' call messag(messg//ftype) stop endif c c make sure that all of re(y), im(y), amp(y), and phase(y) are known do 300 i = 1, ndata if ( ( (yampl(i).eq.zero).and.(yphas(i).eq.zero) ) .and. $ ( (yreal(i).ne.zero).or. (yimag(i).ne.zero) ) ) then yampl(i) = sqrt( yreal(i)**2 + yimag(i)**2 ) yphas(i) = atan2( yimag(i), yreal(i) ) if (i.gt.1) call pijump( yphas(i), yphas(i-1) ) elseif ( (yreal(i).eq.zero).and.(yimag(i).eq.zero) $ .and.(yampl(i).ne.zero) ) then yreal(i) = yampl(i) * cos ( yphas(i) ) yimag(i) = yampl(i) * sin ( yphas(i) ) end if 300 continue c c put data into a single buffer according to data type c convert ftype to the case of this routine. c 'case' controls the the case of this routine call smcase(ftype, 'case') c usually buffer(1) and buffer(2) are xdata(1) and xdata(2) -xdata(1) buffer(1) = xdata(1) buffer(2) = xdata(2) - xdata(1) c xmu: nbuff energy, then nbuff y-values if (ftype.eq.'xmu') then nbuff = 2*ndata do 400 i = 1, ndata buffer(i) = xdata(i) buffer(ndata + i) = yreal(i) 400 continue c chi: kmin, deltak, chi(kmin + i*deltak) elseif (ftype.eq.'chi') then nbuff = ndata + 2 do 500 i = 1, ndata buffer(2 + i) = yreal(i) 500 continue c env: kmin, deltak, phase, amplitude pairs (kmin + i*deltak) elseif ( (ftype.eq.'env').or.(ftype.eq.'rspep') ) then nbuff = 2* (ndata + 1) do 600 i = 1, ndata buffer(2*i+1) = yphas(i) buffer(2*i+2) = yampl(i) 600 continue c rsp: kmin, deltak, real, imaginary pairs (kmin + i*deltak) elseif ( (ftype.eq.'rsp').or.(ftype.eq.'rip') ) then nbuff = 2* (ndata + 1) do 700 i = 1, ndata buffer(2*i+1) = yreal(i) buffer(2*i+2) = yimag(i) 700 continue c other data types not yet supported else call messag('outuwx: not able to decipher ftype ='//ftype) stop end if c c generate skyout for data with hash call hash(buffer, nbuff, docout, idoc, skyout) c check if this record is already in the file, c and decide whether or not to write data and c documentation for the record to the file call gnkey(iounit, skyout, nkey, ier) if ( (iexist.eq.1).and.(nkey.ne.0) ) then skey = ' ' else call putrec(iounit, buffer, nbuff, skyout, 0, ier) call putdoc(iounit, docout, idoc, skyout, 0, ier) skey = skyout end if ftypin = ftype c close file and leave call closrf(iounit, ierr) return c end subroutine outuwx end subroutine getrea(keywrd,string,value) character*(*) keywrd, string character*72 messg integer j, k, id, ie real value logical isnum 400 format(bn,f13.0) 410 format(bn,e19.13) 420 format(bn,d19.13) if (isnum(string)) then call lower(string) ie = index(string, 'e') id = index(string, 'd') if ((id.eq.0).and.(ie.eq.0)) then read(string,400)value elseif (ie.ne.0) then read(string,410)value elseif (id.ne.0) then read(string,420)value endif else j = istrln(string) k = istrln(keywrd) messg = 'Error reading "'//string(:j)//'" as the value for "' $ //keywrd(:k)//'"' call messag(messg) messg = '"'//string(:j)//'" must be a real number.' call messag(messg) stop endif return c end subroutine getrea end c********************************************************************** subroutine getint(keywrd,string,ivalue) character*(*) keywrd, string character*72 messg integer ivalue, j, k 400 format(bn,i10) read(string,400,iostat=ierr)ivalue if (ierr.ne.0) then j = istrln(string) k = istrln(keywrd) messg = 'Error reading '//string(:j)//' as the value for '// $ keywrd(:k) call messag(messg) messg = string(:j)//' must be an integer.' call messag(messg) stop endif return c end subroutine getint end c********************************************************************** subroutine getlog(keywrd,string,lvalue) character*(*) keywrd, string character*72 test*2 logical lvalue test = 'ab' lvalue = .false. call triml(string) call case(test,string) if ( (string(:1).eq.'t') .or. (string(:1).eq.'y') $ .or. (string(:2).eq.'on') ) $ lvalue=.true. return c end subroutine getlog end c********************************************************************** subroutine gettit(keywrd,string,title,ntit) character*(*) keywrd, string, title character*72 messg, toss integer ntit, i, j ntit = ntit + 1 toss = string call case(keywrd,toss) i = index(toss, keywrd) j = istrln(keywrd) title = string(i+j+1:60) call triml(title) if ( (title(:1).eq.'=') .or. (title(:1).eq.',') ) then toss = title(2:) title = toss call triml(title) endif messg = ' title > '//title call messag(messg) return c end subroutine gettit end c subroutine case(test,word) c c returns *word* in the same case as *test* c c note that this is just the reverse of smcase ! c character*(*) test, word c call smcase (word, test) c return c c end subroutine case c end c function nofx(x,array,npts) c c implicit integer(i-n) c implicit real(a-h,o-z) c c implicit double precision(a-h,o-z) c c c c c function nofx c c c c purpose c c given a value x and an array of values, find the index c c corresponding to the array element closest to x c c c c usage c c n = nofx(x,array,npts) c c c c parameters c c x - a given value c c array - array of values, assumed to be stored in order of c c increasing value c c npts - number of elements in array c c c c subroutines and function subprograms required c c none c c c c written 8/11/81 by j.m. tranquada c c c dimension array(npts) c imin = 1 c imax = npts c inc = ( imax - imin ) / 2 c 10 continue c it = imin + inc c xit = array(it) c if ( x .lt. xit ) then c imax = it c else if ( x .gt. xit ) then c imin = it c else c nofx = it c return c endif c inc = ( imax - imin ) / 2 c if ( inc .gt. 0 ) go to 10 c xave = ( array(imin) + array(imin+1) ) / 2. c if ( x .lt. xave ) then c nofx = imin c else c nofx = imin + 1 c endif c return c c end function nofx c end c subroutine openfl(iunit, file, status, iexist, ierr) c c open a file, c if unit <= 0, the first unused unit number greater than 7 will c be assigned. c if status = 'old', the existence of the file is checked. c if the file does not exist iexist is set to -1 c if the file does exist, iexist = iunit. c if any errors are encountered, ierr is set to -1. c c note: iunit, iexist, and ierr may be overwritten by this routine character*(*) file, status, stat*10 integer iunit, iexist, ierr, iulow logical opend, exist data iulow /7/ c c make sure there is a unit number ierr = -3 iexist = 0 if (iunit.le.0) then iunit = iulow 10 continue inquire(unit=iunit, opened = opend) if (opend) then if (iunit.gt.20) return iunit = iunit + 1 go to 10 endif end if c c if status = 'old', check that the file name exists ierr = -2 stat = status call smcase(stat,'a') if (stat.eq.'old') then iexist = -1 inquire(file=file, exist = exist) if (.not.exist) return iexist = iunit end if c c open the file ierr = -1 open(unit=iunit, file=file, status=status, err=100) ierr = 0 100 continue return c end subroutine openfl end subroutine pijump (ph, old) c c removes jumps of 2*pi in phases c ph = current value of phase (may be modified on output, but c only by multiples of 2*pi) c old = previous value of phase ccc implicit double precision (a-h, o-z) parameter (pi = 3.14159 26535 89793 23846 26433) parameter (twopi = 2 * pi) dimension xph(3) xph(1) = ph - old jump = (abs(xph(1))+ pi) / twopi xph(2) = xph(1) - jump*twopi xph(3) = xph(1) + jump*twopi xphmin = min (abs(xph(1)), abs(xph(2)), abs(xph(3))) do 10 i = 1, 3 if (abs (xphmin - abs(xph(i))) .le. 0.01) isave = i 10 continue ph = old + xph(isave) return c end subroutine pijump end c%%% integer function nbrstr(string) c%%% c c%%% c find a number in a string c%%% c given a string that is known to begin with a digit or sign. c%%% c return the position of the end of the number. c%%% c nbrstr : position of end of number c%%% integer istrln, i, ilen, iback c%%% character*(*) string c%%% character*1 digit*10, plus, minus, d, e, decml, s, sp c%%% logical lexp, ldecml c%%% data digit /'1234567890'/ c%%% data minus,plus,d,e,decml /'-','+','d','e','.'/ c%%% c------ c%%% ldecml = .false. c%%% lexp = .false. c%%% ilen = istrln(string) c%%% nbrstr = ilen c%%% if (ilen.gt.1) then c%%% iback = 1 c%%% c find end of number : digits are always ok. c%%% c stop at second d, e, decml, or sign that's not preceded by (d,e) c%%% do 200 i = 2, ilen c%%% sp = string(i-1:i-1) c%%% s = string(i:i) c%%% if (index(digit,s).eq.0) then c%%% if ( ( (s.ne.plus).and.(s.ne.minus).and.(s.ne.d) c%%% $ .and.(s.ne.e).and.(s.ne.decml)) c%%% $ .or.( lexp.and.((s.eq.d).or.(s.eq.e))) c%%% $ .or.( ldecml.and.(s.eq.decml)) c%%% $ .or.( ( (s.eq.plus).or.(s.eq.minus)).and. c%%% $ (sp.ne.d).and.(sp.ne.e)) ) go to 210 c%%% lexp = lexp.or.(s.eq.d).or.(s.eq.e) c%%% ldecml = ldecml.or.(s.eq.decml) c%%% end if c%%% 200 continue c%%% iback = 0 c%%% 210 continue c%%% nbrstr = i - 1 - iback c%%% end if c%%% return c%%% c end function nbrstr c%%% end c%%% subroutine unblnk (string) c%%% c c%%% c remove blanks from a string c%%% integer i, ilen, j c%%% character*(*) string, str*256 c%%% ilen = min(256, max(1, istrln(string))) c%%% j = 0 c%%% str = ' ' c%%% do 10 i = 1, ilen c%%% if (string(i:i).ne.' ') then c%%% j = j+1 c%%% str(j:j) = string(i:i) c%%% end if c%%% 10 continue c%%% string = ' ' c%%% string = str(1:j) c%%% return c%%% c end subroutine unblnk c%%% end subroutine sgetvl(word, xreal, ierr) c reads real number "xreal" from input character string "word". c returns xreal = 0. and (ierr.ne.0) if word cannot be a number character*(*) word real xreal integer ierr logical number, isnum external isnum c xreal = 0. ierr = -999 number = isnum(word) if (number) then ierr = 0 read ( word, '(bn,f30.0)', iostat = ierr) xreal end if return c end subroutine sgetvl end c---------------------------------------------------------------- c%%% subroutine uneql (string) c%%% c replace commas or equal signs with blanks c%%% integer i, ilen c%%% character*(*) string, equal*1 c%%% data equal / '=' / c%%% ilen = max(1, istrln(string)) c%%% do 10 i=1,ilen c%%% if (string(i:i).eq.equal) string(i:i) = ' ' c%%% 10 continue c%%% return c%%% c end subroutine uneql c%%% end c---------------------------------------------------------------- integer function nxtunt(iunit) c this function returns the value of the next unopened unit number equal c to or larger than iunit. it will return neither unit numbers 0, 5, c or 6 nor a negative unit number logical open iun = iunit if (iun.le.0) iun = 1 10 continue if ((iun.eq.5).or.(iun.eq.6)) then iun = 7 goto 10 endif inquire (unit=iun, opened=open) if (open) then iun = iun + 1 goto 10 endif nxtunt = iun return c end integer function nxtunt end logical function isdat(string) c tests if string contains numerical data c returns true if the first (up to eight) words in string can c all be numbers. requires at least two words, and tests only c the first eight columns integer nwords, mwords, i parameter (mwords = 8) character*(30) string*(*), words(mwords), line*(256) logical isnum external isnum c isdat = .false. do 10 i = 1, mwords words(i) = 'no' 10 continue c nwords = mwords line = string call triml(line) call untab(line) call bwords(line, nwords, words) if (nwords.ge.2) then isdat = .true. do 50 i = 1, nwords if (.not. ( isnum( words(i) ) ) ) isdat = .false. 50 continue end if c return end subroutine testrf(flnam, irecl, flform, ier) c c test whether a data file can be interpreted as uwxafs binary c data file or ascii column data file. c c uwxafs binary files use direct access binary files c with word size irecl, which is a machine dependent parameter c c ier = -1 : file not found c ier = -2 : broken uwxafs file? c ier = -3 : not uwxafs file, but can't find data. c ier = -4 : looks like ascii, saw line of minus signs, c but 2nd following line doesn't have data c c copright 1994 university of washington matt newville c ----------------------------------------------------- integer i, irecl, iunit parameter (mwords = 3) character*(*) flnam, flform, line*80 integer*2 indx(4) logical exist, opend, isdat, prevdt, lisdat external isdat c ----------------------------------------------------- flform = 'none' ier = -1 iunit = 7 10 continue inquire(unit=iunit, opened = opend) if (opend) then if (iunit.gt.20) return iunit = iunit + 1 go to 10 endif inquire(file = flnam, exist = exist) if (.not.exist) return ier = -2 c ----------------------------------------------------- c try reading file as a uwxafs binary file c which have patriotic magic numbers embedded in them indx(3) = 0 indx(4) = 0 open(iunit, file= flnam, recl = irecl, err = 20, $ access = 'direct', status = 'old' ) 20 continue read(iunit, rec=1, err = 25) (indx(i), i=1,4) 25 continue if ((indx(3).eq.1776).and.(indx(4).eq.704)) then flform = 'uwxafs' ier = 0 go to 900 end if c ----------------------------------------------------- c try to read file as ascii data file close(iunit) open(iunit, file=flnam, status='old') prevdt = .false. 200 continue ier = -3 read(iunit, '(a)', end = 900, err = 900) line call triml (line) if (line(3:6) .eq. '----') then ier = -4 read(iunit, '(a)', end = 900, err = 900) line read(iunit, '(a)', end = 900, err = 900) line lisdat = isdat(line) if (lisdat ) then flform = 'ascii' ier = 0 end if go to 900 end if c if two lines in a row have all words being numbers, it is data lisdat = isdat(line) if (lisdat.and.prevdt) then flform = 'ascii' ier = 0 go to 900 end if prevdt = lisdat go to 200 c--------------------- 900 continue close(iunit) return c end subroutine testrf end subroutine mcfluo (cent, edge, mel, elem, perc, $ latwei, mutf, mub, muc) c This subroutine will calculate the three total cross sections c mutf, mub and muc needed for calculating self absorption corrections c to fluorescence data. c input variables: c cent - central atom symbol c edge - edge of interest symbol: k, l1, l2, l3 c mel - length of array with atoms symbols and relative amounts c elem - array of characters of length mel c perc - array of reals of length mel c latwei - logical flag = .true. if relative amounts of elements given c in units of atomic weigth. The default is at. % c output variables: c mutf - Total cross section at the fluorescent energy c mub - background absorption xsection c muc - central atom xsection of edge of interest parameter(zero=0.e0, imel=20) parameter( weimin = 1.e-5, alpmin = 0.001, ediff = 0.01 ) dimension perc(imel), abxsec(imel), blxsec(imel), $ flxsec(imel), xsec(10), energy(9), $ weigat(imel) dimension izel(imel) character*2 elem(imel), test, cent, edge character*1 unit character*78 errom1 real mutf, mub, muc, percen, ecent, eabove, ebelow integer incen, izcent logical latwei, erf c Initialize some variables used in this subroutine print*, ' Mcfluo calculating corrections ...' test ='ab' do 5 i=1, imel abxsec(i)=zero blxsec(i)=zero flxsec(i)=zero weigat(i)=zero izel(i)=0 5 continue incen=0 izcent=is2z(cent) if (izcent.eq.0) goto 80 do 10 i= 1,mel call case (test, elem(i)) izel(i)=is2z(elem(i)) if (izel(i).eq.izcent) then incen=i percen=perc(i) endif if (izel(i).eq.0) goto 81 10 continue if (incen.eq.0) goto 79 call smcase(edge, test) if ((edge(1:2) .ne. 'k ').and.(edge(1:2) .ne. 'l1') $ .and.(edge(1:2).ne. 'l2').and.(edge(1:2).ne. 'l3')) $ goto 78 ecent = s2e(cent, edge) eabove = ecent + ediff ebelow = ecent - ediff c initialize mucal variables erf = .false. do 11 i=1,9 energy(i)=zero xsec(i) =zero 11 continue xsec(10)=zero do 30 i=1, mel unit = 'b' c above call mucal(eabove, elem(i), izel(i), unit, xsec, energy, $ erf, ier) abxsec(i)=xsec(4) weigat(i)=xsec(7) c below call mucal(ebelow, elem(i), izel(i), unit, xsec, energy, $ erf, ier) blxsec(i)=xsec(4) if ((i.eq.incen).and.(edge(1:1).eq.'k')) alphfl = energy(6) if ((i.eq.incen).and.(edge(1:1).eq.'l')) alphfl = energy(8) 30 continue if (alphfl.le.alpmin) go to 90 c print*, alphfl do 40 i=1, mel unit = 'b' call mucal (alphfl, elem(i), izel(i), unit, xsec, energy, $ erf, ier) flxsec(i)=xsec(4) 40 continue c If percent is given in atomic weight units divide by atomic weight if (latwei) then do 41 i=1, mel weigat(i) = max( weimin, weigat(i) ) perc(i)=perc(i)/(weigat(i)) 41 continue endif c Calculate mub, the total "background" xsection. c All atoms at eabove except central atom at ebelow c calculate mutf, total xsection at fluorescence energy mub = zero mutf = zero do 50 i=1, mel mub = mub + abxsec(i)*perc(i) mutf = mutf + flxsec(i)*perc(i) 50 continue c correct mub for having taken central atom above edge mub= mub - abxsec(incen)*perc(incen) + blxsec(incen)*perc(incen) c calculate muc, edge of interest xsection muc = (abxsec(incen)-blxsec(incen))*perc(incen) 70 continue return 78 continue errom1 = '** The edge of interest has to be either '// $ ' K, L1, L2 or L3 **' goto 82 79 continue errom1 = '** Your central atom must be one of the '// $ ' elements in the atom list **' goto 82 80 continue errom1 = '** You did not specify a central atom or'// $ ' its symbol is unrecognized **' goto 82 81 continue errom1 = '** The atom symbol "'//elem(i)//'" is not an'// $ ' element **' 82 continue mi = istrln(errom1) call messag (errom1(1:mi)) call messag ('*** Fluo is stopping ***') stop 90 continue call messag ('** no fluorescence energy found **') call messag ('*** Fluo is stopping ***') stop c end subroutine mcfluo end c**************************************************************** subroutine mucal(en,sym,iz,unit,xsec,energy,erf,ier) c--------------------------------------------------------------------- c standardized, made double precision, stylized and adapted for atom c by b ravel, september 1992 c--------------------------------------------------------------------- c this is a routine to calculate x-sections using data from c the may 1969 edition of mcmaster. c written by pathikrit bandyopadhyay at the university of c notre dame. c c this program has data for all the elements from z=1 to 94 c with the following exceptions: c 84 po \ c 85 at | c 87 fr | c 88 ra > mcmaster does not publish data for these elements. c 89 ac | c 91 pa | c 93 np / c c****** please correct data errors as they are found. thanx. ****** c c input: c en: energy at which to calculate the x-section, en<0 means c to skip calculation of absorption c ** symb: name of material c unit: units to be used. 'c' for cm**2/gm, 'b' for barns/atom c ** iz: z number of element c erf: logical print flag (true = print run-time error messages) c **either or both supplied on input, both returned on output c c output: c xsec(1) = photoelectric x-section c xsec(2) = coherent x-section c xsec(3) = incoherent x-section c xsec(4) = total x-section c xsec(5) = conversion factor c xsec(6) = absorption coefficient c xsec(7) = atomic weight c xsec(8) = density c xsec(9) = l2-edge jump c xsec(10) = l3-edge jump c energy(1) = k-edge energy c energy(2) = l1-edge energy c energy(3) = l2-edge energy c energy(4) = l3-edge energy c energy(5) = m-edge energy c energy(6) = k-alpha1 c energy(7) = k-beta1 c energy(8) = l-alpha1 c energy(9) = l-beta1 c ier = error code c c error codes: c ier=1: energy input is zero c ier=2: name does not match z c ier=3: no documentation for given element (z<94) c ier=4: no documentation for given element (z>94) c ier=5: l-edge calculation may be wrong for z<30 as mcmaster c uses l1 only. c ier=6: energy at the middle of edge c ier=7: no name or z supplied c c-------------------------------------------------------------------- implicit integer(i-n) implicit real(a-h,o-z) c implicit double precision(a-h,o-z) parameter(nelem=94, eps=.001) c double precision ka(nelem),kb(nelem),la(nelem),lb(nelem) c double precision l2(nelem),l3(nelem),lj3(nelem),lj2 real ka(nelem),kb(nelem),la(nelem),lb(nelem) real l2(nelem),l3(nelem),lj3(nelem),lj2 logical erf character*1 unit character*2 name(nelem),symb,sym,test dimension ek(nelem),el(nelem),em(nelem),den(nelem), $ atwt(nelem),cf(nelem) dimension ak(0:3,nelem),al(0:3,nelem),am(0:3,nelem), $ an(0:3,nelem) dimension cih(0:3,nelem),coh(0:3,nelem) dimension xsec(10),energy(9) c###the following is the symbol, edge energies, mcmaster c coefficients, flouescent energies, density, conversion c constant, and atomic weight of each supported element c up to plutonium data name( 1), ek( 1),el( 1)/ $ 'h ', 0.140000e-01, 0.000000e+00/ data ak(0, 1),ak(1, 1),ak(2, 1),ak(3, 1)/ $ 0.244964e+01, -0.334953e+01, -0.471370e-01, 0.709962e-02/ data al(0, 1),al(1, 1),al(2, 1),al(3, 1)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data am(0, 1),am(1, 1),am(2, 1),am(3, 1)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data an(0, 1),an(1, 1),an(2, 1),an(3, 1)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den( 1),cf( 1),atwt(1)/1.008, 0.167400e+01,1.008/ data coh(0, 1),coh(1, 1),coh(2, 1),coh(3, 1)/ $ -0.119075e+00, -0.937086e+00, -0.200538e+00, 0.106587e-01/ data cih(0, 1),cih(1, 1),cih(2, 1),cih(3, 1)/ $ -0.215772e+01, 0.132685e+01, -0.305620e+00, 0.185025e-01/ data name( 2), ek( 2),el( 2)/ $ 'he', 0.250000e-01, 0.000000e+00/ data ak(0, 2),ak(1, 2),ak(2, 2),ak(3, 2)/ $ 0.606488e+01, -0.329055e+01, -0.107256e+00, 0.144465e-01/ data al(0, 2),al(1, 2),al(2, 2),al(3, 2)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data am(0, 2),am(1, 2),am(2, 2),am(3, 2)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data an(0, 2),an(1, 2),an(2, 2),an(3, 2)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den( 2),cf( 2),atwt(2)/1.785e-04, 0.664700e+01,4.003/ data coh(0, 2),coh(1, 2),coh(2, 2),coh(3, 2)/ $ 0.104768e+01, -0.851805e-01, -0.403527e+00, 0.269398e-01/ data cih(0, 2),cih(1, 2),cih(2, 2),cih(3, 2)/ $ -0.256357e+01, 0.202536e+01, -0.448710e+00, 0.279691e-01/ data name( 3), ek( 3),el( 3)/ $ 'li', 0.550000e-01, 0.000000e+00/ data ak(0, 3),ak(1, 3),ak(2, 3),ak(3, 3)/ $ 0.775370e+01, -0.281801e+01, -0.241378e+00, 0.262542e-01/ data al(0, 3),al(1, 3),al(2, 3),al(3, 3)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data am(0, 3),am(1, 3),am(2, 3),am(3, 3)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data an(0, 3),an(1, 3),an(2, 3),an(3, 3)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den( 3),cf( 3),atwt(3)/0.534, 0.115200e+02,6.940/ data coh(0, 3),coh(1, 3),coh(2, 3),coh(3, 3)/ $ 0.134366e+01, 0.181557e+00, -0.423981e+00, 0.266190e-01/ data cih(0, 3),cih(1, 3),cih(2, 3),cih(3, 3)/ $ -0.108740e+01, 0.103368e+01, -0.190377e+00, 0.779955e-02/ data name( 4), ek( 4),el( 4)/ $ 'be', 0.112000e+00, 0.000000e+00/ data ak(0, 4),ak(1, 4),ak(2, 4),ak(3, 4)/ $ 0.904511e+01, -0.283487e+01, -0.210021e+00, 0.229526e-01/ data al(0, 4),al(1, 4),al(2, 4),al(3, 4)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data am(0, 4),am(1, 4),am(2, 4),am(3, 4)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data an(0, 4),an(1, 4),an(2, 4),an(3, 4)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den( 4),cf( 4),atwt(4)/1.848, 0.149600e+02,9.012/ data coh(0, 4),coh(1, 4),coh(2, 4),coh(3, 4)/ $ 0.200860e+01, -0.461920e-01, -0.337018e+00, 0.186939e-01/ data cih(0, 4),cih(1, 4),cih(2, 4),cih(3, 4)/ $ -0.690079e+00, 0.946448e+00, -0.171142e+00, 0.651413e-02/ data name( 5), ek( 5),el( 5)/ $ 'b ', 0.188000e+00, 0.000000e+00/ data ak(0, 5),ak(1, 5),ak(2, 5),ak(3, 5)/ $ 0.995057e+01, -0.274174e+01, -0.215138e+00, 0.227845e-01/ data al(0, 5),al(1, 5),al(2, 5),al(3, 5)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data am(0, 5),am(1, 5),am(2, 5),am(3, 5)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data an(0, 5),an(1, 5),an(2, 5),an(3, 5)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den( 5),cf( 5),atwt(5)/2.340, 0.179500e+02,10.811/ data coh(0, 5),coh(1, 5),coh(2, 5),coh(3, 5)/ $ 0.262862e+01, -0.207916e+00, -0.286283e+00, 0.144966e-01/ data cih(0, 5),cih(1, 5),cih(2, 5),cih(3, 5)/ $ -0.791177e+00, 0.121611e+01, -0.239087e+00, 0.117686e-01/ data name( 6), ek( 6),el( 6)/ $ 'c ', 0.284000e+00, 0.000000e+00/ data ak(0, 6),ak(1, 6),ak(2, 6),ak(3, 6)/ $ 0.106879e+02, -0.271400e+01, -0.200530e+00, 0.207248e-01/ data al(0, 6),al(1, 6),al(2, 6),al(3, 6)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data am(0, 6),am(1, 6),am(2, 6),am(3, 6)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data an(0, 6),an(1, 6),an(2, 6),an(3, 6)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den( 6),cf( 6),atwt(6)/2.250, 0.199400e+02,12.01/ data coh(0, 6),coh(1, 6),coh(2, 6),coh(3, 6)/ $ 0.310861e+01, -0.260580e+00, -0.271974e+00, 0.135181e-01/ data cih(0, 6),cih(1, 6),cih(2, 6),cih(3, 6)/ $ -0.982878e+00, 0.146693e+01, -0.293743e+00, 0.156005e-01/ data name( 7), ek( 7),el( 7)/ $ 'n ', 0.402000e+00, 0.000000e+00/ data ak(0, 7),ak(1, 7),ak(2, 7),ak(3, 7)/ $ 0.112765e+02, -0.265400e+01, -0.200445e+00, 0.200765e-01/ data al(0, 7),al(1, 7),al(2, 7),al(3, 7)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den( 7),cf( 7),atwt(7)/1.250e-03, 0.232600e+02,14.008/ data coh(0, 7),coh(1, 7),coh(2, 7),coh(3, 7)/ $ 0.347760e+01, -0.215762e+00, -0.288874e+00, 0.115131e-01/ data cih(0, 7),cih(1, 7),cih(2, 7),cih(3, 7)/ $ -0.123693e+01, 0.174510e+01, -0.354660e+00, 0.198705e-01/ data name( 8), ek( 8),el( 8)/ $ 'o ', 0.537000e+00, 0.000000e+00/ data ak(0, 8),ak(1, 8),ak(2, 8),ak(3, 8)/ $ 0.117130e+02, -0.257229e+01, -0.205893e+00, 0.199244e-01/ data al(0, 8),al(1, 8),al(2, 8),al(3, 8)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den( 8),cf( 8),atwt(8)/1.429e-03, 0.265700e+02,16.0/ data coh(0, 8),coh(1, 8),coh(2, 8),coh(3, 8)/ $ 0.377239e+01, -0.148539e+00, -0.307124e+00, 0.167303e-01/ data cih(0, 8),cih(1, 8),cih(2, 8),cih(3, 8)/ $ -0.173679e+01, 0.217686e+01, -0.449050e+00, 0.264733e-01/ data name( 9), ek( 9),el( 9)/ $ 'f ', 0.686000e+00, 0.000000e+00/ data ak(0, 9),ak(1, 9),ak(2, 9),ak(3, 9)/ $ 0.120963e+02, -0.244148e+01, -0.234461e+00, 0.219954e-01/ data al(0, 9),al(1, 9),al(2, 9),al(3, 9)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den( 9),cf( 9),atwt(9)/1.108, 0.315500e+02,19.0/ data coh(0, 9),coh(1, 9),coh(2, 9),coh(3, 9)/ $ 0.400716e+01, -0.560908e-01, -0.332017e+00, 0.187934e-01/ data cih(0, 9),cih(1, 9),cih(2, 9),cih(3, 9)/ $ -0.187570e+01, 0.232016e+01, -0.475412e+00, 0.280680e-01/ data name(10), ek(10),el(10)/ $ 'ne', 0.867000e+00, 0.000000e+00/ data ak(0,10),ak(1,10),ak(2,10),ak(3,10)/ $ 0.124485e+02, -0.245819e+01, -0.212591e+00, 0.196489e-01/ data al(0,10),al(1,10),al(2,10),al(3,10)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den(10),cf(10),atwt(10)/9.000e-04, 0.335100e+02,20.183/ data coh(0,10),coh(1,10),coh(2,10),coh(3,10)/ $ 0.420151e+01, 0.416247e-01, -0.356754e+00, 0.207585e-01/ data cih(0,10),cih(1,10),cih(2,10),cih(3,10)/ $ -0.175510e+01, 0.224226e+01, -0.447640e+00, 0.255801e-01/ data name(11), ek(11),el(11)/ $ 'na', 0.107200e+01, 0.000000e+00/ data ak(0,11),ak(1,11),ak(2,11),ak(3,11)/ $ 0.126777e+02, -0.224521e+01, -0.274873e+00, 0.250270e-01/ data al(0,11),al(1,11),al(2,11),al(3,11)/ $ 0.102355e+02, -0.255905e+01, -0.119524e+00, 0.000000e+00/ data den(11),cf(11),atwt(11)/0.970, 0.381900e+02,22.997/ data coh(0,11),coh(1,11),coh(2,11),coh(3,11)/ $ 0.426374e+01, 0.134662e+00, -0.370080e+00, 0.214467e-01/ data cih(0,11),cih(1,11),cih(2,11),cih(3,11)/ $ -0.967717e+00, 0.161794e+01, -0.287191e+00, 0.131526e-01/ data ka(11),kb(11),la(11),lb(11)/ $ 0.104100e+01, 0.106700e+01, 0.000000e+00, 0.000000e+00/ data name(12), ek(12),el(12)/ $ 'mg', 0.130500e+01, 0.630000e-01/ data ak(0,12),ak(1,12),ak(2,12),ak(3,12)/ $ 0.128793e+02, -0.212574e+01, -0.299392e+00, 0.267643e-01/ data al(0,12),al(1,12),al(2,12),al(3,12)/ $ 0.105973e+02, -0.289818e+01, 0.234506e+00, 0.000000e+00/ data den(12),cf(12),atwt(12)/1.740, 0.403800e+02,24.320/ data coh(0,12),coh(1,12),coh(2,12),coh(3,12)/ $ 0.439404e+01, 0.137858e+00, -0.359540e+00, 0.202380e-01/ data cih(0,12),cih(1,12),cih(2,12),cih(3,12)/ $ -0.571611e+00, 0.135498e+01, -0.222491e+00, 0.830141e-02/ data ka(12),kb(12),la(12),lb(12)/ $ 0.125400e+01, 0.129700e+01, 0.000000e+00, 0.000000e+00/ data name(13), ek(13),el(13)/ $ 'al', 0.156000e+01, 0.870000e-01/ data ak(0,13),ak(1,13),ak(2,13),ak(3,13)/ $ 0.131738e+02, -0.218203e+01, -0.258960e+00, 0.222840e-01/ data al(0,13),al(1,13),al(2,13),al(3,13)/ $ 0.108711e+02, -0.277860e+01, 0.175853e+00, 0.000000e+00/ data den(13),cf(13),atwt(13)/2.720, 0.447800e+02,26.970/ data coh(0,13),coh(1,13),coh(2,13),coh(3,13)/ $ 0.415995e+01, 0.140549e+00, -0.352441e+00, 0.193692e-01/ data cih(0,13),cih(1,13),cih(2,13),cih(3,13)/ $ -0.439322e+00, 0.130867e+01, -0.211648e+00, 0.754210e-02/ data ka(13),kb(13),la(13),lb(13)/ $ 0.148700e+01, 0.155300e+01, 0.000000e+00, 0.000000e+00/ data name(14), ek(14),el(14)/ $ 'si', 0.183900e+01, 0.118000e+00/ data ak(0,14),ak(1,14),ak(2,14),ak(3,14)/ $ 0.132682e+02, -0.198174e+01, -0.316950e+00, 0.273928e-01/ data al(0,14),al(1,14),al(2,14),al(3,14)/ $ 0.112237e+02, -0.273694e+01, 0.127557e+00, 0.000000e+00/ data den(14),cf(14),atwt(14)/2.330, 0.466300e+02,28.086/ data coh(0,14),coh(1,14),coh(2,14),coh(3,14)/ $ 0.464678e+01, 0.162780e+00, -0.358563e+00, 0.196926e-01/ data cih(0,14),cih(1,14),cih(2,14),cih(3,14)/ $ -0.414971e+00, 0.134868e+01, -0.222315e+00, 0.841959e-02/ data ka(14),kb(14),la(14),lb(14)/ $ 0.174000e+01, 0.183200e+01, 0.000000e+00, 0.000000e+00/ data name(15), ek(15),el(15)/ $ 'p ', 0.214900e+01, 0.153000e+00/ data ak(0,15),ak(1,15),ak(2,15),ak(3,15)/ $ 0.133735e+02, -0.186342e+01, -0.339440e+00, 0.288858e-01/ data al(0,15),al(1,15),al(2,15),al(3,15)/ $ 0.115508e+02, -0.292200e+01, 0.254262e+00, 0.000000e+00/ data den(15),cf(15),atwt(15)/1.820, 0.514300e+02,30.975/ data coh(0,15),coh(1,15),coh(2,15),coh(3,15)/ $ 0.478525e+01, 0.168708e+00, -0.360383e+00, 0.197155e-01/ data cih(0,15),cih(1,15),cih(2,15),cih(3,15)/ $ -0.476903e+00, 0.146032e+01, -0.251331e+00, 0.107202e-01/ data ka(15),kb(15),la(15),lb(15)/ $ 0.201500e+01, 0.213600e+01, 0.000000e+00, 0.000000e+00/ data name(16), ek(16),el(16)/ $ 's ', 0.247200e+01, 0.193000e+00/ data ak(0,16),ak(1,16),ak(2,16),ak(3,16)/ $ 0.137394e+02, -0.204786e+01, -0.273259e+00, 0.229976e-01/ data al(0,16),al(1,16),al(2,16),al(3,16)/ $ 0.118181e+02, -0.264618e+01, -0.968049e-01, 0.000000e+00/ data am(0,16),am(1,16),am(2,16),am(3,16)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data an(0,16),an(1,16),an(2,16),an(3,16)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den(16),cf(16)/2.000, 0.532400e+02/ data em(16),atwt(16)/ 0.170000e-01,32.066/ data coh(0,16),coh(1,16),coh(2,16),coh(3,16)/ $ 0.492707e+01, 0.165746e+00, -0.359424e+00, 0.195505e-01/ data cih(0,16),cih(1,16),cih(2,16),cih(3,16)/ $ -0.656419e+00, 0.165408e+01, -0.298623e+00, 0.142979e-01/ data ka(16),kb(16),la(16),lb(16)/ $ 0.230800e+01, 0.246400e+01, 0.000000e+00, 0.000000e+00/ data name(17), ek(17),el(17)/ $ 'cl', 0.282200e+01, 0.238000e+00/ data ak(0,17),ak(1,17),ak(2,17),ak(3,17)/ $ 0.136188e+02, -0.171937e+01, -0.354154e+00, 0.290841e-01/ data al(0,17),al(1,17),al(2,17),al(3,17)/ $ 0.120031e+02, -0.241694e+01, -0.240897e+00, 0.000000e+00/ data am(0,17),am(1,17),am(2,17),am(3,17)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data an(0,17),an(1,17),an(2,17),an(3,17)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den(17),cf(17)/1.560, 0.588700e+02/ data em(17),atwt(17)/ 0.170000e-01,35.457/ data coh(0,17),coh(1,17),coh(2,17),coh(3,17)/ $ 0.507222e+01, 0.149127e+00, -0.352858e+00, 0.189439e-01/ data cih(0,17),cih(1,17),cih(2,17),cih(3,17)/ $ -0.718627e+00, 0.174294e+01, -0.319429e+00, 0.158429e-01/ data ka(17),kb(17),la(17),lb(17)/ $ 0.262200e+01, 0.281500e+01, 0.000000e+00, 0.000000e+00/ data name(18), ek(18),el(18)/ $ 'ar', 0.320200e+01, 0.287000e+00/ data ak(0,18),ak(1,18),ak(2,18),ak(3,18)/ $ 0.139491e+02, -0.182276e+01, -0.328827e+00, 0.274382e-01/ data al(0,18),al(1,18),al(2,18),al(3,18)/ $ 0.122960e+02, -0.263279e+01, -0.736600e-01, 0.000000e+00/ data am(0,18),am(1,18),am(2,18),am(3,18)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data an(0,18),an(1,18),an(2,18),an(3,18)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den(18),cf(18)/1.784e-03, 0.663200e+02/ data em(18),atwt(18)/ 0.270000e-01,39.944/ data coh(0,18),coh(1,18),coh(2,18),coh(3,18)/ $ 0.521079e+01, 0.135618e+00, -0.347214e+00, 0.184333e-01/ data cih(0,18),cih(1,18),cih(2,18),cih(3,18)/ $ -0.682105e+00, 0.174279e+01, -0.317646e+00, 0.156467e-01/ data ka(18),kb(18),la(18),lb(18)/ $ 0.295700e+01, 0.319200e+01, 0.000000e+00, 0.000000e+00/ data name(19), ek(19),el(19)/ $ 'k ', 0.360700e+01, 0.341000e+00/ data ak(0,19),ak(1,19),ak(2,19),ak(3,19)/ $ 0.137976e+02, -0.154015e+01, -0.394528e+00, 0.323561e-01/ data al(0,19),al(1,19),al(2,19),al(3,19)/ $ 0.124878e+02, -0.253656e+01, -0.104892e+00, 0.000000e+00/ data am(0,19),am(1,19),am(2,19),am(3,19)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data an(0,19),an(1,19),an(2,19),an(3,19)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den(19),cf(19)/0.862, 0.649300e+02/ data em(19),atwt(19)/ 0.340000e-01,39.102/ data coh(0,19),coh(1,19),coh(2,19),coh(3,19)/ $ 0.525587e+01, 0.188040e+00, -0.359623e+00, 0.193085e-01/ data cih(0,19),cih(1,19),cih(2,19),cih(3,19)/ $ -0.344007e+00, 0.149236e+01, -0.254135e+00, 0.107684e-01/ data ka(19),kb(19),la(19),lb(19)/ $ 0.331300e+01, 0.358900e+01, 0.000000e+00, 0.000000e+00/ data name(20), ek(20),el(20)/ $ 'ca', 0.403800e+01, 0.400000e+00/ data ak(0,20),ak(1,20),ak(2,20),ak(3,20)/ $ 0.142950e+02, -0.188644e+01, -0.283647e+00, 0.226263e-01/ data al(0,20),al(1,20),al(2,20),al(3,20)/ $ 0.127044e+02, -0.255011e+01, -0.943195e-01, 0.000000e+00/ data am(0,20),am(1,20),am(2,20),am(3,20)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data an(0,20),an(1,20),an(2,20),an(3,20)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den(20),cf(20)/1.550, 0.665500e+02/ data em(20),atwt(20)/ 0.440000e-01,40.080/ data coh(0,20),coh(1,20),coh(2,20),coh(3,20)/ $ 0.532375e+01, 0.206685e+00, -0.361664e+00, 0.193328e-01/ data cih(0,20),cih(1,20),cih(2,20),cih(3,20)/ $ -0.982420e-01, 0.132829e+01, -0.213747e+00, 0.773065e-02/ data ka(20),kb(20),la(20),lb(20)/ $ 0.369100e+01, 0.401200e+01, 0.000000e+00, 0.000000e+00/ data name(21), ek(21),el(21)/ $ 'sc', 0.449300e+01, 0.463000e+00/ data ak(0,21),ak(1,21),ak(2,21),ak(3,21)/ $ 0.139664e+02, -0.140872e+01, -0.414365e+00, 0.334355e-01/ data al(0,21),al(1,21),al(2,21),al(3,21)/ $ 0.128949e+02, -0.240609e+01, -0.177791e+00, 0.000000e+00/ data am(0,21),am(1,21),am(2,21),am(3,21)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data an(0,21),an(1,21),an(2,21),an(3,21)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den(21), cf(21)/ 0.299200e+01, 0.746500e+02/ data em(21),atwt(21)/ 0.540000e-01,44.960/ data coh(0,21),coh(1,21),coh(2,21),coh(3,21)/ $ 0.543942e+01, 0.200174e+00, -0.359064e+00, 0.191027e-01/ data cih(0,21),cih(1,21),cih(2,21),cih(3,21)/ $ -0.159831e+00, 0.139055e+01, -0.225849e+00, 0.851954e-02/ data ka(21),kb(21),la(21),lb(21)/ $ 0.409000e+01, 0.446000e+01, 0.000000e+00, 0.000000e+00/ data name(22), ek(22),el(22)/ $ 'ti', 0.496500e+01, 0.531000e+00/ data ak(0,22),ak(1,22),ak(2,22),ak(3,22)/ $ 0.143506e+02, -0.166322e+01, -0.331539e+00, 0.262065e-01/ data al(0,22),al(1,22),al(2,22),al(3,22)/ $ 0.131075e+02, -0.253576e+01, -0.957177e-01, 0.000000e+00/ data am(0,22),am(1,22),am(2,22),am(3,22)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data an(0,22),an(1,22),an(2,22),an(3,22)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den(22), cf(22)/ 0.454000e+01, 0.795300e+02/ data em(22),atwt(22)/ 0.590000e-01,47.9/ data coh(0,22),coh(1,22),coh(2,22),coh(3,22)/ $ 0.555039e+01, 0.197697e+00, -0.357694e+00, 0.189866e-01/ data cih(0,22),cih(1,22),cih(2,22),cih(3,22)/ $ -0.230573e+00, 0.145848e+01, -0.239160e+00, 0.938528e-02/ data ka(22),kb(22),la(22),lb(22)/ $ 0.451000e+01, 0.493100e+01, 0.000000e+00, 0.000000e+00/ data name(23), ek(23),el(23)/ $ 'v ', 0.546500e+01, 0.604000e+00/ data ak(0,23),ak(1,23),ak(2,23),ak(3,23)/ $ 0.147601e+02, -0.188867e+01, -0.271861e+00, 0.215792e-01/ data al(0,23),al(1,23),al(2,23),al(3,23)/ $ 0.132514e+02, -0.249765e+01, -0.106383e+00, 0.000000e+00/ data am(0,23),am(1,23),am(2,23),am(3,23)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data an(0,23),an(1,23),an(2,23),an(3,23)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den(23), cf(23)/ 0.611000e+01, 0.845900e+02/ data em(23),atwt(23)/ 0.670000e-01,50.942/ data coh(0,23),coh(1,23),coh(2,23),coh(3,23)/ $ 0.565514e+01, 0.199533e+00, -0.357487e+00, 0.189691e-01/ data cih(0,23),cih(1,23),cih(2,23),cih(3,23)/ $ -0.308103e+00, 0.152879e+01, -0.252768e+00, 0.102571e-01/ data ka(23),kb(23),la(23),lb(23)/ $ 0.495200e+01, 0.542700e+01, 0.000000e+00, 0.000000e+00/ data name(24), ek(24),el(24)/ $ 'cr', 0.598900e+01, 0.682000e+00/ data ak(0,24),ak(1,24),ak(2,24),ak(3,24)/ $ 0.148019e+02, -0.182430e+01, -0.279116e+00, 0.217324e-01/ data al(0,24),al(1,24),al(2,24),al(3,24)/ $ 0.134236e+02, -0.251532e+01, -0.101999e+00, 0.000000e+00/ data am(0,24),am(1,24),am(2,24),am(3,24)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data an(0,24),an(1,24),an(2,24),an(3,24)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den(24), cf(24)/ 0.719000e+01, 0.863400e+02/ data em(24),atwt(24)/ 0.740000e-01,51.996/ data coh(0,24),coh(1,24),coh(2,24),coh(3,24)/ $ 0.577399e+01, 0.203858e+00, -0.359699e+00, 0.192225e-01/ data cih(0,24),cih(1,24),cih(2,24),cih(3,24)/ $ -0.387641e+00, 0.159727e+01, -0.266240e+00, 0.111523e-01/ data ka(24),kb(24),la(24),lb(24)/ $ 0.541400e+01, 0.594600e+01, 0.000000e+00, 0.000000e+00/ data name(25), ek(25),el(25)/ $ 'mn', 0.654000e+01, 0.754000e+00/ data ak(0,25),ak(1,25),ak(2,25),ak(3,25)/ $ 0.148965e+02, -0.179872e+01, -0.283664e+00, 0.222095e-01/ data al(0,25),al(1,25),al(2,25),al(3,25)/ $ 0.135761e+02, -0.249761e+01, -0.105493e+00, 0.000000e+00/ data am(0,25),am(1,25),am(2,25),am(3,25)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data an(0,25),an(1,25),an(2,25),an(3,25)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den(25), cf(25)/ 0.742000e+01, 0.912200e+02/ data em(25),atwt(25)/ 0.840000e-01,54.940/ data coh(0,25),coh(1,25),coh(2,25),coh(3,25)/ $ 0.584604e+01, 0.213814e+00, -0.359718e+00, 0.191459e-01/ data cih(0,25),cih(1,25),cih(2,25),cih(3,25)/ $ -0.247059e+00, 0.149722e+01, -0.238781e+00, 0.893208e-02/ data ka(25),kb(25),la(25),lb(25)/ $ 0.589800e+01, 0.649000e+01, 0.000000e+00, 0.000000e+00/ data name(26), ek(26),el(26)/ $ 'fe', 0.711200e+01, 0.842000e+00/ data ak(0,26),ak(1,26),ak(2,26),ak(3,26)/ $ 0.143456e+02, -0.123491e+01, -0.423491e+00, 0.321661e-01/ data al(0,26),al(1,26),al(2,26),al(3,26)/ $ 0.136696e+02, -0.239195e+01, -0.137680e+00, 0.000000e+00/ data den(26),cf(26)/7.860,9.274e+01/ data em(26),atwt(26)/0.094,55.850/ data coh(0,26),coh(1,26),coh(2,26),coh(3,26)/ $ 0.593292e+01, 0.225038e+00, -0.361748e+00, 0.193024e-01/ data cih(0,26),cih(1,26),cih(2,26),cih(3,26)/ $ -0.342379e+00, 0.157245e+01, -0.253198e+00, 0.985822e-02/ data ka(26),kb(26),la(26),lb(26)/ $ 0.640300e+01, 0.705700e+01, 0.000000e+00, 0.000000e+00/ data name(27), ek(27),el(27)/ $ 'co', 0.770900e+01, 0.929000e+00/ data ak(0,27),ak(1,27),ak(2,27),ak(3,27)/ $ 0.147047e+02, -0.138933e+01, -0.386631e+00, 0.303286e-01/ data al(0,27),al(1,27),al(2,27),al(3,27)/ $ 0.138699e+02, -0.250669e+01, -0.869945e-01, 0.000000e+00/ data am(0,27),am(1,27),am(2,27),am(3,27)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data an(0,27),an(1,27),an(2,27),an(3,27)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den(27), cf(27)/ 0.890000e+01, 0.978500e+02/ data em(27),atwt(27)/ 0.101000e+00,58.933/ data coh(0,27),coh(1,27),coh(2,27),coh(3,27)/ $ 0.601478e+01, 0.237959e+00, -0.364056e+00, 0.194754e-01/ data cih(0,27),cih(1,27),cih(2,27),cih(3,27)/ $ -0.428804e+00, 0.164129e+01, -0.266013e+00, 0.106512e-01/ data ka(27),kb(27),la(27),lb(27)/ $ 0.693000e+01, 0.764900e+01, 0.000000e+00, 0.000000e+00/ data name(28), ek(28),el(28)/ $ 'ni', 0.833300e+01, 0.101200e+01/ data ak(0,28),ak(1,28),ak(2,28),ak(3,28)/ $ 0.142388e+02, -0.967736e+00, -0.478070e+00, 0.366138e-01/ data al(0,28),al(1,28),al(2,28),al(3,28)/ $ 0.139848e+02, -0.248080e+01, -0.888115e-01, 0.000000e+00/ data am(0,28),am(1,28),am(2,28),am(3,28)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data an(0,28),an(1,28),an(2,28),an(3,28)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den(28), cf(28)/ 0.890000e+01, 0.974500e+02/ data em(28),atwt(28)/ 0.113000e+00,58.690/ data coh(0,28),coh(1,28),coh(2,28),coh(3,28)/ $ 0.609024e+01, 0.222770e+00, -0.366568e+00, 0.196586e-01/ data cih(0,28),cih(1,28),cih(2,28),cih(3,28)/ $ -0.504360e+00, 0.170040e+01, -0.276443e+00, 0.112628e-01/ data ka(28),kb(28),la(28),lb(28)/ $ 0.747700e+01, 0.826400e+01, 0.000000e+00, 0.000000e+00/ data name(29), ek(29),el(29)/ $ 'cu', 0.897900e+01, 0.110000e+01/ data ak(0,29),ak(1,29),ak(2,29),ak(3,29)/ $ 0.145808e+02, -0.118375e+01, -0.413850e+00, 0.312088e-01/ data al(0,29),al(1,29),al(2,29),al(3,29)/ $ 0.142439e+02, -0.258677e+01, -0.667398e-01, 0.000000e+00/ data am(0,29),am(1,29),am(2,29),am(3,29)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data an(0,29),an(1,29),an(2,29),an(3,29)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den(29), cf(29)/ 0.894000e+01, 0.105500e+03/ data em(29),atwt(29)/ 0.120000e+00,63.540/ data coh(0,29),coh(1,29),coh(2,29),coh(3,29)/ $ 0.617739e+01, 0.273123e+00, -0.372360e+00, 0.201638e-01/ data cih(0,29),cih(1,29),cih(2,29),cih(3,29)/ $ -0.570210e+00, 0.175042e+01, -0.2845550e+00, 0.116930e-01/ data ka(29),kb(29),la(29),lb(29)/ $ 0.804700e+01, 0.890400e+01, 0.000000e+00, 0.000000e+00/ data name(30), ek(30),el(30)/ $ 'zn', 0.965900e+01, 0.119600e+01/ data ak(0,30),ak(1,30),ak(2,30),ak(3,30)/ $ 0.144118e+02, -0.933083e+00, -0.477357e+00, 0.362829e-01/ data al(0,30),al(1,30),al(2,30),al(3,30)/ $ 0.143221e+02, -0.262384e+01, -0.264926e-01, 0.000000e+00/ data am(0,30),am(1,30),am(2,30),am(3,30)/ $ 0.120597e+02, -0.110258e+01, 0.000000e+00, 0.000000e+00/ data an(0,30),an(1,30),an(2,30),an(3,30)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den(30), cf(30)/ 0.714000e+01, 0.108600e+03/ data em(30),atwt(30)/ 0.139000e+00,65.380/ data coh(0,30),coh(1,30),coh(2,30),coh(3,30)/ $ 0.623402e+01, 0.284312e+00, -0.372143e+00, 0.200525e-01/ data cih(0,30),cih(1,30),cih(2,30),cih(3,30)/ $ -0.420535e+00, 0.163400e+01, -0.253646e+00, 0.927233e-02/ data ka(30),kb(30),la(30),lb(30)/ $ 0.863800e+01, 0.957100e+01, 0.100900e+01, 0.103200e+01/ data name(31), ek(31),el(31)/ $ 'ga', 0.103670e+02, 0.130200e+01/ data ak(0,31),ak(1,31),ak(2,31),ak(3,31)/ $ 0.136182e+02, -0.318459e+00, -0.611348e+00, 0.458138e-01/ data al(0,31),al(1,31),al(2,31),al(3,31)/ $ 0.144792e+02, -0.254469e+01, -0.757204e-01, 0.000000e+00/ data am(0,31),am(1,31),am(2,31),am(3,31)/ $ 0.122646e+02, -0.268965e+01, 0.000000e+00, 0.000000e+00/ c this is pretty interesting so far, isn't it! data an(0,31),an(1,31),an(2,31),an(3,31)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data coh(0,31),coh(1,31),coh(2,31),coh(3,31)/ $ 0.628298e+01, 0.291334e+00, -0.369391e+00, 0.197029e-01/ data cih(0,31),cih(1,31),cih(2,31),cih(3,31)/ $ -0.358218e+00, 0.160050e+01, -0.244908e+00, 0.861898e-02/ data den(31), cf(31)/ 0.590300e+01, 0.115800e+03/ data em(31),atwt(31)/ 0.158000e+00,69.720/ data ka(31),kb(31),la(31),lb(31)/ $ 0.925100e+01, 0.102630e+02, 0.109600e+01, 0.112200e+01/ data name(32), ek(32),el(32)/ $ 'ge', 0.111040e+02, 0.141400e+01/ data ak(0,32),ak(1,32),ak(2,32),ak(3,32)/ $ 0.139288e+02, -0.479613e+00, -0.572897e+00, 0.431277e-01/ data al(0,32),al(1,32),al(2,32),al(3,32)/ $ 0.146813e+02, -0.269285e+01, -0.208355e-01, 0.000000e+00/ data am(0,32),am(1,32),am(2,32),am(3,32)/ $ 0.124133e+02, -0.253085e+01, 0.000000e+00, 0.000000e+00/ data an(0,32),an(1,32),an(2,32),an(3,32)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data coh(0,32),coh(1,32),coh(2,32),coh(3,32)/ $ 0.633896e+01, 0.291512e+00, -0.365643e+00, 0.192896e-01/ data cih(0,32),cih(1,32),cih(2,32),cih(3,32)/ $ -0.334383e+00, 0.160327e+01, -0.245555e+00, 0.871239e-02/ data den(32), cf(32)/ 0.532300e+01, 0.120500e+03/ data em(32),atwt(32)/ 0.181000e+00,72.590/ data ka(32),kb(32),la(32),lb(32)/ $ 0.988500e+01, 0.109810e+02, 0.118600e+01, 0.121600e+01/ data name(33), ek(33),el(33)/ $ 'as', 0.118680e+02, 0.153000e+01/ data ak(0,33),ak(1,33),ak(2,33),ak(3,33)/ $ 0.134722e+02, -0.773513e-01, -0.660456e+00, 0.492177e-01/ data al(0,33),al(1,33),al(2,33),al(3,33)/ $ 0.146431e+02, -0.248397e+01, -0.796180e-01, 0.000000e+00/ data am(0,33),am(1,33),am(2,33),am(3,33)/ $ 0.125392e+02, -0.241380e+01, 0.000000e+00, 0.000000e+00/ data an(0,33),an(1,33),an(2,33),an(3,33)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data den(33), cf(33)/ 0.573000e+01, 0.124400e+03/ data em(33),atwt(33)/ 0.206000e+00,74.92/ data coh(0,33),coh(1,33),coh(2,33),coh(3,33)/ $ 0.639750e+01, 0.288866e+00, -0.361747e+00, 0.188788e-01/ data cih(0,33),cih(1,33),cih(2,33),cih(3,33)/ $ -0.339189e+00, 0.162535e+01, -0.250783e+00, 0.909103e-02/ data ka(33),kb(33),la(33),lb(33)/ $ 0.105430e+02, 0.117250e+02, 0.128200e+01, 0.131700e+01/ data name(34), ek(34),el(34)/ $ 'se', 0.126580e+02, 0.165300e+01/ data ak(0,34),ak(1,34),ak(2,34),ak(3,34)/ $ 0.130756e+02, 0.183235e+00, -0.694264e+00, 0.502280e-01/ data al(0,34),al(1,34),al(2,34),al(3,34)/ $ 0.147048e+02, -0.238853e+01, -0.105877e+00, 0.000000e+00/ data am(0,34),am(1,34),am(2,34),am(3,34)/ $ 0.126773e+02, -0.239750e+01, 0.000000e+00, 0.000000e+00/ data an(0,34),an(1,34),an(2,34),an(3,34)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(34),den(34),cf(34)/ 0.789600e+02, 0.479000e+01, $ 0.131100e+03/ data em(34)/ 0.232000e+00/ data coh(0,34),coh(1,34),coh(2,34),coh(3,34)/ $ 0.645637e+01, 0.286737e+00, -0.358794e+00, 0.185618e-01/ data cih(0,34),cih(1,34),cih(2,34),cih(3,34)/ $ -0.432927e+00, 0.172833e+01, -0.277138e+00, 0.111735e-01/ data ka(34),kb(34),la(34),lb(34)/ $ 0.112210e+02, 0.124950e+02, 0.141900e+01, 0.137900e+01/ data name(35), ek(35),el(35)/ $ 'br', 0.134740e+02, 0.178200e+01/ data ak(0,35),ak(1,35),ak(2,35),ak(3,35)/ $ 0.132273e+02, 0.137130e+00, -0.683203e+00, 0.495424e-01/ data al(0,35),al(1,35),al(2,35),al(3,35)/ $ 0.148136e+02, -0.242347e+01, -0.914590e-01, 0.000000e+00/ data am(0,35),am(1,35),am(2,35),am(3,35)/ $ 0.127612e+02, -0.237730e+01, 0.000000e+00, 0.000000e+00/ data an(0,35),an(1,35),an(2,35),an(3,35)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(35),den(35),cf(35)/ 0.799200e+02, 0.312000e+01, $ 0.132700e+03/ data em(35)/ 0.257000e+00/ data coh(0,35),coh(1,35),coh(2,35),coh(3,35)/ $ 0.651444e+01, 0.286324e+00, -0.357027e+00, 0.183557e-01/ data cih(0,35),cih(1,35),cih(2,35),cih(3,35)/ $ -0.448001e+00, 0.176082e+01, -0.285099e+00, 0.117865e-01/ data ka(35),kb(35),la(35),lb(35)/ $ 0.119230e+02, 0.132900e+02, 0.148000e+01, 0.152600e+01/ data name(36), ek(36),el(36),em(36)/ $ 'kr', 0.143220e+02, 0.192000e+01, 0.288000e+00/ data ak(0,36),ak(1,36),ak(2,36),ak(3,36)/ $ 0.135927e+02, -0.305214e-01, -0.651340e+00, 0.477616e-01/ data al(0,36),al(1,36),al(2,36),al(3,36)/ $ 0.149190e+02, -0.242418e+01, -0.876447e-01, 0.000000e+00/ data am(0,36),am(1,36),am(2,36),am(3,36)/ $ 0.128898e+02, -0.226021e+01, 0.000000e+00, 0.000000e+00/ data an(0,36),an(1,36),an(2,36),an(3,36)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(36),den(36),cf(36)/ $ 0.838000e+02, 0.374000e-02, 0.139100e+03/ data coh(0,36),coh(1,36),coh(2,36),coh(3,36)/ $ 0.657113e+01, 0.287711e+00, -0.356311e+00, 0.182470e-01/ data cih(0,36),cih(1,36),cih(2,36),cih(3,36)/ $ -0.391810e+00, 0.173010e+01, -0.276824e+00, 0.111280e-01/ data ka(36),kb(36),la(36),lb(36)/ $ 0.126480e+02, 0.141120e+02, 0.158700e+01, 0.163800e+01/ data name(37), ek(37),el(37),em(37)/ $ 'rb', 0.152000e+02, 0.206500e+01, 0.322000e+00/ data ak(0,37),ak(1,37),ak(2,37),ak(3,37)/ $ 0.130204e+02, 0.382736e+00, -0.732427e+00, 0.529874e-01/ data al(0,37),al(1,37),al(2,37),al(3,37)/ $ 0.149985e+02, -0.239108e+01, -0.959473e-01, 0.000000e+00/ data am(0,37),am(1,37),am(2,37),am(3,37)/ $ 0.130286e+02, -0.238693e+01, 0.000000e+00, 0.000000e+00/ data an(0,37),an(1,37),an(2,37),an(3,37)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(37),den(37),cf(37)/ $ 0.854800e+02, 0.153200e+01, 0.141900e+03/ data coh(0,37),coh(1,37),coh(2,37),coh(3,37)/ $ 0.659750e+01, 0.302389e+00, -0.356755e+00, 0.181706e-01/ data cih(0,37),cih(1,37),cih(2,37),cih(3,37)/ $ -0.128039e+00, 0.153044e+01, -0.227403e+00, 0.739033e-02/ data ka(37),kb(37),la(37),lb(37)/ $ 0.133940e+02, 0.149600e+02, 0.169400e+01, 0.175200e+01/ data name(38), ek(38),el(38),em(38)/ $ 'sr', 0.161050e+02, 0.221600e+01, 0.358000e+00/ data ak(0,38),ak(1,38),ak(2,38),ak(3,38)/ $ 0.135888e+02, 0.220194e-02, -0.638940e+00, 0.460070e-01/ data al(0,38),al(1,38),al(2,38),al(3,38)/ $ 0.150114e+02, -0.228169e+01, -0.126485e+00, 0.000000e+00/ data am(0,38),am(1,38),am(2,38),am(3,38)/ $ 0.131565e+02, -0.236655e+01, 0.000000e+00, 0.000000e+00/ data an(0,38),an(1,38),an(2,38),an(3,38)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(38),den(38),cf(38)/ $ 0.876200e+02, 0.254000e+01, 0.145500e+03/ data coh(0,38),coh(1,38),coh(2,38),coh(3,38)/ $ 0.662203e+01, 0.324559e+00, -0.361651e+00, 0.184800e-01/ data cih(0,38),cih(1,38),cih(2,38),cih(3,38)/ $ 0.799161e-01, 0.138397e+01, -0.192225e+00, 0.478611e-02/ data ka(38),kb(38),la(38),lb(38)/ $ 0.141640e+02, 0.158340e+02, 0.180600e+01, 0.187200e+01/ data name(39), ek(39),el(39),em(39)/ $ 'y ', 0.170800e+02, 0.237300e+01, 0.395000e+00/ data ak(0,39),ak(1,39),ak(2,39),ak(3,39)/ $ 0.134674e+02, 0.191023e+00, -0.686616e+00, 0.497356e-01/ data al(0,39),al(1,39),al(2,39),al(3,39)/ $ 0.151822e+02, -0.238946e+01, -0.881174e-01, 0.000000e+00/ data am(0,39),am(1,39),am(2,39),am(3,39)/ $ 0.132775e+02, -0.243174e+01, 0.000000e+00, 0.000000e+00/ data an(0,39),an(1,39),an(2,39),an(3,39)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(39),den(39),cf(39)/ $ 0.889050e+02, 0.440500e+01, 0.147600e+03/ data coh(0,39),coh(1,39),coh(2,39),coh(3,39)/ $ 0.667096e+01, 0.325070e+00, -0.360613e+00, 0.183326e-01/ data cih(0,39),cih(1,39),cih(2,39),cih(3,39)/ $ 0.629057e-01, 0.141577e+01, -0.199713e+00, 0.533312e-02/ data ka(39),kb(39),la(39),lb(39)/ $ 0.149570e+02, 0.167360e+02, 0.192200e+01, 0.199600e+01/ data name(40), ek(40),el(40),em(40)/ $ 'zr', 0.179980e+02, 0.253200e+01, 0.431000e+00/ data ak(0,40),ak(1,40),ak(2,40),ak(3,40)/ $ 0.127538e+02, 0.697409e+00, -0.789307e+00, 0.564531e-01/ data al(0,40),al(1,40),al(2,40),al(3,40)/ $ 0.152906e+02, -0.238703e+01, -0.912292e-01, 0.000000e+00/ data am(0,40),am(1,40),am(2,40),am(3,40)/ $ 0.134508e+02, -0.250201e+01, 0.000000e+00, 0.000000e+00/ data an(0,40),an(1,40),an(2,40),an(3,40)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(40),den(40),cf(40)/ $ 0.912200e+02, 0.653000e+01, 0.151500e+03/ data coh(0,40),coh(1,40),coh(2,40),coh(3,40)/ $ 0.672275e+01, 0.323964e+00, -0.359463e+00, 0.181890e-01/ data cih(0,40),cih(1,40),cih(2,40),cih(3,40)/ $ 0.366697e-01, 0.145207e+01, -0.208122e+00, 0.595139e-02/ data ka(40),kb(40),la(40),lb(40)/ $ 0.157740e+02, 0.176660e+02, 0.204200e+01, 0.212400e+01/ data name(41), ek(41),el(41),em(41)/ $ 'nb', 0.189860e+02, 0.269800e+01, 0.468000e+00/ data ak(0,41),ak(1,41),ak(2,41),ak(3,41)/ $ 0.133843e+02, 0.281028e+00, -0.686607e+00, 0.486607e-01/ data al(0,41),al(1,41),al(2,41),al(3,41)/ $ 0.152088e+02, -0.220278e+01, -0.136759e+00, 0.000000e+00/ data am(0,41),am(1,41),am(2,41),am(3,41)/ $ 0.135434e+02, -0.250135e+01, 0.000000e+00, 0.000000e+00/ data an(0,41),an(1,41),an(2,41),an(3,41)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(41),den(41),cf(41)/ $ 0.929060e+02, 0.857000e+01, 0.154300e+03/ data coh(0,41),coh(1,41),coh(2,41),coh(3,41)/ $ 0.679013e+01, 0.311282e+00, -0.355233e+00, 0.178231e-01/ data cih(0,41),cih(1,41),cih(2,41),cih(3,41)/ $ 0.202289e-03, 0.149347e+01, -0.217419e+00, 0.662245e-02/ data ka(41),kb(41),la(41),lb(41)/ $ 0.166140e+02, 0.186210e+02, 0.216600e+01, 0.225700e+01/ data name(42), ek(42),el(42),em(42)/ $ 'mo', 0.199990e+02, 0.286600e+01, 0.505000e+00/ data ak(0,42),ak(1,42),ak(2,42),ak(3,42)/ $ 0.139853e+02, -0.117426e+00, -0.591094e+00, 0.417843e-01/ data al(0,42),al(1,42),al(2,42),al(3,42)/ $ 0.153494e+02, -0.226640e+01, -0.116881e+00, 0.000000e+00/ data am(0,42),am(1,42),am(2,42),am(3,42)/ $ 0.136568e+02, -0.248480e+01, 0.000000e+00, 0.000000e+00/ data an(0,42),an(1,42),an(2,42),an(3,42)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(42),den(42),cf(42)/ $ 0.959500e+02, 0.102200e+02, 0.159300e+03/ data coh(0,42),coh(1,42),coh(2,42),coh(3,42)/ $ 0.684600e+01, 0.302790e+00, -0.351131e+00, 0.174403e-01/ data cih(0,42),cih(1,42),cih(2,42),cih(3,42)/ $ -0.562860e-01, 0.155778e+01, -0.233341e+00, 0.785506e-02/ data ka(42),kb(42),la(42),lb(42)/ $ 0.174780e+02, 0.196070e+02, 0.229300e+01, 0.239500e+01/ data name(43), ek(43),el(43),em(43)/ $ 'tc', 0.210450e+02, 0.304300e+01, 0.544000e+00/ data ak(0,43),ak(1,43),ak(2,43),ak(3,43)/ $ 0.128214e+02, 0.751993e+00, -0.787006e+00, 0.558668e-01/ data al(0,43),al(1,43),al(2,43),al(3,43)/ $ 0.155086e+02, -0.233733e+01, -0.987857e-01, 0.000000e+00/ data am(0,43),am(1,43),am(2,43),am(3,43)/ $ 0.137498e+02, -0.244730e+01, 0.000000e+00, 0.000000e+00/ data an(0,43),an(1,43),an(2,43),an(3,43)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(43),den(43),cf(43)/ $ 0.990000e+02, 0.115000e+02, 0.164400e+03/ data coh(0,43),coh(1,43),coh(2,43),coh(3,43)/ $ 0.687599e+01, 0.326165e+00, -0.358969e+00, 0.180482e-01/ data cih(0,43),cih(1,43),cih(2,43),cih(3,43)/ $ 0.757616e-01, 0.144950e+01, -0.204890e+00, 0.564745e-02/ data ka(43),kb(43),la(43),lb(43)/ $ 0.184100e+02, 0.205850e+02, 0.242400e+01, 0.253800e+01/ data name(44), ek(44),el(44),em(44)/ $ 'ru', 0.221170e+02, 0.322400e+01, 0.585000e+00/ data ak(0,44),ak(1,44),ak(2,44),ak(3,44)/ $ 0.126658e+02, 0.885020e+00, -0.811144e+00, 0.573759e-01/ data al(0,44),al(1,44),al(2,44),al(3,44)/ $ 0.154734e+02, -0.223080e+01, -0.119454e+00, 0.000000e+00/ data am(0,44),am(1,44),am(2,44),am(3,44)/ $ 0.138782e+02, -0.248066e+01, 0.000000e+00, 0.000000e+00/ data an(0,44),an(1,44),an(2,44),an(3,44)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(44),den(44),cf(44)/ $ 0.101070e+03, 0.124100e+02, 0.167800e+03/ data coh(0,44),coh(1,44),coh(2,44),coh(3,44)/ $ 0.693136e+01, 0.334794e+00, -0.363497e+00, 0.184429e-01/ data cih(0,44),cih(1,44),cih(2,44),cih(3,44)/ $ -0.424981e-01, 0.154639e+01, -0.226470e+00, 0.718375e-02/ data ka(44),kb(44),la(44),lb(44)/ $ 0.192780e+02, 0.216550e+02, 0.255800e+01, 0.268300e+01/ data name(45), ek(45),el(45),em(45)/ $ 'rh', 0.232200e+02, 0.341200e+01, 0.627000e+00/ data ak(0,45),ak(1,45),ak(2,45),ak(3,45)/ $ 0.121760e+02, 0.119682e+01, -0.866697e+00, 0.606931e-01/ data al(0,45),al(1,45),al(2,45),al(3,45)/ $ 0.155757e+02, -0.224976e+01, -0.113377e+00, 0.000000e+00/ data am(0,45),am(1,45),am(2,45),am(3,45)/ $ 0.140312e+02, -0.261303e+01, 0.000000e+00, 0.000000e+00/ data an(0,45),an(1,45),an(2,45),an(3,45)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(45),den(45),cf(45)/ $ 0.102910e+03, 0.124400e+02, 0.170900e+03/ data coh(0,45),coh(1,45),coh(2,45),coh(3,45)/ $ 0.697547e+01, 0.346394e+00, -0.367794e+00, 0.187885e-01/ data cih(0,45),cih(1,45),cih(2,45),cih(3,45)/ $ -0.160399e+00, 0.164861e+01, -0.250238e+00, 0.893818e-02/ data ka(45),kb(45),la(45),lb(45)/ $ 0.202140e+02, 0.227210e+02, 0.269600e+01, 0.283400e+01/ data name(46), ek(46),el(46),em(46)/ $ 'pd', 0.243500e+02, 0.360500e+01, 0.670000e+00/ data ak(0,46),ak(1,46),ak(2,46),ak(3,46)/ $ 0.139389e+02, 0.164528e+00, -0.662117e+00, 0.476289e-01/ data al(0,46),al(1,46),al(2,46),al(3,46)/ $ 0.155649e+02, -0.217229e+01, -0.127652e+00, 0.000000e+00/ data am(0,46),am(1,46),am(2,46),am(3,46)/ $ 0.141392e+02, -0.257206e+01, 0.000000e+00, 0.000000e+00/ data an(0,46),an(1,46),an(2,46),an(3,46)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(46),den(46),cf(46)/ $ 0.106400e+03, 0.121600e+02, 0.176700e+03/ data coh(0,46),coh(1,46),coh(2,46),coh(3,46)/ $ 0.703216e+01, 0.349838e+00, -0.370099e+00, 0.189983e-01/ data cih(0,46),cih(1,46),cih(2,46),cih(3,46)/ $ -0.267564e+00, 0.173740e+01, -0.269883e+00, 0.103248e-01/ data ka(46),kb(46),la(46),lb(46)/ $ 0.211750e+02, 0.238160e+02, 0.283800e+01, 0.299000e+01/ data name(47), ek(47),el(47),em(47)/ $ 'ag', 0.255140e+02, 0.380600e+01, 0.717000e+00/ data ak(0,47),ak(1,47),ak(2,47),ak(3,47)/ $ 0.133926e+02, 0.441380e+00, -0.693711e+00, 0.482085e-01/ data al(0,47),al(1,47),al(2,47),al(3,47)/ $ 0.156869e+02, -0.222636e+01, -0.112223e+00, 0.000000e+00/ data am(0,47),am(1,47),am(2,47),am(3,47)/ $ 0.141673e+02, -0.248078e+01, 0.000000e+00, 0.000000e+00/ data an(0,47),an(1,47),an(2,47),an(3,47)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(47),den(47),cf(47)/ $ 0.107880e+03, 0.105000e+02, 0.179100e+03/ data coh(0,47),coh(1,47),coh(2,47),coh(3,47)/ $ 0.706446e+01, 0.363456e+00, -0.373597e+00, 0.192478e-01/ data cih(0,47),cih(1,47),cih(2,47),cih(3,47)/ $ -0.166475e+00, 0.165794e+01, -0.248740e+00, 0.866218e-02/ data ka(47),kb(47),la(47),lb(47)/ $ 0.221620e+02, 0.249420e+02, 0.298400e+01, 0.315100e+01/ data name(48), ek(48),el(48),em(48)/ $ 'cd', 0.267110e+02, 0.401800e+01, 0.770000e+00/ data ak(0,48),ak(1,48),ak(2,48),ak(3,48)/ $ 0.125254e+02, 0.107714e+01, -0.831424e+00, 0.579120e-01/ data al(0,48),al(1,48),al(2,48),al(3,48)/ $ 0.159668e+02, -0.238363e+01, -0.801104e-01, 0.000000e+00/ data am(0,48),am(1,48),am(2,48),am(3,48)/ $ 0.143497e+02, -0.252756e+01, 0.000000e+00, 0.000000e+00/ data an(0,48),an(1,48),an(2,48),an(3,48)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(48),den(48),cf(48)/ $ 0.112410e+03, 0.865000e+01, 0.186600e+03/ data coh(0,48),coh(1,48),coh(2,48),coh(3,48)/ $ 0.709856e+01, 0.372199e+00, -0.375345e+00, 0.193481e-01/ data cih(0,48),cih(1,48),cih(2,48),cih(3,48)/ $ -0.516701e-01, 0.157426e+01, -0.227646e+00, 0.705650e-02/ data ka(48),kb(48),la(48),lb(48)/ $ 0.231720e+02, 0.260930e+02, 0.313300e+01, 0.331600e+01/ data name(49), ek(49),el(49),em(49)/ $ 'in', 0.279400e+02, 0.423800e+01, 0.825000e+00/ data ak(0,49),ak(1,49),ak(2,49),ak(3,49)/ $ 0.118198e+02, 0.145768e+01, -0.888529e+00, 0.605982e-01/ data al(0,49),al(1,49),al(2,49),al(3,49)/ $ 0.162101e+02, -0.251838e+01, -0.540061e-01, 0.000000e+00/ data am(0,49),am(1,49),am(2,49),am(3,49)/ $ 0.144115e+02, -0.249401e+01, 0.000000e+00, 0.000000e+00/ data an(0,49),an(1,49),an(2,49),an(3,49)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(49),den(49),cf(49)/ $ 0.114820e+03, 0.728000e+01, 0.190700e+03/ data coh(0,49),coh(1,49),coh(2,49),coh(3,49)/ $ 0.712708e+01, 0.382082e+00, -0.376855e+00, 0.194151e-01/ data cih(0,49),cih(1,49),cih(2,49),cih(3,49)/ $ -0.817283e-02, 0.155865e+01, -0.224492e+00, 0.685776e-02/ data ka(49),kb(49),la(49),lb(49)/ $ 0.242070e+02, 0.272740e+02, 0.328700e+01, 0.348700e+01/ data name(50), ek(50),el(50),em(50)/ $ 'sn', 0.292000e+02, 0.446500e+01, 0.884000e+00/ data ak(0,50),ak(1,50),ak(2,50),ak(3,50)/ $ 0.130323e+02, 0.790788e+00, -0.762349e+00, 0.527872e-01/ data al(0,50),al(1,50),al(2,50),al(3,50)/ $ 0.158638e+02, -0.219019e+01, -0.113539e+00, 0.000000e+00/ data am(0,50),am(1,50),am(2,50),am(3,50)/ $ 0.145572e+02, -0.256792e+01, 0.000000e+00, 0.000000e+00/ data an(0,50),an(1,50),an(2,50),an(3,50)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(50),den(50),cf(50)/ $ 0.118690e+03, 0.576000e+01, 0.197100e+03/ data coh(0,50),coh(1,50),coh(2,50),coh(3,50)/ $ 0.716085e+01, 0.385512e+00, -0.376481e+00, 0.193305e-01/ data cih(0,50),cih(1,50),cih(2,50),cih(3,50)/ $ 0.142151e-01, 0.155754e+01, -0.224736e+00, 0.691395e-02/ data ka(50),kb(50),la(50),lb(50)/ $ 0.252700e+02, 0.284830e+02, 0.344400e+01, 0.366200e+01/ data name(51), ek(51),el(51),em(51)/ $ 'sb', 0.304910e+02, 0.469800e+01, 0.944000e+00/ data ak(0,51),ak(1,51),ak(2,51),ak(3,51)/ $ 0.906990e+01, 0.328791e+01, -0.126203e+01, 0.853470e-01/ data al(0,51),al(1,51),al(2,51),al(3,51)/ $ 0.157557e+02, -0.204460e+01, -0.140745e+00, 0.000000e+00/ data am(0,51),am(1,51),am(2,51),am(3,51)/ $ 0.146268e+02, -0.255562e+01, 0.000000e+00, 0.000000e+00/ data an(0,51),an(1,51),an(2,51),an(3,51)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(51),den(51),cf(51)/ $ 0.121760e+03, 0.669100e+01, 0.202200e+03/ data coh(0,51),coh(1,51),coh(2,51),coh(3,51)/ $ 0.719665e+01, 0.385543e+00, -0.375054e+00, 0.191608e-01/ data cih(0,51),cih(1,51),cih(2,51),cih(3,51)/ $ 0.156362e-01, 0.157175e+01, -0.228753e+00, 0.726386e-02/ data ka(51),kb(51),la(51),lb(51)/ $ 0.263570e+02, 0.297230e+02, 0.360500e+01, 0.384300e+01/ data name(52), ek(52),el(52),em(52)/ $ 'te', 0.318130e+02, 0.493900e+01, 0.100600e+01/ data ak(0,52),ak(1,52),ak(2,52),ak(3,52)/ $ 0.116656e+02, 0.171052e+01, -0.948281e+00, 0.653213e-01/ data al(0,52),al(1,52),al(2,52),al(3,52)/ $ 0.161087e+02, -0.227876e+01, -0.929405e-01, 0.000000e+00/ data am(0,52),am(1,52),am(2,52),am(3,52)/ $ 0.147125e+02, -0.254324e+01, 0.000000e+00, 0.000000e+00/ data an(0,52),an(1,52),an(2,52),an(3,52)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(52),den(52),cf(52)/ $ 0.127600e+03, 0.624000e+01, 0.211900e+03/ data coh(0,52),coh(1,52),coh(2,52),coh(3,52)/ $ 0.723460e+01, 0.382493e+00, -0.372715e+00, 0.189194e-01/ data cih(0,52),cih(1,52),cih(2,52),cih(3,52)/ $ -0.407579e-01, 0.164267e+01, -0.247890e+00, 0.880567e-02/ data ka(52),kb(52),la(52),lb(52)/ $ 0.274710e+02, 0.309930e+02, 0.376900e+01, 0.402900e+01/ data name(53), ek(53),el(53),em(53)/ $ 'i ', 0.331690e+02, 0.518800e+01, 0.107200e+01/ data ak(0,53),ak(1,53),ak(2,53),ak(3,53)/ $ 0.121075e+02, 0.143635e+01, -0.882038e+00, 0.603575e-01/ data al(0,53),al(1,53),al(2,53),al(3,53)/ $ 0.164086e+02, -0.248214e+01, -0.507179e-01, 0.000000e+00/ data am(0,53),am(1,53),am(2,53),am(3,53)/ $ 0.147496e+02, -0.248179e+01, 0.000000e+00, 0.000000e+00/ data an(0,53),an(1,53),an(2,53),an(3,53)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(53),den(53),cf(53)/ $ 0.126910e+03, 0.494000e+01, 0.210700e+03/ data coh(0,53),coh(1,53),coh(2,53),coh(3,53)/ $ 0.727415e+01, 0.377223e+00, -0.369728e+00, 0.186280e-01/ data cih(0,53),cih(1,53),cih(2,53),cih(3,53)/ $ -0.404420e-01, 0.165596e+01, -0.251067e+00, 0.904874e-02/ data ka(53),kb(53),la(53),lb(53)/ $ 0.286100e+02, 0.322920e+02, 0.393700e+01, 0.422000e+01/ data name(54), ek(54),el(54),em(54)/ $ 'xe', 0.345820e+02, 0.545200e+01, 0.114300e+01/ data ak(0,54),ak(1,54),ak(2,54),ak(3,54)/ $ 0.110857e+02, 0.208357e+01, -0.101209e+01, 0.690310e-01/ data al(0,54),al(1,54),al(2,54),al(3,54)/ $ 0.163098e+02, -0.231679e+01, -0.854498e-01, 0.000000e+00/ data am(0,54),am(1,54),am(2,54),am(3,54)/ $ 0.147603e+02, -0.245068e+01, 0.000000e+00, 0.000000e+00/ data an(0,54),an(1,54),an(2,54),an(3,54)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(54),den(54),cf(54)/ $ 0.131300e+03, 0.590000e-02, 0.218000e+03/ data coh(0,54),coh(1,54),coh(2,54),coh(3,54)/ $ 0.731469e+01, 0.370315e+00, -0.366280e+00, 0.183025e-01/ data cih(0,54),cih(1,54),cih(2,54),cih(3,54)/ $ -0.282407e-02, 0.164039e+01, -0.247642e+00, 0.882144e-02/ data ka(54),kb(54),la(54),lb(54)/ $ 0.298020e+02, 0.336440e+02, 0.411100e+01, 0.442200e+01/ data name(55), ek(55),el(55),em(55)/ $ 'cs', 0.359850e+02, 0.571300e+01, 0.121800e+01/ data ak(0,55),ak(1,55),ak(2,55),ak(3,55)/ $ 0.113750e+02, 0.194161e+01, -0.983232e+00, 0.671986e-01/ data al(0,55),al(1,55),al(2,55),al(3,55)/ $ 0.165418e+02, -0.246363e+01, -0.542849e-01, 0.000000e+00/ data am(0,55),am(1,55),am(2,55),am(3,55)/ $ 0.149713e+02, -0.253145e+01, 0.000000e+00, 0.000000e+00/ data an(0,55),an(1,55),an(2,55),an(3,55)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(55),den(55),cf(55)/ $ 0.132910e+03, 0.187300e+01, 0.220700e+03/ data coh(0,55),coh(1,55),coh(2,55),coh(3,55)/ $ 0.733490e+01, 0.376825e+00, -0.365713e+00, 0.181843e-01/ data cih(0,55),cih(1,55),cih(2,55),cih(3,55)/ $ 0.184861e+00, 0.150030e+01, -0.213333e+00, 0.624264e-02/ data ka(55),kb(55),la(55),lb(55)/ $ 0.309700e+02, 0.349840e+02, 0.428600e+01, 0.462000e+01/ data name(56), ek(56),el(56),em(56)/ $ 'ba', 0.374410e+02, 0.598700e+01, 0.129300e+01/ data ak(0,56),ak(1,56),ak(2,56),ak(3,56)/ $ 0.102250e+02, 0.267835e+01, -0.112648e+01, 0.762669e-01/ data al(0,56),al(1,56),al(2,56),al(3,56)/ $ 0.166217e+02, -0.248972e+01, -0.449623e-01, 0.000000e+00/ data am(0,56),am(1,56),am(2,56),am(3,56)/ $ 0.150844e+02, -0.256341e+01, 0.000000e+00, 0.000000e+00/ data an(0,56),an(1,56),an(2,56),an(3,56)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(56),den(56),cf(56)/ $ 0.137360e+03, 0.350000e+01, 0.228100e+03/ data coh(0,56),coh(1,56),coh(2,56),coh(3,56)/ $ 0.735812e+01, 0.379361e+00, -0.364099e+00, 0.179817e-01/ data cih(0,56),cih(1,56),cih(2,56),cih(3,56)/ $ 0.344376e+00, 0.138742e+01, -0.186356e+00, 0.424917e-02/ data ka(56),kb(56),la(56),lb(56)/ $ 0.321910e+02, 0.363760e+02, 0.446700e+01, 0.482800e+01/ data name(57), ek(57),el(57),em(57)/ $ 'la', 0.389250e+02, 0.626700e+01, 0.136300e+01/ data ak(0,57),ak(1,57),ak(2,57),ak(3,57)/ $ 0.109780e+02, 0.238140e+01, -0.103549e+01, 0.702339e-01/ data al(0,57),al(1,57),al(2,57),al(3,57)/ $ 0.163134e+02, -0.220156e+01, -0.980569e-01, 0.000000e+00/ data am(0,57),am(1,57),am(2,57),am(3,57)/ $ 0.151863e+02, -0.258287e+01, 0.000000e+00, 0.000000e+00/ data an(0,57),an(1,57),an(2,57),an(3,57)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(57),den(57),cf(57)/ $ 0.138920e+03, 0.615000e+01, 0.230700e+03/ data coh(0,57),coh(1,57),coh(2,57),coh(3,57)/ $ 0.739532e+01, 0.369895e+00, -0.359376e+00, 0.175406e-01/ data cih(0,57),cih(1,57),cih(2,57),cih(3,57)/ $ 0.409104e+00, 0.133075e+01, -0.170883e+00, 0.304111e-02/ data ka(57),kb(57),la(57),lb(57)/ $ 0.334400e+02, 0.377990e+02, 0.465100e+01, 0.504300e+01/ data name(58), ek(58),el(58),em(58)/ $ 'ce', 0.404440e+02, 0.654900e+01, 0.143400e+01/ data ak(0,58),ak(1,58),ak(2,58),ak(3,58)/ $ 0.102725e+02, 0.274562e+01, -0.114174e+01, 0.774162e-01/ data al(0,58),al(1,58),al(2,58),al(3,58)/ $ 0.165862e+02, -0.236288e+01, -0.654708e-01, 0.000000e+00/ data am(0,58),am(1,58),am(2,58),am(3,58)/ $ 0.152693e+02, -0.258174e+01, 0.000000e+00, 0.000000e+00/ data an(0,58),an(1,58),an(2,58),an(3,58)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(58),den(58),cf(58)/ $ 0.140130e+03, 0.667000e+01, 0.232700e+03/ data coh(0,58),coh(1,58),coh(2,58),coh(3,58)/ $ 0.744255e+01, 0.371328e+00, -0.359642e+00, 0.175852e-01/ data cih(0,58),cih(1,58),cih(2,58),cih(3,58)/ $ 0.439881e+00, 0.130925e+01, -0.164548e+00, 0.252641e-02/ data ka(58),kb(58),la(58),lb(58)/ $ 0.347170e+02, 0.392550e+02, 0.484000e+01, 0.526200e+01/ data name(59), ek(59),el(59),em(59)/ $ 'pr', 0.419910e+02, 0.683500e+01, 0.150800e+01/ data ak(0,59),ak(1,59),ak(2,59),ak(3,59)/ $ 0.110156e+02, 0.222056e+01, -0.102216e+01, 0.690465e-01/ data al(0,59),al(1,59),al(2,59),al(3,59)/ $ 0.167118e+02, -0.240326e+01, -0.612619e-01, 0.000000e+00/ data am(0,59),am(1,59),am(2,59),am(3,59)/ $ 0.153379e+02, -0.257086e+01, 0.000000e+00, 0.000000e+00/ data an(0,59),an(1,59),an(2,59),an(3,59)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(59),den(59),cf(59)/ $ 0.140920e+03, 0.676900e+01, 0.234000e+03/ data coh(0,59),coh(1,59),coh(2,59),coh(3,59)/ $ 0.748347e+01, 0.368431e+00, -0.357689e+00, 0.174099e-01/ data cih(0,59),cih(1,59),cih(2,59),cih(3,59)/ $ 0.449124e+00, 0.130351e+01, -0.161841e+00, 0.227394e-02/ data ka(59),kb(59),la(59),lb(59)/ $ 0.360230e+02, 0.407460e+02, 0.503400e+01, 0.548900e+01/ data name(60), ek(60),el(60),em(60)/ $ 'nd', 0.435690e+02, 0.712600e+01, 0.157500e+01/ data ak(0,60),ak(1,60),ak(2,60),ak(3,60)/ $ 0.117632e+02, 0.179481e+01, -0.936661e+00, 0.635332e-01/ data al(0,60),al(1,60),al(2,60),al(3,60)/ $ 0.165964e+02, -0.226007e+01, -0.872426e-01, 0.000000e+00/ data am(0,60),am(1,60),am(2,60),am(3,60)/ $ 0.154335e+02, -0.259006e+01, 0.000000e+00, 0.000000e+00/ data an(0,60),an(1,60),an(2,60),an(3,60)/ $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00/ data atwt(60),den(60),cf(60)/ $ 0.144270e+03, 0.696000e+01, 0.239600e+03/ data coh(0,60),coh(1,60),coh(2,60),coh(3,60)/ $ 0.752334e+01, 0.366462e+00, -0.356048e+00, 0.172620e-01/ data cih(0,60),cih(1,60),cih(2,60),cih(3,60)/ $ 0.437283e+00, 0.131370e+01, -0.162866e+00, 0.229377e-02/ data ka(60),kb(60),la(60),lb(60)/ $ 0.373590e+02, 0.422690e+02, 0.543100e+01, 0.595600e+01/ data name(61), ek(61),el(61),em(61)/ $ 'pm', 0.451840e+02, 0.742800e+01, 0.165100e+01/ data ak(0,61),ak(1,61),ak(2,61),ak(3,61)/ $ 0.113864e+02, 0.205593e+01, -0.988180e+00, 0.669106e-01/ data al(0,61),al(1,61),al(2,61),al(3,61)/ $ 0.168337e+02, -0.238810e+01, -0.645041e-01, 0.000000e+00/ data am(0,61),am(1,61),am(2,61),am(3,61)/ $ 0.155131e+02, -0.259623e+01, 0.000000e+00, 0.000000e+00/ data an(0,61),an(1,61),an(2,61),an(3,61)/ $ 0.155131e+02, -0.259623e+01, 0.000000e+00, 0.000000e+00/ data atwt(61),den(61),cf(61)/ $ 0.147000e+03, 0.678200e+01, 0.244100e+03/ data coh(0,61),coh(1,61),coh(2,61),coh(3,61)/ $ 0.756222e+01, 0.365055e+00, -0.354511e+00, 0.171214e-01/ data cih(0,61),cih(1,61),cih(2,61),cih(3,61)/ $ 0.405823e+00, 0.133837e+01, -0.167229e+00, 0.255570e-02/ data ka(61),kb(61),la(61),lb(61)/ $ 0.386490e+02, 0.439450e+02, 0.563600e+01, 0.620600e+01/ data name(62), ek(62),el(62),em(62)/ $ 'sm', 0.468350e+02, 0.773700e+01, 0.172900e+01/ data ak(0,62),ak(1,62),ak(2,62),ak(3,62)/ $ 0.119223e+02, 0.179546e+01, -0.942902e+00, 0.644202e-01/ data al(0,62),al(1,62),al(2,62),al(3,62)/ $ 0.168725e+02, -0.239051e+01, -0.601080e-01, 0.000000e+00/ data am(0,62),am(1,62),am(2,62),am(3,62)/ $ 0.156000e+02, -0.261328e+01, 0.000000e+00, 0.000000e+00/ data an(0,62),an(1,62),an(2,62),an(3,62)/ $ 0.156006e+02, -0.261328e+01, 0.000000e+00, 0.000000e+00/ data atwt(62),den(62),cf(62)/ $ 0.150350e+03, 0.753600e+01, 0.249600e+03/ data coh(0,62),coh(1,62),coh(2,62),coh(3,62)/ $ 0.760020e+01, 0.364134e+00, -0.353086e+00, 0.169894e-01/ data cih(0,62),cih(1,62),cih(2,62),cih(3,62)/ $ 0.355383e+00, 0.137733e+01, -0.174941e+00, 0.306213e-02/ data ka(62),kb(62),la(62),lb(62)/ $ 0.401240e+02, 0.454000e+02, 0.584600e+01, 0.645600e+01/ data name(63), ek(63),el(63),em(63)/ $ 'eu', 0.485200e+02, 0.805200e+01, 0.180000e+01/ data ak(0,63),ak(1,63),ak(2,63),ak(3,63)/ $ 0.116168e+02, 0.197533e+01, -0.970901e+00, 0.658459e-01/ data al(0,63),al(1,63),al(2,63),al(3,63)/ $ 0.170692e+02, -0.248046e+01, -0.447055e-01, 0.000000e+00/ data am(0,63),am(1,63),am(2,63),am(3,63)/ $ 0.157063e+02, -0.263481e+01, 0.000000e+00, 0.000000e+00/ data an(0,63),an(1,63),an(2,63),an(3,63)/ $ 0.157063e+02, -0.263481e+01, 0.000000e+00, 0.000000e+00/ data atwt(63),den(63),cf(63)/ $ 0.152000e+03, 0.525900e+01, 0.252400e+03/ data coh(0,63),coh(1,63),coh(2,63),coh(3,63)/ $ 0.763711e+01, 0.363957e+00, -0.351900e+00, 0.168783e-01/ data cih(0,63),cih(1,63),cih(2,63),cih(3,63)/ $ 0.280316e+00, 0.144016e+01, -0.188641e+00, 0.401226e-02/ data ka(63),kb(63),la(63),lb(63)/ $ 0.415290e+02, 0.470270e+02, 0.605900e+01, 0.671400e+01/ data name(64), ek(64),el(64)/ $ 'gd', 0.502400e+02, 0.837600e+01/ data em(64),atwt(64)/1.882,157.260/ data ak(0,64),ak(1,64),ak(2,64),ak(3,64)/ $ 0.991968e+01, 0.303111e+01, -0.117520e+01, 0.786750e-01/ data al(0,64),al(1,64),al(2,64),al(3,64)/ $ 0.171159e+02, -0.247838e+01, -0.437107e-01, 0.000000e+00/ data am(0, 64),am(1, 64),am(2, 64),am(3, 64)/ $ 1.57159e+01, -2.60843e+00, 0.000000e+00, 0.000000e+00/ data an(0, 64),an(1, 64),an(2, 64),an(3, 64)/ $ 1.57159e+01, -2.60843e+00, 0.000000e+00, 0.000000e+00/ data den(64),cf(64)/7.950,2.611e+02/ data coh(0,64),coh(1,64),coh(2,64),coh(3,64)/ $ 0.766938e+01, 0.359752e+00, -0.348899e+00, 0.165890e-01/ data cih(0,64),cih(1,64),cih(2,64),cih(3,64)/ $ 0.273133e+00, 0.143842e+01, -0.186137e+00, 0.375240e-02/ data ka(64),kb(64),la(64),lb(64)/ $ 0.429830e+02, 0.487180e+02, 0.627500e+01, 0.697900e+01/ data name(65), ek(65),el(65),em(65)/ $ 'tb', 0.519960e+02, 0.870800e+01, 0.196700e+01/ data ak(0,65),ak(1,65),ak(2,65),ak(3,65)/ $ 0.113818e+02, 0.214447e+01, -0.999222e+00, 0.675569e-01/ data al(0,65),al(1,65),al(2,65),al(3,65)/ $ 0.171499e+02, -0.245507e+01, -0.471370e-01, 0.000000e+00/ data am(0,65),am(1,65),am(2,65),am(3,65)/ $ 0.158412e+02, -0.264040e+01, 0.000000e+00, 0.000000e+00/ data an(0,65),an(1,65),an(2,65),an(3,65)/ $ 0.158415e+02, -0.264040e+01, 0.000000e+00, 0.000000e+00/ data atwt(65),den(65),cf(65)/ $ 0.158930e+03, 0.827200e+01, 0.263900e+03/ data coh(0,65),coh(1,65),coh(2,65),coh(3,65)/ $ 0.770798e+01, 0.365345e+00, -0.350031e+00, 0.166927e-01/ data cih(0,65),cih(1,65),cih(2,65),cih(3,65)/ $ 0.257539e+00, 0.145064e+01, -0.187591e+00, 0.379932e-02/ data ka(65),kb(65),la(65),lb(65)/ $ 0.444700e+02, 0.503910e+02, 0.649500e+01, 0.724900e+01/ data name(66), ek(66),el(66),em(66)/ $ 'dy', 0.537890e+02, 0.904700e+01, 0.204600e+01/ data ak(0,66),ak(1,66),ak(2,66),ak(3,66)/ $ 0.114845e+02, 0.210451e+01, -0.989870e+00, 0.669382e-01/ data al(0,66),al(1,66),al(2,66),al(3,66)/ $ 0.173446e+02, -0.254821e+01, -0.317606e-01, 0.000000e+00/ data am(0,66),am(1,66),am(2,66),am(3,66)/ $ 0.159225e+02, -0.265289e+01, 0.000000e+00, 0.000000e+00/ data an(0,66),an(1,66),an(2,66),an(3,66)/ $ 0.159225e+02, -0.265289e+01, 0.000000e+00, 0.000000e+00/ data atwt(66),den(66),cf(66)/ $ 0.162510e+03, 0.853600e+01, 0.269800e+03/ data coh(0,66),coh(1,66),coh(2,66),coh(3,66)/ $ 0.774188e+01, 0.367107e+00, -0.349433e+00, 0.166273e-01/ data cih(0,66),cih(1,66),cih(2,66),cih(3,66)/ $ 0.242685e+00, 0.146266e+01, -0.189102e+00, 0.385628e-02/ data ka(66),kb(66),la(66),lb(66)/ $ 0.459850e+02, 0.521780e+02, 0.672000e+01, 0.752800e+01/ data name(67), ek(67),el(67),em(67)/ $ 'ho', 0.556180e+02, 0.939500e+01, 0.212700e+01/ data ak(0,67),ak(1,67),ak(2,67),ak(3,67)/ $ 0.875203e+01, 0.371822e+01, -0.129273e+01, 0.855026e-01/ data al(0,67),al(1,67),al(2,67),al(3,67)/ $ 0.176583e+02, -0.272523e+01, -0.819409e-03, 0.000000e+00/ data am(0,67),am(1,67),am(2,67),am(3,67)/ $ 0.160140e+02, -0.267903e+01, 0.000000e+00, 0.000000e+00/ data an(0,67),an(1,67),an(2,67),an(3,67)/ $ 0.160140e+02, -0.267903e+01, 0.000000e+00, 0.000000e+00/ data atwt(67),den(67),cf(67)/ $ 0.164940e+03, 0.880300e+01, 0.273900e+03/ data coh(0,67),coh(1,67),coh(2,67),coh(3,67)/ $ 0.777470e+01, 0.369722e+00, -0.349132e+00, 0.165862e-01/ data cih(0,67),cih(1,67),cih(2,67),cih(3,67)/ $ 0.228493e+00, 0.147438e+01, -0.190559e+00, 0.390903e-02/ data ka(67),kb(67),la(67),lb(67)/ $ 0.475280e+02, 0.539340e+02, 0.694800e+01, 0.781000e+01/ data name(68), ek(68),el(68),em(68)/ $ 'er', 0.574860e+02, 0.975200e+01, 0.221200e+01/ data ak(0,68),ak(1,68),ak(2,68),ak(3,68)/ $ 0.120195e+02, 0.184815e+01, -0.939582e+00, 0.638106e-01/ data al(0,68),al(1,68),al(2,68),al(3,68)/ $ 0.177988e+02, -0.274671e+01, -0.287580e-02, 0.000000e+00/ data am(0,68),am(1,68),am(2,68),am(3,68)/ $ 0.160672e+02, -0.267580e+01, 0.000000e+00, 0.000000e+00/ data an(0,68),an(1,68),an(2,68),an(3,68)/ $ 0.160672e+02, -0.267580e+01, 0.000000e+00, 0.000000e+00/ data atwt(68),den(68),cf(68)/ $ 0.167270e+03, 0.905100e+01, 0.277700e+03/ data coh(0,68),coh(1,68),coh(2,68),coh(3,68)/ $ 0.780643e+01, 0.373226e+00, -0.349147e+00, 0.165710e-01/ data cih(0,68),cih(1,68),cih(2,68),cih(3,68)/ $ 0.225233e+00, 0.148545e+01, -0.191908e+00, 0.395645e-02/ data ka(68),kb(68),la(68),lb(68)/ $ 0.490990e+02, 0.556900e+02, 0.718100e+01, 0.810300e+01/ data name(69), ek(69),el(69),em(69)/ $ 'tm', 0.593900e+02, 0.101160e+02, 0.230700e+01/ data ak(0,69),ak(1,69),ak(2,69),ak(3,69)/ $ 0.125613e+02, 0.157523e+01, -0.890467e+00, 0.609779e-01/ data al(0,69),al(1,69),al(2,69),al(3,69)/ $ 0.174250e+02, -0.251103e+01, -0.329450e-01, 0.000000e+00/ data am(0,69),am(1,69),am(2,69),am(3,69)/ $ 0.161269e+02, -0.267886e+01, 0.000000e+00, 0.000000e+00/ data an(0,69),an(1,69),an(2,69),an(3,69)/ $ 0.161269e+02, -0.267886e+01, 0.000000e+00, 0.000000e+00/ data atwt(69),den(69),cf(69)/ $ 0.168940e+03, 0.933200e+01, 0.280500e+03/ data coh(0,69),coh(1,69),coh(2,69),coh(3,69)/ $ 0.783711e+01, 0.377547e+00, -0.349441e+00, 0.165780e-01/ data cih(0,69),cih(1,69),cih(2,69),cih(3,69)/ $ 0.202656e+00, 0.149625e+01, -0.193234e+00, 0.400233e-02/ data ka(69),kb(69),la(69),lb(69)/ $ 0.507300e+02, 0.575760e+02, 0.718100e+01, 0.810300e+01/ data name(70), ek(70),el(70),em(70)/ $ 'yb', 0.613320e+02, 0.104880e+02, 0.239800e+01/ data ak(0,70),ak(1,70),ak(2,70),ak(3,70)/ $ 0.742791e+01, 0.428955e+01, -0.135167e+01, 0.866136e-01/ data al(0,70),al(1,70),al(2,70),al(3,70)/ $ 0.169795e+02, -0.222577e+01, -0.732557e-01, 0.000000e+00/ data am(0,70),am(1,70),am(2,70),am(3,70)/ $ 0.161794e+02, -0.267715e+01, 0.000000e+00, 0.000000e+00/ data an(0,70),an(1,70),an(2,70),an(3,70)/ $ 0.139111e+02, -0.240380e+01, 0.000000e+00, 0.000000e+00/ data atwt(70),den(70),cf(70)/ $ 0.173040e+03, 0.697700e+01, 0.287300e+03/ data coh(0,70),coh(1,70),coh(2,70),coh(3,70)/ $ 0.786662e+01, 0.382933e+00, -0.350126e+00, 0.166173e-01/ data cih(0,70),cih(1,70),cih(2,70),cih(3,70)/ $ 0.202248e+00, 0.148804e+01, -0.189143e+00, 0.362264e-02/ data ka(70),kb(70),la(70),lb(70)/ $ 0.523600e+02, 0.593520e+02, 0.741400e+01, 0.840100e+01/ data name(71), ek(71),el(71),em(71)/ $ 'lu', 0.633140e+02, 0.108700e+02, 0.249200e+01/ data ak(0,71),ak(1,71),ak(2,71),ak(3,71)/ $ 0.126387e+02, 0.155476e+01, -0.881094e+00, 0.602036e-01/ data al(0,71),al(1,71),al(2,71),al(3,71)/ $ 0.172638e+02, -0.237189e+01, -0.495994e-01, 0.000000e+00/ data am(0,71),am(1,71),am(2,71),am(3,71)/ $ 0.162289e+02, -0.267128e+01, 0.000000e+00, 0.000000e+00/ data an(0,71),an(1,71),an(2,71),an(3,71)/ $ 0.139813e+02, -0.240841e+01, 0.000000e+00, 0.000000e+00/ data atwt(71),den(71),cf(71)/ $ 0.174990e+03, 0.984200e+01, 0.290600e+03/ data coh(0,71),coh(1,71),coh(2,71),coh(3,71)/ $ 0.789137e+01, 0.386034e+00, -0.349756e+00, 0.165480e-01/ data cih(0,71),cih(1,71),cih(2,71),cih(3,71)/ $ 0.197176e+00, 0.150264e+01, -0.192474e+00, 0.385751e-02/ data ka(71),kb(71),la(71),lb(71)/ $ 0.540630e+02, 0.612820e+02, 0.765400e+01, 0.870800e+01/ data name(72), ek(72),el(72),em(72)/ $ 'hf', 0.653510e+02, 0.112720e+02, 0.260200e+01/ data ak(0,72),ak(1,72),ak(2,72),ak(3,72)/ $ 0.758160e+01, 0.447037e+01, -0.142808e+01, 0.939044e-01/ data al(0,72),al(1,72),al(2,72),al(3,72)/ $ 0.164329e+02, -0.182851e+01, -0.132268e+00, 0.000000e+00/ data am(0,72),am(1,72),am(2,72),am(3,72)/ $ 0.162758e+02, -0.266220e+01, 0.000000e+00, 0.000000e+00/ data an(0,72),an(1,72),an(2,72),an(3,72)/ $ 0.140548e+02, -0.242829e+01, 0.000000e+00, 0.000000e+00/ data atwt(72),den(72),cf(72)/ $ 0.178500e+03, 0.133000e+02, 0.296400e+03/ data coh(0,72),coh(1,72),coh(2,72),coh(3,72)/ $ 0.791803e+01, 0.387021e+00, -0.348881e+00, 0.164406e-01/ data cih(0,72),cih(1,72),cih(2,72),cih(3,72)/ $ 0.199469e+00, 0.150233e+01, -0.191385e+00, 0.374011e-02/ data ka(72),kb(72),la(72),lb(72)/ $ 0.557570e+02, 0.632090e+02, 0.789800e+01, 0.902100e+01/ data name(73), ek(73),el(73),em(73)/ $ 'ta', 0.674140e+02, 0.116800e+02, 0.270300e+01/ data ak(0,73),ak(1,73),ak(2,73),ak(3,73)/ $ 0.865271e+01, 0.373117e+01, -0.126359e+01, 0.823539e-01/ data al(0,73),al(1,73),al(2,73),al(3,73)/ $ 0.172411e+02, -0.230313e+01, -0.591006e-01, 0.000000e+00/ data am(0,73),am(1,73),am(2,73),am(3,73)/ $ 0.163038e+02, -0.266148e+01, 0.000000e+00, 0.000000e+00/ data an(0,73),an(1,73),an(2,73),an(3,73)/ $ 0.141313e+02, -0.247214e+01, 0.000000e+00, 0.000000e+00/ data atwt(73),den(73),cf(73)/ $ 0.180950e+03, 0.166000e+02, 0.300500e+03/ data coh(0,73),coh(1,73),coh(2,73),coh(3,73)/ $ 0.794534e+01, 0.387299e+00, -0.347926e+00, 0.163230e-01/ data cih(0,73),cih(1,73),cih(2,73),cih(3,73)/ $ 0.196871e+00, 0.150623e+01, -0.191396e+00, 0.370889e-02/ data ka(73),kb(73),la(73),lb(73)/ $ 0.575240e+02, 0.652100e+02, 0.814500e+01, 0.934100e+01/ data name(74), ek(74),el(74),em(74)/ $ 'w ', 0.695240e+02, 0.120980e+02, 0.281800e+01/ data ak(0,74),ak(1,74),ak(2,74),ak(3,74)/ $ 0.757541e+01, 0.428874e+01, -0.134998e+01, 0.865200e-01/ data al(0,74),al(1,74),al(2,74),al(3,74)/ $ 0.172533e+02, -0.223874e+01, -0.727338e-01, 0.000000e+00/ data am(0,74),am(1,74),am(2,74),am(3,74)/ $ 0.162613e+02, -0.260672e+01, 0.000000e+00, 0.000000e+00/ data an(0,74),an(1,74),an(2,74),an(3,74)/ $ 0.141536e+02, -0.232580e+01, 0.000000e+00, 0.000000e+00/ data atwt(74),den(74),cf(74)/ $ 0.183920e+03, 0.193000e+02, 0.305400e+03/ data coh(0,74),coh(1,74),coh(2,74),coh(3,74)/ $ 0.797266e+01, 0.387704e+00, -0.347155e+00, 0.162372e-01/ data cih(0,74),cih(1,74),cih(2,74),cih(3,74)/ $ 0.191015e+00, 0.151240e+01, -0.191220e+00, 0.371450e-02/ data ka(74),kb(74),la(74),lb(74)/ $ 0.593100e+02, 0.672330e+02, 0.839600e+01, 0.967000e+01/ data name(75), ek(75),el(75),em(75)/ $ 're', 0.716760e+02, 0.125250e+02, 0.293100e+01/ data ak(0,75),ak(1,75),ak(2,75),ak(3,75)/ $ 0.136944e+01, 0.779444e+01, -0.199822e+01, 0.126225e+00/ data al(0,75),al(1,75),al(2,75),al(3,75)/ $ 0.178750e+02, -0.261051e+01, -0.136093e-01, 0.000000e+00/ data am(0,75),am(1,75),am(2,75),am(3,75)/ $ 0.163564e+02, -0.262453e+01, 0.000000e+00, 0.000000e+00/ data an(0,75),an(1,75),an(2,75),an(3,75)/ $ 0.142392e+02, -0.235326e+01, 0.000000e+00, 0.000000e+00/ data atwt(75),den(75),cf(75)/ $ 0.186200e+03, 0.210200e+02, 0.309200e+03/ data coh(0,75),coh(1,75),coh(2,75),coh(3,75)/ $ 0.799940e+01, 0.388739e+00, -0.346726e+00, 0.161751e-01/ data cih(0,75),cih(1,75),cih(2,75),cih(3,75)/ $ 0.189644e+00, 0.150867e+01, -0.189570e+00, 0.349584e-02/ data ka(75),kb(75),la(75),lb(75)/ $ 0.611310e+02, 0.692980e+02, 0.865100e+01, 0.100080e+02/ data name(76), ek(76),el(76),em(76)/ $ 'os', 0.738720e+02, 0.129640e+02, 0.305000e+01/ data ak(0,76),ak(1,76),ak(2,76),ak(3,76)/ $ 0.137534e+02, 0.102122e+01, -0.777126e+00, 0.538811e-01/ data al(0,76),al(1,76),al(2,76),al(3,76)/ $ 0.173525e+02, -0.228550e+01, -0.588047e-01, 0.000000e+00/ data am(0,76),am(1,76),am(2,76),am(3,76)/ $ 0.164233e+02, -0.263163e+01, 0.000000e+00, 0.000000e+00/ data an(0,76),an(1,76),an(2,76),an(3,76)/ $ 0.142795e+02, -0.221971e+01, 0.000000e+00, 0.000000e+00/ data atwt(76),den(76),cf(76)/ $ 0.190200e+03, 0.225000e+02, 0.315800e+03/ data coh(0,76),coh(1,76),coh(2,76),coh(3,76)/ $ 0.802574e+01, 0.390458e+00, -0.346658e+00, 0.161455e-01/ data cih(0,76),cih(1,76),cih(2,76),cih(3,76)/ $ 0.116448e+00, 0.157615e+01, -0.205332e+00, 0.466731e-02/ data ka(76),kb(76),la(76),lb(76)/ $ 0.629910e+02, 0.714040e+02, 0.891000e+01, 0.103540e+02/ data name(77), ek(77),el(77),em(77)/ $ 'ir', 0.761120e+02, 0.134240e+02, 0.317200e+01/ data ak(0,77),ak(1,77),ak(2,77),ak(3,77)/ $ 0.125506e+02, 0.163090e+01, -0.875676e+00, 0.592011e-01/ data al(0,77),al(1,77),al(2,77),al(3,77)/ $ 0.165270e+02, -0.176315e+01, -0.135232e+00, 0.000000e+00/ data am(0,77),am(1,77),am(2,77),am(3,77)/ $ 0.165144e+02, -0.264832e+01, 0.000000e+00, 0.000000e+00/ data an(0,77),an(1,77),an(2,77),an(3,77)/ $ 0.143422e+02, -0.240183e+01, 0.000000e+00, 0.000000e+00/ data atwt(77),den(77),cf(77)/ $ 0.192200e+03, 0.224200e+02, 0.319100e+03/ data coh(0,77),coh(1,77),coh(2,77),coh(3,77)/ $ 0.805150e+01, 0.393143e+00, -0.347052e+00, 0.161570e-01/ data cih(0,77),cih(1,77),cih(2,77),cih(3,77)/ $ 0.719908e-01, 0.161204e+01, -0.213186e+00, 0.520497e-02/ data ka(77),kb(77),la(77),lb(77)/ $ 0.648860e+02, 0.735490e+02, 0.917300e+01, 0.107060e+02/ data name(78), ek(78),el(78),em(78)/ $ 'pt', 0.783950e+02, 0.138920e+02, 0.329700e+01/ data ak(0,78),ak(1,78),ak(2,78),ak(3,78)/ $ 0.127882e+02, 0.163605e+01, -0.898523e+00, 0.618550e-01/ data al(0,78),al(1,78),al(2,78),al(3,78)/ $ 0.173636e+02, -0.221120e+01, -0.730934e-01, 0.000000e+00/ data am(0,78),am(1,78),am(2,78),am(3,78)/ $ 0.167024e+02, -0.271631e+01, 0.000000e+00, 0.000000e+00/ data an(0,78),an(1,78),an(2,78),an(3,78)/ $ 0.143785e+02, -0.234834e+01, 0.000000e+00, 0.000000e+00/ data atwt(78),den(78),cf(78)/ $ 0.195090e+03, 0.213700e+02, 0.323900e+03/ data coh(0,78),coh(1,78),coh(2,78),coh(3,78)/ $ 0.808084e+01, 0.395790e+00, -0.348032e+00, 0.162345e-01/ data cih(0,78),cih(1,78),cih(2,78),cih(3,78)/ $ 0.420186e-01, 0.163611e+01, -0.217964e+00, 0.552670e-02/ data ka(78),kb(78),la(78),lb(78)/ $ 0.668200e+02, 0.757360e+02, 0.944100e+01, 0.110690e+02/ data name(79), ek(79),el(79),em(79)/ $ 'au', 0.807230e+02, 0.143530e+02, 0.342500e+01/ data ak(0,79),ak(1,79),ak(2,79),ak(3,79)/ $ 0.496352e+01, 0.579212e+01, -0.161842e+01, 0.102911e+00/ data al(0,79),al(1,79),al(2,79),al(3,79)/ $ 0.174240e+02, -0.223911e+01, -0.663720e-01, 0.000000e+00/ data am(0,79),am(1,79),am(2,79),am(3,79)/ $ 0.164734e+02, -0.257834e+01, 0.000000e+00, 0.000000e+00/ data an(0,79),an(1,79),an(2,79),an(3,79)/ $ 0.144398e+02, -0.232838e+01, 0.000000e+00, 0.000000e+00/ data atwt(79),den(79),cf(79)/ $ 0.197200e+03, 0.193700e+02, 0.327400e+03/ data coh(0,79),coh(1,79),coh(2,79),coh(3,79)/ $ 0.810524e+01, 0.400576e+00, -0.349340e+00, 0.163264e-01/ data cih(0,79),cih(1,79),cih(2,79),cih(3,79)/ $ 0.156916e-01, 0.165406e+01, -0.220982e+00, 0.570751e-02/ data ka(79),kb(79),la(79),lb(79)/ $ 0.687790e+02, 0.779680e+02, 0.971100e+01, 0.114390e+02/ data name(80), ek(80),el(80),em(80)/ $ 'hg', 0.831030e+02, 0.148460e+02, 0.356200e+01/ data ak(0,80),ak(1,80),ak(2,80),ak(3,80)/ $ 0.197594e+02, -0.197990e+01, -0.276981e+00, 0.268856e-01/ data al(0,80),al(1,80),al(2,80),al(3,80)/ $ 0.171857e+02, -0.208470e+01, -0.853294e-01, 0.000000e+00/ data am(0,80),am(1,80),am(2,80),am(3,80)/ $ 0.165903e+02, -0.260670e+01, 0.000000e+00, 0.000000e+00/ data an(0,80),an(1,80),an(2,80),an(3,80)/ $ 0.145195e+02, -0.233016e+01, 0.000000e+00, 0.000000e+00/ data atwt(80),den(80),cf(80)/ $ 0.200610e+03, 0.135460e+02, 0.333100e+03/ data coh(0,80),coh(1,80),coh(2,80),coh(3,80)/ $ 0.812542e+01, 0.405858e+00, -0.350329e+00, 0.163772e-01/ data cih(0,80),cih(1,80),cih(2,80),cih(3,80)/ $ 0.114587e+00, 0.158076e+01, -0.202968e+00, 0.435692e-02/ data ka(80),kb(80),la(80),lb(80)/ $ 0.708210e+02, 0.802580e+02, 0.998700e+01, 0.118230e+02/ data name(81), ek(81),el(81),em(81)/ $ 'tl', 0.855280e+02, 0.153440e+02, 0.370000e+01/ data ak(0,81),ak(1,81),ak(2,81),ak(3,81)/ $ 0.152879e+02, 0.273664e+00, -0.638890e+00, 0.457495e-01/ data al(0,81),al(1,81),al(2,81),al(3,81)/ $ 0.177379e+02, -0.237745e+01, -0.433223e-01, 0.000000e+00/ data am(0,81),am(1,81),am(2,81),am(3,81)/ $ 0.166564e+02, -0.261593e+01, 0.000000e+00, 0.000000e+00/ data an(0,81),an(1,81),an(2,81),an(3,81)/ $ 0.145473e+02, -0.226773e+01, 0.000000e+00, 0.000000e+00/ data atwt(81),den(81),cf(81)/ $ 0.204390e+03, 0.118600e+02, 0.339400e+03/ data coh(0,81),coh(1,81),coh(2,81),coh(3,81)/ $ 0.814399e+01, 0.408692e+00, -0.349802e+00, 0.162880e-01/ data cih(0,81),cih(1,81),cih(2,81),cih(3,81)/ $ 0.147052e+00, 0.156695e+01, -0.200347e+00, 0.420901e-02/ data ka(81),kb(81),la(81),lb(81)/ $ 0.728600e+02, 0.825580e+02, 0.102660e+02, 0.122100e+02/ data name(82), ek(82),el(82)/ $ 'pb', 0.880060e+02,0.158600e+02/ data ak(0,82),ak(1,82),ak(2,82),ak(3,82)/ $ 0.863374e+01, 0.369400e+01, -0.121312e+01, 0.774601e-01/ data al(0,82),al(1,82),al(2,82),al(3,82)/ $ 0.177963e+02, -0.237691e+01, -0.455883e-01, 0.000000e+00/ data am(0,82),am(1,82),am(2,82),am(3,82)/ $ 0.167131e+02, -0.261538e+01, 0.000000e+00, 0.000000e+00/ data an(0,82),an(1,82),an(2,82),an(3,82)/ $ 0.145771e+02, -0.225279e+01, 0.000000e+00, 0.000000e+00/ data coh(0,82),coh(1,82),coh(2,82),coh(3,82)/ $ 8.15996e+00, 4.18031e-01, -3.52330e-01, 1.64660e-02/ data cih(0,82),cih(1,82),cih(2,82),cih(3,82)/ $ 1.82167e-01, 1.54661e+00, -1.95973e-01, 3.90772e-03/ data den(82), cf(82)/ 0.113400e+02, 0.344100e+03/ data ka(82),kb(82),la(82),lb(82)/ $ 0.749570e+02, 0.849220e+02, 0.105490e+02, 0.126110e+02/ data em(82),atwt(82)/ 0.385000e+01,207.210/ data name(83), ek(83),el(83),em(83)/ $ 'bi', 0.905270e+02, 0.163850e+02, 0.399900e+01/ data ak(0,83),ak(1,83),ak(2,83),ak(3,83)/ $ 0.944293e+01, 0.344965e+01, -0.119886e+01, 0.783484e-01/ data al(0,83),al(1,83),al(2,83),al(3,83)/ $ 0.175348e+02, -0.223353e+01, -0.596161e-01, 0.000000e+00/ data am(0,83),am(1,83),am(2,83),am(3,83)/ $ 0.167078e+02, -0.258648e+01, 0.000000e+00, 0.000000e+00/ data an(0,83),an(1,83),an(2,83),an(3,83)/ $ 0.146832e+02, -0.230940e+01, 0.000000e+00, 0.000000e+00/ data atwt(83),den(83),cf(83)/ $ 0.209000e+03, 0.980000e+01, 0.347000e+03/ data coh(0,83),coh(1,83),coh(2,83),coh(3,83)/ $ 0.817489e+01, 0.427916e+00, -0.355068e+00, 0.166601e-01/ data cih(0,83),cih(1,83),cih(2,83),cih(3,83)/ $ 0.189860e+00, 0.156125e+01, -0.200932e+00, 0.436768e-02/ data ka(83),kb(83),la(83),lb(83)/ $ 0.770970e+02, 0.873350e+02, 0.108360e+02, 0.130210e+02/ data name(86), ek(86),el(86),em(86)/ $ 'rn', 0.984170e+02, 0.180550e+02, 0.447800e+01/ data ak(0,86),ak(1,86),ak(2,86),ak(3,86)/ $ 0.151782e+02, 0.349020e+00, -0.637638e+00, 0.451377e-01/ data al(0,86),al(1,86),al(2,86),al(3,86)/ $ 0.175028e+02, -0.214876e+01, -0.724638e-01, 0.000000e+00/ data am(0,86),am(1,86),am(2,86),am(3,86)/ $ 0.169000e+02, -0.260945e+01, 0.000000e+00, 0.000000e+00/ data an(0,86),an(1,86),an(2,86),an(3,86)/ $ 0.147243e+02, -0.212905e+01, 0.000000e+00, 0.000000e+00/ data atwt(86),den(86),cf(86)/ $ 0.222000e+03, 0.973000e-02, 0.368600e+03/ data coh(0,86),coh(1,86),coh(2,86),coh(3,86)/ $ 0.822553e+01, 0.451478e+00, -0.362056e+00, 0.171556e-01/ data cih(0,86),cih(1,86),cih(2,86),cih(3,86)/ $ 0.196619e+00, 0.160080e+01, -0.213800e+00, 0.551717e-02/ data ka(86),kb(86),la(86),lb(86)/ $ 0.838000e+02, 0.948770e+02, 0.117240e+02, 0.143160e+02/ data name(90), ek(90),el(90),em(90)/ $ 'th', 0.109649e+03, 0.204700e+02, 0.518200e+01/ data ak(0,90),ak(1,90),ak(2,90),ak(3,90)/ $ 0.134336e+02, 0.134805e+01, -0.813280e+00, 0.555664e-01/ data al(0,90),al(1,90),al(2,90),al(3,90)/ $ 0.185481e+02, -0.261281e+01, -0.790574e-02, 0.000000e+00/ data am(0,90),am(1,90),am(2,90),am(3,90)/ $ 0.170483e+02, -0.258569e+01, 0.000000e+00, 0.000000e+00/ data an(0,90),an(1,90),an(2,90),an(3,90)/ $ 0.147730e+02, -0.191192e+01, 0.000000e+00, 0.000000e+00/ data atwt(90),den(90),cf(90)/ $ 0.232000e+03, 0.117000e+02, 0.385200e+03/ data coh(0,90),coh(1,90),coh(2,90),coh(3,90)/ $ 0.827843e+01, 0.479056e+00, -0.367657e+00, 0.174621e-01/ data cih(0,90),cih(1,90),cih(2,90),cih(3,90)/ $ 0.170890e+00, 0.165561e+01, -0.229702e+00, 0.692516e-02/ data ka(90),kb(90),la(90),lb(90)/ $ 0.933340e+02, 0.105592e+03, 0.129660e+02, 0.162000e+02/ data name(92), ek(92),el(92),em(92)/ $ 'u ', 0.115603e+03, 0.217560e+02, 0.554900e+01/ data ak(0,92),ak(1,92),ak(2,92),ak(3,92)/ $ 0.137951e+02, 0.123983e+01, -0.801545e+00, 0.553596e-01/ data al(0,92),al(1,92),al(2,92),al(3,92)/ $ 0.175258e+02, -0.207237e+06, -0.723932e-01, 0.000000e+00/ data am(0,92),am(1,92),am(2,92),am(3,92)/ $ 0.170353e+02, -0.256903e+01, 0.000000e+00, 0.000000e+00/ data an(0,92),an(1,92),an(2,92),an(3,92)/ $ 0.149036e+02, -0.212148e+01, 0.000000e+00, 0.000000e+00/ data atwt(92),den(92),cf(92)/ $ 0.238070e+03, 0.190500e+02, 0.395300e+03/ data coh(0,92),coh(1,92),coh(2,92),coh(3,92)/ $ 0.833010e+01, 0.478314e+00, -0.367250e+00, 0.174129e-01/ data cih(0,92),cih(1,92),cih(2,92),cih(3,92)/ $ 0.108277e+00, 0.174158e+01, -0.254104e+00, 0.895056e-02/ data ka(92),kb(92),la(92),lb(92)/ $ 0.984280e+02, 0.111289e+03, 0.136130e+02, 0.172180e+02/ data name(94), ek(94),el(94),em(94)/ $ 'pu', 0.121760e+03, 0.230950e+02, 0.591400e+01/ data ak(0,94),ak(1,94),ak(2,94),ak(3,94)/ $ 0.182787e+02, -0.117371e+01, -0.368344e+00, 0.298738e-01/ data al(0,94),al(1,94),al(2,94),al(3,94)/ $ 0.175519e+02, -0.202162e+01, -0.822940e-01, 0.000000e+00/ data am(0,94),am(1,94),am(2,94),am(3,94)/ $ 0.172953e+02, -0.262164e+01, 0.000000e+00, 0.000000e+00/ data an(0,94),an(1,94),an(2,94),an(3,94)/ $ 0.148535e+02, -0.187733e+01, 0.000000e+00, 0.000000e+00/ data atwt(94),den(94),cf(94)/ $ 0.239100e+03, 0.197000e+02, 0.397000e+03/ data coh(0,94),coh(1,94),coh(2,94),coh(3,94)/ $ 0.838174e+01, 0.477085e+00, -0.366556e+00, 0.173422e-01/ data cih(0,94),cih(1,94),cih(2,94),cih(3,94)/ $ 0.388791e-01, 0.182229e+01, -0.276099e+00, 0.107392e-01/ data ka(94),kb(94),la(94),lb(94)/ $ 0.103653e+03, 0.117146e+03, 0.142790e+02, 0.182780e+02/ data name(84),name(85),name(87),name(88),name(89),name(91) $ ,name(93)/' ',' ',' ',' ',' ',' ',' '/ c###the following is the l2 and l3 edge energies and l3 edge jumps c for each supported element between nickel and plutonium. c lj2 is the l2 edge jump. data lj2 / 0.141000e+01 / data l2(28),l3(28),lj3(28)/ $ 0.872000e+00, 0.885000e+00, 0.277200e+01/ data l2(29),l3(29),lj3(29)/ $ 0.952000e+00, 0.932000e+00, 0.287400e+01/ data l2(30),l3(30),lj3(30)/ $ 0.104400e+01, 0.102100e+01, 0.568400e+01/ data l2(31),l3(31),lj3(31)/ $ 0.114200e+01, 0.111500e+01, 0.567100e+01/ data l2(32),l3(32),lj3(32)/ $ 0.124900e+01, 0.121800e+01, 0.570400e+01/ data l2(33),l3(33),lj3(33)/ $ 0.136000e+01, 0.132500e+01, 0.487500e+01/ data l2(34),l3(34),lj3(34)/ $ 0.147700e+01, 0.143600e+01, 0.458700e+01/ data l2(35),l3(35),lj3(35)/ $ 0.159600e+01, 0.155000e+01, 0.455700e+01/ data l2(36),l3(36),lj3(36)/ $ 0.172600e+01, 0.167500e+01, 0.417000e+01/ data l2(37),l3(37),lj3(37)/ $ 0.186300e+01, 0.180500e+01, 0.422300e+01/ data l2(38),l3(38),lj3(38)/ $ 0.200700e+01, 0.194000e+01, 0.390600e+01/ data l2(39),l3(39),lj3(39)/ $ 0.215600e+01, 0.208000e+01, 0.403600e+01/ data l2(40),l3(40),lj3(40)/ $ 0.230700e+01, 0.222300e+01, 0.397600e+01/ data l2(41),l3(41),lj3(41)/ $ 0.246500e+01, 0.237100e+01, 0.377400e+01/ data l2(42),l3(42),lj3(42)/ $ 0.262500e+01, 0.252000e+01, 0.367500e+01/ data l2(43),l3(43),lj3(43)/ $ 0.279300e+01, 0.267700e+01, 0.359100e+01/ data l2(44),l3(44),lj3(44)/ $ 0.296700e+01, 0.283800e+01, 0.343100e+01/ data l2(45),l3(45),lj3(45)/ $ 0.314600e+01, 0.300300e+01, 0.372100e+01/ data l2(46),l3(46),lj3(46)/ $ 0.333000e+01, 0.317300e+01, 0.340200e+01/ data l2(47),l3(47),lj3(47)/ $ 0.352400e+01, 0.335100e+01, 0.322300e+01/ data l2(48),l3(48),lj3(48)/ $ 0.372700e+01, 0.353700e+01, 0.324900e+01/ data l2(49),l3(49),lj3(49)/ $ 0.393800e+01, 0.373000e+01, 0.325500e+01/ data l2(50),l3(50),lj3(50)/ $ 0.415600e+01, 0.392900e+01, 0.306000e+01/ data l2(51),l3(51),lj3(51)/ $ 0.438100e+01, 0.413200e+01, 0.293900e+01/ data l2(52),l3(52),lj3(52)/ $ 0.461200e+01, 0.434100e+01, 0.297900e+01/ data l2(53),l3(53),lj3(53)/ $ 0.485200e+01, 0.455700e+01, 0.285600e+01/ data l2(54),l3(54),lj3(54)/ $ 0.510000e+01, 0.478100e+01, 0.287900e+01/ data l2(55),l3(55),lj3(55)/ $ 0.535900e+01, 0.501200e+01, 0.284700e+01/ data l2(56),l3(56),lj3(56)/ $ 0.562400e+01, 0.524700e+01, 0.283900e+01/ data l2(57),l3(57),lj3(57)/ $ 0.589100e+01, 0.548300e+01, 0.271700e+01/ data l2(58),l3(58),lj3(58)/ $ 0.616500e+01, 0.572400e+01, 0.273700e+01/ data l2(59),l3(59),lj3(59)/ $ 0.644100e+01, 0.596500e+01, 0.269500e+01/ data l2(60),l3(60),lj3(60)/ $ 0.672200e+01, 0.620800e+01, 0.266200e+01/ data l2(61),l3(61),lj3(61)/ $ 0.701300e+01, 0.646000e+01, 0.270200e+01/ data l2(62),l3(62),lj3(62)/ $ 0.731200e+01, 0.671700e+01, 0.268000e+01/ data l2(63),l3(63),lj3(63)/ $ 0.761800e+01, 0.697700e+01, 0.272300e+01/ data l2(64),l3(64),lj3(64)/ $ 0.793100e+01, 0.724300e+01, 0.270100e+01/ data l2(65),l3(65),lj3(65)/ $ 0.825200e+01, 0.751500e+01, 0.271300e+01/ data l2(66),l3(66),lj3(66)/ $ 0.858100e+01, 0.779000e+01, 0.904700e+01/ data l2(67),l3(67),lj3(67)/ $ 0.891900e+01, 0.807100e+01, 0.286300e+01/ data l2(68),l3(68),lj3(68)/ $ 0.926500e+01, 0.835800e+01, 0.293300e+01/ data l2(69),l3(69),lj3(69)/ $ 0.961800e+01, 0.864800e+01, 0.275800e+01/ data l2(70),l3(70),lj3(70)/ $ 0.997800e+01, 0.894300e+01, 0.257300e+01/ data l2(71),l3(71),lj3(71)/ $ 0.103490e+02, 0.924400e+01, 0.262000e+01/ data l2(72),l3(72),lj3(72)/ $ 0.107390e+02, 0.956100e+01, 0.241500e+01/ data l2(73),l3(73),lj3(73)/ $ 0.111360e+02, 0.988100e+01, 0.260000e+01/ data l2(74),l3(74),lj3(74)/ $ 0.115420e+02, 0.102040e+02, 0.261700e+01/ data l2(75),l3(75),lj3(75)/ $ 0.119570e+02, 0.105340e+02, 0.267500e+01/ data l2(76),l3(76),lj3(76)/ $ 0.123840e+02, 0.108710e+02, 0.252900e+01/ data l2(77),l3(77),lj3(77)/ $ 0.128240e+02, 0.112150e+02, 0.238700e+01/ data l2(78),l3(78),lj3(78)/ $ 0.138920e+02, 0.115640e+02, 0.263200e+01/ data l2(79),l3(79),lj3(79)/ $ 0.137330e+02, 0.119180e+02, 0.243900e+01/ data l2(80),l3(80),lj3(80)/ $ 0.142090e+02, 0.122840e+02, 0.240000e+01/ data l2(81),l3(81),lj3(81)/ $ 0.146980e+02, 0.126570e+02, 0.249800e+01/ data l2(82),l3(82),lj3(82)/ $ 0.151980e+02, 0.130550e+02, 0.246600e+01/ data l2(83),l3(83),lj3(83)/ $ 0.157080e+02, 0.134180e+02, 0.233800e+01/ data l2(86),l3(86),lj3(86)/ $ 0.173340e+02, 0.146120e+02, 0.234400e+01/ data l2(90),l3(90),lj3(90)/ $ 0.196920e+02, 0.163000e+02, 0.238800e+01/ data l2(92),l3(92),lj3(92)/ $ 0.209470e+02, 0.171670e+02, 0.229200e+01/ data l2(94),l3(94),lj3(94)/ $ 0.222630e+02, 0.180530e+02, 0.225100e+01/ c----------------------------------------------------------------- c no data for element numbers 84, 85, 87, 88, 89, 91, and 93 c----------------------------------------------------------------- c%%%% test = 'ab' c----------------------------------------------------------------- c begin by error checking name and z c----------------------------------------------------------------- c check for null or undocumented elements symb = sym call case(test,symb) if (symb.eq.'nu') then symb = 'h' iz = 1 endif if ( (iz.eq.84).or.(iz.eq.85).or.(iz.eq.87).or.(iz.eq.88).or. $ (iz.eq.89).or.(iz.eq.91).or.(iz.eq.93) ) then ier = 3 if (erf) call messag('No data for Z=84,85,87-89,91,93') goto 999 elseif ( (symb.eq.'po').or.(symb.eq.'at').or.(symb.eq.'fr').or. $ (symb.eq.'ra').or.(symb.eq.'ac').or.(symb.eq.'pa').or. $ (symb.eq.'np') ) then ier = 3 if (erf) $ call messag('No data for Po, At, Fr, Ra, Ac, Pa, or Np') goto 999 endif c check for missing element if ((symb.eq.' ').and.(iz.eq.0)) then ier = 7 if (erf) call messag('No atomic symbol or Z specified '// $ 'on input') goto 999 endif c search for symbol or z is too big or get symbol from z if (symb.ne.' ') then goto 10 elseif (iz.gt.nelem) then ier = 4 if (erf) call messag('No data for z > 94.') goto 999 else symb = name(iz) goto 20 endif 10 continue c search for z from symbol j = 1 15 continue if ((j.le.nelem).and.(symb.eq.name(j))) then iz = j elseif (j.gt.nelem) then ier = 2 if (erf) call messag('Unknown atomic symbol') goto 999 else j = j+1 goto 15 endif c------------------------------------------------------------- c checking for zero energy input, if en<0 skip over abs. calc. c------------------------------------------------------------- 20 continue if (en.lt.0) goto 160 if (abs(en).lt.eps) then ier = 1 if (erf) call messag('Can not calculate absorption '// $ 'for zero input energy') goto 999 else e = en endif c----------------------------------------------------------------- c check for middle of edge input c----------------------------------------------------------------- if ((e.lt.ek(iz)+eps).and.(e.gt.ek(iz)-eps)) then goto 24 elseif ((e.lt.el(iz)+eps).and.(e.gt.el(iz)-eps)) then goto 24 elseif ((e.lt.em(iz)+eps).and.(e.gt.em(iz)-eps)) then goto 24 else goto 27 endif 24 continue ier = 6 if (erf) call messag('Energy at the middle of edge, '// $ 'using pre-edge fit results may be wrong') c----------------------------------------------------------------- c initialize c----------------------------------------------------------------- 27 continue bsum = 0 sum = 0 chs = 0 csum = 0 cis = 0 cisum = 0 c----------------------------------------------------------------- c select correct range c----------------------------------------------------------------- if (e.ge.ek(iz)) then goto 90 elseif ( (e.lt.ek(iz)).and.(e.ge.el(iz)) ) then goto 30 elseif ( (e.lt.el(iz)).and.(e.ge.em(iz)) ) then goto 50 elseif (e.lt.em(iz)) then goto 70 endif c----------------------------------------------------------------- c start calculation, polynomial exapnsions in log(e) c for photoelectric, coherent, and incoherent x-sections c n is the z number of the atom. c----------------------------------------------------------------- c between k and l edges 30 continue do 40 i=0,3 sum = sum + al(i,iz)*(log(e))**i 40 continue goto 120 c between l and m edges 50 continue do 60 i=0,3 sum = sum + am(i,iz)*(log(e))**i 60 continue goto 110 c below m edge 70 continue do 80 i=0,3 sum = sum + an(i,iz)*(log(e))**i 80 continue goto 120 c above k-edge 90 continue do 100 i=0,3 sum = sum + ak(i,iz)*(log(e))**i 100 continue goto 120 110 if (iz.le.29) then ier = 5 if (erf) call messag('Warning: McMaster uses L1 edge '// $ 'results. May be wrong for Z<30') bax = exp(sum) goto 130 endif 120 continue bax = exp(sum) c---------------------------------------------------------------- c correct for l-edges since mcmaster uses l1-edge c use edge jumps for correct x-sections c---------------------------------------------------------------- if(e.ge.l3(iz).and.e.lt.l2(iz)) bax = bax*lj3(iz) if(e.ge.l2(iz).and.e.lt.el(iz)) bax = bax*lj3(iz)*lj2 c---------------------------------------------------------------- c calculate coherent and incoherent x-sections c---------------------------------------------------------------- 130 continue do 140 i=0,3 chs = chs + coh(i,iz)*(log(e))**i 140 continue bcox = exp(chs) do 150 i=0,3 cis = cis + cih(i,iz)*(log(e))**i 150 continue binx = exp(cis) c sum for total x-section btox = bax+bcox+binx c--------------------------------------------------------------- c calculation done in barns/atom, convert to cm^2/gm if desired c--------------------------------------------------------------- call case(test,unit) if (unit.eq.'c') then xsec(1) = bax/cf(iz) xsec(2) = bcox/cf(iz) xsec(3) = binx/cf(iz) xsec(4) = btox/cf(iz) else xsec(1) = bax xsec(2) = bcox xsec(3) = binx xsec(4) = btox endif xsec(5) = cf(iz) xsec(6) = btox*den(iz)/cf(iz) 160 continue xsec(7) = atwt(iz) xsec(8) = den(iz) if (iz.ge.29) xsec(9) = lj2 xsec(10) = lj3(iz) c---------------------------------------------------------------- c fill the energy array c---------------------------------------------------------------- energy(1) = ek(iz) energy(2) = el(iz) energy(3) = l2(iz) energy(4) = l3(iz) energy(5) = em(iz) energy(6) = ka(iz) energy(7) = kb(iz) energy(8) = la(iz) energy(9) = lb(iz) 999 continue return c end subroutine mucal end c*********************************************************************** c%%% c%%% information from the periodic table c%%% c%%% contents: c%%% s z2s(isym,elem): returns z number from symbol c%%% f is2z(sym): returns symbol from z number c%%% f s2e(elem,edge): returns energy from symbol and edge c%%% s case(test,word): returns word in same case as test c%%% s upper(word): returns word in upper case c%%% s lower(word): returns word in lower case c%%% s fixsym(elem): returns char*2 as 1st upper, 2nd lower c%%% f istrln(word): returns length (not dimension) of word c%%% c%%% s = subroutine, f = function c%%% c%%% upper, lower, istrln written by steven zabinsky c%%% all others written by bruce ravel and matthew newville c%%% c%%% all executable statements in this source code must be in the same case c%%% that is, the executable statements must be uniformly upper case c%%% or uniformly lower case c%%% subroutine z2s(isym,elem) c--------------------------------------------------------------------------- c copyright 1993 university of washington matt newville and bruce ravel c--------------------------------------------------------------------------- c returns atomic symbol from z number: default is ' ' character*2 elem,symbol(103) c data (symbol(i),i=1,103) /'h','he','li','be','b','c','n','o', $'f','ne','na','mg','al','si','p','s','cl','ar','k','ca','sc', $'ti','v','cr','mn','fe','co','ni','cu','zn','ga','ge','as','se', $'br','kr','rb','sr','y','zr','nb','mo','tc','ru','rh','pd','ag', $'cd','in','sn','sb','te','i','xe','cs','ba','la','ce','pr','nd', $'pm','sm','eu','gd','tb','dy','ho','er','tm','yb','lu','hf','ta', $'w','re','os','ir','pt','au','hg','tl','pb','bi','po','at','rn', $'fr','ra','ac','th','pa','u','np','pu','am','cm','bk','cf','es', $'fm','md','no','lr'/ elem = ' ' if ((isym.ge.1).and.(isym.le.103)) elem = symbol(isym) return end c************************************************************************ function s2e(elem,edge) c double precision function s2e(elem,edge) c implicit double precision (a-h,o-z) c this function returns the energy of a given element and edeg c input: c elem: 2 character atomic symbol (case insensitive) c edge: k, l1, l2, l3 (case insensitive) c output: c s2e: energy value parameter (nelem=94) character*2 elem, edge, ed, el, test dimension ek(nelem), el1(nelem), el2(nelem), el3(nelem) data (ek(i),i=1,60) / $ 0.140000e-01, 0.250000e-01, 0.550000e-01, 0.112000e+00, $ 0.188000e+00, 0.284000e+00, 0.402000e+00, 0.537000e+00, $ 0.686000e+00, 0.867000e+00, 0.107200e+01, 0.130500e+01, $ 0.156000e+01, 0.183900e+01, 0.214900e+01, 0.247200e+01, $ 0.282200e+01, 0.320200e+01, 0.360700e+01, 0.403800e+01, $ 0.449300e+01, 0.496500e+01, 0.546500e+01, 0.598900e+01, $ 0.654000e+01, 0.711200e+01, 0.770900e+01, 0.833300e+01, $ 0.897900e+01, 0.965900e+01, 0.103670e+02, 0.111040e+02, $ 0.118680e+02, 0.126580e+02, 0.134740e+02, 0.143220e+02, $ 0.152000e+02, 0.161050e+02, 0.170800e+02, 0.179980e+02, $ 0.189860e+02, 0.199990e+02, 0.210450e+02, 0.221170e+02, $ 0.232200e+02, 0.243500e+02, 0.255140e+02, 0.267110e+02, $ 0.279400e+02, 0.292000e+02, 0.304910e+02, 0.318130e+02, $ 0.331690e+02, 0.345820e+02, 0.359850e+02, 0.374410e+02, $ 0.389250e+02, 0.404440e+02, 0.419910e+02, 0.435690e+02/ data (ek(i),i=61,nelem) / $ 0.451840e+02, 0.468350e+02, 0.485200e+02, 0.502400e+02, $ 0.519960e+02, 0.537890e+02, 0.556180e+02, 0.574860e+02, $ 0.593900e+02, 0.613320e+02, 0.633140e+02, 0.653510e+02, $ 0.674140e+02, 0.695240e+02, 0.716760e+02, 0.738720e+02, $ 0.761120e+02, 0.783950e+02, 0.807230e+02, 0.831030e+02, $ 0.855280e+02, 0.880060e+02, 0.905270e+02, 0.931050e+02, $ 0.957300e+02, 0.984170e+02, 0.101137e+03, 0.103922e+02, $ 0.106755e+03, 0.109651e+03, 0.112601e+03, 0.115603e+03, $ 0., 0.121760e+03/ data (el1(i),i=1,60) / $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, $ 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.630000e-01, $ 0.870000e-01, 0.118000e+00, 0.153000e+00, 0.193000e+00, $ 0.238000e+00, 0.287000e+00, 0.341000e+00, 0.400000e+00, $ 0.463000e+00, 0.531000e+00, 0.604000e+00, 0.682000e+00, $ 0.754000e+00, 0.842000e+00, 0.929000e+00, 0.101200e+01, $ 0.110000e+01, 0.119600e+01, 0.130200e+01, 0.141400e+01, $ 0.153000e+01, 0.165300e+01, 0.178200e+01, 0.192000e+01, $ 0.206500e+01, 0.221600e+01, 0.237300e+01, 0.253200e+01, $ 0.269800e+01, 0.286600e+01, 0.304300e+01, 0.322400e+01, $ 0.341200e+01, 0.360500e+01, 0.380600e+01, 0.401800e+01, $ 0.423800e+01, 0.446500e+01, 0.469800e+01, 0.493900e+01, $ 0.518800e+01, 0.545200e+01, 0.571300e+01, 0.598700e+01, $ 0.626700e+01, 0.654900e+01, 0.683500e+01, 0.712600e+01/ data (el1(i),i=61,nelem) / $ 0.742800e+01, 0.773700e+01, 0.805200e+01, 0.837600e+01, $ 0.870800e+01, 0.904700e+01, 0.939500e+01, 0.975200e+01, $ 0.101160e+02, 0.104880e+02, 0.108700e+02, 0.112720e+02, $ 0.116800e+02, 0.120980e+02, 0.125250e+02, 0.129640e+02, $ 0.134240e+02, 0.138920e+02, 0.143530e+02, 0.148460e+02, $ 0.153440e+02, 0.158600e+02, 0.163850e+02, 0.169390e+02, $ 0.174930e+02, 0.180490e+02, 0.186390e+02, 0.192370e+02, $ 0.198400e+02, 0.204700e+02, 0.211050e+02, 0.217560e+02, $ 0., 0.230950e+02/ data (el2(i), i=1,63) / 27*0., $ 0.872000e+00, 0.952000e+00, 0.104400e+01, 0.114200e+01, $ 0.124900e+01, 0.136000e+01, 0.147700e+01, 0.159600e+01, $ 0.172600e+01, 0.186300e+01, 0.200700e+01, 0.215600e+01, $ 0.230700e+01, 0.246500e+01, 0.262500e+01, 0.279300e+01, $ 0.296700e+01, 0.314600e+01, 0.333000e+01, 0.352400e+01, $ 0.372700e+01, 0.393800e+01, 0.415600e+01, 0.438100e+01, $ 0.461200e+01, 0.485200e+01, 0.510000e+01, 0.535900e+01, $ 0.562400e+01, 0.589100e+01, 0.616500e+01, 0.644100e+01, $ 0.672200e+01, 0.701300e+01, 0.731200e+01, 0.761800e+01/ data (el2(i), i=64,nelem) / $ 0.793100e+01, 0.825200e+01, 0.858100e+01, 0.891900e+01, $ 0.926500e+01, 0.961800e+01, 0.997800e+01, 0.103490e+02, $ 0.107390e+02, 0.111360e+02, 0.115420e+02, 0.119570e+02, $ 0.123840e+02, 0.128240e+02, 0.138920e+02, 0.137330e+02, $ 0.142090e+02, 0.146980e+02, 0.151980e+02, 0.157080e+02, $ 0.162440e+02, 0.167850e+02, 0.173370e+02, 0.179070e+02, $ 0.184840e+02, 0.190830e+02, 0.196920e+02, 0.203140e+02, $ 0.209470e+02, 0, 0.222630e+02/ data (el3(i), i=1,63) / 27*0., $ 0.885000e+00, 0.932000e+00, 0.102100e+01, 0.111500e+01, $ 0.121800e+01, 0.132500e+01, 0.143600e+01, 0.155000e+01, $ 0.167500e+01, 0.180500e+01, 0.194000e+01, 0.208000e+01, $ 0.222300e+01, 0.237100e+01, 0.252000e+01, 0.267700e+01, $ 0.283800e+01, 0.300300e+01, 0.317300e+01, 0.335100e+01, $ 0.353700e+01, 0.373000e+01, 0.392900e+01, 0.413200e+01, $ 0.434100e+01, 0.455700e+01, 0.478100e+01, 0.501200e+01, $ 0.524700e+01, 0.548300e+01, 0.572400e+01, 0.596500e+01, $ 0.620800e+01, 0.646000e+01, 0.671700e+01, 0.697700e+01/ data (el3(i), i=64,nelem) / $ 0.724300e+01, 0.751500e+01, 0.779000e+01, 0.807100e+01, $ 0.835800e+01, 0.864800e+01, 0.894300e+01, 0.924400e+01, $ 0.956100e+01, 0.988100e+01, 0.102040e+02, 0.105340e+02, $ 0.108710e+02, 0.112150e+02, 0.115640e+02, 0.119180e+02, $ 0.122840e+02, 0.126570e+02, 0.130550e+02, 0.134180e+02, $ 0.138140e+02, 0.142140e+02, 0.146120e+02, 0.150310e+02, $ 0.154440e+02, 0.158710e+02, 0.163000e+02, 0.167330e+02, $ 0.171670e+02, 0., 0.180530e+02/ i=is2z(elem) ed=edge el=elem c check case c ab must be the same case as k,l1,l2,l3 below for this to work! test = 'ab' call case(test,ed) call case(test,el) if (el.eq.'nu') then s2e = 0. goto 99 endif if (ed.eq.'k') then s2e = ek(i) elseif (ed.eq.'l1') then s2e = el1(i) elseif (ed.eq.'l2') then s2e = el2(i) elseif (ed.eq.'l3') then s2e = el3(i) endif 99 continue return c end function s2e end c********************************************************************** integer function is2z(sym) c--------------------------------------------------------------------------- c copyright 1993 university of washington matt newville and bruce ravel c--------------------------------------------------------------------------- c returns z number given atomic symbol: default is 0 character*2 symbol(103), sym, sy, test c data (symbol(i),i=1,103) /'h','he','li','be','b','c','n','o', $'f','ne','na','mg','al','si','p','s','cl','ar','k','ca','sc', $'ti','v','cr','mn','fe','co','ni','cu','zn','ga','ge','as','se', $'br','kr','rb','sr','y','zr','nb','mo','tc','ru','rh','pd','ag', $'cd','in','sn','sb','te','i','xe','cs','ba','la','ce','pr','nd', $'pm','sm','eu','gd','tb','dy','ho','er','tm','yb','lu','hf','ta', $'w','re','os','ir','pt','au','hg','tl','pb','bi','po','at','rn', $'fr','ra','ac','th','pa','u','np','pu','am','cm','bk','cf','es', $'fm','md','no','lr'/ c test case of this routine : note that 'ab' must be the same c case as all of the symbols above sy = sym test = 'ab' call case(test,sy) is2z = 0 do 110 i=1,103 if (sy.eq.symbol(i)) then is2z = i goto 120 end if 110 continue 120 return c end integer function is2z end c******************************************************************** subroutine fixsym(sym) c returns a word with the first letter capitalized and the remaining c letters in lower case c this is very useful for writing out atomic symbols character*2 toss, sym*(*) toss = sym call upper(toss(1:1)) call lower(toss(2:2)) ii = istrln(sym) if (ii.gt.2) then call lower(sym(3:ii)) sym = toss(1:1)//toss(2:2)//sym(3:ii) else sym = toss(1:1)//toss(2:2) endif return c end subroutine fixsym end function dfridr(func, x, h, err) integer ntab real dfridr, err, h, x, func, con, con2, big, safe parameter(con=1.4, con2=con*con, big=1.e30, ntab=10, safe=2.) external func integer i, j real errt, fac, hh, a(ntab,ntab) if (h.eq.0) then print*,'h must be non-zero in dfridr' stop endif hh = h a(1,1) = (func(x+hh) - func(x-hh)) / (2.0*hh) err = big do 120 i=2,ntab hh = hh/con a(1,i) = (func(x+hh) - func(x-hh)) / (2.0*hh) fac = con2 do 110 j=2,i a(j,i) = (a(j-1,i)*fac - a(j-1,i-1)) / (fac-1.0) fac = con2*fac errt = max( abs(a(j,i)-a(j-1,i)), abs(a(j,i)-a(j-1,i-1)) ) if (errt.lt.err) then err = errt dfridr = a(j,i) endif 110 continue if (abs(a(i,i)-a(i-1,i-1)).ge.safe*err) return 120 continue return c end function dfridr end function xmuint(x) c interpolate to find the value of xmu at any energy parameter(nptx=2**11) dimension energy(nptx), xmuout(nptx) common /idata/ nxmu save /idata/ common /rdata/ e0, energy, xmuout save /rdata/ nterms = 3 call interp(energy,xmuout,nxmu,nterms,x,y) xmuint = y return c end function xmuint end function istrln (string) c c returns index of last non-blank character. c returns zero if string is null or all blank. character*(*) string c-- if null string or blank string, return length zero. istrln = 0 if (string(1:1).eq.char(0)) return if (string.eq.' ') return c c-- find rightmost non-blank character. ilen = len (string) do 20 i = ilen, 1, -1 if (string (i:i) .ne. ' ') goto 30 20 continue 30 istrln = i return c end function istrln end subroutine triml (string) c removes leading blanks. character*(*) string jlen = istrln(string) c c-- all blank and null strings are special cases. if (jlen .eq. 0) return c-- find first non-blank char do 10 i = 1, jlen if (string (i:i) .ne. ' ') goto 20 10 continue 20 continue c-- if i is greater than jlen, no non-blanks were found. if (i .gt. jlen) return c-- remove the leading blanks. string = string (i:) return c end subroutine triml end subroutine untab (string) c replace tabs with blanks : tab is ascii dependent integer itab , i, ilen parameter (itab = 9) character*(*) string, tab*1 tab = char(itab) ilen = max(1, istrln(string)) 10 continue i = index(string(:ilen), tab ) if (i .ne. 0) then string(i:i) = ' ' go to 10 end if return c end subroutine untab end subroutine uncomm (string) c replace commas with blanks integer i, ilen character*(*) string, comma*1 data comma / ',' / ilen = max(1, istrln(string)) 10 continue i = index(string(:ilen), comma ) if (i .ne. 0) then string(i:i) = ' ' go to 10 end if return c end subroutine uncomm end subroutine upper (str) c changes a-z to upper case. ascii specific c- for ascii: ichar(upper case 'a') = 65 c- ichar(lower case 'a') = 97 character*(*) str integer iupa, iloa, iloz, idif data iupa, iloa / 65, 97/ idif = iloa - iupa iloz = iloa + 25 jlen = max(1, istrln (str) ) do 10 i = 1, jlen ic = ichar (str(i:i)) if ((ic.ge.iloa).and.(ic.le.iloz)) str(i:i) = char(ic-idif) 10 continue return c end subroutine upper end subroutine lower (str) c changes a-z to lower case. ascii specific c- for ascii: ichar(upper case 'a') = 65 c- ichar(lower case 'a') = 97 character*(*) str integer iupa, iloa, iupz, idif data iupa, iloa / 65, 97/ idif = iloa - iupa iupz = iupa + 25 jlen = max(1, istrln (str) ) do 10 i = 1, jlen ic = ichar (str(i:i)) if ((ic.ge.iupa).and.(ic.le.iupz)) str(i:i) = char(ic+idif) 10 continue return c end subroutine lower end subroutine smcase (str, contrl) c convert case of string *str*to be the same case c as the first letter of string *contrl* c if contrl(1:1) is not a letter, *str* will be made lower case. character*(*) str, contrl, s1*1, t1*1 s1 = contrl(1:1) t1 = s1 call lower(t1) if (t1.eq.s1) call lower(str) if (t1.ne.s1) call upper(str) return c end subroutine smcase end subroutine case(test,word) c returns *word* in the same case as *test* c note that this is just the reverse of smcase ! character*(*) test, word call smcase (word, test) return c end subroutine case end logical function isnum (string) c returns true if string can be a number, else returns false c recognizes e and d exponentials, but is not foolproof c to be a number, a string must contain: c - only characters in 'de.+-, 1234567890' (case is checked) c - no more than one 'd' or 'e' c - no more than one '.' character*(*) string, str*70, number*20 integer iexp, idec, i, ilen, ier, j, istrln real x external istrln c note: layout of *number* is important: don't change this!! data number /'de.+-, 1234567890'/ c- isnum = .false. iexp = 0 idec = 0 str = string ilen = max(1, istrln (str) ) call smcase(str, number ) do 100 i = 1, ilen j = index(number,str(i:i) ) if (j.le.0) go to 200 if((j.eq.1).or.(j.eq.2)) iexp = iexp + 1 if (j.eq.3) idec = idec + 1 100 continue c every character in the string was found in *number* c so the string probably is a number isnum = .true. c but let's do a few more tests: c number of exponential and decimal markers if (iexp.ge.2) isnum = .false. if (idec.ge.2) isnum = .false. c read with iostat (this may report an error, but not always) read(str,150,iostat=ier) x 150 format (bn,f70.0) if (ier.ne.0) isnum = .false. c all tests done 200 continue return c end logical function isnum end subroutine bwords (s, nwords, words) c c breaks string into words. words are seperated by one or more c blanks, or a comma or equal sign and zero or more blanks. c c args i/o description c ---- --- ----------- c s i char*(*) string to be broken up c nwords i/o input: maximum number of words to get c output: number of words found c words(nwords) o char*(*) words(nwords) c contains words found. words(j), where j is c greater then nwords found, are undefined on c output. c c written by: steven zabinsky, september 1984 c c************************** deo soli gloria ************************** c-- no floating point numbers in this routine. implicit integer (a-z) character*(*) s, words(nwords) character blank, comma, equal parameter (blank = ' ', comma = ',', equal = '=') c-- betw .true. if between words c comfnd .true. if between words and a comma or equal has c already been found logical betw, comfnd c-- maximum number of words allowed wordsx = nwords c-- slen is last non-blank character in string slen = istrln (s) c-- all blank string is special case if (slen .eq. 0) then nwords = 0 return endif c-- begc is beginning character of a word begc = 1 nwords = 0 betw = .true. comfnd = .true. do 10 i = 1, slen if (s(i:i) .eq. blank) then if (.not. betw) then nwords = nwords + 1 words (nwords) = s (begc : i-1) betw = .true. comfnd = .false. endif elseif ((s(i:i).eq.comma).or.(s(i:i).eq.equal)) then if (.not. betw) then nwords = nwords + 1 words (nwords) = s(begc : i-1) betw = .true. elseif (comfnd) then nwords = nwords + 1 words (nwords) = blank endif comfnd = .true. else if (betw) then betw = .false. begc = i endif endif if (nwords .ge. wordsx) return 10 continue c if (.not. betw .and. nwords .lt. wordsx) then nwords = nwords + 1 words (nwords) = s (begc :slen) endif return c end subroutine bwords end subroutine strclp(str,str1,str2,strout) c c a rather complex way of clipping a string: c strout = the part of str that begins with str2. c str1 and str2 are subsrtings of str, (str1 coming before str2), c and even if they are similar, strout begins with str2 c for example: c 1. str = "title title my title" with str1 = str2 = "title" c gives strout = "title my title" c 2. str = "id 1 1st path label" with str1 = "1", str2 = "1st" c gives strout = "1st path label" c character*(*) str, str1, str2, strout integer i1, i2, ibeg, iend, istrln, ilen external istrln ilen = len(strout) i1 = max(1, istrln(str1)) i2 = max(1, istrln(str2)) i1e = index(str,str1(1:i1)) + i1 ibeg = index(str(i1e:),str2(1:i2) ) + i1e - 1 iend = min(ilen+ibeg, istrln(str) ) strout = str(ibeg:iend) return c end subroutine strclp end subroutine cabort(messg, abortf) c c copyright university of washington 1981 c part of uwxafs binary (rdf) file handling c conditional abort; if abortf is .true. c character*(*) messg logical abortf call messag(messg) if (abortf) then call messag('* uwxafs data file handling abort *') stop endif return c end subroutine cabort end subroutine getdoc(iounit,doc,nl,skey,nkey,ntl,ier) c c copyright university of washington 1981 c part of uwxafs binary (rdf) file handling c c get documentation lines from data file c c input parameters c iounit : i/o unit number (integer,input) c doc : document array (character,in-out) c nl : no of lines to get (integer,input) c skey : symbolic key of record (character,input) c nkey : numeric key of record (integer,input) c ntc : number of lines actually got(integer,output) c ier : error code (integer,output) c 1 - unit not declared c 2 - nl negative or zero c 3 - nkey .gt. indxl .or. nkey .lt. 1 c 4 - nkey is not on file c 5 - skey and nkey don't match c 6 - skey is not on file c c if skey is blank, nkey is used. if nkey is 0, skey is used. c if both are given, they should match. c c if doc is equal to 'doc' then data are transferred c to internal buffer. nl is ignored in this case. c see routines getline, addline and prntdoc for c the manipulation of this buffer c implicit integer(a-z) parameter (indxl = 191, nu = 2 ) c character*80 fname(nu), doctmp*100 character*2048 cindx(nu) integer*2 indx(4,0:indxl,nu) logical abortf, safe, rewrt, modify(nu) c common /uwdata/ unit(nu),modify,abortf,safe,rewrt,nldoc,indx common /uwdocs/ cindx,fname c save /uwdata/,/uwdocs/ c parameter (maxl=20, maxchr=100 ) character*(maxchr) dbuf(maxl) common /uwdbuf/ dbuf save /uwdbuf/ c character*(*) skey,doc(*) c call gunit(iounit,u) if(u.eq.0) then call cabort('getdoc: unit not declared',abortf) ier=1 return endif c c convert form to the case of this routine. c 'case' controls the the case of this routine doctmp = doc(1) call smcase(doctmp, 'case') c if ((nl.le.0).and.(doctmp.ne.'doc')) then call cabort('getdoc: no document lines to get',abortf) ier=2 return endif c if(nkey.ne.0) then if(nkey.lt.0.or.nkey.gt.indxl) then call cabort('getdoc: nkey out of bounds',abortf) ier=3 return elseif(indx(1,nkey,u).eq.0) then call cabort('getdoc: nkey does not exist',abortf) ier=4 return elseif(skey.ne.' '.and. $ skey.ne.cindx(u)(nkey*10+1:nkey*10+10)) then call cabort('getdoc: skey mismatch',abortf) ier=5 return endif iord=nkey c skey is not given else call gnkey(iounit,skey,iord,ier) if(iord.eq.0) then call cabort('getdoc: skey not found',abortf) ier=6 return endif endif c pru=indx(3,iord,u) c c to document buffer if (doctmp.eq.'doc') then nldoc=indx(4,iord,u) nblk=(nldoc+4)/5 ntl=nldoc do 10 i=1,nblk read(iounit,rec=pru+i-1) (dbuf(j),j=i*5-4,i*5) 10 continue c to doc else ntl=indx(4,iord,u) ntl=min(ntl,nl) nblk=(ntl+4)/5 do 20 i=1,nblk lblk=min(i*5,ntl) read(iounit,rec=pru+i-1) (doc(j),j=i*5-4,lblk) 20 continue endif ier=0 return c end subroutine getdoc end subroutine getrec(iounit,array,nw,skey,nkey,ntw,ier) c c copyright university of washington 1981 c part of uwxafs binary (rdf) file handling c c get data record from a data file c iounit : i/o unit number (integer,input) c array : array to receive data (any type,output) c nw : maximum no of words to receive(integer,input) c skey : symbolic key of data record to get(character,input) c nkey : numeric key of data record to get(integer,input) c ntw : actual no of words received(integer,output) c ier : error code(integer,output) c 1 - unit not declared c 2 - nw .lt. 0 c 3 - nkey .gt. indxl .or. .lt. 0 c 4 - nkey not on file c 5 - two keys do not match c 6 - skey not on file c either nkey(with skey=' ') or skey(with nkey=0) may be c given. if both are given, they should match. c implicit integer(a-z) c parameter (indxl=191, nu=2, iword=128 ) c character*(*) skey character*80 fname(nu) character*2048 cindx(nu) integer*2 indx(4,0:indxl,nu) logical abortf,safe,rewrt,modify(nu) c common /uwdata/ unit(nu),modify,abortf,safe,rewrt,nldoc,indx common /uwdocs/ cindx,fname c save /uwdata/,/uwdocs/ c real array(nw) call gunit(iounit,u) if(u.eq.0) then call cabort('getrec: iounit not declared',abortf) ier=1 return endif c if(nw.le.0) then call cabort('getrec: no data to get',abortf) ier=2 return endif c last=indx(2,0,u) if(nkey.ne.0) then if(nkey.lt.0.or.nkey.gt.indxl) then call cabort('getrec: nkey out of bounds',abortf) ier=3 return elseif(nkey.gt.last) then call cabort('getrec: nkey does not exist',abortf) ier=4 return elseif((skey.ne.' ').and.(skey.ne.cindx(u)(nkey*10+1: $ nkey*10+10))) then call cabort('getrec: skey mismatch',abortf) ier=5 return endif iord=nkey c c skey is given c else call gnkey(iounit,skey,iord,iii) if(iord.eq.0) then call cabort('getrec: skey not found',abortf) ier=6 return endif endif c pru=indx(1,iord,u) ntw=indx(2,iord,u) if(ntw.gt.nw) then ntw=nw call messag('getrec: field shorter than data') endif c c iword is no of words in a block nblk = ( ntw + iword - 1) / iword do 10 i = 1, nblk l = min(i*iword, ntw) read(iounit,rec=pru+i-1)(array(j),j=(i-1)*iword+1,l) 10 continue ier=0 return c end subroutine getrec end subroutine gunit(iounit,u) c c copyright university of washington 1981 c part of uwxafs binary (rdf) file handling c c match unit no with iounit no c iounit : fortran i/o unit no(integer,input) c u : corresponding index(1,2.. upto nu, in order of c call to openrf routine)(integer,output) c relation unit(u)=iounit implicit integer(a-z) c parameter (indxl=191, nu=2 ) c character*80 fname(nu) character*2048 cindx(nu) integer*2 indx(4,0:indxl,nu) logical abortf,safe,rewrt,modify(nu) c common /uwdata/ unit(nu),modify,abortf,safe,rewrt,nldoc,indx common /uwdocs/ cindx,fname c save /uwdata/,/uwdocs/ c u=0 do 10 n=1,nu if(unit(n).eq.iounit) then u = n go to 20 end if 10 continue 20 continue return c end subroutine gunit end subroutine hash(array, narray, doc, ndoc, skey) c c copyright university of washington 1981 c part of uwxafs binary (rdf) file handling c c generate 5 character string composed of alphanumerics c c since this routine is not called too often, execution speed c is not as important as getting a reasonably pseudo-random skey. c *** note: random numbers done as in numerical recipes c----------------------------------------------------------------------- parameter(maxpts = 4096) parameter(imod = 6655, imul = 936, iadd = 1399) parameter(jmod = 14406, jmul = 967, jadd = 3041) parameter(kmod = 7875, kmul = 211, kadd = 1663) parameter(lmod = 6075, lmul = 106, ladd = 1283) parameter(ihalf= 3333, imost= 5555, jhalf= 7201, jthrd= 4821) c double precision work(maxpts), sum(4) character*(*) doc(*) character*(5) skey real array(*) integer ikey(5), iran(131) logical loop c c initialize nwork = 2*narray skey = 'skey0' if (narray.le.0) return ichr0 = ichar('0') ichra = ichar('a') do 15 i = 1, 5 ikey(i) = 0 15 continue do 18 i = 1, 4 sum(i) = 0.0 18 continue c get measure of the magnitude of the data two different ways c ihalf ~= imod/2 , so that roughly half the values c are used to make the partial sum aver = 0.1 part = 0.2 do 30 i = 1, narray aver = abs( array(i) / narray ) + aver irndm = mod(imul*(i + narray) + iadd, imod) if ( (i.gt.2) .and. (irndm.gt.ihalf) ) then part = abs( array(i) / narray ) + part if (irndm.gt.imost) then part = abs( array(i-2) / narray ) + part else part = abs( array(i-1) / narray ) + part end if end if 30 continue c create work array such that all values of work are of the order 1. c - most values are scaled to one of the two different measures of c magnitude from above. c - a few values get scaled to be on the order of 1 without reference c to these values (should prevent against two arrays differing only c by a constant factor from having the same skey) c - couldn't resist the golden mean and fine structure constant. if (abs(aver).le.0.001) aver = 0.01 if (abs(part).le.0.001) part = 0.01 do 150 i = 1, narray loop = .true. work(i) = abs( array(i) / aver ) work(i+narray) = abs( array(i) / part ) j = i*kmul + kadd + narray jrndm = mod(jmul*j + jadd, jmod) if ( jrndm.lt.jthrd ) then work(i) = work(i) * 3.14159265 elseif ( jrndm.gt.jhalf ) then j = mod(i*j + jrndm, narray - 1 ) + 1 work(i+narray) = abs( (array(i)+array(j)) / 2.71828183) 100 continue if(loop.and.(work(i+narray).le.(0.0072974))) then work(i+narray) = work(i+narray) * 1.618033989 loop = .false. go to 100 elseif(loop.and.(work(i+narray).ge.( 137.036) )) then work(i+narray) = work(i+narray) * 0.618033989 loop = .false. go to 100 endif end if 150 continue c c generate iran: a list of 131 pseudo-random integers that c do not depend on the data at all do 200 i = 1, 131 j = i*imul + iadd jtmp = mod(jmul*j + jadd, jmod) + 1 ktmp = mod(kmul*jtmp + kadd, kmod) + 3 iran(i) = mod(lmul*ktmp + ladd, lmod) + 7 200 continue c c collect 4 different sums from the work array: c each value in the work array is multiplied by a pseudo-randomly c selected integer between 20 and 120, and 4 different sums are made. c since each value in work() is on the order of 1, each of the sums c should be of the order 100*narray. with narray being something c between 100 and 1000, each of the different sums will be a random c number on the order of 50 000. this value mod 36 should be a good c random number. c do 500 i = 1, nwork c get some random numbers, and make them bigger than 50 i1 = mod (i * imul + iadd, imod ) + 53 i2 = mod (i * jmul + jadd, jmod ) + 67 i3 = mod (i * kmul + kadd, kmod ) + 31 i4 = mod (i * lmul + ladd, lmod ) + 79 c use these to make random numbers between [1, 130 ] j1 = mod( jmul*i1 + kadd, 109) + 3 j2 = mod( kmul*i2 + ladd, 119) + 5 j3 = mod( lmul*i3 + iadd, 111) + 7 j4 = mod( imul*i4 + jadd, 123) + 1 c use these for the iran array of random numbers to get a set c of numbers between [20 and 150] k1 = mod( jmul*( i4 + iran(j1)) + kadd, 73 ) + 43 k2 = mod( kmul*( i2 + iran(j2)) + ladd, 111 ) + 37 k3 = mod( lmul*( i1 + iran(j3)) + iadd, 91 ) + 29 k4 = mod( imul*( i3 + iran(j4)) + jadd, 121 ) + 19 c do "randomly weighted" sum of work array sum(1) = sum(1) + work(i) * k1 sum(2) = sum(2) + work(i) * k2 sum(3) = sum(3) + work(i) * k3 sum(4) = sum(4) + work(i) * k4 500 continue c turn the sums to integers between 1 and 36 for ikey(1) - ikey(4) do 900 i = 1, 4 880 continue if (abs(sum(i)).ge.100 000 000) then sum(i) = sum(i) / 1.618033989 go to 880 end if isum = int( sum(i) ) ikey(i) = mod(isum, 36) 900 continue c ikey(5) : sum from document array isum = 0 im = mod(iran(16) * ndoc + iran(61), 353) + 27 ia = mod(iran(77) * ndoc + iran(52), 347) + 19 do 2000 i = 1, ndoc call triml( doc(i) ) jlen = max(1, istrln( doc(i))) do 1800 j = 1, jlen kseed = mod( (j + 2*i) * imul + jadd , 127) + 1 k = mod( iran(kseed) * im + ia , 13) + 3 isum = isum + k * ichar( doc(i)(j:j) ) 1800 continue 2000 continue ikey(5) = mod(isum, 36) c c map integers 1 to 36 to numerals and letters c ascii assumed but not required. the numerals must be c ordered 0 - 9 and the letters must be ordered a - z. do 4000 i = 1, 5 if (ikey(i).le.9) then ikey(i) = ikey(i) + ichr0 else ikey(i) = ikey(i) - 10 + ichra end if 4000 continue c write skey from ikey do 5000 i = 1, 5 skey(i:i) = char( ikey(i) ) 5000 continue call upper(skey) return c end subroutine hash end subroutine openrf(iounit,lfn,aflag,sflag,ftype,irecl,ier) c c copyright university of washington 1981 c part of uwxafs binary (rdf) file handling c c openrf uses direct access binary files with word size irecl. c irecl = 128 on vax, 512 otherwise? c c structure of random data file c c block size=512 byte=128 word c c block 1-3 : indx(4,0:indxl,u) c block 4-7 : cindx(u)*2048 c c first data block is 8 c c indx(1,0,u) address of eof (non exisiting) c indx(2,0,u)=no of entries c indx(3,0,u)=1776 for identification c indx(4,0,u)=704 same purpose c c indx(1,n,u)=address of data n c indx(2,n,u)=no of words for data n c indx(3,n,u)=address of doc n c indx(4,n,u)=no of lines for doc n c c cindx(u)(1:10)=ftype c cindx(u)(n*10+1:n*10+10)=skey for nkey n c--------------------------- implicit integer(a-z) parameter (indxl = 191, nu = 2) character*(*) ftype,lfn,aflag,sflag character*80 fname(nu),fn, aflg*10, sflg*10 character*2048 cindx(nu) integer*2 indx(4,0:indxl,nu) logical abortf,safe,rewrt,exist logical clor, modify(nu) c common /uwdata/ unit(nu),modify,abortf,safe,rewrt,nldoc,indx common /uwdocs/ cindx,fname c save /uwdata/,/uwdocs/ c ccc data unit,nldoc/nu*0,0/ ccc data abortf,safe/.true.,.false./ c clor =.false. c find out the case of this routine. c 'case' controls the the case of this routine sflg = sflag aflg = aflag call smcase(aflg, 'case') call smcase(sflg, 'case') c abortf = (aflg.ne.'noabort') safe = safe.or.(sflg.eq.'safe') c c if lfn is blank, lfn=for0un where un is fortran i/o no if(lfn.eq.' ') then fn = 'for0' write(fn(5:6),1,err=9999)iounit 1 format(i2.2) else fn=lfn endif inquire (file=fn,exist=exist) c call gunit(iounit,u) if(u.eq.0) then c assign iounit no to unit(n), a table do 100 n=1,nu if(unit(n).eq.0) then unit(n)=iounit u=n fname(u)=fn go to 110 endif 100 continue c call cabort('openrf: max no of files exceeded',abortf) ier=1 return 110 continue c else call messag('openrf: unit reopened') endif c if(ftype.eq.' ') then c no modify permit modify(u)=.true. else modify(u)=.false. endif c if(exist) then c file exists if(.not.modify(u)) then c can modify the file open(iounit, file=fn, recl=irecl, access='direct', $ status='old', iostat=iosb, err=9999) else c cannot modify the file open(iounit, file=fn, recl=irecl, access='direct', $ status='old') cccccc $ ,readonly) endif c read in existing index do 10 i=1,3 read(iounit,rec=i)((indx(k,l,u),k=1,4),l=i*64-64,i*64-1) 10 continue do 20 i=1,4 read(iounit,rec=i+3)cindx(u)(i*512-511:i*512) 20 continue c if(indx(3,0,u).ne.1776 .or. indx(4,0,u).ne.704) then call cabort('openrf: wrong file',abortf) ier=4 return endif ier=0 c if(ftype.ne.' '.and.ftype.ne.cindx(u)(1:10))then call cabort('openrf: wrong file type',abortf) ier=3 return endif return c c new file c else if(ftype.eq.' ') then call cabort('openrf: ftype needed for file creation',abortf) ier=5 return endif c create a new file open(iounit,file=fn,recl=irecl,access='direct', x status='new') cindx(u)=ftype c index initialization do 30 i=1,indxl do 30 j=1,4 indx(j,i,u)=0 30 continue indx(1,0,u)=8 indx(2,0,u)=0 indx(3,0,u)=1776 indx(4,0,u)=704 ier=0 go to 99 endif c entry wrindx(iounit) c write out index call gunit(iounit,u) 99 continue do 101 i=1,3 write(iounit,rec=i,err=9999) $ ((indx(j,k,u),j=1,4),k=i*64-64,i*64-1) 101 continue do 111 i=1,4 write(iounit,rec=i+3,err=9999)cindx(u)(i*512-511:i*512) 111 continue if(clor) go to 77 9998 return 9999 ier=iosb go to 9998 c entry closrf(iounit,ier) c call gunit(iounit,u) if(u.eq.0) then call cabort('openrf: unit not declared',abortf) ier=1 return endif c reset i/o unit table unit(u)=0 ier=0 c if modified, rewrite index clor=.false. if(modify(u)) then clor=.true. go to 99 endif 77 continue close(iounit) clor=.false. return c entry rwrtrf(iounit,ier) c rewrite interlock. safeguard. call gunit(iounit,u) if(u.eq.0) then call cabort('openrf: unit not declared',abortf) ier=1 return endif c rewrt=.true. return c end subroutine openrf end subroutine putdoc(iounit,doc,nl,skey,nkey,ier) c c copyright university of washington 1981 c part of uwxafs binary (rdf) file handling c c put documentation lines to a random access data file c iounit : fortran i/o unit no(integer,input) c doc : document lines(character,input) c nl : no of lines to be written(integer,input) c skey : symbolic key associated with the record(char,in) c nkey : numeric key associated with the record(integer,in) c ier : error code (integer,output) c 1 - unit not declared c 2 - nl .le. 0 c 3 - c 4 - skey not given c 5 - rewrite interlock not cleared c 6 - skey not found c if doc='doc' , internal buffer is used c symbolic key must be given c for new record, nkey should be zero. c implicit integer(a-z) c parameter (indxl=191, nu=2 ) c character*(*) skey character*80 fname(nu), doctmp*100 character*2048 cindx(nu) integer*2 indx(4,0:indxl,nu) logical abortf,safe,rewrt,modify(nu) c common /uwdata/ unit(nu),modify,abortf,safe,rewrt,nldoc,indx common /uwdocs/ cindx,fname c save /uwdata/,/uwdocs/ c parameter(maxl=20, maxchr=100 ) character*(maxchr) dbuf(maxl) common /uwdbuf/ dbuf save /uwdbuf/ c character*(*) doc(*) c call gunit(iounit,u) if(u.eq.0) then call cabort('putdoc: unit not declared',abortf) ier=1 return endif c c convert form to the case of this routine. c 'case' controls the the case of this routine doctmp = doc(1) call smcase(doctmp, 'case') c if ((nl.le.0).and.(doctmp.ne.'doc')) then call cabort('putdoc: no documents to write',abortf) ier=2 return endif c if(skey.eq.' ') then call cabort('putdoc: skey must be given',abortf) ier=4 return endif c last=indx(2,0,u) c if(nkey.ne.0) then c existing record if(.not.rewrt) then call cabort('putdoc: rewrite interlock not cleared',abortf) ier=5 return endif rewrt=.false. iord=nkey else c new record call gnkey(iounit,skey,iord,ier) if(iord.eq.0) then call cabort('putdoc: skey not found',abortf) ier=6 return endif endif c c write at the end, always pru=indx(1,0,u) if (doctmp.eq.'doc') then c write internal buffer indx(4,iord,u)=nldoc nblk=(nldoc+4)/5 do 10 i=1,nblk lblk=min(i*5,nldoc) write(iounit, rec=pru+i-1) (dbuf(j),j=i*5-4,lblk) 10 continue else c write doc indx(4,iord,u)=nl nblk=(nl+4)/5 do 20 i=1,nblk lblk=min(i*5,nl) write(iounit, rec=pru+i-1)(doc(j),j=i*5-4,lblk) 20 continue endif c adjust index indx(3,iord,u)=pru indx(1,0,u)=pru+nblk c if (safe) write out new index if(safe) then call wrindx(iounit) else c otherwise, mark it c new index will be written when closrf is called modify(u)=.true. endif ier=0 return c end subroutine putdoc end subroutine putrec(iounit,array,nw,skey,nkey,ier) c c copyright university of washington 1981 c part of uwxafs binary (rdf) file handling c c put data record to a random access data file c c iounit : i/o unit no(integer,input) c array : data to be put(any type,input) c nw : no of words in array(integer,input) c skey : symbolic key for the record(character,input) c nkey : numeric key for the record(integer,input) c ier : error code(integer,output) c 1 - unit not declared c 2 - nw .lt. 0 c 3 - file protection violated c 4 - skey is blank c 5 - rewrite interlock not cleared c 6 - record length is different for rewrite c 7 - nkey does not exist c 8 - index full c c skey must be given always and nkey should be zero for new record c non zero nkey means rewrite. - nw should be equal to current one c implicit integer(a-z) c parameter (indxl=191, nu=2, iword = 128 ) c character*(*) skey character*80 fname(nu),ctmp character*2048 cindx(nu) integer*2 indx(4,0:indxl,nu) logical abortf,safe,rewrt,modify(nu) c common /uwdata/ unit(nu),modify,abortf,safe,rewrt,nldoc,indx common /uwdocs/ cindx,fname c save /uwdata/,/uwdocs/ c real array(nw) c call gunit(iounit,u) if(u.eq.0) then call cabort('putrec: unit not declared',abortf) ier=1 return endif c if(nw.lt.0) then call cabort('putrec: no data to write',abortf) ier=2 return endif c if(modify(u)) then call cabort('putrec: '//fname(u)// $ ' has no write permission',abortf) ier=3 return endif c if(skey.eq.' ') then call cabort('putrec: skey not given',abortf) ier=4 return endif c last=indx(2,0,u) if(nkey.ne.0) then c rewrite if(.not.rewrt) then call cabort('putrec: rewrite interlock not cleared',abortf) ier=5 return endif rewrt=.false. c rewrite if(nkey.gt.last) then call cabort('putrec: old nkey does not exist',abortf) ier=7 return endif if(nw.ne.indx(2,nkey,u)) then call cabort('putrec: old/new record length dont match',abortf) ier=6 return endif pru=indx(1,nkey,u) c c new record else if(last.ge.indxl) then call cabort('putrec: index full',abortf) ier=8 return endif pru=indx(1,0,u) endif c nblk = (nw + iword - 1)/iword do 10 i = 1 , nblk l = min(i*iword, nw) write(iounit, rec=i+pru-1)(array(j),j=(i-1)*iword+1,l) 10 continue c c new index for rewrite if(nkey.ne.0) then call messag('putrec: symbolic key overwritten'// $ cindx(u)(nkey*10+1:nkey*10+10)) cindx(u)(nkey*10+1:nkey*10+10)=skey else c new index for new record cindx(u)(last*10+11:last*10+20)=skey indx(1,last+1,u)=indx(1,0,u) indx(2,last+1,u)=nw endif c new no of entries and eof block no indx(1,0,u)=indx(1,0,u)+nblk indx(2,0,u)=last+1 ier=0 c check duplicate key ctmp=skey ntmp=nkey 100 continue do 500 n=1,last if(ctmp.ne.cindx(u)(n*10+1:n*10+10)) go to 500 if(n.ne.ntmp) go to 510 500 continue c no duplicate key go to 600 c put an asterisk, if one is not there already 510 continue do 520 i=6,10 if(cindx(u)(n*10+i:n*10+i).eq.'*') go to 520 cindx(u)(n*10+i:n*10+i)='*' c check again if new skey is duplicate ctmp=cindx(u)(n*10+1:n*10+10) ntmp=n go to 100 520 continue call cabort('putrec: same skey occured five times',abortf) 600 continue if(safe) then call wrindx(iounit) else modify(u)=.true. endif return c end subroutine putrec end subroutine rfmisc c c copyright university of washington 1981 c part of uwxafs binary (rdf) file handling c c miscellaneous routines for handling uwexafs files. most of c the entries here are to find out what's inside a file. c copyright university of washington 1981 c implicit integer(a-z) c parameter (indxl=191, nu=2) c character*(*) skey,ftype,lfn character*80 fname(nu) character*2048 cindx(nu) integer*2 indx(4,0:indxl,nu) logical abortf,safe,rewrt,modify(nu) c common /uwdata/ unit(nu),modify,abortf,safe,rewrt,nldoc,indx common /uwdocs/ cindx,fname c save /uwdata/,/uwdocs/ c entry gskey(iounit,nkey,skey,ier) c c get symbolic key from a numeric key c call gunit(iounit,u) if(u.eq.0) then call cabort('gskey: unit not declared',abortf) ier=1 return endif if(nkey.lt.0.or.nkey.gt.indxl) then call cabort('gskey: nkey out of range',abortf) ier=2 return endif c if nkey does not exist, skey=' ' skey=cindx(u)(nkey*10+1:nkey*10+10) ier=0 return c entry gnkey(iounit,skey,nkey,ier) c c get a numeric key from a symbolic key c call gunit(iounit,u) if(u.eq.0) then call cabort('gnkey: unit not declared',abortf) ier=1 return endif c last=indx(2,0,u) do 300 n=1,last if(skey.eq.cindx(u)(n*10+1:n*10+10)) go to 310 300 continue c skey not found nkey=0 return c 310 continue nkey=n ier=0 return c entry gftype(iounit,ftype,ier) c c get file-type from a i/o unit no c call gunit(iounit,u) if(u.eq.0) then call cabort('gftype: iounit not declared',abortf) ier=1 return endif c ftype=cindx(u)(1:10) ier=0 return c entry glfn(iounit,lfn,ier) c c get filename from i/o unit no c call gunit(iounit,u) if(u.eq.0) then call cabort('gfln: iounit not declared',abortf) ier=1 return endif c lfn=fname(u) ier=0 return c entry gflen(iounit,flen,ier) c c get the length of a file from i/o unit no c call gunit(iounit,u) if(u.eq.0) then call cabort('gflen: iounit not declared',abortf) ier=1 return endif c flen=indx(1,0,u)-1 return c entry gnie(iounit,nie,ier) c c get number of entries from i/o unit no c call gunit(iounit,u) if(u.eq.0) then call cabort('gnie: iounit not declared ',abortf) ier=1 return endif c nie=indx(2,0,u) ier=0 return c entry grlen(iounit,nkey,rlen,ier) c c get the length of a data record (in words) from a numeric key c call gunit(iounit,u) if(u.eq.0) then call cabort('grlen: iounit not declared ',abortf) ier=1 return endif c if(nkey.lt.0.or.nkey.gt.indxl) then call cabort('grlen: nkey out of range',abortf) ier=2 return endif c rlen=indx(2,nkey,u) ier=0 return c entry gdlen(iounit,nkey,dlen,ier) c c get the length of document (in lines) from a numeric key c call gunit(iounit,u) if(u.eq.0) then call cabort('gdlen: unit not declared',abortf) ier=1 return endif c if(nkey.lt.0.or.nkey.gt.indxl) then call cabort('gdlen: nkey out of range',abortf) ier=2 return endif c dlen=indx(4,nkey,u) ier=0 return c c end subroutine rfmisc end block data uwbdat c c copyright university of washington 1981 c part of uwxafs binary (rdf) file handling c c block data statements for uwxafs c implicit integer(a-z) parameter (indxl = 191, nu = 2) integer*2 indx(4,0:indxl,nu) logical abortf,safe,rewrt logical modify(nu) c common /uwdata/ unit(nu),modify,abortf,safe,rewrt,nldoc,indx save /uwdata/ data unit,nldoc/nu*0,0/ data abortf,safe/.true.,.false./ c end block data end subroutine addocd(ndoc, doc, commnt, maxder, emader) parameter(ndocx=25) character*78 commnt character*100 doc(ndocx) real maxder c If file comment is given, need to make room for 2 new lines; if not c just need one document line. Copy line i to line i+2 (or i+1). The c first line of doc in not altered. nnd = 2 if (commnt.ne.' ') then do 10 i=min(ndoc, ndocx-nnd), nnd, -1 doc(i+2) = doc(i) 10 continue write(doc(nnd),400)commnt(2:) 400 format('fluo: ',a) write(doc(1+nnd),410)maxder, emader 410 format('fluo: Max. deriv = ',f7.4', energy max. der. = ', $ f8.2) else do 20 i=min(ndoc, ndocx-nnd), nnd, -1 doc(i+1) = doc(i) 20 continue write(doc(nnd),410)maxder, emader endif return c end subroutine addocd end subroutine interp(x,y,npts,nterms,xin,yout) c interpolation routine from bevington's book double precision deltax,delta,a,prod,sum dimension x(*),y(*) dimension delta(10),a(10) c search for appropriate value of x(1) do 19 i=1,npts if(xin-x(i)) 13,17,19 13 i1=i-nterms/2 if(i1) 15,15,21 15 i1=1 go to 21 17 yout=y(i) 18 go to 61 19 continue i1=npts-nterms+1 21 i2=i1+nterms-1 if(npts-i2) 23,31,31 23 i2=npts i1=i2-nterms+1 25 if(i1) 26,26,31 26 i1=1 nterms=i2-i1+1 c evaluate deviations delta 31 denom=x(i1+1)-x(i1) deltax=(xin-x(i1))/denom do 35 i=1,nterms ix=i1+i-1 35 delta(i)=(x(ix)-x(i1))/denom c accumulate coefficients a a(1)=y(i1) 41 do 50 k=2,nterms prod=1. sum=0. imax=k-1 ixmax=i1+imax do 49 i=1,imax j=k-i prod=prod*(delta(k)-delta(j)) 49 sum=sum-a(j)/prod 50 a(k)=sum+y(ixmax)/prod c accumulate sum of expansion sum=a(1) do 57 j=2,nterms prod=1. imax=j-1 do 56 i=1,imax 56 prod=prod*(deltax-delta(i)) 57 sum=sum+a(j)*prod yout=sum 61 return c end subroutine interp end subroutine polyft(xfit1,xfit2,xdata,ydata,ndata,nterms,aout) c c get coefficients for polynomial fit : c ydata = aout(1) + aout(2)*xdata + aout(3) *xdata^2 + ... c the fit is done between xdata = [xfit1, xfit2] c c arguments: c xfit1 lower bound of fitting range (single precision) (in) c xfit2 upper bound of fitting range (single precision) (in) c xdata array of abscissa values for data (single precision) (in) c ydata array of ordinate values for data (single precision) (in) c ndata length of data arrays (in) c nterms number of terms in polynomial (in) c aout coefficients of fitted polynomial (single precision) (out) c c copyright 1992 university of washington : matt newville c c requires functions nofx and determ. c note that double and single precision are mixed here. c most internal, working arrays use dp (as does routine determ) c c see bevington pg 104 for details c parameter (max= 5, max2m1 = 2*max-1, zero = 0.) real xdata(ndata), ydata(ndata), aout(nterms) double precision sumx(max2m1), sumy(max) double precision array(max,max), ain(max), delta, determ external determ, nofx c c initialize internal arrays c nmax = 2 * nterms - 1 do 100 i=1, nmax sumx(i) = zero 100 continue do 120 i = 1, nterms ain(i) = zero sumy(i) = zero do 110 j = 1, nterms array(i,j) = zero 110 continue 120 continue c c find points closest to endpoints of fitting range c nfit1 = nofx(xfit1,xdata,ndata) nfit2 = nofx(xfit2,xdata,ndata) if (nfit1.gt.nfit2) then ntemp = nfit1 nfit1 = nfit2 nfit2 = ntemp end if if(nfit1.eq.nfit2) go to 300 c c collect sums of data, sum of squares of data, etc. c do 200 i = nfit1, nfit2 xi = xdata(i) yi = ydata(i) xterm = 1.0 do 180 n=1, nmax sumx(n) = sumx(n) + xterm xterm = xterm * xi 180 continue yterm = yi do 190 n=1,nterms sumy(n) = sumy(n) + yterm yterm = yterm * xi 190 continue 200 continue c c construct matrices and evaluate coefficients c do 220 j=1,nterms do 210 k=1,nterms array(j,k) = sumx(j + k - 1) 210 continue 220 continue c c take determinant, get coefficients c delta = determ(array,nterms,max) if (delta.ne.zero) then do 260 l=1,nterms do 250 j=1,nterms do 240 k=1,nterms array(j,k) = sumx(j+k-1) 240 continue array(j,l) = sumy(j) 250 continue ain(l) = determ(array,nterms,max)/delta 260 continue end if c c convert coefficients to single precision, leave c 300 continue do 400 i = 1, nterms aout(i) = sngl(ain(i)) 400 continue return c end subroutine polyft end function nofx(x,array,npts) c c function nofx c c notes: x and array are real c c purpose c given a value x and an array of values, find the index c corresponding to the array element closest to x c usage c n = nofx(x,array,npts) c c parameters c x - a given value c array - array of values, assumed to be stored in order of c increasing value c npts - number of elements in array c c subroutines and function subprograms required c none c c written 8/11/81 by j.m. tranquada c real array(npts), x imin = 1 imax = npts inc = ( imax - imin ) / 2 10 continue it = imin + inc xit = array(it) if ( x .lt. xit ) then imax = it else if ( x .gt. xit ) then imin = it else nofx = it return endif inc = ( imax - imin ) / 2 if ( inc .gt. 0 ) go to 10 xave = ( array(imin) + array(imin+1) ) / 2. if ( x .lt. xave ) then nofx = imin else nofx = imin + 1 endif return c end function nofx end double precision function determ(array,nord,nrows) c c calculate determinate of a square matrix c (from bevington "data reduction and error analysis c for the physical sciences" pg 294) c array: matrix to be analyzed c nterms: order of matrix c nrows: first dimension of matrix in calling routine c double precision array(nrows,nrows) determ = 1. do 150 k=1,nord c c if (array(k,k).ne.0) go to 130 do 100 j=k,nord if (array(k,j).ne.0) go to 110 100 continue determ = 0. go to 160 c 110 do 120 i=k,nord save = array(i,j) array(i,j) = array(i,k) 120 array(i,k) = save determ = -determ c c 130 determ = determ*array(k,k) if (k.ge.nord) go to 150 k1 = k+1 do 140 i=k1,nord do 140 j=k1,nord 140 array(i,j) = array(i,j)-array(i,k)*array(k,j)/array(k,k) 150 continue 160 return c end double precision function determ end subroutine messag(messg) c c write message to standard ouput with (1x,a) format c character*(*) messg write(*,10) messg 10 format(1x,a) return c end subroutine messag end