program Gausfit
c
c  Gaussian deblending
c
c  Input files:
c data read from fort.4, line fit parameters from fort.3
c
c  Output files:
c fort.11- input data- medianized, gaussian smoothing, SNR, standard deviation 
c          from median
c fort.12- input data, initial gaussians guess, residual of guess subtractions
c fort.13- final gaussian components
c fort.14- final fit, sum of gaussians in fort.13
c fort.15- standard deviations of final fit from the data
c fort.16- best fit reflected wrt point of maximum flux, reflected-observed
c fort.17- data reflected as in fort.16
c fort.21- points deleted from the fit due to large deviations
c fort.22- points used for the fit
c fort.23- velocity profile: velocity wrt rest, wrt line peak, data, fit, 
c          reflection ratio, kinematic ratio, observed ratio
c fort.26- template fitting: velocity wrt rest, data, error
c
c
      parameter (ne=6000,me=30)
      integer lista(me)
      real data(ne),sig(ne),x(ne),xmd(ne),xg(ne,2)
     &  ,wav(ne),fc(ne),fl(ne),a(me),covar(me,me),ygb(ne)
     & ,alpha(me,me),fgaus(ne,10),flz(ne),yga(ne)
     & ,flz1(ne),wav1(ne),sig1(ne),fc1(ne)
     & ,wl(5),dw(ne),wavz(ne),fcz(ne),wrk(ne,3),yga1(ne) 
     & ,wavz1(ne),wrat(ne),xsig(ne), xag(ne,2),fcl(me)
      external fgauss 
      common /log/ wav,flz,xg,a,sig
      common /err/ wav1,flz1,sig1
      data data,xg/ne*0.,ne*0.,ne*0./
      data a/me*0./
c
c   xs1 - standard deviation limit used for median filtering
c   xs2 - as above using the fitted model
c
      xs1=2.5
cc      xs2=3.
      xs2=100.
      ismooth=1
      read(3,*) wline, ma,wcon1,wcon2,iblue,nfix
      mft=ma-nfix
      read(3,*)(lista(i),i=1,ma)
      read(3,*)(a(lista(i)),i=ma,ma-nfix+1,-1)
      print*, wline,ma,wcon1,wcon2,iblue, nfix
c
c  line width, i.e. FWHM/2.3548, should be put in units of A
c  preset values should be different from 0.
c
      write(*,*)(a(lista(i)),i=ma,ma-nfix+1,-1)
      if(ma.gt.9) read(3,*) (wl(i),i=1,ma/3-3)
       write(*,*) (wl(i),i=1,ma/3-3)
c
c  ismooth=1  use data after gaussian filtering, 0- unsmoothed data
c   iblue=1-use blue wing only, 0-use red wing only, 2-use both wings  
c
cc      iblue=1 
      read(4,*) nd,z,ext
      print*,nd,' data points  z =',z,'  E(B-V) =',ext
c
c  in optical data do not read errors, assume S/N=50 as typical in UV after
c  smoothing on resolution element size (1 in optical)
c  all optical data files are shorter than 3000 lines
c
cc      if(nd.gt.3000) then
      if(z.gt..01) then
       read(4,*) (wrk(j,1),wrk(j,2),wrk(j,3),j=1,nd)
      else
c
c  optical data files are in rest frame, and do not include flux error
c  Assume an arbitrary S/N=50
c
       read(4,*) (wrk(j,1),wrk(j,2),j=1,nd)
       do 42 j=1,nd
 42     wrk(j,3)=wrk(j,2)/50.
c
c   subtract fe2 from optical spectrum, put in xx.opt1 files
c
cc       call fe2(wrk,nd,z)
cc      stop 'done subtracting Fe II'
      endif
c
c  cut out only the interesting region
c
      nc=0
      wc1=.99*wcon1*(1.+z)
      wc2=1.01*wcon2*(1.+z)
      do 32 i=1,nd
       if(wrk(i,1).lt.wc1.or.wrk(i,1).gt.wc2) goto 32 
       nc=nc+1
       wav(nc)=wrk(i,1)/(1.+z)
       wavz(nc)=wav(nc)
         data(nc)=wrk(i,2)
         sig(nc)=wrk(i,3)
 32   continue 
      nd=nc
      if(nc.eq.0) stop 'no data points in requested region' 
c
c  dereddning the data using E(B-V)=ext, note that error vertor sig() is corrected
c
      if(ext.ne.0.) call dred(z,wav,data,sig,nd,ext)
      ndz=nd
      if(wav(1)*(1.+z).lt.3270.) ns=10
      if(wav(1)*(1.+z).ge.3270.) ns=1
c-------------- note that median box size is 2*ns+1
      xav=0.
      do 2 i=1,nd
       np=0
       do 3 j=-ns,ns
        if(i+j.lt.1.or.i+j.gt.nd) goto 3
         np=np+1
         x(np)=data(i+j)
         xav=xav+x(np)
         xsig(np)=sig(i+j)
 3     continue
       xav=xav/np
       call median(x,np,xmed)
       xmd(i)=xmed
 2    continue
      ngaus=4
      sigt=0.
      do 4 i=1,nd
       xg(i,1)=0.
       xg(i,2)=0.
       wt=0.
