C Checked and the sigma problem debugged June 10, 2002
C Added better support for the issue of "too few stars" in icfit 6/27/03
	program fitem
C written Phil Massey 5 April 2002
	character*1 filt
	dimension ifa(8),va(8)
	character*80 inline
	open(unit=1,file='tempmatchem.dat',status='old')
	open(unit=2,file='tempcoefs.dat',status='new')
	open(unit=11,file='tempchip',status='old')
	read(11,*) inum
	close(unit=11)
	open(unit=9,
     &  file='/tofu/home/massey/survey/kpnocoefs.dat',status='old')
	read(1,2) filt
2	format(a1)
	rewind(1)
4	read(9,3,end=6) inline
3	format(a80)
	if(inline(1:1).eq.filt) go to 5
	go to 4
6	stop 'problem with filter name down in fitem'
5	read(inline(2:80),*)(ifa(ik),va(ik),ik=1,8)
	iflag=ifa(inum)
	value=va(inum)
	if(filt.eq.'U') then
	read(9,3) inline
	read(inline(2:80),*)(ifa(ik),va(ik),ik=1,8)
	iflag2=ifa(inum)
	value2=va(inum)
	end if
	close(unit=9)
	if(filt.eq.'B') then
	call bfit(a,sigmaa,b,sigmab,nfit,r,iflag,value,sd)
	else if (filt.eq.'V') then
	call vfit(a,sigmaa,b,sigmab,nfit,r,iflag,value,sd)
	else if (filt.eq.'R') then
	call rfit(a,sigmaa,b,sigmab,nfit,r,iflag,value,sd)
Cwrite(6,*) a,b
	else if (filt.eq.'I') then
	call ifit(a,sigmaa,b,sigmab,nfit,r,iflag,value,sd)
	else if (filt.eq.'U') then
	call u1fit(a,sigmaa,b,sigmab,nfit,r,iflag,value,sd)
		if(iflag2.eq.0) then
		rewind(1)
		value=a
C	write(6,123) value
123	format('About to call u2fit with value=',f6.2)
		call u2fit(value,uprob,uprobe,nfit2,ru,sdu)
		else
		uprob=value2
		sd=0.
		uprobe=0.0
		end if
	end if
	close(unit=1)
	write(6,1) filt,a,sigmaa,b,sigmab,nfit,r,iflag,sd
	write(2,1) filt,a,sigmaa,b,sigmab,nfit,r,iflag,sd
1	format(a1,1x,4f10.3,i4,f6.2,i2,f8.2)
11	format(a1,1x,2f10.3,i4,f6.2,i2,f8.2)
	if(filt.eq.'U') then
	write(6,11) filt,uprob,uprobe,nfit2,ru,iflag2,sdu
	write(2,11) filt,uprob,uprobe,nfit2,ru,iflag2,sdu
	end if
	stop
	end
	subroutine bfit(a,sigmaa,b,sigmab,numt,r,iflag,value,sd)
C maximum number of stars
	parameter(maxstars=250000)
C temporary arrays
	dimension tx(maxstars),ty(maxstars),tsig(maxstars)
	ik=0
1	read(1,206,end=99) dist,bi,bie,v,ve,bmv,bmve
206	format(2x,f6.1,12f7.3)
	if(bmv.gt.90.) go to 1
	if(v.gt.90.) go to 1
	if(bmv.lt.-.4) go to 1
	if(bmv.gt.2.8) go to 1
	stand=v+bmv
	standerr2=bmve**2-ve**2
	if(standerr2.le.25.0e-6) then
	stande=0.005
	else
	stande=sqrt(standerr2)
	end if
	total=sqrt(stande**2+bie**2)
	dy=bi-stand
	ik=ik+1
	if(ik.gt.maxstars) stop 'too many stars in fitem'
C Fill up the temporary arrays:
	tx(ik)=bmv
	ty(ik)=dy
	tsig(ik)=total
	go to 1
