       subroutine binemis(lun11,lpri,
     $       t,vturbi,epi,ncn2,dpthc,
     $       idat1,rdat1,kdat1,nptrs,
     $       npar,npnxt,npfi,
     $       nplin,nlsvn,
     $       eliml,elimh,elum,zrems,zremsz,ilsv,
     $       zrtmp,ewsv,elsv,nlsv)
c
C
      implicit none
c
      include './PARAM'
      integer nbtpp
      parameter (nbtpp=10*ncn)
c
c
c     passed parameters
      integer lun11
c     master data 
      integer idat1(nidat1),nptrs(nptt,ndat2)
      real rdat1(nrdat1)
      character*1 kdat1(nkdat1)
c     pointers to master data 
      integer npar(ndat2),npnxt(ndat2)
      integer npfi(ntyp,nni)
c     pointers to line data 
      integer nplin(nnnl)
c     line luminosities
      real elum(3,nnnl)
c     energy bins
      real epi(ncn)
c     continuum lum
      real zrems(4,ncn),zremsz(ncn)
c     continuum optical depths
      real dpthc(2,ncn)
c     temporaries for dread
      real rdat(20000)
      integer idat(20000)
      real zrtmp(6,ncn)
      character*1 kdat(20000),kdat2(20000)
      integer kl, nilin, nkdt,nidt,lcon,lrtyp,ltyp,ml
      integer nlsvn, ln, ll, numcon, lnn, nbtmp, nrdt
      integer  verbose
      integer lup,ndtmp,mllz,nkdt2,iion,nitmp,
     $        ltyp2,lrtyp2,lcon2,nrdt2,nidt2,idat2(20000),iltmp,mlpar
      real eliml, elimh, elmmtpp, dele, etmp, elin, aasmall
      real dele2,egam,profile,deler,delea,rdat2(20000),vturbi,aatmp,bbb
      real optpp,e00,deleepi,etptst,tst,sume,zrsum1,zrsum2,deletpp,
     $        deleused,tmpe,zrtp2,zrtp1
      integer ml1,mlmin,mlmax,ij,mlm,ldir,mloff,ml1m,ml1mn,mlc,ncut,
     $        ml1mx,ml2
c     arrays containing printout info
      integer ilsv(nnnl)
      real ewsv(nnnl),elsv(nnnl)
      real etpp(nbtpp),zrtpp2(2,nbtpp),zrtmps(2,nbtpp)
      integer ldon(2)
c
      integer lpri,ncn2,nlsv,nelin
      real t,delet,deleturb,deleth,e0,vth,vturb
c
c     externally defined functions
      integer nbinc
      real voigte
c
      verbose=lpri
c
c     open and prepare the fits file for spectral data
      if(verbose.gt.0) write (lun11,*)'in binemis:'

c     build spectra data tables
      numcon=ncn2
      bbb=vturbi
      do ll=1,ncn2
        zrtmp(4,ll)=0.
        zrtmp(5,ll)=0.
        zrtmps(1,ll)=0.
        zrtmps(2,ll)=0.
        zrtpp2(1,ll)=0.
        zrtpp2(2,ll)=0.
        enddo
      do  lnn=1,nlsvn
        ln=lnn
        ml=nplin(ln)
        call dread(ltyp,lrtyp,lcon,
     $          nrdt,rdat,nidt,idat,nkdt,kdat,ml-1,
     $          idat1,rdat1,kdat1,nptrs,0,lun11)
        elin=abs(rdat(1))
        egam=rdat(3)
        lup=idat(2)
        nilin=npar(ml)
        call dread(ltyp,lrtyp,lcon,
     $               nrdt,rdat,nidt,idat,nkdt,kdat,nilin-1,
     $               idat1,rdat1,kdat1,nptrs,0,lun11)
        nelin=npar(nilin)
        nilin=idat(3)
c       get nuclear mass       
        ml=nelin
        call dread(ltyp,lrtyp,lcon,
     $          nrdt,rdat,nidt,idat,nkdt,kdat,ml-1,
     $          idat1,rdat1,kdat1,nptrs,0,lun11)
        aatmp=rdat(2)
        elmmtpp=(elum(2,ln)+elum(1,ln))/2.
        if (lpri.ne.0)
     $     write (lun11,*)ln,elin,elmmtpp,nilin,nelin,egam,
     $              lup,aatmp
        if (((ln.gt.0).and.(ln.le.nnnl)) 
     $    .and.((elin.gt.eliml).and.(elin.lt.elimh)) 
     $    .and.(elmmtpp.gt.1.e-24).and.(aatmp.gt.1.e-24)
     $    .and.((nilin.gt.0).and.(nilin.le.nni)))
     $    then

c         line parameters
          etmp=12398.54/elin
          nbtmp=nbinc(etmp,epi,ncn2)