c
c  do a Gaussian smoothing (for UV data only)
c
      if(wav(i)*(1.+z).lt.3300.) then
       if(i.eq.1)print*,'data smoothed with sigma=1.75pix Gaussian'
       do 5 j=-ngaus,ngaus
        if(i+j.lt.1.or.i+j.gt.nd) goto 5
        si=j/1.75
        w=.39894*exp(-si**2/2.) 
        wt=wt+w 
        xg(i,1)=xg(i,1)+data(i+j)*w
        xg(i,2)=xg(i,2)+(sig(i+j)*w)**2
 5    continue
       xg(i,1)=xg(i,1)/wt
       xg(i,2)=sqrt(xg(i,2))/wt
      else
       if(i.eq.1) print*,'data is not smoothed'
       xg(i,1)=data(i)
       xg(i,2)=sig(i)
      endif
       xag(i,1)=data(i)
       xag(i,2)=sig(i)
       if(xmd(i).eq.0) goto 4
       diff=(xg(i,1)-xmd(i))/xg(i,2)
       snr=xg(i,1)/xg(i,2)
c--- comment following line for fig.1 and 2
       write(11,*) wav(i),xmd(i),xg(i,1),snr,diff
       d1=abs(wav(i)-1216.)
       d2=abs(wav(i)-1549.)
       d3=abs(wav(i)-1909.)
       d4=abs(wav(i)-2799.)
       if(diff.lt.-3.) goto 4 
c--- following lines used to produce fig.1 and 2 
       if(d1.lt.5..or.d2.lt.5..or.d3.lt.5..or.d4.lt.5.)then
cc         write(11,*) wav(i),xg(i,1),xg(i,1)
cc         write(11,'(f7.2,1p,e10.2)') wav(i),xg(i,1)
        else
cc         write(11,*) wav(i),xmd(i),xg(i,1)
c         write(11,'(f7.2,1p,e11.3)') wav(i),xmd(i)
cc         write(11,'(f7.2,1p,e11.3)') wav(i),xg(i,1)
       endif
 4    continue 
cc       stop
c
c  set continuum level fc
c
      do 6 i=1,nd-1
       d1=(wcon1-wav(i))/(wav(i+1)-wav(i)) 
       d2=(wcon2-wav(i))/(wav(i+1)-wav(i))
       if(d1.ge.0.and.d1.lt.1.)then
         fcon1=(1.-d1)*xmd(i)+d1*xmd(i+1)
c
c err1,err2= the average 1 sigma error of the median continuum flux 
c
         err1=0.
         do 61 j=-ns,ns
 61       err1=err1+xg(i+j,2)
         err1=err1/(1.+2.*ns)
       endif  
       if(d2.ge.0.and.d2.lt.1.)then
         fcon2=(1.-d2)*xmd(i)+d2*xmd(i+1)
         err2=0.
         do 62 j=-ns,ns
 62       err2=err2+xg(i+j,2)
         err2=err2/(1.+2.*ns)
       endif 
 6     continue
       print*,'continuum windows:'
       print*,wcon1,fcon1,err1
       print*,wcon2,fcon2,err2
       wcc1=wcon1
       wcc2=wcon2
c
c--  o6int=1. extrapolate continuum instead of interpolate (for O VI in some ob)
c
       o6int=0.
       if(z.eq.0.23405.or.z.eq.0.15754) then
        if(wline.eq.1033.816) o6int=1.
       endif 
       if(o6int.eq.1.) then
         print*,'type in wc1 fc1 err1 wc2 fc2 err2'
         read*,wcc1,fcon1,err1,wcc2,fcon2,err2
       endif 
c
c  determine whether to take observed continuum or observed + sigma, observed -sigma
c
       print*,'type in iadd  -1 0 1'
       read*,iadd 
ccc       iadd=0
       fcon1=fcon1+iadd*err1
       fcon2=fcon2+iadd*err2
      ac=log(fcon2/fcon1)/log(wcc2/wcc1) 
      print*,'alpha_nu = ',-ac-2.
cc      stop
      flmx=0.
      imax=0
c
c  set continuum level at line centers
c
       fcl(1)=fcon1*exp( ac*log(wline/wcc1) )
       if(ma.ge.12)fcl(2)=fcon1*exp( ac*log(wl(1)/wcc1) )  
       if(ma.ge.15)fcl(3)=fcon1*exp( ac*log(wl(2)/wcc1) )  
       if(ma.ge.18)fcl(4)=fcon1*exp( ac*log(wl(3)/wcc1) )  
c
c  subtract continuum  fc  to get fl, set sig1, 
c  delete points which deviate by more than -xs1 sigma from the median
c  also delete points outside wcon1-wcon2
c
       itc=0
       do 7 i=1,nd
        if(o6int.eq.1.) then
         xfc=log(fcon1)+ac*log(wav(i)/wcc1)
        else
         xfc=log(fcon1)+ac*log(wav(i)/wcon1)
        endif
        fcz(i)=exp(xfc)
        if(wav(i).lt.wcon1.or.wav(i).gt.wcon2) goto 7
c
c  looking for maximum deviation within +-3 points (region of point rejection)
c
        devmax=10000.
        do 41 j=-3,3
         if(i+j.lt.1.or.i+j.gt.nd) goto 41
         if(ismooth.eq.0) then
          devi=(data(i+j)-xmd(i+j))/sig(i+j)
         else
          devi=(xg(i+j,1)-xmd(i+j))/xg(i+j,2)
         endif
         if(devi.lt.devmax) devmax=devi
 41     continue
c-- following used in Ly alpha fit for 0405, to prevent spurious absorption 
        if(abs(wav(i)-1215.).lt.5.and.z.eq..5728) goto 1215
c----  do not reject points if they are in the vally inside [O III]     
        if(abs(wav(i)-4977.).gt.5.) then
         if(devmax.lt.-xs1) then
            write(21,*) wav(i),xg(i,1),devmax  
            goto 7
         endif
        endif
 1215  continue