99	numt=ik
C OK,now do the fitting:
	call icfit(tx,ty,tsig,numt,a,sigmaa,b,sigmab,r,iflag,value,sd)
	return
	end
	subroutine vfit(a,sigmaa,b,sigmab,numt,r,iflag,value,sd)
C maximum number of stars
	parameter(maxstars=250000)
C temporary arrays
	dimension tx(maxstars),ty(maxstars),tsig(maxstars)
	ik=0
1	read(1,206,end=99) dist,bi,bie,v,ve,bmv,bmve
206	format(2x,f6.1,12f7.3)
	if(bmv.gt.90.) go to 1
	if(v.gt.90.) go to 1
	if(bmv.lt.-.4) go to 1
	if(bmv.gt.2.8) go to 1
	stand=v
	standerr2=ve**2
	if(standerr2.le.25.0e-6) then
	stande=0.005
	else
	stande=sqrt(standerr2)
	end if
	total=sqrt(stande**2+vie**2)
	dy=bi-stand
	ik=ik+1
	if(ik.gt.maxstars) stop 'too many stars in fitem'
C Fill up the temporary arrays:
	tx(ik)=bmv
	ty(ik)=dy
	tsig(ik)=total
	go to 1
99	numt=ik
C OK,now do the fitting:
	call icfit(tx,ty,tsig,numt,a,sigmaa,b,sigmab,r,iflag,value,sd)
	return
	end
	subroutine rfit(a,sigmaa,b,sigmab,numt,r,iflag,value,sd)
C maximum number of stars
	parameter(maxstars=250000)
C temporary arrays
	dimension tx(maxstars),ty(maxstars),tsig(maxstars)
	ik=0
1	read(1,206,end=99) dist,bi,bie,v,ve,bmv,bmve,umb,umbe,vmr,vmre
206	format(2x,f6.1,12f7.3)
	if(vmr.gt.90.) go to 1
	if(v.gt.90.) go to 1
	if(vmr.gt.2.0) go to 1
	if(vmr.lt.-.2) go to 1
	stand=v-vmr
	standerr2=vmre**2-ve**2
	if(standerr2.le.25.0e-6) then
	stande=0.005
	else
	stande=sqrt(standerr2)
	end if
	total=sqrt(stande**2+bie**2)
	dy=bi-stand
	ik=ik+1
	if(ik.gt.maxstars) stop 'too many stars in fitem'
C Fill up the temporary arrays:
	tx(ik)=vmr
	ty(ik)=dy
	tsig(ik)=total
	go to 1
99	numt=ik
C OK,now do the fitting:
	do 100 i=1,numt
Cwrite(19,*) tx(i),ty(i),tsig(i)
100	continue
	call icfit(tx,ty,tsig,numt,a,sigmaa,b,sigmab,r,iflag,value,sd)
Cwrite(6,*) a,b
	return
	end
	subroutine ifit(a,sigmaa,b,sigmab,numt,r,iflag,value,sd)
C maximum number of stars
	parameter(maxstars=250000)
C temporary arrays
	dimension tx(maxstars),ty(maxstars),tsig(maxstars)
	ik=0
1	read(1,206,end=99) dist,bi,bie,v,ve,bmv,bmve,umb,umbe,vmr,vmre,
     &  rmi,rmie
206	format(2x,f6.1,12f7.3)
	if(v.gt.90.) go to 1
	if(rmi.gt.90.) go to 1
	if(vmr.gt.90.) go to 1
	if(vmr.gt.2.) go to 1
	if(vmr.lt.-.2) go to 1
	if(rmi.lt.-0.5) go to 1
	if(rmi.gt.3.0) go to 1
	stand=v-vmr-rmi
C this may not be right
	standerr2=rmie**2-vmre**2-ve**2
	if(standerr2.le.25.0e-6) then
	stande=0.005
	else
	stande=sqrt(standerr2)
	end if
	total=sqrt(stande**2+bie**2)
	dy=bi-stand
	ik=ik+1
	if(ik.gt.maxstars) stop 'too many stars in fitem'
C Fill up the temporary arrays:
	tx(ik)=rmi
	ty(ik)=dy
	tsig(ik)=total
	go to 1