c
c         find associated type 86 data
          iion=0
          idat(nidt)=0
          if (lpri.ne.0)
     $      write (lun11,*)'searching for ion'
          do while ((idat(nidt).ne.nilin).and.(iion.lt.nni))
            iion=iion+1
            nitmp=npfi(13,iion)
            call dread(ltyp,lrtyp,lcon,
     $               nrdt,rdat,nidt,idat,nkdt,kdat,nitmp-1,
     $               idat1,rdat1,kdat1,nptrs,0,lun11)
            if (lpri.ne.0)
     $        write (lun11,*)iion,idat(nidt),nilin,nitmp
            enddo
          ndtmp=npfi(41,iion)
          delea=0.
          if (ndtmp.gt.0) then
            if (lpri.ne.0)
     $        write (lun11,*)'  found ion',lup,ndtmp
            mllz=npar(ndtmp)
            call dread(ltyp2,lrtyp2,lcon2,
     $         nrdt2,rdat2,nidt2,idat2,nkdt2,kdat2,ndtmp-1,
     $         idat1,rdat1,kdat1,nptrs,0,lun11)
            iltmp=idat2(2)
            mlpar=mllz
            do while ((ndtmp.ne.0).and.(lup.ne.iltmp)
     $         .and.(mlpar.eq.mllz)) 
              call dread(ltyp2,lrtyp2,lcon2,
     $          nrdt2,rdat2,nidt2,idat2,nkdt2,kdat2,ndtmp-1,
     $          idat1,rdat1,kdat1,nptrs,0,lun11)
              iltmp=idat2(2)
              if (lpri.ne.0)
     $           write (lun11,*)'   ',nidt2,iltmp,ndtmp
              ndtmp=npnxt(ndtmp)     
              mlpar=0
              if (ndtmp.ne.0) mlpar=npar(ndtmp)
              enddo
            endif
          if (lup.eq.iltmp) then
            delea=rdat2(3)*(4.14e-15)
            egam=rdat2(4)
            endif
c
c         cheat for narrow line plot
c         delea=0.        
c
c         a list of all the deles
c           delea=auger natural width in eV
c           deleturb=turbulent width
c           deleth=thermal Doppler width
c           dele=thermal+turbulent width
c           deler=radiative natural width
c           dele2=total width, natural + Doppler
c           deletpp=goal of resolution of internal grid=dele/8
c           deleepi=xstar grid spacing
c           deleused=spacing of internal grid=deleepi/int(deleepi/depetpp)
c           delet=energy offset from line center in units of dele (local)
c
c         thermal width quantities
          vth=(1.2e+1)*sqrt(t/aatmp)
          vturb=max(bbb,vth)
          e0=(12398.54)/amax1(elin,1.e-24)
          deleturb=e0*(vturb/3.e+5)
          deleth=e0*(vth/3.e+5)
          dele=deleth+deleturb
          deler=egam*(4.14e-15)
          dele2=delea+deler+dele
          aasmall=(delea+deler)/(1.e-36+dele)/12.56


          ml1=nbtmp
          if (lpri.ne.0) write (lun11,*)'calculating profile',
     $                     ml1,epi(ml1),epi(ml1+1),optpp
c
c         calculate profile on temporary grid
          e00=epi(ml1)
          deleepi=epi(ml1+1)-epi(ml1)
          deletpp=dele/8.
          ncut=int(deleepi/deletpp)
          ncut=max(ncut,2)
          deleused=deleepi/float(ncut)
          if (lpri.ne.0) write (lun11,*)'ncut=',ncut,deleused,deletpp,
     $                                    deleepi
