	real a(2048*2048)
	real a101(101*101)
	real b(101*101,100),scale(100),white(100)
	real d(100),s(100)
	character*30 f1,f2,f3
	if(iargc().lt.1)then
	write(*,*)
	write(*,*)'       Usage: cob101 file.fit'
	write(*,*)
	stop
	endif
	call getarg(1,f1)
	k=index(f1,'.')
	if(k.eq.0)f1=f1(1:lnblnk(f1))//'.fit'
	if(k.ne.0)f1=f1(1:k)//'fit'
	inquire(file=f1,exist=i)
	if(i.eq.0)then
	  write(*,*)f1,'file not found !'
	  stop
	endif
	k=index(f1,'.')
	f2=f1(1:k)//'lst6'
	f3=f1(1:k)//'con'
	open(3,file=f3,status='unknown')
	inquire(file=f2,exist=i)
	if(i.eq.0)then
	  write(*,*)f2,'file not found !'
	  stop
	endif
	call readfits(f1,a,n1,n2)
	call readlst(f2,a,b,scale,white,m)
	  gain=3.3
	  s61=6.
	  s62=6.
	  if(m.gt.30)m=30
	  write(*,*)(scale(i),i=1,m)
	  write(*,*)(white(i),i=1,m)
	do 215 i=2,m
215	scale(i)=scale(i)/scale(1)
	scale(1)=1.
	do 216 i=1,m
	d(i)=sqrt(white(i)/gain)
216	s(i)=s62*d(i)/scale(i)
	sigma=xmedian(d,m)
	k=0
	do 240 i=1,101*101
	do 241 j=1,m
241	d(j)=b(i,j)
	k=k+1
240	a101(k)=allc(d,white,scale,m,sigma,s,s61)		
c	call disp(a101,m)
	call phot(a101,m,ratio)
	call phot6(f2,a,white,m,ratio)
	close(3)
	end

	function allc(d,white,scale,m,sigma,s,s61)
	real d(1),white(1),scale(1),s(1)
	real d1(100)			! normal value
        real d2(100)                     ! bad working unit
        bad=-999.                       ! s(i)=6*sigma(i)/scale(i)
c bad pixel
        allc=bad
	do 5 i=1,m
5	if(d(i).gt.29000.or.d(i).lt.0)d(i)=bad
        do 20 i=1,m
20      if(d(i).gt.bad)d1(i)=(d(i)-white(i))/scale(i)  !ne
        n=0
        do 30 i=1,m
        if(d(i).le.bad)goto 30                    ! eq
        n=n+1
        d2(n)=d1(i)
30      continue
        x=xmedian(d2,n)
        do 40 i=1,m
40      if(d(i).le.bad)d1(i)=x                    ! eq
c cosmos ray
        do 110 i=1,m
110     d2(i)=d1(i)
        x=xmedian(d2,m)
c	  write(*,'(8f6.1)')(d2(i),i=1,m),x,sigma
c        if(x.gt.sigma)then              ! star
          if(d2(1).lt.s61*x)goto 200
          do 120 i=1,m
120       if(d1(i).gt.s61*x1)d1(i)=x
c        else                            ! background
c          do 130 i=1,m
c130       if( d1(i).eq.d2(1) .and. d1(i).gt.s(i) )d1(i)=x
c130       if(d1(i).gt.s(i) )d1(i)=x
c        endif
c combine
200     x=0.
        do 210 i=1,m
        y= d1(i)*scale(i)+white(i)
210     x=x+y
220     allc=x
				! s(i)=6*sigma(i)/scale(i)
	end	

        function xmedian(a,n)
        real a(n)
        do 10 i=1,n-1
        do 10 j=i+1,n
        if(a(i).gt.a(j))goto 10
        x=a(i)   
        a(i)=a(j)
        a(j)=x  
10      continue
        m=n/2
        xmedian=a(m+1)
        if(m+m.ne.n)return
	xmedian=(a(m)+a(m+1))/2
        end

	subroutine phot6(f1,a,w,m,ratio)
	character*30 f1
	real a(2048,2048)
	real u(100)
	real w(100)
	open(1,file=f1,status='old')
	k=0
	do 20 i=1,m
21	read(1,*)j,x2,y2,z2,z1
 	if(x2.lt.51.or.x2.gt.1997.or.y2.lt.51.or.y2.gt.1997)goto 21
	call phot6_5(a,w(i),x2,y2,adu)
	xmag=25.-2.5*log10(adu*ratio)
	write(*,'(i6,3f8.1,2f10.3)')j,x2,y2,z2,xmag,z1
	write(3,'(i6,3f8.1,2f10.3)')j,x2,y2,z2,xmag,z1
	k=k+1
	u(k)=xmag-z1
20	continue
	close(1)
        zz=sigma(u,k)
        write(*,*)'total, sigma: ',k,zz
        write(3,*)'total, sigma: ',k,zz

	end

        function sigma(b,k)
        real b(k)
        x=0
        do 20 j=1,k
20      x=x+b(j)
        x=x/k
        write(*,*)' average: ',x
        write(3,*)' average: ',x
        xx=0
        do 30 j=1,k
30      xx=xx+(b(j)-x)**2
        sigma=sqrt(xx/k)
        end


	subroutine readlst(f1,a,b,s,w,m)
	character*30 f1
	real a(2048,2048)
	real b(101*101,100)
	real s(100),w(100)
	integer*2 map(101*101)
	open(1,file=f1,status='old')
	kk=0
	do 20 i=1,100