99	numt=ik
C OK,now do the fitting:
	call icfit(tx,ty,tsig,numt,a,sigmaa,b,sigmab,r,iflag,value,sd)
	return
	end
	subroutine u1fit(a,sigmaa,b,sigmab,numt,r,iflag,value,sd)
C maximum number of stars
	parameter(maxstars=250000)
C temporary arrays
	dimension tx(maxstars),ty(maxstars),tsig(maxstars)
	ik=0
1	read(1,206,end=99) dist,bi,bie,v,ve,bmv,bmve,umb,umbe,vmr,vmre,
     &  rmi,rmie
206	format(2x,f6.1,12f7.3)
	if(umb.gt.90.) go to 1
	if(v.gt.90.) go to 1
	if(umb.lt.0.) go to 1
	if(umb.gt.2.4) go to 1
	if(bmv.gt.90.) go to 1
	if(bmv.gt.2.8) go to 1
	if(bmv.lt.-0.4) go to 1
	stand=v+umb+bmv
C this may not be right
	standerr2=umbe**2-bmve**2-ve**2
	if(standerr2.le.25.0e-6) then
	stande=0.005
	else
	stande=sqrt(standerr2)
	end if
	total=sqrt(stande**2+bie**2)
	dy=bi-stand
	ik=ik+1
	if(ik.gt.maxstars) stop 'too many stars in fitem'
C Fill up the temporary arrays:
	tx(ik)=umb
	ty(ik)=dy
	tsig(ik)=total
	go to 1
99	numt=ik
C OK,now do the fitting:
	call icfit(tx,ty,tsig,numt,a,sigmaa,b,sigmab,r,iflag,value,sd)
	return
	end
	subroutine u2fit(value,uprob,uprobe,numt,r,sd)
C maximum number of stars
	parameter(maxstars=250000)
C temporary arrays
	dimension tx(maxstars),ty(maxstars),tsig(maxstars)
	iflag=-1
	write(6,11) value
11	format('In u2:a=',f8.2)
	ik=0
1	read(1,206,end=99) dist,bi,bie,v,ve,bmv,bmve,umb,umbe
206	format(2x,f6.1,12f7.3)
	if(v.gt.90.) go to 1
	if(umb.gt.90.) go to 1
	if(umb.gt.0.) go to 1
	if(umb.lt.-1.3) go to 1
	if(bmv.gt.90.) go to 1
	if(bmv.gt.2.8) go to 1
	if(bmv.lt.-.4) go to 1
	stand=v+umb+bmv
C this may not be right
	standerr2=umbe**2-bmve**2-ve**2
	if(standerr2.le.25.0e-6) then
	stande=0.005
	else
	stande=sqrt(standerr2)
	end if
	total=sqrt(stande**2+bie**2)
	dy=bi-stand
	ik=ik+1
	if(ik.gt.maxstars) stop 'too many stars in fitem'
C Fill up the temporary arrays:
	tx(ik)=umb
	ty(ik)=dy
	tsig(ik)=total
	go to 1
99	numt=ik
	if(numt.eq.0) then
	write(6,1001) 
1001	format('No data points')
	a=0.00
	sigmaa=0.0
	b=0.00
	sigmab=0.0
	r=0.0
	return
	end if
C OK,now do the fitting:
	call icfit(tx,ty,tsig,numt,a,sigmaa,b,sigmab,r,iflag,value,sd)
	write(6,1002) numt
1002	format('Number of data points:', i5)
	uprob=b
	uprobe=sigmab
	return
	end
	subroutine icfit(tx,ty,tsig,numt,a,sigmaa,b,sigmab,r,iflag,
     &   value,sd)
C added support for too few stars 6/27/03
	parameter(maxstars=250000)
	dimension tx(*),ty(*),tsig(*),
     &             x(maxstars),y(maxstars),sig(maxstars)
C if zero stars for get it!
	a=99.9
	sigmaa=99.9
	b=99.9
	sigmab=99.9
	r=99.9
	sd=99.9
	if(numt.lt.3) return
	call median(ty,numt,amed)
