c                             PTrap17j    G.Audi  WJ.Huang     v 27 jan 2017
c
c  Conversion of Frequency Ratios to Linear Equations
c  including electron and molecular binding energies
c
      real*8 xzero,mel,mref,smref,mrefk,rap,srap,coef
      real*8 prov,membre,sigmem,m118,sm118,meb,eref,e118,y4
      integer*4 q118,qref,znum,val,nbion,i,qel
      character txref*4,tx118*4,rev*2,ael*4
      character txiref*4,txi118*4
      character*30 filea,fileb,filec
      dimension znum(450),ael(450),qel(450),meb(450)
c
c        mel : electron mass in micro-u
c        mref,smref,m118,sm118 : Masses and uncert. for reference (ref) and mesured (118)
c
      filea='ptkl.equat                      '                 ! output file
      fileb='ptkl.freq                       '                 ! input file
      filec='ionization.data                 '                 ! electron+mol. file
      open(unit=1,file=filea,form='formatted',status='new')
      open(unit=3,file=fileb,form='formatted',status='old')
      open(unit=10,file=filec,form='formatted',status='old')
      mel   = 548.57990907d0                                   ! mass of electron in micro-u
      xzero = 9.314940038d-1                                   ! conversion factor micro-u to keV

      do i=1,450                                               ! read electron+mol. file
        read(10,*,iostat=val) znum(i),ael(i),qel(i),meb(i)
        meb(i) = (meb(i)*1.d-3)/xzero                          ! in micro-u
        if(val.lt.0) then
          nbion = i - 1                                        ! nbion: number of lines in elec+mol file
          exit
        end if
      end do
      close(10)

   12 read(3,1001,err=99) iaref,txref,qref,mref,rev,smref      ! read ref.name, mass(micro-u) and charge
 1001 format(i4,a4,i4,f19.8,a2,f13.8)
      mrefk = mref  * xzero                                    ! ref. mass in keV
c
   15 read(3,1001,end=90,err=99) ia118,tx118,q118,rap,rev,srap ! read frequency ratio
      if(tx118.eq.'NEW ') go to 12                             ! reset reference
      if(rev.eq.'  ') then                                     ! if reversed freq. ratio: rev=-1
        rap  = rap  / 1.d+6
        srap = srap / 1.d+6
      else
        rap  = 1.d+6 / rap
        srap = rap*rap * srap/1.d+6
      endif
      coef = 1000.*ia118/iaref
      coef = anint(coef) / 1000.                               ! calculate 3-digit coefficient
c
      prov = (ia118*1.d+0)/iaref - rap*q118/qref               ! start calculating the equation value
      membre = mref*(rap*q118/qref-coef) + mel*q118*(1-rap)
     *         - iaref*1.d+6*prov                              ! value (in micro-u) for the equation
      sigmem = srap * iaref * 1.d+6 * q118/qref                ! its uncertainty
c
c       Correct for electron and molecular binding energies
c
      eref = 0.
      e118 = 0.
      do i=1,nbion                                             ! loop over the nbion lines in elec+mol file
         txiref=txref
         if(txiref(3:3).eq.'m' .or. txiref(3:3).eq.'n' .or. txiref(3:3).eq.'x') then
           txiref(3:4)=txiref(4:4)//' '                        ! isomers m n or x : treated as gs
           if(txiref(2:2).eq.'x') txiref(2:3)=txiref(3:3)//' '
         endif
         txi118=tx118
         if(txi118(3:3).eq.'m' .or. txi118(3:3).eq.'n' .or. txi118(3:3).eq.'x') then
           txi118(3:4)=txi118(4:4)//' '
           if(txi118(2:2).eq.'x') txi118(2:3)=txi118(3:3)//' '
         endif
         if(trim(txiref)==trim(ael(i)).and.qref==qel(i)) eref = meb(i)
         if(trim(txi118)==trim(ael(i)).and.q118==qel(i)) e118 = meb(i)
      end do

      if(eref.ne.0. .and. e118.ne.0.) then
        y4 = eref*rap*q118/qref - e118                         ! calculate electronic and molecular correction
      else
        y4 = 0.
        write (1,*) "WARNING : el.+mol. binding en. not found"
      end if
      membre = membre + y4
c
      write (1,1020) ia118,tx118,iaref,txref,coef,membre,sigmem
 1020 format(5x,i6,a4,'-',i4,a4,'*',f6.3,' =', f15.5,' (',f9.5,')')
      m118  = membre + coef*mref
      sm118 = sqrt(sigmem**2 + (coef*smref)**2)
      write (1,1030) ia118,tx118,m118,sm118
 1030 format(13x,i4,a4,' =',f15.6,' +/-',f11.6,' micro-u')
      m118  = m118  * xzero
      sm118 = sm118 * xzero
      write (1,1032) m118,sm118
 1032 format(13x,8x,' =',f15.6,' +/-',f11.6,' keV',/)
c
      go to 15
c
   90 write (1,1990)
 1990 format(1H0,'Normal End of Freq.Ratios to Equations Conversion')
      stop
   99 write(1,1999)
 1999 format(1H0,'Error in File Reading')
      stop
      end