c         apply binned opacity
          delet=(e00-etmp)/dele
          if (aasmall.gt.1.e-6) then
              profile=voigte(abs(delet),aasmall)
            else
              profile=exp(-delet*delet)/1.772
            endif 
          ml2=nbtpp/2
          etpp(ml2)=e00
          zrtpp2(1,ml2)=elum(1,ln)*profile
          zrtpp2(2,ml2)=elum(2,ln)*profile
          mlc=0
          ldir=1
          ldon(1)=0
          ldon(2)=0
          mlmin=nbtpp
          mlmax=1
          do while ((ldon(1)*ldon(2).eq.0).and.(mlc.lt.nbtpp/2)) 
            mlc=mlc+1
            do ij=1,2
              ldir=-ldir
              if (ldon(ij).ne.1) then
                mlm=ml2+ldir*mlc
                etptst=e00+float(ldir*mlc)*deleused
                if ((mlm.lt.nbtpp).and.(mlm.gt.1)
     $            .and.(etptst.gt.0.).and.(etptst.lt.epi(ncn2))) then
                  mlmin=min(mlm,mlmin)
                  mlmax=max(mlm,mlmax)
                  etpp(mlm)=e00+float(ldir*mlc)*deleused
                  delet=(etpp(mlm)-etmp)/dele
                  if (aasmall.gt.1.e-9) then
                      profile=voigte(abs(delet),aasmall)
                    else
                      profile=exp(-delet*delet)/1.772
                    endif 
                  profile=profile/dele/(1.602197e-12)
                  zrtpp2(1,mlm)=elum(1,ln)*profile
                  zrtpp2(2,mlm)=elum(2,ln)*profile
                  if (lpri.ne.0) write (lun11,*) mlm,etpp(mlm),ij,
     $               deleused,delet,mlmin,mlmax,
     $               mlc,mloff,mod(mloff,ncut),profile,zrtpp2(1,mlm),tst
                  endif
                if ((mlm.le.1).or.(mlm.ge.nbtpp)
     $                 .or.(etptst.le.0.).or.(etptst.ge.epi(ncn2))
     $                 .or.(mlc.gt.nbtpp).or.(abs(delet).gt.50.)) 
     $                  ldon(ij)=1
                endif
              enddo
            enddo
          sume=0.
          zrsum1=0.
          zrsum2=0.
          if (lpri.ne.0) write (lun11,*)'renormalizing profile'
          mlmin=max(mlmin,2)        
          mlmax=min(mlmax,nbtpp)        
          ml1m=nbinc(etpp(mlmin),epi,ncn2)
          ml1mn=ml1m
          do mlm=mlmin,mlmax
            mloff=mlm-ml2
            if (((etpp(mlm).gt.epi(ml1m)).or.(mlm.eq.mlmax))
     $          .and.(mlm.ne.mlmin)) then 
                if (mlm.eq.mlmax) ml1m=max(1,ml1m-1)
                tmpe=abs(etpp(mlm)-etpp(mlm-1))
                if (mlm.gt.mlmin+1) then
                  sume=sume+tmpe
                  zrsum1=zrsum1+(zrtpp2(1,mlm)+zrtpp2(1,mlm-1))*tmpe/2.
                  zrsum2=zrsum2+(zrtpp2(2,mlm)+zrtpp2(2,mlm-1))*tmpe/2.
                  endif
                if (sume.gt.1.d-4) then
                  zrtp2=zrsum2/sume
                  zrtp1=zrsum1/sume
                  endif
                do while ((etpp(mlm).gt.epi(ml1m)).and.(ml1m.lt.ncn2)) 
                  zrtmps(1,ml1m)=zrtp1
                  zrtmps(2,ml1m)=zrtp2
                  ml1m=ml1m+1
                  enddo
                if (lpri.ne.0) write (lun11,*)mlm,mloff,ml1m,
     $               epi(ml1m),epi(ml1m+1),etpp(mlm),zrtmps(2,ml1m),
     $               zrtp1,zrtp2
                zrsum2=0.
                zrsum1=0.
                sume=0.
              else
                tmpe=abs(etpp(mlm)-etpp(mlm-1))
                sume=sume+tmpe
                zrsum1=zrsum1+(zrtpp2(1,mlm)+zrtpp2(1,mlm-1))*tmpe/2.
                zrsum2=zrsum2+(zrtpp2(2,mlm)+zrtpp2(2,mlm-1))*tmpe/2.
                sume=sume+tmpe
                if (lpri.ne.0) write (lun11,*)mlm,mloff,ml1m,
     $                epi(ml1m),etpp(mlm),zrtpp2(2,mlm),zrsum2,sume
              endif
            enddo
          ml1mx=nbinc(etpp(mlmax),epi,ncn2)
          ml1mx=max(ml1mx,ml1mn+1)
          zrsum1=0.
          zrsum2=0.
          do ml1m=ml1mn+1,ml1mx
            zrsum2=zrsum2+(zrtmps(2,ml1m)+zrtmps(2,ml1m-1))
     $                   *(epi(ml1m)-epi(ml1m-1))/2.
            zrsum1=zrsum1+(zrtmps(1,ml1m)+zrtmps(1,ml1m-1))
     $                   *(epi(ml1m)-epi(ml1m-1))/2.
            enddo
          zrsum1=zrsum1*(1.602197e-12)/(1.e-34+elum(1,ln))
          zrsum2=zrsum2*(1.602197e-12)/(1.e-34+elum(2,ln))
          do ml1m=ml1mn,ml1mx
            zrtmps(2,ml1m)=zrtmps(2,ml1m)/(1.e-34+zrsum2)
            zrtmp(5,ml1m)=zrtmp(5,ml1m)+zrtmps(2,ml1m)
            zrtmps(1,ml1m)=zrtmps(1,ml1m)/(1.e-34+zrsum1)
            zrtmp(4,ml1m)=zrtmp(4,ml1m)+zrtmps(1,ml1m)
            enddo
          endif
        enddo


      if (lpri.ne.0) write (lun11,*)'after first binemis loop'
      do kl=1,numcon
         if (lpri.ne.0) write (lun11,*)kl,epi(kl),zrems(2,kl),
     $          zrems(3,kl),zrtmp(4,kl),zrtmp(5,kl),dpthc(1,kl),
     $          zremsz(kl)
         zrtmp(4,kl)=zrtmp(4,kl)+zrems(2,kl)
         zrtmp(5,kl)=zrtmp(5,kl)+zrems(3,kl)
         zrems(2,kl)=zrtmp(4,kl)
         zrems(3,kl)=zrtmp(5,kl)
         zrtmp(3,kl)=zremsz(kl)*exp(-dpthc(1,kl))
         zrtmp(2,kl)=zremsz(kl)
         zrtmp(1,kl)=epi(kl)
         zrtmp(6,kl)=zrems(4,kl)
         enddo
c

      return
      end
