parameter (n=90000)	   ! n= *.ap star number
	logical logi

	character*60 f1,f2,f3
	character*80 head(36)

	real a_b(4096*4032)
   	real aa(n,17),a(n,17),b(n,4)
	real result(200,6)
	real seeing(200)
	integer*2 dummy(n)

        if(iargc().lt.1)then
        write(*,*)
        write(*,*)
        write(*,*)'     **********   DAO_AUTOPSF  3*3 (2004.5)*********'
        write(*,*)
        write(*,*)
        write(*,*)'     Usage:  j_auto file'
        write(*,*)
        write(*,*)
        stop
        endif
        call getarg(1,f1)
        k=index(f1,'.')
        if(k.eq.0)f1=f1(1:lnblnk(f1))//'.fit'
        inquire(file=f1,exist=logi)
        if(.not.logi)stop '.fit file not found !'
	call readfits(f1,a_b,n1,n2)
c	write(*,*)n1,n2
 
        open(1,file=f1,status='old',access='direct',recl=2880)
        read(1,rec=1)head
        ipos=indexpos(head,'NAXIS1  ')
        read(head(ipos)(65:74),'(2i5)')ix1,ix2
        ipos=indexpos(head,'NAXIS2  ')
        read(head(ipos)(65:74),'(2i5)')iy1,iy2
        close(1)
	if(ix2.eq.0)then
	  ix1=1
	  ix2=n1
	  iy1=1
	  iy2=n2
	endif

        k=index(f1,'.')
        f2=f1(1:k)//'ap'
        inquire(file=f2,exist=logi)
        if(.not.logi)stop '.ap file not found !'
        f3=f1(1:k)//'lst'

	open(1,file=f2,status='old')
        do 19 i=1,3
19      read(1,'(a)')j
        do 29 i=1,n
        read(1,*,end=28)(aa(i,j),j=1,15)
        read(1,*,end=28)aa(i,16),aa(i,17)
29	continue
28      mn=i-1
	close(1)

c get mag_200
	im=0
	do 910 i=1,mn
	do 920 j=4,9			
	if(aa(i,j).gt.90.)goto 910
920	continue
	im=im+1
	dummy(im)=aa(i,4)*100.
910	continue
	call sort1(dummy,im)
	xm=dummy(im)
	if(im.gt.50)xm=dummy(50)
	xm=xm*0.01
	write(*,*)'xm ',xm

	im=0
	i15=9		! 4-11
        i30=25          ! 19-30  not used
	x05=0.
	r2=225

1001	m=0
        do 27 i=1,mn
	m=m+1
	do 26 j=1,17
26	a(m,j)=aa(i,j)
27	continue
	write(*,*)'all_star: ',m

	do 5 i=1,m
	do 5 j=1,4
5	b(i,j)=a(i,j)
	mall=m

	l=0
	do 10 i=1,m
c delete >95 point
	do 20 j=4,i15			! 15
	if(a(i,j).gt.90)goto 10
20	continue
c delete edge point 20 pixel
	if(a(i,2).lt.20. .or. a(i,2).gt.4044)goto 10
	if(a(i,3).lt.20. .or. a(i,3).gt.4044)goto 10

	l=l+1	
	do 90 j=1,17
90	a(l,j)=a(i,j)
10	continue
	m=l

	write(*,*)'af 4,11>90; >20 edge_star  ',m

c get sky,s_sky
	xs=0
	do 110 i=1,m
110	xs=xs+a(i,16)
	xm2=xm+x05
	xs=xs/m
	write(*,*)'mag,sky: ',xm2,xs

	l=0
	xs2=xs*1.1

	do 210 i=1,m
	if(a(i,4).gt.xm2)goto 210
	if(a(i,16).gt.xs2)goto 210
	l=l+1
	do 290 j=1,17
290	a(l,j)=a(i,j)
210	continue
	m=l
 	write(*,*)'af xm2, sky>1.1: ',m
c      if(i.eq.i)goto 10090
c delete galacy
	l=0
	do 310 i=1,m
	a4=a(i,4)
	af=a4
	do 320 j=5,i15
	if(a(i,j)-af .gt. 0.1)goto 310
	if(a4-a(i,j) .gt. 1.0)goto 310
	af=a(i,j)