c--  eliminating blue noisy end part of the spectrum (used for PG 0953 +414)
        if(wav(i).lt.980..and.z.eq.0.23405) then
          goto 7
        endif
c
c  delete all points close to geocoronal Ly alpha
c   
        if(abs(wav(i)*(1.+z)-1216.).lt.9.*(1.+z)) then
           write(21,*) wav(i),xg(i,1),devmax   
           goto 7
        endif 
c-- in H beta fitting ignores region near [O III]
c-- but do not reject points in vally inside [O III]
          dis1=abs(wav(i)-4959.) 
          dis2=abs(wav(i)-5007.) 
          if(dis1.lt.20..or.dis2.lt.20..and.wline.eq.4861.32) then
           if(abs(wav(i)-4977.).gt.5.) then
             write(21,*) wav(i),xg(i,1)
             goto 7
           endif
          endif
        itc=itc+1 
        fc(itc)=fcz(i)
        if(ismooth.eq.0) then
          fl(itc)=data(i)-fc(itc) 
          sig1(itc)=sig(i)
         else
          fl(itc)=xg(i,1)-fc(itc)
          sig1(itc)=xg(i,2)
        endif 
        if(xg(i,1).eq.0) fl(itc)=data(i)-fc(itc)
        flz(itc)=fl(itc)
        if(fl(itc).gt.flmx) then
          flmx=fl(itc)
          imax=itc
        endif
        if(sig1(itc).eq.0) sig1(itc)=sig(i) 
c
c wav1, flz1, sig1, fc1 are the arrays cleaned of possible absorption lines
c
        flz1(itc)=flz(itc)
        wav1(itc)=wav(i)
        fc1(itc)=fc(itc)
c--  prints data points remaining after median rejection
cc       write(20,*) wav1(itc),fc(itc),fc(itc)+fl(itc)
 7     continue
        print*,itc, ' points left after median rejection'
        nd=itc
        flmxz=flmx
       dw(1)=wav1(2)-wav1(1)
       dw(nd)=wav1(nd)-wav1(nd-1)
       do 17 i=1,nd
        wav(i)=wav1(i)
        sig(i)=sig1(i) 
        if(i.gt.1.and.i.lt.nd) dw(i)=(wav1(i+1)-wav1(i-1))/2.
 17     continue
c
c  calculate first three moments x1,x2,x3 for up to 3 components.  
c
       itr=0
cc        print*,'Initial guess for gaussian components'
 20    itr=itr+1
        x1=0.
        x2=0.
        x3=0.
        idis1=imax-1
        idis2=nd-imax
        idmin=idis1
        if(idis2.lt.idis1) idmin=idis2 
c
c  force 3rd component to be around the peak
c
        if(itr.eq.3.and.abs(wline-wmx).gt.10.)wmx=wline
c-- in H beta fitting ignore the higher O III peak    
       if(itr.eq.2.and.wav(i)*(1.+z).gt.3300.)wmx=wline
       do 8 i=1,nd
        if(itr.le.3) then
c
c  use a symmetric range within idmin of imax (to have the gaussian peak at imax)
c
          if(abs(i-imax).gt.idmin) goto 8
          iref=2*imax-i
          if(iref.lt.1) iref=1
          if(iref.gt.nd) iref=nd
          flx=fl(i)  
          if(i.gt.imax.and.iblue.eq.1)flx=fl(iref)
          if(i.lt.imax.and.iblue.eq.0)flx=fl(iref)
cc          if(flx.gt..2*flmx.and.itr.eq.1) flx=0.2*flmx
          if(flx.gt..1*flmx.and.itr.eq.1) flx=0.1*flmx
          if(flx.gt..7*flmx.and.itr.eq.2) flx=0.7*flmx
         else
          flx=flz(i) 
        endif 
c
c        third (narrow) component is fit only with data near center.
c         range of fit for UV data is smaller than in optical
c
        if(flx.lt.0.) flx=0. 
        if(wav(i)*(1.+z).lt.3300.) then
         if(abs(wav(i)-wmx).gt.20.and.itr.eq.2) goto 8 
         if(abs(wav(i)-wmx).gt.5.and.itr.eq.3) goto 8 
         if(abs(wav(i)-wl(1)).gt.10.and.itr.eq.4) goto 8 
         if(abs(wav(i)-wl(2)).gt.10.and.itr.eq.5) goto 8 
         if(abs(wav(i)-wl(3)).gt.10.and.itr.eq.6) goto 8 
ccc         if(abs(wav(i)-wl(4)).gt.4.and.itr.eq.7) goto 8 
        else 
         if(abs(wav(i)-wmx).gt.60.and.itr.eq.2) goto 8 
         if(abs(wav(i)-wmx).gt.5.and.itr.eq.3) goto 8 
         if(abs(wav(i)-wl(1)).gt.10.and.itr.eq.4) goto 8 
         if(abs(wav(i)-wl(2)).gt.20.and.itr.eq.5) goto 8 
         if(abs(wav(i)-wl(3)).gt.20.and.itr.eq.6) goto 8 
        endif
cc        write(17,*) wav(i),flx
        x1=x1+flx*dw(i)
        x2=x2+wav(i)*flx*dw(i)
        x3=x3+(wav(i)**2)*flx*dw(i)
 8     continue
        wcen=x2/x1
        delta=sqrt(abs(x3/x1-wcen**2))
cc        print*,x1,wcen,delta
        flmx=0.
        wmx=0. 
        do 9 i=1,nd
         dev=((wav(i)-wcen)/delta)**2/2.
         fgaus(i,itr)=x1*.3989/delta*exp(-dev)