21	read(1,*,end=201)j,x2,y2
	kk=kk+1
 	if(x2.lt.51.or.x2.gt.1997)goto 21
 	if(y2.lt.51.or.y2.gt.1997)goto 21
	i1=nint(x2)
	i2=nint(y2)
	k=0
	k1=0
	do 203 j1=i1-50,i1+50
	do 203 j2=i2-50,i2+50
	k=k+1
	b(k,i)=a(j1,j2)
	mmm=(j1-i1)**2+(j2-i2)**2
c        if(mmm.lt.900.and.mmm.gt.400)then
        if(mmm.lt.625)then
	  k1=k1+1
	  map(k1)=a(j1,j2)
	endif
203	continue
	call whitexblack(map,1,k1,w(i),dummy)
	s(i)=0
	do 204 j1=i1-5,i1+5
	do 204 j2=i2-5,i2+5
204	s(i)=s(i)+a(j1,j2)-w(i)
20	continue
201	m=i-1
	write(*,*)'number: ',m,kk
	write(3,*)'number: ',m,kk
	close(1)
	end

	subroutine disp(a,m)
	real a(101*101)
	integer*2 map(101*101)
	do 10 i=1,101*101
10	map(i)=a(i)/m
	call whitexblack(map,101,101,white,sigma)
	call pgbegin(0,'/xw',2,2)
        call pgenv(1.,101.,1.,101.,1,-1)
	black=white+20.*sigma
	call pgxgray(map,101,101,1,101,1,101,white,black,irat)
	call pgend
	end

	subroutine phot6_5(a,w,x2,y2,adu)
	real a(2048,2048)
	xc=x2
	yc=y2
	lx=xc-9
	mx=xc+9
	ly=yc-9
	my=yc+9
	par=5.
	adu=0.
      DO 2130 J=LY,MY
      DYSQ=(J-YC)**2
         DO 2125 I=LX,MX
            RSQ=DYSQ+(I-XC)**2
            DATUM=a(I,J)-w
            R=SQRT(RSQ)-0.5
               IF (R .GT. par ) GO TO 2120
               FRACTN=AMAX1(0.0, AMIN1(1.0,par-R))
               adu = adu+FRACTN*DATUM
 2120       CONTINUE
 2125    CONTINUE
 2130 CONTINUE
	end

	subroutine phot(a,m,ratio)
	parameter (naper=18)
	integer*2 map(101*101)
	real a(101,101)
	real apmag(naper)
	real ap(naper)
	integer par(naper)

	k=0
	do 10 i=51-20,51+20
	do 10 j=51-20,51+20
        mmm=(51-i)**2+(51-j)**2                !
c	if(mmm.lt.900.and.mmm.gt.400)then
	if(mmm.lt.625)then
 	  K=k+1
	  map(k)=amin0(a(i,j)/m,32700.)+0.5
	endif
10	continue
	call whitexblack(map,1,k,white,sigma)
	       write(*,*)white
	xc=51
	yc=51
	n1=naper+3
	lx=xc-n1
	mx=xc+n1
	ly=yc-n1
	my=yc+n1

	DO 2010 I=1,NAPER
	par(i)=i+2
2010    APMAG(I) = 0
C
      DO 2130 J=LY,MY
      DYSQ=(J-YC)**2
C
         DO 2125 I=LX,MX
            RSQ=DYSQ+(I-XC)**2
            DATUM=a(I,J)-white*m
C
C The inclusion of partial pixels inside the aperture is done as
C follows:  if the distance of the center of the current pixel from the
C centroid of the star [radius vector r(i,j)] is exactly equal to the
C radius of the aperture [R(k)], then one-half of the counts in the
C pixel are included.  If r(i,j) < R(k)-0.5, then the entire pixel is
C included, while if r(i,j) > R(k)+0.5, the pixel is wholly excluded.
C In between, viz. for  R(k)-0.5 < r(i,j) < R(k)+0.5, the fraction of
C the counts included varies linearly.  Therefore a circular aperture
C is approximated by an irregular (not even convex) polygon.
C
C If this pixel falls completely outside the LARGEST aperture, go on
C to the next pixel.  Notice that APMXSQ has actually been defined
C as (R(k)+0.5)**2 for the largest value of R(k), in accordance with
C the formula used for the partial pixels.
C
            R=SQRT(RSQ)-0.5
C
            DO 2120 K=1,NAPER
               IF (R .GT. PAR(K)) GO TO 2120
               FRACTN=AMAX1(0.0, AMIN1(1.0,PAR(K)-R))
               APMAG(K) = APMAG(K)+FRACTN*DATUM
 2120       CONTINUE
 2125    CONTINUE
C
 2130 CONTINUE


      DO 2220 I=1,NAPER
	ap(i)=apmag(i)
      APMAG(I)=25.-2.5*LOG10(APMAG(I))
 2220 CONTINUE

	do 99 i=1,naper
	write(*,'(i2,f12.2,f12.5)')par(i),ap(i),apmag(i)
99	write(3,'(i2,f12.2,f12.5)')par(i),ap(i),apmag(i)
	do 101 i=1,naper
	if(apmag(i+1).gt.apmag(i))goto 102
101	continue
102	if(i.gt.13)i=13
	  i=13
	write(*,*)'100% aperture_radius is:',i+2,' pixel'
	write(3,*)'100% aperture_radius is:',i+2,' pixel'
	ratio=ap(i)/ap(3)
	write(*,*)' while using 5 pixel, adjust % is:',ratio
	write(3,*)' while using 5 pixel, adjust % is:',ratio
	end
