	character*80 f1,head(108)
        real*8 a8(8),adcoef(2,3)
        real a(2048,2048)
	integer mag(25)

        pi=4.*atan(1.)
        cx=pi/12.
        cy=pi/180.

	if(iargc().lt.2)then
	  write(*,*)
	  write(*,*)' ********* sub_gscno in_fits out_fits scale shift'
	  write(*,*)'                                   2000.7'
	  write(*,*)' 1. 1st run gscno_z, to  produce gscno_z.tmp'
	  write(*,*)' 2. 2nd edit sub_gscno.dat (mag,radius)'
	  write(*,*)' 3. scale ~1.,   0.5~2.0'
	  write(*,*)' 4. shift ~0,,   -2.~2.0'  
          write(*,*)
          write(*,*)' RADIUS(mag) = TABLE(gscno_mag)*scale+shift'
          write(*,*)
	  stop
	endif

	scale=1.
	shift=0.
	if(iargc().ge.3)then
	  call getarg(3,f1)
	  read(f1(1:),*)scale
	endif
	if(iargc().ge.4)then
	  call getarg(4,f1)
	  read(f1(1:),*)shift
	endif
	open(1,file='sub_gscno.dat',status='old')
5	read(1,*,end=6)i,x
	mag(i)=x*scale+0.5+shift
	goto 5
6	close(1)
	write(*,*)mag

	call getarg(1,f1)         
	k=lnblnk(f1)+1
	f1(k:k)=char(0)
	write(*,'(2a)')'read: ',f1

	call gethead(f1,head)
        call geta8(head,a8)
	call getsky(head,sky)
        call xytoad(a8,adcoef)

	call readfits(f1,a,n1,n2)

	open(1,file='gscno_z.tmp',status='old')
15	read(1,'(a)',end=20)f1
	if(f1(2:2).eq.' ')goto 15
        i= itoh(f1(2:),alpha,delta)
          alp=alpha*cx
          del=delta*cy
	a87=a8(7)
	a88=a8(8)
          call standc(a87,a88,alp,del,xi,xn)
          alp=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3)
          del=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)
          ia=2049.-alp+.5
          id=2049.-del+.5
c        write(*,*)f1(2:30),ia,id
        if(ia.gt.2045.or.ia.lt.4)goto 15
        if(id.gt.2045.or.id.lt.4)goto 15
	read(f1(24:),*)x
	i5=x+0.5
	i5=mag(i5)
	if(i5.eq.0)goto 15
	kk=i5*i5
	do 25 i=-i5,i5
	do 25 j=-i5,i5
	k=i*i+j*j
	if(ia+i.gt.2048.or.ia+i.lt.1)goto 25
	if(id+j.gt.2048.or.id+j.lt.1)goto 25
	if(k.le.kk)a(ia+i,id+j)=sky
25	continue
	goto 15
20	close(1)
	nhead=indexpos(head,'END     ')
        if(nhead.lt.72)nhead=72
        if(nhead.gt.72)then
          nhead=72
	  do 10 i=1,80
10	  head(72)(i:i)=char(32)
	  head(72)(1:3)="END"
	endif
	call getarg(2,f1)         
	k=lnblnk(f1)+1
	f1(k:k)=char(0)
	write(*,'(2a)')'out: ',f1
	call swap4(a,n1*n2*4)
	call writefits(f1,head,80*nhead,1)
	call writefits(f1,a,n1*n2*4,0)
c do shift end
	end

        function itoh(cc,alpha,delta)
        character*(*) cc
        character*20 dd
        byte bell
        bell=7
        itoh=1
        k1=index(cc,' ')
        k2=index(cc,',')
        k=min(k1,k2)
        if(k1.eq.0)k=k2
        if(k2.eq.0)k=k1
        if(k.lt.2.or.k.gt.13)goto 100
        dd=cc(1:k)
        do 10 i=1,9
10      if(dd(i:i).eq.':')dd(i:i)=','
        x=0
        y=0
        z=0
        read(dd(1:),*,end=15,err=100)x,y,z
15      alpha=x+y/60.+z/3600.
        do 16 i=k+1,k+6
        if(cc(i:i).ne.' ')goto 17
16      continue
        goto 100
17      dd=cc(i:)
        k1=index(dd,' ')
        k2=index(dd,',')
        k=min(k1,k2)
        if(k1.eq.0)k=k2
        if(k2.eq.0)k=k1
        do 20 i=1,k
20      if(dd(i:i).eq.':')dd(i:i)=','
        x=0
        y=0
        z=0
        read(dd(1:K),*,end=30,err=100)x,y,z
30      if(x.lt.0.)x=-x
        delta=x+y/60.+z/3600.
        if(index(cc,'-').ne.0)delta=-delta
100     end


	subroutine gethead(f1,head)
	character*80 f1,head(108)
	open(52,file=f1,status='old',access='direct',recl=8640)
	read(52,rec=1)head
	close (52)
	end
	
        function indexpos(head,f1)
        character*80 head(108),f1*8
        do 10 indexpos=1,108
10      if(head(indexpos)(1:8).eq.f1)return
        end

        subroutine geta8(head,a8)
        real*8 a8(8)
        character*80 head(72)
        ipos=indexpos(head,'A81     ')-1
        do 10 i=1,8
10      read(head(i+ipos)(12:),*)a8(i)
        end

        subroutine getsky(head,sky)
        character*80 head(72)
        ipos=indexpos(head,'VOLT2   ')
        read(head(ipos)(20:),*)sky
	end
