program composite
c
c  This program averages spectra to produce a composite.
c  Averaging is either of log f_nu, or of dlog f_nu/dnu,
c  i.e. the local spectral slope (set islope=1).
c  nob- number of objects to be averaged.
c  nm- number of wavelengths for re-sampling of spectra.
c  Note that a large oversampling is required for islope=1  (nm~10,000)
c  It is also important to smooth the spectra before making a
c  slope composite, otherwise the noise in the composite is too high.
c  Input: unit 1, Output: unit 2
c
       parameter (nob=1,nm=10000,nn=nob*nm)
       real wav(nob,6000),fl(nob,6000),wavm(nm),fln(nob,nm)
     & ,flcmp(nm),x(nob),nspec(nm),sn(nob,6000),fav(nob)
     & ,snn(nob,nm),z(nob),ext(nob),xlum(nob)
       integer np(nob)
       common /red/wav,fl
       data fln,fav/nn*0.,nob*0./
c 
c  islope=1 use slope averaging, islope=0 use log(flux) averaging
c 
        islope=1
c
c  iaverage=1- calculate average and dispersion (for sigma clipping)
c  iaverage=0- calculate median (a generally noisier option)
c
        iaverage=1
c
c  xl1, xl2 wavelength range for the composite
c
       xl1=1150.
       xl2=2200.
c
c  Load all spectra from unit 1
c  all objects are in the same file
c  first line for each object includes np = number of data points,
c  z- the redshift, ext- E(B-V)
c  data in format lambda_observed(A), f_lambda, flux error_lambda
c
       do 1 i=1,nob
        read(1,*) np(i),z(i),ext(i)
        do 1 j=1,np(i)
cc         read(1,*) wav(i,j),fl(i,j),err
         read(1,*) wav(i,j),x,x,x,x,fl(i,j),x,x
         sn(i,j)=fl(i,j)/err
         if(fl(i,j).lt.0.) fl(i,j)=-fl(i,j)/100.
c
c  uncomment to check input spectrum
c
         write(3,*) wav(i,j)/(1.+z(i)),fl(i,j)
 1     continue
c
c  dred- dereddens all spectra
c
       call dred(np,ext,nob)
       write(2,'(2g14.4)')(wav(i,j)/(1.+z(i)),fl(i,j),j=1,1454)
       stop
c
c  smooth- smooths the spectra with a Gaussian equal to the wavelength resolution
c
       call smooth(fl,sn,np)
c
c  calculate absolute luminosity (nuL_nu) at 1350A, 
c
c       do 12 i=1,nob
        f1350=0.
        do 14 j=1,np(i)
c
c  change to rest frame wavelength
c
         wav(i,j)=wav(i,j)/(1.+z(i))
         fl(i,j)=alog10( fl(i,j) )
         if(j.eq.1) goto 14
          if(wav(i,j).gt.1350.and.wav(i,j-1).le.1350) f1350=fl(i,j)
 14     continue
c
c  uncomment following line for q_0=0.5, otherwise q_0=0.
c
cc       zc = 2*( 1+z(i)-sqrt(1+z(i)) )
       zc=z(i)*(1.+z(i)/2.)
c
c  H_0=75 ---> D=1.232E28  4piD**2 = 10**57.280
c  H_0=50 ---> D=1.848E28  4piD**2 = 10**57.633
c
       ac = 57.633 + alog10( zc*zc )
       xlum(i)=f1350+alog10(1350.*(1.+z(i)))+ac
cc       if(f1350.eq.0) print*,i,' f_1350 not available'
cc       print*,i,z(i),xlum(i),f1350
 12    continue
c
c  set wavelengths for resampling
c
       dl=(xl2/xl1)**(1./(nm-1.))
       do 2 i=1,nm
 2      wavm(i)=xl1*dl**(i-1.) 
c
c  resample all spectra at wavm
c
       do 3 i=1,nob
        kz=1
        do 4 j=1,nm
         do 5  k=kz+1,np(i) 
          if(wav(i,k).lt.wavm(j)) goto 5
          if(wav(i,k-1).gt.wavm(j)) goto 4
          if(wav(i,k).ge.wavm(j).and.wav(i,k-1).le.wavm(j)) then
            kz=k-1
c 
c  for islope=1 fln is local slope, for islope=0 fln is local flux
c 
            if(islope.eq.0) then    
             a=(fl(i,k)-fl(i,k-1))/(wav(i,k)-wav(i,k-1))
             b=fl(i,k)-a*wav(i,k)
             fln(i,j)=a*wavm(j)+b
            else
             fratio= fl(i,k)-fl(i,k-1)
             wratio= wav(i,k)/wav(i,k-1)
             fln(i,j)=fratio/alog(wratio)
            endif
c
c  do not include regions with a large gap in wavelength
c
             if(wav(i,k)/wav(i,k-1).gt.1.01) fln(i,j)=0.
c
c  snn- the S/N of the uniformly resampled spectrum
c
            snn(i,j)=(sn(i,k)+sn(i,k-1))/2.