c--- fl() = the residual after gaussians subtractions 
         fl(i)=fl(i)-fgaus(i,itr)
         if(fl(i).gt.flmx) then
           flmx=fl(i)
           wmx=wav(i)
         endif
         write(12,*) wav(i),flz(i)+fc(i),fc(i)+fgaus(i,itr)
     &  ,fc(i)+fl(i)
 9      continue
        nit=3*(itr-1)
c
c  change a(i) which are read from fort.3 (i.e. be fixed) to proper units
c
        if(a(nit+3).ne.0.) a(nit+3)=a(nit+3)*sqrt(2.)
        if(a(nit+1).ne.0.) a(nit+1)=a(nit+1)*.3989*sqrt(2.)/a(nit+3)
c
c  calculating a(i) which were not read (these were set to 0 in the data state') 
c
        if(a(nit+1).eq.0.) a(nit+1)=x1*.3989/delta
        if(a(nit+2).eq.0.) a(nit+2)=wcen
        if(a(nit+3).eq.0.) a(nit+3)=delta*sqrt(2.)
        if(itr.lt.ma/3) goto 20
        do 15 i=1,nd
         ygb(i)=0.
         do 16 j=1,ma/3
 16       ygb(i)=ygb(i)+fgaus(i,j)  
c--   prints initial guess for the shape
ccc         write(13,*) wav(i),fc(i)+ygb(i)
 15     continue   
c
c  Initial decomposition completed, now iterate for best fit
c  Set order of fitting parameters
c
       ndd=nd 
       nm=0
      nk=20
      dv=(vmax-vmin)/float(nk-1.)
      do 10 k=1,nk
       niter1=0
       vct=vmax-dv*(k-1.)
       zshift=vct/3.e5
 21     alamda=-5.
        chisq=0.
        niter1=niter1+1
        chibef=0.
       do 11 i=1,30
        i3=i/3
        chibef=chisq 
        call mrqmin(wav1,flz1,sig1,ndd,a,ma,lista,mft,covar,alpha,ma,
     &   chisq,fgauss,alamda)
c--  stop iterating if chisq didn't change for i>15
        if(chisq.eq.chibef.and.i.gt.15) goto 111
 11     continue
c
c  adding all components to get the fitted model yga()
c 
 111   print*,'ch^2=',chisq,' after',i,' iterations'
       print*,'lambda_line  f_cont    f_line    EW   v shift   FWHM'
       flinetot=0.
       do 12 j=1,ma-2,3
        delta=a(j+2)/sqrt(2.)
        xflx=a(j)*abs(delta)/.3989
cc        print*,xflx,a(j+1),delta
        if(j.lt.10) then
         wll=wline
         fcc=fcl(1)
        else
         wll=wl((j-7)/3)
         fcc=fcl((j-4)/3)
        endif
        fline=xflx*(1.+z)
        ew=xflx/fcc
        iveldev=3.e5*((a(j+1)-wll)/wll)
        iveldis=delta*2.3548*3.e5/a(j+1) 
        flinetot=flinetot+fline
        write(6,127) wll,fcc,fline,ew,iveldev,iveldis
 127    format(f10.2,1p,2e10.2,0p,f6.1,2i8)
        do 12 i=1,ndd
         if(j.eq.1)yga(i)=0.
          arg=(wav1(i)-a(j+1))/a(j+2)
          yg=a(j)*exp(-arg**2)
          yga(i)=yga(i)+yg
 12     continue
c
c  store  points used in this iteration
c
       do 26 i=1,ndd
        flz(i)=flz1(i)
        wav(i)=wav1(i)
        sig(i)=sig1(i)
        fc(i)=fc1(i)
 26    continue
c
c  deleting data points which deviate by more than xs2 sigma (below) the
c  the model fit
c  calculate chi**2 before and after mrqmin, ignoring points which deviate
c  by more than xs2 sigma
c
        chis1=0.
        phot=0. 
       if(niter1.eq.1) nt1=0
        do 14 i=1,ndd  
         dch1a=(flz(i)-yga(i))/sig(i)
         dch1=dch1a**2
         chis1=chis1+dch1 
         write(15,*) wav(i),dch1a
         velo=abs( (wav(i)-wline)/wline*3.e5 )
         if(velo.lt.6000.) phot=phot+(flz(i)/sig(i))**2.
         if(niter1.eq.2) then
           goto 14
         endif
c
c  looking for maximum deviation within +-3 points (region of point rejection)
c
        chmax=0.
        do 24 j=-3,3
         if(i+j.lt.1.or.i+j.gt.ndd) goto 24
         dchi=(flz(i+j)-yga(i+j))/sig(i+j)
cc         if(dchi**2.gt.chmax**2) chmax=dchi
         if(dchi.lt.chmax) chmax=dchi
 24     continue
         dchi=(flz(i)-yga(i))/sig(i)
cc         if(chmax**2.gt.xs2**2) then
         if(chmax.lt.-xs2) then
           write(21,*) wav(i),fc(i)+flz(i),dchi   
           goto 14
         endif 
         nt1=nt1+1
c
c          store new values
c
         flz1(nt1)=flz(i)
         wav1(nt1)=wav(i)
         sig1(nt1)=sig(i)
         fc1(nt1)=fc(i) 
cc---  change to .lt.2 if iterating again after first fit (if points are rejected)
        if(niter1.lt.1) goto 14
         write(22,*) wav1(nt1),fc1(i)+flz1(nt1),dchi 
 14     continue   
cc        print*, ' iterated fit=',chis1/(nt1-ma+nfix),nt1-ma+nfix
        ndd=nt1
cc---  change as above.
       if(niter1.lt.1) goto 21
c
c  calculating model using all original frequencies
c
       fmmx=0.
       immx=0
c
c   to eliminate 1175A feature from Lya red wing template (PG1241+176, 
c   PG1634+706) reset value of ma to ma-3, so that the template is without
c   the 1175A component.
c
cc      ma=12
      do 27 j=1,ma-2,3
        delta=a(j+2)/sqrt(2.)
        xflx=a(j)*abs(delta)/.3989
       do 27 i=1,ndz
         if(j.eq.1)yga(i)=0.
         arg=(wavz(i)-a(j+1))/a(j+2)
         yg=a(j)*exp(-arg**2)
         yga(i)=yga(i)+yg
         write(13,*) wavz(i),fcz(i)+yg
c--  print first 3 components, i.e. principal line only
cc         if(j.eq.7) write(14,*) wavz(i),fcz(i)+yga(i)
         if(j.eq.ma-2) write(14,*) wavz(i),fcz(i)+yga(i)
         if(yga(i).gt.fmmx.and.j.eq.ma-2)then
          fmmx=yga(i)
          immx=i
         endif
 27     continue
        goto 40
 10    continue
c
c  approximate the accurate position of the peak by passing a parabola around fmmx
c  then calculate reflected spectrum around fmmx
c
 40    ap1=(yga(immx)-yga(immx-1))/(wavz(immx)-wavz(immx-1))
       ap2=(yga(immx+1)-yga(immx-1))/(wavz(immx+1)-wavz(immx-1))
       ap=(ap1-ap2)/(wavz(immx)-wavz(immx+1))
       bp=ap1-ap*(wavz(immx)+wavz(immx-1)) 
       wavmx=-bp/(2.*ap)
       do 29 j=1,ma-2,3
        delta=a(j+2)/sqrt(2.)
        xflx=a(j)*abs(delta)/.3989
       do 29 i=1,ndz
        wavsym=2.*wavmx-wavz(i)
        wavz1(i)=wavsym
        arg1=(wavsym-a(j+1))/a(j+2)
        yg1=a(j)*exp(-arg1**2)
        yga1(i)=yga1(i)+yg1
        if(j.eq.ma-2) write(16,*) wavz(i),fcz(i)+yga1(i)
     &   ,fcz(i)+yga(i)-yga1(i)
cc        if(j.eq.ma-2) write(17,*) wavsym,xg(i,1)
 29    continue
       do 36 i=1,ndz
        do 36 j=ndz,2,-1
          d=(wavz(i)-wavz1(j-1))/(wavz1(j)-wavz1(j-1))
          if(d.le.0.or.d.gt.1) goto 36
          fint=(1.-d)*(xg(j-1,1)-fcz(j-1))+d*(xg(j,1)-fcz(j))
          fcin=(1.-d)*fcz(j-1)+d*fcz(j)
          vel=(wavz(i)-wline)/wline*3.e5
c--  wrat()  the wings ratio reflected through wavsym
          wrat(i)=(xg(i,1)-fcz(i))/fint
          write(17,*) wavz(i),fint+fcz(i), xg(i,1)-fint
 36     continue          
c
c  print line profile in velocity units normalized to the continuum flux at the
c  midpoint, use all observed points (including rejected ones)
c
        fnor2=0.
        v1=0.
        do 28 i=1,ndz
c-- vel=velocity wrt expected rest wavelength, vel1=wrt observed peak
         vel=(wavz(i)-wline)/wline*3.e5
         vel1=(wavz(i)-wavmx)/wavmx*3.e5
c--  rat- expected ratio given special relatativity kinematics
         rat=(1.-vel1/3.e5)**3/(1.+vel1/3.e5)**3
         fprev=fnor2
         fnor1=(xg(i,1)-fcz(i))/fmmx
         fnor2=yga(i)/fmmx
         fnor3=yga1(i)/fmmx
c
c find FWHM of model fit
c
         sign=(fnor2-0.5)*(fprev-0.5)
         if(sign.le.0.) then
          aa=(vel-(wavz(i-1)-wline)/wline*3.e5)/(fnor2-fprev)
          bb=vel-aa*fnor2
          if(v1.eq.0.) v1=aa*0.5+bb
          if(v1.ne.0.) v2=aa*0.5+bb
         endif
         write(23,'(7g12.4)') vel,vel1,fnor1,fnor2,fnor2/fnor3
     &     ,rat,wrat(i)
cc         write(23,'(f9.1,g12.4)') vel1,fnor2/fnor3
         write(28,'(f9.1,g12.4)') vel1,fnor2
         write(27,'(2g12.4)') vel,xg(i,1)-fcz(i)
 28     continue
       print*,'parameters of fitted model'
       ivshift=(wavmx-wline)/wline*3.e5
       ifwhm=v2-v1
       ewtot=flinetot/fcl(1)/(1.+z)
       write(6,127) wavmx,fcl(1),flinetot,ewtot,ivshift,ifwhm
cc       print*,'profile normalized to',fmmx,' FWHM=',v2-v1
cc       print*,'model maximum is at',wavmx
c
c  print velocity profile including only points used for the fit (i.e. not absorbed)
c
       do 37 i=1,ndd
        vel=(wav1(i)-wline)/wline*3.e5
        write(26,'(3g12.4)') vel,flz1(i),sig1(i)
 37    continue
cc        q=gammq(.5*(ndd-ma+nfix),.5*chisq)
cc        print*,'Pr for chi^2>',chisq,' for' 
cc     &   ,ndd-ma+nfix, ' d.o.f is ',q
      end
      subroutine median (x,n,xmed)
      real x(n)
      call sort(n,x)
      n2=n/2
      if(2*n2.eq.n) then
        xmed=.5*(x(n2)+x(n2+1))
       else
        xmed=x(n2+1)
       endif
       return
       end 
      subroutine sort(n,ra)
      real ra(n)
      l=n/2+1
      ir=n
10    continue
        if(l.gt.1)then
          l=l-1
          rra=ra(l)
        else
          rra=ra(ir)
          ra(ir)=ra(1)
          ir=ir-1
          if(ir.eq.1)then
            ra(1)=rra
            return
          endif
        endif
        i=l
        j=l+l
20      if(j.le.ir)then
          if(j.lt.ir)then
            if(ra(j).lt.ra(j+1))j=j+1
          endif
          if(rra.lt.ra(j))then
            ra(i)=ra(j)
            i=j
            j=j+j
          else
            j=ir+1
          endif
        go to 20
        endif
        ra(i)=rra
      go to 10
      end
      SUBROUTINE MRQMIN(X,Y,SIG,NDATA,A,MA,LISTA,MFIT,
     *    COVAR,ALPHA,NCA,CHISQ,FUNCS,ALAMDA)
      PARAMETER (MMAX=20)
      DIMENSION X(NDATA),Y(NDATA),SIG(NDATA),A(MA),LISTA(MA),
     *  COVAR(NCA,NCA),ALPHA(NCA,NCA),ATRY(MMAX),BETA(MMAX),DA(MMAX)
      IF(ALAMDA.LT.0.)THEN
        KK=MFIT+1
        DO 12 J=1,MA
          IHIT=0
          DO 11 K=1,MFIT
            IF(LISTA(K).EQ.J)IHIT=IHIT+1
11        CONTINUE
          IF (IHIT.EQ.0) THEN
            LISTA(KK)=J
            KK=KK+1
          ELSE IF (IHIT.GT.1) THEN
            PAUSE 'Improper permutation in LISTA'
          ENDIF
12      CONTINUE
        IF (KK.NE.(MA+1)) PAUSE 'Improper permutation in LISTA'
        ALAMDA=0.001
        CALL MRQCOF(X,Y,SIG,NDATA,A,MA,LISTA,MFIT,ALPHA,BETA,NCA,CHISQ,F
     *UNCS)
        OCHISQ=CHISQ
        DO 13 J=1,MA
          ATRY(J)=A(J)
13      CONTINUE
      ENDIF
      DO 15 J=1,MFIT
        DO 14 K=1,MFIT
          COVAR(J,K)=ALPHA(J,K)
14      CONTINUE
        COVAR(J,J)=ALPHA(J,J)*(1.+ALAMDA)
        DA(J)=BETA(J)
15    CONTINUE
      CALL GAUSSJ(COVAR,MFIT,NCA,DA,1,1)
      IF(ALAMDA.EQ.0.)THEN
        CALL COVSRT(COVAR,NCA,MA,LISTA,MFIT)
        RETURN
      ENDIF
      DO 16 J=1,MFIT
        ATRY(LISTA(J))=A(LISTA(J))+DA(J)
16    CONTINUE
c
c-- following condition prevents negative gaussians
c
      if(atry(1).lt.0.) atry(1)=atry(1)/100.
      if(atry(4).lt.0.) atry(4)=atry(4)/100.
      CALL MRQCOF(X,Y,SIG,NDATA,ATRY,MA,LISTA,MFIT,COVAR,DA,NCA,CHISQ,FU
     *NCS)
      IF(CHISQ.LT.OCHISQ)THEN
        ALAMDA=0.1*ALAMDA
        OCHISQ=CHISQ
        DO 18 J=1,MFIT
          DO 17 K=1,MFIT
            ALPHA(J,K)=COVAR(J,K)
17        CONTINUE
          BETA(J)=DA(J)
          A(LISTA(J))=ATRY(LISTA(J))
18      CONTINUE
      if(a(1).lt.0.) a(1)=a(1)/100.
      if(a(4).lt.0.) a(4)=a(4)/100.
      ELSE
        ALAMDA=10.*ALAMDA
        CHISQ=OCHISQ
      ENDIF
      RETURN
      END
      SUBROUTINE MRQCOF(X,Y,SIG,NDATA,A,MA,LISTA,MFIT,ALPHA,BETA,NALP,CH
     *ISQ,FUNCS)
      PARAMETER (MMAX=20)
      external funcs
      DIMENSION X(NDATA),Y(NDATA),SIG(NDATA),ALPHA(NALP,NALP),BETA(MA),
     *    DYDA(MMAX),LISTA(MFIT),A(MA)
      DO 12 J=1,MFIT
        DO 11 K=1,J
          ALPHA(J,K)=0.
11      CONTINUE
        BETA(J)=0.
12    CONTINUE
      CHISQ=0.
      DO 15 I=1,NDATA
        CALL FUNCS(X(I),A,YMOD,DYDA,MA)
        SIG2I=1./(SIG(I)*SIG(I))
        DY=Y(I)-YMOD
        DO 14 J=1,MFIT
          WT=DYDA(LISTA(J))*SIG2I
          DO 13 K=1,J
            ALPHA(J,K)=ALPHA(J,K)+WT*DYDA(LISTA(K))
13        CONTINUE
          BETA(J)=BETA(J)+DY*WT
14      CONTINUE
        CHISQ=CHISQ+DY*DY*SIG2I
15    CONTINUE
      DO 17 J=2,MFIT
        DO 16 K=1,J-1
          ALPHA(K,J)=ALPHA(J,K)
16      CONTINUE
17    CONTINUE
      RETURN
      END
      SUBROUTINE COVSRT(COVAR,NCVM,MA,LISTA,MFIT)
      DIMENSION COVAR(NCVM,NCVM),LISTA(MFIT)
      DO 12 J=1,MA-1
        DO 11 I=J+1,MA
          COVAR(I,J)=0.
11      CONTINUE
12    CONTINUE
      DO 14 I=1,MFIT-1
        DO 13 J=I+1,MFIT
          IF(LISTA(J).GT.LISTA(I)) THEN
            COVAR(LISTA(J),LISTA(I))=COVAR(I,J)
          ELSE
            COVAR(LISTA(I),LISTA(J))=COVAR(I,J)
          ENDIF
13      CONTINUE
14    CONTINUE
      SWAP=COVAR(1,1)
      DO 15 J=1,MA
        COVAR(1,J)=COVAR(J,J)
        COVAR(J,J)=0.
15    CONTINUE
      COVAR(LISTA(1),LISTA(1))=SWAP
      DO 16 J=2,MFIT
        COVAR(LISTA(J),LISTA(J))=COVAR(1,J)
16    CONTINUE
      DO 18 J=2,MA
        DO 17 I=1,J-1
          COVAR(I,J)=COVAR(J,I)
17      CONTINUE
18    CONTINUE
      RETURN
      END
      SUBROUTINE FGAUSS(X,A,Y,DYDA,NA)
      DIMENSION A(NA),DYDA(NA)
      Y=0.
      DO 11 I=1,NA-1,3
        ARG=(X-A(I+1))/A(I+2)
        EX=EXP(-ARG**2)
        FAC=A(I)*EX*2.*ARG
        Y=Y+A(I)*EX
        DYDA(I)=EX
        DYDA(I+1)=FAC/A(I+2)
        DYDA(I+2)=FAC*ARG/A(I+2)
11    CONTINUE
      RETURN
      END
      SUBROUTINE GAUSSJ(A,N,NP,B,M,MP)
      PARAMETER (NMAX=50)
      DIMENSION A(NP,NP),B(NP,MP),IPIV(NMAX),INDXR(NMAX),INDXC(NMAX)
      DO 11 J=1,N
        IPIV(J)=0
11    CONTINUE
      DO 22 I=1,N
        BIG=0.
        DO 13 J=1,N
          IF(IPIV(J).NE.1)THEN
            DO 12 K=1,N
              IF (IPIV(K).EQ.0) THEN
                IF (ABS(A(J,K)).GE.BIG)THEN
                  BIG=ABS(A(J,K))
                  IROW=J
                  ICOL=K
                ENDIF
              ELSE IF (IPIV(K).GT.1) THEN
                PAUSE 'Singular matrix'
              ENDIF
12          CONTINUE
          ENDIF
13      CONTINUE
        IPIV(ICOL)=IPIV(ICOL)+1
        IF (IROW.NE.ICOL) THEN
          DO 14 L=1,N
            DUM=A(IROW,L)
            A(IROW,L)=A(ICOL,L)
            A(ICOL,L)=DUM
14        CONTINUE
          DO 15 L=1,M
            DUM=B(IROW,L)
            B(IROW,L)=B(ICOL,L)
            B(ICOL,L)=DUM
15        CONTINUE
        ENDIF
        INDXR(I)=IROW
        INDXC(I)=ICOL
        IF (A(ICOL,ICOL).EQ.0.) PAUSE 'Singular matrix.'
        PIVINV=1./A(ICOL,ICOL)
        A(ICOL,ICOL)=1.
        DO 16 L=1,N
          A(ICOL,L)=A(ICOL,L)*PIVINV
16      CONTINUE
        DO 17 L=1,M
          B(ICOL,L)=B(ICOL,L)*PIVINV
17      CONTINUE
        DO 21 LL=1,N
          IF(LL.NE.ICOL)THEN
            DUM=A(LL,ICOL)
            A(LL,ICOL)=0.
            DO 18 L=1,N
              A(LL,L)=A(LL,L)-A(ICOL,L)*DUM
18          CONTINUE
            DO 19 L=1,M
              B(LL,L)=B(LL,L)-B(ICOL,L)*DUM
19          CONTINUE
          ENDIF
21      CONTINUE
22    CONTINUE
      DO 24 L=N,1,-1
        IF(INDXR(L).NE.INDXC(L))THEN
          DO 23 K=1,N
            DUM=A(K,INDXR(L))
            A(K,INDXR(L))=A(K,INDXC(L))
            A(K,INDXC(L))=DUM
23        CONTINUE
        ENDIF
24    CONTINUE
      RETURN
      END
      FUNCTION GAMMQ(A,X)
      IF(X.LT.0..OR.A.LE.0.)PAUSE
      IF(X.LT.A+1.)THEN
        CALL GSER(GAMSER,A,X,GLN)
        GAMMQ=1.-GAMSER
      ELSE
        CALL GCF(GAMMCF,A,X,GLN)
        GAMMQ=GAMMCF
      ENDIF
      RETURN
      END
      SUBROUTINE GSER(GAMSER,A,X,GLN)
      PARAMETER (ITMAX=100,EPS=3.E-7)
      GLN=GAMMLN(A)
      IF(X.LE.0.)THEN
        IF(X.LT.0.)PAUSE
        GAMSER=0.
        RETURN
      ENDIF
      AP=A
      SUM=1./A
      DEL=SUM
      DO 11 N=1,ITMAX
        AP=AP+1.
        DEL=DEL*X/AP
        SUM=SUM+DEL
        IF(ABS(DEL).LT.ABS(SUM)*EPS)GO TO 1
11    CONTINUE
      PAUSE 'A too large, ITMAX too small'
1     GAMSER=SUM*EXP(-X+A*LOG(X)-GLN)
      RETURN
      END
      SUBROUTINE GCF(GAMMCF,A,X,GLN)
      PARAMETER (ITMAX=100,EPS=3.E-7)
      GLN=GAMMLN(A)
      GOLD=0.
      A0=1.
      A1=X
      B0=0.
      B1=1.
      FAC=1.
      DO 11 N=1,ITMAX
        AN=FLOAT(N)
        ANA=AN-A
        A0=(A1+A0*ANA)*FAC
        B0=(B1+B0*ANA)*FAC
        ANF=AN*FAC
        A1=X*A0+ANF*A1
        B1=X*B0+ANF*B1
        IF(A1.NE.0.)THEN
          FAC=1./A1
          G=B1*FAC
          IF(ABS((G-GOLD)/G).LT.EPS)GO TO 1
          GOLD=G
        ENDIF
11    CONTINUE
      PAUSE 'A too large, ITMAX too small'
1     GAMMCF=EXP(-X+A*ALOG(X)-GLN)*G
      RETURN
      END
      FUNCTION GAMMLN(XX)
      REAL*8 COF(6),STP,HALF,ONE,FPF,X,TMP,SER
      DATA COF,STP/76.18009173D0,-86.50532033D0,24.01409822D0,
     *    -1.231739516D0,.120858003D-2,-.536382D-5,2.50662827465D0/
      DATA HALF,ONE,FPF/0.5D0,1.0D0,5.5D0/
      X=XX-ONE
      TMP=X+FPF
      TMP=(X+HALF)*LOG(TMP)-TMP
      SER=ONE
      DO 11 J=1,6
        X=X+ONE
        SER=SER+COF(J)/X
11    CONTINUE
      GAMMLN=TMP+LOG(STP*SER)
      RETURN
      END
       subroutine dred(z,w,f,sig,nd,ext)
       parameter(np=5904)
       real w(np),f(np),sig(np),e(26) 
c
c  Seaton 1979 (MN, 187, 73p) extinction curve
c
      data e/.14,.27,.41,.55,.69,.83,.97,1.18,1.36,1.44,1.84,2.04,
     1  2.24, 2.44, 2.66, 2.88, 3.14, 3.36, 3.56, 3.77, 3.96, 4.15,
     2 4.26, 4.4, 4.52, 4.64/
       do 1 i=1,nd
        xw=10000./(w(i)*(1.+z))
       if(xw.le.2.7)then
        do 19 j=1,25
         d1=xw-(.2+(j-1)*.1)
         d2=xw-(.2+j*.1)
         if(d1*d2.le.0)then
           d1=d1/.1
           d2=-d2/.1
           ax=.4*ext*(d1*e(j+1)+d2*e(j))
           goto 12
         endif
 19     continue
        endif
        if(xw.gt.2.7.and.xw.le.3.65)then
          if(xw.lt.3) then
            ax=.4*ext*(1.56+1.048*xw+1.01/((xw-4.6)**2+.28))
          else
            ax=.4*ext*(1.56+1.048*xw+1.01/((xw-4.6)**2+.28))
          endif
         else
           if(xw.gt.3.65.and.xw.le.7.14)then
              ax=.4*ext*(2.29+.848*xw+1.01/((xw-4.6)**2+.28))
             else
               if(xw.gt.7.14.and.xw.le.10)then
                 ax=.4*ext*(16.17-3.2*xw+.297*xw*xw)
                else
                 ax=0.
                endif
           endif
        endif
  12   f(i)=f(i)*10**ax
       sig(i)=sig(i)*10**ax
  1    continue
       return
         end
      subroutine fe2(wrk,nd,z)
       parameter (ne=6000,me=30)
       real wrk(ne,3),wv(ne),dat(ne),fe(1834,2), feg(1834)
       read(55,*) (fe(i,1),fe(i,2),i=1,1834)
c
c   scale and broaden fe 
c   xs=v/c
c
       print*,'type in a and xs'
       read*,a,xs 
       xs=xs/3.e5 
       do 1 i=1,1834 
        feg(i)=0.
        wt=0. 
        sig=xs*fe(i,1) 
        ngaus=3.*sig/(fe(i+1,1)-fe(i,1)) 
       do 2 j=-ngaus,ngaus
        if(i+j.lt.1.or.i+j.gt.1834) goto 2
        sg2=((fe(i,1)-fe(i+j,1))/sig)**2
        w=.39894*exp(-sg2/2.) 
        wt=wt+w 
        feg(i)=feg(i)+fe(i+j,2)*w
 2    continue
       feg(i)=feg(i)/wt*a
cc       write(27,*) fe(i,1),fe(i,2),feg(i)
 1    continue
      do 3 i=1,nd
       wv(i)=wrk(i,1)/(1.+z)
       if(wv(i).lt.4250..or.wv(i).gt.6999.5) goto 3
       xn=(wv(i)-4250.)/1.5
       n=xn
       fegint=feg(n+1)*(xn-n)+feg(n)*(n+1.-xn)
       dat(i)=wrk(i,2)-fegint
cc       write(27,*) wv(i),wrk(i,2),dat(i),fegint
       write(27,*) wrk(i,1),dat(i)
 3    continue
      return
      end