	character*80 inline,inline2
	character*1 filt
	parameter (maxnum=30)
	dimension x(8,5,maxnum),xe(8,5,maxnum),xu(8,maxnum),xue(8,maxnum)
        dimension xs(maxnum),sig(maxnum)
	dimension num(8,5)
	do 101 ichip=1,8
	do 101 ifilt=1,5
	num(ichip,ifilt)=0
101	continue
	open(unit=1,file='fitresults',status='old')
	open(unit=2,file='calcoefsout',status='new')
1	read(1,2,end=99) inline
2	format(a80)
	do 3 j=2,30
	if(inline(j:j).eq.' ') go to 4
3	continue
	stop 'problem readinf fitresults'
4	read(inline(j-1:j-1),*) ichip
	read(inline(j+1:j+1),5) filt
5	format(a1)
	if(filt.eq.'U') ifilt=1
	if(filt.eq.'B') ifilt=2
	if(filt.eq.'V') ifilt=3
	if(filt.eq.'R') ifilt=4
	if(filt.eq.'I') ifilt=5
	num(ichip,ifilt)=num(ichip,ifilt)+1
	if(num(ichip,ifilt).gt.maxnum) stop 'too many frames'
	read(inline(j+2:80),*) zeropt,zeroerr,color,colerr,nfit,rfit,
     &   iflag,sd
	x(ichip,ifilt,num(ichip,ifilt))=color
	xe(ichip,ifilt,num(ichip,ifilt))=colerr
	if(filt.eq.'U') then
	read(1,2) inline2
	read(inline2(j+2:80),*) u2,u2err,nfit2,rfit2,iflag2,sd2
	xu(ichip,num(ichip,1))=u2
	xue(ichip,num(ichip,1))=u2err
	end if
	go to 1
99	continue
	do 201 ichip=1,8
	do 201 ifilt=1,5
	ns=num(ichip,ifilt)
	if(ns.lt.3) go to 201
	do 203 inum=1,ns
	xs(inum)=x(ichip,ifilt,inum)
	sig(inum)=xe(ichip,ifilt,inum)
203	continue
	call xfit(xs,sig,ns,1,ave,sigmam,sigma)
	call median(xs,ns,value)
	write(2,204) ichip,ifilt,ave,sigmam,value
204	format('Chip=',i2,' Filt=',i2,' Color term:',f8.3,'+/-',f8.3,
     &  ' Median=',f8.3)
	if(ifilt.eq.1) then
	do 205 inum=1,ns
	xs(inum)=xu(ichip,inum)
	sig(inum)=xue(ichip,inum)
205	continue
	call xfit(xs,sig,ns,1,ave,sigmam,sigma)
	call median(xs,ns,value)
	write(2,204) ichip,ifilt,ave,sigmam,value
	end if
201	continue
	stop
	end
c subroutine xfit.f
c
c source
c   Bevington, page 76.
c
c purpose
c   calculate the mean and estimated errors for a set of data points
c
c usage
c   call xfit (x, sigmax, npts, mode, xmean, sigmam, sigma)
c
c description of parameters
c   x      - array of data points
c   sigmax - array of standard deviations for data points
c   npts   - number of data points
c   mode   - determines method of weighting
c            +1 (instrumental) weight(i) = 1./sigmax(i)**2
c             0 (no weighting) weight(i) = 1.
c            -1 (statistical)  weight(i) = 1.
c   xmean  - weighted mean
c   sigmam - standard deviation of mean
c   sigma  - standard deviation of data
c
c subroutines and function subprograms required
c   none
c
	subroutine xfit (x,sigmax,npts,mode,xmean,sigmam,sigma)
	double precision sum,sumx,weight,free
	dimension x(*),sigmax(*)
c
c accumulate weighted sums
c
11	sum=0.
	sumx=0.
	sigma=0.
	sigmam=0.
20	do 32 i=1,npts
21	if (mode) 22,22,24
22	weight=1.
	goto 31
24	weight=1./sigmax(i)**2
31	sum=sum+weight
32	sumx=sumx+weight*x(i)
c
c evaluate mean and standard deviations
c
41	xmean=sumx/sum
51	do 52 i=1,npts
52	sigma=sigma+(x(i)-xmean)**2
	free=npts-1
54	sigma=dsqrt(sigma/free)
61	if (mode) 62,64,66
62	sigmam=dsqrt(xmean/sum)
	goto 70
64	sigmam=sigma/dsqrt(sum)
	goto 70
66	sigmam=dsqrt(1./sum)
70	return
	end
        subroutine median(data,n,amed)
        dimension data(*),indx(20000)
        call sort(n,data,indx)
        if( ((n/2)*2).eq.n) then
        amed=(data(indx(n/2))+data(indx(n/2+1)))/2.
        else
        amed=data(indx(n/2+1))
        end if
        return
        end
        subroutine sort(n,arrin,indx)
        dimension arrin(*),indx(*)
        do 11 j=1,n
        indx(j)=j
11      continue
        if(n.eq.1) return
        l=n/2+1
        ir=n
10      continue
                if(l.gt.1) then
                        l=l-1
                        indxt=indx(l)
                        q=arrin(indxt)
                ELSE
                        indxt=indx(ir)
                        q=arrin(indxt)
                        indx(ir)=indx(1)
                        ir=ir-1
                        if(ir.eq.1) then
                                indx(1)=indxt
                                return
                        endif

                endif
        i=l
        j=l+l
20      if(j.le.ir) then
                if(j.lt.ir) then
                        if(arrin(indx(j)).lt.arrin(indx(j+1)))j=j+1
                endif
                if(q.lt.arrin(indx(j)))then
                        indx(i)=indx(j)
                        i=j
                        j=j+j
                else
                        j=ir+1
                endif
        go to 20
        endif
        indx(i)=indxt
        go to 10
        end