c
c   to produce composites for various subsets of objects
c
cc            if(z(i).le..6) fln(i,j)=0.
cc           if(xlum(i).gt.45.96) fln(i,j)=0.
cc           if(xlum(i).lt.46.65) fln(i,j)=0.
cc            if(xlum(i).le.45.96.or.xlum(i).ge.46.65) fln(i,j)=0.
cc            if(snn(i,j).lt.5.) fln(i,j)=0.
c
c  delete points near Mg II and Ly alpha Galactic absorption, 
c
            wob=wavm(j)*(1.+z(i)) 
            if(wob.gt.2795..and.wob.lt.2822.) fln(i,j)=0.
            if(wob.gt.1210.7.and.wob.lt.1220.7) fln(i,j)=0.
            if(wob.gt.2595..and.wob.lt.2605.) fln(i,j)=0.
            if(wob.gt.2581..and.wob.lt.2591.) fln(i,j)=0.
            goto 4
          endif
 5       continue
 4      continue
 3     continue 
c--  uncomment to print all resampled spectra (or slopes) 
cc      write(2,'(2f12.4)')((wavm(j),fln(i,j),j=1,nm),
cc     &    i=1,nob)
cc      stop
      fint=1.
c
c  DO 6 runs to the end of the main program (each wavm treated independently) 
c
      do 6 j=1,nm
        ndat=0.
        weight=0. 
        inoz=0
        xav=0.
        sig2=0.
c
c  ndat- # of objects with data at wavm(j), xav- average flux, sig- dispersion
c
        do 7 i=1,nob
         if(fln(i,j).ne.0.) then
          inoz=i
          ndat=ndat+1
          x(ndat)=fln(i,j)
          xav=xav+fln(i,j)*snn(i,J)
          sig2=sig2+fln(i,j)**2*snn(i,j)
          weight=weight+snn(i,j)
         endif
 7      continue
        nspec(j)=ndat
        if(ndat.gt.1.and.iaverage.eq.1) then
         xav=xav/weight
         sig2=sig2-weight*xav**2
         sig=sqrt(sig2*ndat/weight/(ndat-1.))
c
c  recalculate average excluding points xsig units of sigma away
c
         xsig=2.
         ndat=0
         xavnew=0.
         weight=0. 
         do 11 i=1,nob
          if(fln(i,j).eq.0.) goto 11
          dev=( (fln(i,j)-xav)/sig )**2
          if(dev.gt.xsig**2) goto 11
          ndat=ndat+1
          xavnew=xavnew+fln(i,j)*snn(i,j)
          weight=weight+snn(i,j)
 11      continue
         xav=xavnew/weight
       endif
       xmed=0. 
c
c  if iaverage=0 use median flux rather than average.
c
       if(ndat.ge.2.and.iaverage.eq.0) then 
         call median(x,ndat,xmed)
         flcmp(j)=xmed
       endif
       if(ndat.ge.2.and.iaverage.eq.1) flcmp(j)=xav
       if(islope.eq.0) then
         if(ndat.eq.1) flcmp(j)=fln(inoz,j)
         if(ndat.eq.0) flcmp(j)=0.
       endif
c
c  if islope=1, Integrate slope spectrum to get back flux spectrum.
c
       if(islope.eq.1) then       
        if(j.gt.1)fint=fint+flcmp(j)*alog( wavm(j)/wavm(j-1) )
        write(2,'(f7.2,f9.3,i4)') wavm(j),fint, ndat
       endif
       if(islope.eq.0)write(2,'(f7.2,f9.3,i4)')wavm(j),flcmp(j),ndat
 6    continue 
      end
c
c  Subroutines for finding median, deredening, and smoothing.
c
      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 dred(np,ext)
       parameter (nob=1)
       real e(26),wav(nob,6000),fl(nob,6000),ext(nob) 
       integer np(nob)
       common /red/ wav,fl
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,nob
       do 1 j=1,np(i)
        xw=10000./wav(i,j)
       if(xw.le.2.7)then
        do 19 k=1,25
         d1=xw-(.2+(k-1)*.1)
         d2=xw-(.2+k*.1)
         if(d1*d2.le.0)then
           d1=d1/.1
           d2=-d2/.1
           ax=.4*ext(i)*(d1*e(k+1)+d2*e(k))
           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(i)*(1.56+1.048*xw+1.01/((xw-4.6)**2+.28))
          else
            ax=.4*ext(i)*(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(i)*(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(i)*(16.17-3.2*xw+.297*xw*xw)
                else
                 ax=0.
                endif
           endif
        endif
  12   fl(i,j)=fl(i,j)*10**ax
  1    continue
       return
       end
      subroutine smooth(fl,sn,np)
c
c  Gaussian smoothing (for UV data only), smoothing is
c  with a sigma of 1.75 pixels (FWHM=4 pixels)
c
      parameter (nob=1)
      real fl(nob,6000),sn(nob,6000),xg(6000,2)
      integer np(nob)
      ngaus=4
      sigt=0.
      do 1 i=1,nob
       do 2 j=1,np(i)
       xg(j,1)=0.
       xg(j,2)=0.
       wt=0.
       do 3 k=-ngaus,ngaus
        if(j+k.lt.1.or.j+k.gt.np(i)) goto 3
        si=float(k)/1.75
        w=.39894*exp(-si**2/2.) 
        wt=wt+w 
        xg(j,1)=xg(j,1)+fl(i,j+k)*w
        xg(j,2)=xg(j,2)+sn(i,j+k)*w
 3    continue
       xg(j,1)=xg(j,1)/wt
       xg(j,2)=xg(j,2)/wt
 2    continue
      do 4 j=1,np(i)
       if(xg(j,1).ne.0) fl(i,j)=xg(j,1)
       if(xg(j,2).ne.0) sn(i,j)=xg(j,2)
 4    continue
 1    continue
      return
      end