	parameter (n=40000)	   ! n= *.ap star number
	parameter (ed=31)	   ! edge
	character*60 f1,f2
   	real a(n,16)
	integer a1(n)
	real a2(n),a3(n),a12(n),a16(n),aa(n)
	real a2048(2048*2048)
	logical logi

	id=iargc()
        if(id.lt.2)then
        write(*,*)
        write(*,*)'       Usage: j_gc2 file radius [!]'
        write(*,*)
	write(*,*)'        input file: *.fit, *.ap'
	write(*,*)'        output file: *.gc, *.lst'
        write(*,*)
        stop
        endif

	iradius=5
	if(id.ge.2)then
	call getarg(2,f1)
	read(f1(1:),*) iradius
	if(iradius.lt.4)iradius=4
	if(iradius.ge.9)iradius=10
	endif

        call getarg(1,f1)
        k=index(f1,'.')-1
	if(k.lt.0)k=lnblnk(f1)
        f2=f1(1:k)//'.ap'
        inquire(file=f2,exist=logi)
        if(.not.logi)stop '.ap file not found !'

	x03=0.3
102	open(1,file=f2,status='old')
        do 10 i=1,3
10      read(1,'(a)')j
        do 20 i=1,n
        read(1,*,end=21)(a(i,j),j=1,15)
        read(1,'(f13.3)',err=20)a(i,16)
20	continue
21      m=i-1
	close(1)

c       find_out brightest star's magnitude
	x=99.
	do 30 i=1,m
30	if(a(i,4).lt.x)x=a(i,4)

c delete faint star
	y=x+5.
	do 40 i=1,m
	if(a(i,4).gt.y)goto 41
	goto 40
41	a(i,1)=0.
40	continue

c delete nearby_star, 20_pixel
	do 50 i=1,m-1
	if(a(i,1).eq.0)goto 50
	do 51 j=i+1,m
	if(a(j,1).eq.0)goto 51
	r=(a(i,2)-a(j,2))**2+(a(i,3)-a(j,3))**2
	if(r.lt.400)then
	  a(i,1)=-1
	  a(j,1)=-1
	endif
51	continue
50	continue

c delete edge star
        do 60 i=1,m
	if(a(i,1).le.0)goto 60
        if(a(i,2).lt.ed-1.or.a(i,2).gt.2048-ed)goto 61
        if(a(i,3).lt.ed-1.or.a(i,3).gt.2048-ed)goto 61
        goto 60
61      a(i,1)=-2.
60      continue


c delete galaxic & staturs star & non star
	do 70 i=1,m
	if(a(i,1).le.0)goto 70
	do 72 j=5,13                          ! between 3--16 
	if(a(i,j).gt.a(i,j-1)+0.005)goto 71
72	continue
	if(a(i,5)-a(i,8).gt.x03)then          ! staturs & non_star
c  50415 if seeing bad, change 0.3 ---> 0.4
	  goto 71
	endif
	goto 70
71	a(i,1)=-3
70	continue	

c delete more star
	y=x+5.
82	k=0
	do 80 i=1,m
	if(a(i,1).le.0)goto 80
	if(a(i,4).gt.y)goto 81
	k=k+1
	goto 80
81	a(i,1)=0
80	continue	
	if(k.gt.100)then
	  y=y-0.5
	  goto 82
	endif

	if(k.lt.10 .and. x03.eq.0.3)then
	  x03=0.4
	  goto 102
	endif
	if(k.lt.8)then
          write(*,*)' bad fits_image. left star ',k
	  stop
	endif

c sort star
	k=0
	do 90 i=1,m
 	if (a(i,1).gt.0.)then
          k=k+1
	  a1(k)=a(i,1)+0.1
	  a2(k)=a(i,2)
	  a3(k)=a(i,3)
	  a12(k)=a(i,12)
	  a16(k)=a(i,16)
	  aa(k)=a(i,6)-(a(i,13)+a(i,12))*.5
	endif
90	continue

	call sort(a1,a2,a3,a12,a16,k,aa)         ! no, x,y, mag_(pixel 14), sky

	if(id.gt.2)then
	do 99 i=1,k
