FUNCTION filter(i,x,y,n,z)
 
c...... Returns flux through i filter. The SED is y(x) with n data
c...... points. z is the redshift to be applied to the SED

      INCLUDE 'filter.dec'
      REAL rlin,xlam(mxwf),ff(mxwf),x(1),y(1)
      EXTERNAL rlin,trapz1
      filter=0.
      IF (i.gt.mxft) RETURN
      if (iread.eq.0) CALL FILTER0
      l=0
      DO m=1,np(i)
         xlam(m)=r(m,i,1)
         ff(m)=r(m,i,2)*rlin(xlam(m)/(1.+z),X,Y,n,l)
      END DO
      filter=trapz1(xlam,ff,np(i))/(1.+z)
      RETURN
      END

      FUNCTION color(i,j,x,y,n,z)

c......	Returns color for filter pair (i,j) for SED in (x,y) with n
c...... data points observed at redshift z. The zeropoint is not added
c...... and must be added in the calling program.

      INCLUDE 'filter.dec'	
      DIMENSION a(2),iz(2),x(1),y(1)
      iz(1)=i
      iz(2)=j
      DO k=1,2
         a(k)=filter(iz(k),x,y,n,z)
         if (a(k).le.0.) GO TO 1
      END DO
      color=-2.5*LOG10(a(1)/a(2))
      RETURN
c     used to be 1.e11
 1    color=111.	
      RETURN
      END

      SUBROUTINE FILTER0

c ..... Reads filter response functions
c ..... rebinning if n. points > 200

      INCLUDE 'filter.dec'
      INTEGER mxx
      PARAMETER (mxx=10000)
      REAL wh(mxx),res(mxx),rlin
      CHARACTER filtresfile*80,line*80
      EXTERNAL rlin
      COMMON/filterres/filtresfile
      COMMON/compose/icomp
      
      CLOSE(81)
      OPEN(81,FILE=filtresfile,STATUS='old')
      i=1
 101  READ(81,'(a)',ERR=102,END=102) line
      IF (line.eq.' ') GO TO 101
      READ(line,*) np(i)
      nf=i
      IF (np(i).gt.mxx) THEN
         WRITE(*,*) 'Too many points for filter  ', i
         WRITE(*,*) 'Maximum number of points is ', mxx
         STOP
      END IF
      DO j=1,np(i)
         READ(81,*,ERR=1) l,wh(j),res(j)
      END DO
      IF (np(i).gt.200) THEN    ! rebinning
         dlambda=wh(np(i))-wh(1)
         ii=1
 103     dx=ii*5.
         new=INT(dlambda/dx)+1
         IF (new.le.200) GO TO 104
         ii=ii+1
         GO TO 103
 104     DO j=1,new
            r(j,i,1)=wh(1)+j*dx
            ll=0
            r(j,i,2)=rlin(r(j,i,1),wh,res,np(i),ll)
         END DO
         np(i)=new
      ELSE
         DO j=1,np(i)
            r(j,i,1)=wh(j)
            r(j,i,2)=res(j)
         END DO
      END IF
      i=i+1
      GO TO 101
 102  CLOSE(81)
      iread=1
      icomp=i-1
      RETURN
 1    STOP 'Program exits because of error reading file FILTER.RES'
      END