320	continue
	l=l+1
	do 390 j=1,17
390	a(l,j)=a(i,j)
310	continue
	m=l
 	write(*,*)'af galaxy: ',m
c save data for if loop
	l=30000
	do 950 i=1,m
	do 950 j=1,17
950	a(i+l,j)=a(i,j)
	mm=m

501	write(*,*)'r2=',r2
c501	continue
c delete neibur star of remain   
	l=5000
	do 410 i=1,m
	x=a(i,2)
	y=a(i,3)
 	do 420 j=1,m
	if(i.eq.j)goto 420
	z=(a(j,2)-x)**2+(a(j,3)-y)**2
	if(z.lt.r2)goto 410
420	continue
	l=l+1
	do 490 j=1,17
490	a(l,j)=a(i,j)
410	continue
	l=l-5000
	do 491 i=1,l
	do 491 j=1,17
491	a(i,j)=a(i+5000,j)
	m=l
 	write(*,*)'af_neibur: ',m
c delete neibur stars
	l=0
c	do 505 i=1,m
c	k=nint(a(i,1))
c	b(k,2)=0.
c505	b(k,3)=0.
	do 510 i=1,m
	x=a(i,2)
	y=a(i,3)
	do 520 j=1,mall
	z=(b(j,2)-x)**2+(b(j,3)-y)**2
	if(z.lt.0.0001)goto 520
 	if(z.lt.r2)goto 510
520	continue
	l=l+1
	do 590 j=1,17
590	a(l,j)=a(i,j)
510	continue
	m=l
	write(*,*)'af_neibur2: ',m

997	do 600 i=1,m
	y=see(a_b,n1,a(i,2),a(i,3),a(i,16))
	seeing(i)=y
600	a(i,5)=y
	x=xmedian(seeing,m)
	seex=x
	seeh=seex*0.75
	l=0
	do 610 i=1,m
 	if(a(i,5).gt.seex)goto 610
	if(a(i,5).lt.seeh)goto 610
	l=l+1
	do 690 j=1,17
690	a(l,j)=a(i,j)
610	continue
	m=l
	write(*,*)'af seeing: ',m

	x=0.
	do 6601 l=1,m
6601	x=x+a(l,16)
	x16=x/m
	x=0.
	do 661 l=1,m
	ix=nint(a(l,2))
	iy=nint(a(l,3))
	z=0
        i70=70
	do 662 i=ix-i70,ix+i70
	if(i.lt.1)goto 662	
	if(i.gt.n1)goto 662	
	do 663 j=iy-i70,iy+i70
	if(j.lt.1)goto 663	
	if(j.gt.n2)goto 663	
	y=(i-ix)**2+(j-iy)**2
	y=sqrt(y)
	if(y.lt.15)goto 663
	k=(j-1)*n1+i
	y=a_b(k)-x16
	if(y.gt.0)z=z+1
663	continue
662	continue
	x=x+z
661	a(l,7)=z
	x=x/m*1.2	
	l=0
	do 664 i=1,m
	if(a(i,7).gt.x)goto 664
	l=l+1
	do 665 j=1,17
665	a(l,j)=a(i,j)
664	continue
	m=l
	write(*,*)'af bright_star near: ',m

c delete bad wing
998	do 700 i=1,m
	y=wing(a_b,n1,n2,a(i,2),a(i,3),seex,r2)
700	a(i,6)=y*(a(i,4)-xm+10.)*0.01
	x12=1.1
704	m3=0
705	x=0
	m3=m3+1
	do 701 i=1,m
701	x=x+a(i,6)
	x=x/m*x12
	l=0
	do 710 i=1,m
	if(a(i,6).gt.x)goto 710
	l=l+1
	do 790 j=1,17
790	a(l,j)=a(i,j)
710	continue
	m=l
	write(*,*)'af bad_wing: ',m
	if(x05.lt.0.9 .and. m.le.20)then
	  x05=x05+0.5
	  i15=i15-1
	  r2=100.
	  goto 1001
	endif
10090	continue
	do 800 i=1,m
	im=im+1
	do 801 j=1,4