99	write(*,'(i6,4f8.1)')a1(i),a2(i),a3(i),a12(i),a16(i)
	write(*,*)k
	endif

	write(*,'(10f8.3)')(a12(i),i=1,k)
	write(*,'(10f8.3)')(aa(i),i=1,k)

	if(k.gt.20)k=20
	nstar=ndel3s(a2,a3,a16,aa,k)
	write(*,*)'nstar: ',nstar
	
	if(nstar.lt.10)then
          write(*,*)' bad fits_image.'
	  stop
	endif

        k=index(f1,'.')-1
	if(k.lt.0)k=lnblnk(f1)
        f2=f1(1:k)//'.fit'
        inquire(file=f2,exist=logi)
        if(.not.logi)stop '.fit file not found !'

	call readfits(f2,a2048,n1,n2)


	f2=f1(1:k)//'.gc'
	open(3,file=f2,status='unknown')
	ratio=graw(a2048,a2,a3,a16,nstar,id,iradius)
	close(3)
c**************************************	
	f2=f1(1:k)//'.lst'
	open(3,file=f2,status='unknown')

	do 100 i=1,m
100	a(i,1)=1
c delete faint star
	y=x+7.
	do 140 i=1,m
	if(a(i,iradius+1).gt.y)goto 141
	goto 140
141	a(i,1)=0.
140	continue

c delete nearby_star, 10_pixel
c	do 150 i=1,m-1
c	if(a(i,1).eq.0)goto 150
c	do 151 j=i+1,m
c	if(a(j,1).eq.0)goto 151
c	r=(a(i,2)-a(j,2))**2+(a(i,3)-a(j,3))**2
c	if(r.lt.100)then
c	  a(i,1)=-1
c	  a(j,1)=-1
c	endif
c151	continue
c150	continue

c delete galaxic & staturs star & non star
c	iss=iradius+2
c	if(iss.gt.10)iss=10
c	do 170 i=1,m
c	if(a(i,1).le.0)goto 170
c	do 172 j=5,iss                          ! between 3--8 
c	if(a(i,j).gt.a(i,j-1)+0.01)goto 171
c172	continue
c	if(a(i,5)-a(i,8).gt.0.5)then          ! staturs & non_star
c	  goto 171
c	endif
c	goto 170
c171	a(i,1)=-3
c170	continue	

c delete more star
c	y=x+8.
c182	k=0
c	do 180 i=1,m
c	if(a(i,1).le.0)goto 180
c	if(a(i,4).gt.y)goto 181
c	k=k+1
c	goto 180
c181	a(i,1)=0
c180	continue	
c	if(k.gt.500)then
c	  y=y-0.5
c	  goto 182
c	endif

	dm=2.5*alog10(ratio)
	write(*,*)'dm ',dm
c sort star
	k=0
	kk=iradius
	if(kk.lt.10)kk=kk+1
	do 190 i=1,m
 	if (a(i,1).gt.0.)then
          k=k+1
	  a1(k)=a(i,1)+0.1
	  a2(k)=a(i,2)
	  a3(k)=a(i,3)
	  a12(k)=a(i,kk)-dm                    ! radius=6  (i,7),  =8~9, =10~10
	  a16(k)=a(i,16)
	endif
190	continue

	call sort(a1,a2,a3,a12,a16,k,aa)         ! no, x,y, mag_(pixel 5), sky

	do 199 i=1,k
          if(a16(i).gt.9999.)a16(i)=9999.
199	write(3,'(i6,2f8.1,2f8.3)')i,a2(i),a3(i),a12(i),a16(i)
	write(*,*)'nstar: ',k
	close(3)

        end

	function ndel3s(a2,a3,a16,aa,m)
	real a2(1),a3(1),a16(1),aa(1)
	real b2(90),b3(90),b16(90),ba(90)
	nstar=m
11	k=0
	x=0
	do 10 i=1,nstar
10	x=x+aa(i)
	x=x/nstar
	y=0
	do 20 i=1,nstar
20	y=y+(aa(i)-x)**2
	y=sqrt(y/nstar)*3
	if(nstar.ne.m.and.y.gt.0.04)y=0.04
	do 30 i=1,nstar
	if(abs(aa(i)-x).gt.y)then
	  a2(i)=0
	  k=k+1
	endif
30	continue
	if(k.ne.0)then
	  k=0
	  do 40 i=1,nstar
	  if(a2(i).ne.0)then
	    k=k+1
	    b2(k)=a2(i)
	    b3(k)=a3(i)
	    b16(k)=a16(i)
	    ba(k)=aa(i)
	  endif
40	  continue
	  nstar=k
	  do 50 i=1,nstar
	  a2(i)=b2(i)
	  a3(i)=b3(i)
	  a16(i)=b16(i)