C find the sd:
	sum2=0.
	do 105 i=1,numt
	write(3,*) tx(i),ty(i),tsig(i)
	sum2=sum2+(ty(i)-amed)**2
105	continue
	sd=sqrt(sum2/numt)
	write(6,101) numt,amed,sd
101	format('Number of stars initial',i5,
     &   ' with median inst-stand=',f8.2, 'and sd=',f8.2)
C Now fill up the arrays and exclude any stars that are wildly deviant.
C My assumption here is that the largest color coefficent will be 0.3
C Over a range of 3 mags or about 1 mag is VERY GENEROUS
	ik=0
	do 102 i=1,numt
C Use 1 mag for the first cut!
	if(abs(ty(i)-amed).gt.0.5) go to 102
	ik=ik+1
	if(ik.gt.maxstars) stop 'too many stars in fitem/icfit'
	x(ik)=tx(i)
	y(ik)=ty(i)
	sig(ik)=tsig(i)
102	continue
	num=ik
	a=99.9
	sigmaa=99.9
	b=99.9
	sigmab=99.9
	r=99.9
	sd=99.9
	if(num.lt.3) return
	call linfit(x,y,sig,num,1,a,sigmaa,b,sigmab,r,iflag,value)
100	format('Num=',i5, ' ZEROPT=',f7.3,'+/-',f7.3,' COLORTRM=',
     &   f7.3,'+/-',f7.3,' R=',f6.2,' SD=',f6.2)
C Now make 3 passes
C First, find the sd:
	sum2=0.
	do 205 i=1,num
	sum2=sum2+(y(i)-a-b*x(i))**2
205	continue
	sd=sqrt(sum2/num)
	write(6,100) num,a,sigmaa,b,sigmab,r,sd
	ik=0
	do 202 i=1,num
	predict=a+b*x(i)
	if(abs(y(i)-predict).gt.3.*sd) go to 202
	ik=ik+1
	if(ik.gt.maxstars) stop 'too many stars in fitem/icfit'
	ty(ik)=y(i)
	tx(ik)=x(i)
	tsig(ik)=sig(i)
202	continue
	numt=ik
	a=99.9
	sigmaa=99.9
	b=99.9
	sigmab=99.9
	r=99.9
	sd=99.9
	if(numt.lt.3) return
	call linfit(tx,ty,tsig,numt,1,a,sigmaa,b,sigmab,r,iflag,value)
	sum2=0
	do 305 i=1,numt
	sum2=sum2+(ty(i)-a-b*tx(i))**2
305	continue
	sd=sqrt(sum2/numt)
	write(6,100) numt,a,sigmaa,b,sigmab,r,sd
C Pass 2:
	ik=0
	do 302 i=1,numt
	predict=a+b*tx(i)
	if(abs(ty(i)-predict).gt.2.5*sd) go to 302
	ik=ik+1
	if(ik.gt.maxstars) stop 'too many stars in fitem/icfit'
	y(ik)=ty(i)
	x(ik)=tx(i)
	sig(ik)=tsig(i)
302	continue
	num=ik
	a=99.9
	sigmaa=99.9
	b=99.9
	sigmab=99.9
	r=99.9
	sd=99.9
	if(num.lt.3) return
	call linfit(x,y,sig,num,1,a,sigmaa,b,sigmab,r,iflag,value)
	sum2=0.
	do 405 i=1,num
	sum2=sum2+(y(i)-a-b*x(i))**2
405	continue
	sd=sqrt(sum2/num)
	write(6,100) num,a,sigmaa,b,sigmab,r,sd
C Pass 3 and final
	ik=0
	do 402 i=1,num
	predict=a+b*x(i)
	if(abs(y(i)-predict).gt.2.5*sd) go to 402
	ik=ik+1
	if(ik.gt.maxstars) stop 'too many stars in fitem/icfit'
	ty(ik)=y(i)
	tx(ik)=x(i)
	tsig(ik)=sig(i)
	write(4,*) tx(ik),ty(ik),tsig(ik)
