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

	id=iargc()
        if(id.lt.1)then
        write(*,*)
        write(*,*)'       Usage: j_gc1 file [!]'
        write(*,*)
	write(*,*)'        input file: *.fit, *.ap'
        write(*,*)
        stop
        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 !'
	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,*,end=21)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+8.
	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.0.3)then          ! staturs & non_star
	  goto 71
	endif
	goto 70
71	a(i,1)=-3
70	continue	

c delete more star
	y=x+8.
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.20)then
          write(*,*)' bad fits_image.'
	  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)
	endif
90	continue

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

	if(id.gt.1)then
	do 99 i=1,k
99	write(*,'(i6,4f8.1)')a1(i),a2(i),a3(i),a12(i),a16(i)
	write(*,*)k
	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')
	nstar=10
	ratio=graw(a2048,a2,a3,a16,nstar,id)
	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+8.
	do 140 i=1,m
	if(a(i,4).gt.y)goto 141
	goto 140
141	a(i,1)=0.
140	continue

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

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

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

	dm=2.5*alog10(ratio)
	write(*,*)'dm ',dm
c sort star
	k=0
	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,6)-dm
	  a16(k)=a(i,16)
	endif
190	continue

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

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

        end

	subroutine sort(a,b,c,d,e,n)
	integer a(1)
	real b(1),c(1),d(1),e(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
	endif
10	continue
	end

	function graw(a2048,ax,ay,as,m,id)
	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.1)call disp(b,c,m)
	call phot(c,cs,m,ratio,id)
	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)
	call whitexblack(map,101,101,white,sigma)
	
	k=icheck(b(1,i),as(ii),sigma)
        if(k.gt.3)then
	  write(*,*)'delete No: ',ii,k
	  goto 31
	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)
	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(3)
	write(3,*)' while using 5 pixel, adjust % is:',ratio
	write(*,*)' while using 5 pixel, adjust % is:',ratio
	end