50	  aa(i)=ba(i)
	  goto 11
	endif
	ndel3s=nstar
	end

	subroutine sort(a,b,c,d,e,n,f)
	integer a(1)
	real b(1),c(1),d(1),e(1),f(1)

	do 10 i=1,n-1
	do 10 j=i+1,n
	if(d(i).gt.d(j))then
	  k=a(i)
	  a(i)=a(j)
	  a(j)=k
          x=b(i)
	  b(i)=b(j)
	  b(j)=x
          x=c(i)
	  c(i)=c(j)
	  c(j)=x
          x=d(i)
	  d(i)=d(j)
	  d(j)=x
          x=e(i)
	  e(i)=e(j)
	  e(j)=x
          x=f(i)
	  f(i)=f(j)
	  f(j)=x
	endif
10	continue
	end

	function graw(a2048,ax,ay,as,m,id,iradius)
	real ax(1),ay(1),as(1)
	real a2048(2048*2048)
	real b(101*101,20)
	real c(101*101)
	call readlst(a2048,b,ax,ay,as,m)
	call addall(b,c,as,cs,m)	
	if(id.gt.2)call disp(b,c,m)
	call phot(c,cs,m,ratio,id,iradius)
	graw=ratio
	end
	
	subroutine addall(b,c,as,cs,m)
	real b(101*101,20),c(101*101)
	real as(1)
	do 20 j=1,101*101
	x=0
	do 10 i=1,m
10	x=x+b(j,i)
20	c(j)=x
	cs=0
	do 30 i=1,m
30	cs=cs+as(i)
	end

	subroutine readlst(a2048,b,ax,ay,as,m)
	real a2048(2048,2048)
	real b(101*101,20)
	real ax(1),ay(1),as(1)
	integer*2 map(101*101)
	ii=0
	do 30 i=1,m
31	ii=ii+1
	i1=nint(ax(ii))
	i2=nint(ay(ii))
	k=0
	do 20 j1=i1-50,i1+50
	do 20 j2=i2-50,i2+50
	k=k+1
	b(k,i)=a2048(j1,j2)
20      map(k)=a2048(j1,j2)
c	call whitexblack(map,101,101,white,sigma)
	
c	k=icheck(b(1,i),as(ii),sigma)
c       if(k.gt.1)then
c	  write(*,*)'delete No: ',ii,k
c	  goto 31
c	endif
30	continue
	end

	function icheck(a,sky,sigma)
	real a(-50:50,-50:50)
	k=0
	do 10 i=-15,15
	do 10 j=-15,15
	ir=i*i+j*j
	if(ir.lt.7*7)goto 10
	if(ir.gt.15*15)goto 10
	if(a(i,j).gt.sky+1000.*sigma/ir)k=k+1
10	continue
	icheck=k
	end
	
	subroutine disp(b,c,m)
	real b(101*101,20),c(101*101)
	integer*2 map(101*101)
	call pgbegin(0,'?',4,4)
	do 100 j=1,m
	do 10 i=1,101*101
10	map(i)=b(i,j)
        call pgenv(1.,101.,1.,101.,1,-1)
	call whitexblack(map,101,101,white,sigma)
	black=white+20.*sigma
100	call pgxgray(map,101,101,1,101,1,101,white,black,1)
	do 20 i=1,101*101
20	map(i)=c(i)/m
        call pgenv(1.,101.,1.,101.,1,-1)
        call pgenv(1.,101.,1.,101.,1,-1)
	call whitexblack(map,101,101,white,sigma)
	black=white+20.*sigma
	call pgxgray(map,101,101,1,101,1,101,white,black,1)
	call pgend
	end

	subroutine phot(a,cs,m,ratio,id,iradius)
	parameter (naper=18)
	real a(101,101)
	real apmag(naper)
	real ap(naper)
	integer par(naper)
	white=cs
	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
      DO 2130 J=LY,MY
      DYSQ=(J-YC)**2
         DO 2125 I=LX,MX
            RSQ=DYSQ+(I-XC)**2
            DATUM=a(I,J)-white
            R=SQRT(RSQ)-0.5
            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
 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(3,'(i2,f12.2,f12.5)')par(i),ap(i),apmag(i)
99	write(*,'(i2,f12.2,f12.5)')par(i),ap(i),apmag(i)
	write(3,*)'100% aperture_radius is 15 pixel'
	write(*,*)'100% aperture_radius is 15 pixel'
	ratio=ap(13)/ap(iradius-2)
	write(3,*)' while using',iradius,' pixel, adjust % is:',ratio
	write(*,*)' while using',iradius,' pixel, adjust % is:',ratio
	end