C write these in case we need to understand a failure
402	continue
	numt=ik
	a=99.9
	sigmaa=99.9
	b=99.9
	sigmab=99.9
	r=99.9
	sd=99.9
	if(numt.lt.3) return
	call linfit(tx,ty,tsig,numt,1,a,sigmaa,b,sigmab,r,iflag,value)
	sum2=0.
	do 505 i=1,numt
	sum2=sum2+(ty(i)-a-b*tx(i))**2
505	continue
	sd=sqrt(sum2/numt)
	write(6,100) numt,a,sigmaa,b,sigmab,r,sd
	return
	end
c subroutine linfit.f
C MODIFIED!!!!
c
c source
c   Bevington, pages 104-105.
c
c purpose
c   make a least-squares fit to a data set with a straight line
c
c usage
c   call linfit (x, y, sigmay, npts, mode, a, sigmaa, b, sigmab, r,
C    iflag, value
C   if iflag=0 proceed normally
C   if iflag=1 fix b=value
C   if iflag=-1 fix a value
c
c description of parameters
c   x      - array of data points for independent variable
c   y      - array of data points for dependent variable
c   sigmay - array of standard deviations for y data points
c   npts   - number of pairs of data points
c   mode   - determines method of weighting least-squares fit
c            +1 (instrumental) weight(i) = 1./sigmay(i)**2
c             0 (no weighting) weight(i) = 1.
c            -1 (statistical)  weight(i) = 1./y(i)
c   a      - y intercept of fitted straight line
c   sigmaa - standard deviation of a
c   b      - slope of fitted straight line
c   sigmab - standard deviation of b
c   r      - linear correlation coefficient
c
c subroutines and function subprograms required
c   none
c
      subroutine linfit (x,y,sigmay,npts,mode,a,sigmaa,b,
     &  sigmab,r,iflag,value)
      double precision sum,sumx,sumy,sumx2,sumxy,sumy2
      double precision xi,yi,weight,delta,varnce
      dimension x(*),y(*),sigmay(*)
c
c accumulate weighted sums
c
11    sum=0.
      sumx=0.
      sumy=0.
      sumx2=0.
      sumxy=0.
      sumy2=0.
21    do 50 i=1,npts
      xi=x(i)
      yi=y(i)
      if (mode) 31,36,38
31    if (yi) 34,36,32
32    weight=1./yi
      goto 41
34    weight=1./(-yi)
      goto 41
36    weight=1.
      goto 41
38    weight=1./sigmay(i)**2
41    sum=sum+weight
      sumx=sumx+weight*xi
      sumy=sumy+weight*yi
      sumx2=sumx2+weight*xi*xi
      sumxy=sumxy+weight*xi*yi
      sumy2=sumy2+weight*yi*yi
50    continue
c
c calculate coefficients and standard deviations
c
51    delta=sum*sumx2-sumx*sumx
	if(iflag.eq.1) then
	a=(sumy-value*sumx)/sum
	b=value
	else if (iflag.eq.0) then
        a=(sumx2*sumy-sumx*sumxy)/delta
        b=(sumxy*sum-sumx*sumy)/delta
	else if (iflag.eq.-1) then
	a=value
	b=(sumy-a*sum)/sumx
	end if
61    if (mode) 62,64,62
62    varnce=1.
      goto 67
64    c=npts-2
      varnce=(sumy2+a*a*sum+b*b*sumx2
     *-2.*(a*sumy+b*sumxy-a*b*sumx))/c
67    sigmaa=dsqrt(varnce*sumx2/delta)
68    sigmab=dsqrt(varnce*sum/delta)
Cwrite(6,*) sigmab,varnce,sum,delta
	if(iflag.eq.1) sigmab=0.
71    r=(sum*sumxy-sumx*sumy)/
     *dsqrt(delta*(sum*sumy2-sumy*sumy))
      return
      end
	subroutine median(data,n,amed)
	parameter(maxstars=250000)
	dimension data(*),indx(maxstars)
	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