801	result(im,j)=a(i,j)
	result(im,5)=a(i,16)
800	result(im,6)=a(i,5)

999	call sort6(result,im,4)

        open(3,file=f3,status='UNKNOWN')
        if(im.gt.30)im=30
	do 900 i=1,im
	write(3,2)nint(result(i,1)),(result(i,j),j=2,6)
900	write(*,2)nint(result(i,1)),(result(i,j),j=2,6)
2	format(i6,7f8.1)
	close(3)
	write(*,*)im
	end

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

	function wing(a,n1,n2,x,y,see,r2)
	real a(n1,n2)
	real ux(961),uy(961)
	j=0
	z=0
	k=20
	i1=nint(x)
	i2=nint(y)
	do 10 j1=i1-15,i1+15
	do 10 j2=i2-15,i2+15
	j=j+1
	uy(j)=a(j1,j2)
10      ux(j)=sqrt((j1-x)**2+(j2-y)**2)
	call sort(ux,uy,961)
	do 20 i=1,961
	if(ux(i).gt.see)goto 30
20	continue
30	do 40 j=i,i+r2
	z=z+abs(uy(j)-uy(j+1))/ux(j)
	k=k-1
	if(k.eq.0)k=1
40	continue
	wing=z
	end

	subroutine sort1(a,n)
	integer*2 a(1),b
	do 10 i=1,n-1
	do 10 j=i+1,n
	if(a(i).gt.a(j))then
	  b=a(i)
	  a(i)=a(j)
	  a(j)=b
	endif
10	continue
	end
	

      subroutine sort6(a,n,m)
      real a(200,6)
	do 10 i=1,n-1
	do 10 j=i+1,n
	if(a(i,m).lt.a(j,m))goto 10
        do 20 k=1,6
        x=a(i,k)
	a(i,k)=a(j,k)
20	a(j,k)=x	
10	continue
      end

      subroutine sort(ia,ib,n)
      real ia(n),ib(n),it,iu
      int=2
   10 int=2*int
      if(int.lt.n)goto 10
      int=min0(n,(3*int)/4-1)
   20 int=int/2
      ifin=n-int
      do 70 ii=1,ifin
      i=ii
      j=i+int
      if(ia(i).le.ia(j))goto 70
      it=ia(j)
      iu=ib(j)
   40 ia(j)=ia(i)
      ib(j)=ib(i)
      j=i
      i=i-int
      if(i.le.0)goto 60
      if(ia(i).gt.it)goto 40
   60 ia(j)=it
      ib(j)=iu
   70 continue
      if(int.gt.1)goto 20
      end

        function see(map,n1,x,y,white)
c (a81+a84)/2 *180*3600/pi = arcsec / pixel
c schmidt is 1.709"/pixel
c bok is 0.45"/pixel
c seeing= 2.355*sigma*1.709
        real map(n1,n1)
        integer b(113)        ! (29-1)*4+1
        ix=x+0.5
        iy=y+0.5
        j=1
        do 10 i=ix-14,ix+14
        b(j)=map(i,iy)-white
10      j=j+4
        do 30 j=1,112,4
30      b(j+2)=(b(j)+b(j+4))*0.5+0.5
        do 40 j=1,112,2
40      b(j+1)=(b(j)+b(j+2))*0.5+0.5
        call histat(b,mode,maxh,xmean,xpeak,fsigma1,gsigma,113)
        j=1
        do 110 i=iy-14,iy+14
        b(j)=map(ix,i)-white
110     j=j+4
        do 130 j=1,112,4
130     b(j+2)=(b(j)+b(j+4))*0.5+0.5
        do 140 j=1,112,2
140     b(j+1)=(b(j)+b(j+2))*0.5+0.5
        call histat(b,mode,maxh,xmean,xpeak,fsigma2,gsigma,113)
c       write(*,'(a,f9.2,a)')'seeing:',(fsigma1+fsigma2)/16.*2.355*1.709,'"'
        see=(fsigma1+fsigma2)/8.*2.355*0.45        ! for BOK
        end
 
        function indexpos(head,f1)
        character*80 head(36),f1*8
        do 10 indexpos=1,36
10      if(head(indexpos)(1:8).eq.f1)return
        end