c
	parameter (n1=4096)
	parameter (n2=4032)
	character*69 f1,aa
	integer a(90000)
	real b(90000,8)
	character*80 head(72)
	logical logi
        real a1(n1*n2),b1(n1*n2)
	
	k=iargc()
	if(k.lt.1)then
	  write(*,*)
	  write(*,*)'***** Usage: u_4to1 file (2010,10)********'
	  write(*,*)
	  write(*,*)'   for 4096*4032 ccd, 90,000 stars'
	  write(*,*)
	  stop
	endif
	call getarg(1,f1)
        k=index(f1,'.')
	if(k.ne.0)f1=f1(1:k-1)//'  '

	k=lnblnk(f1)
ccccccccc for als
        call system('cp '//f1(1:lnblnk(f1))//'1.als '
     c                   //f1(1:lnblnk(f1))//'1.tmp')

        call j_23(f1,2,1,a,b)      ! produce 2~3.tmp
        call j_23(f1,3,1,a,b)      ! according overlap, moidfy mag
        call j_4(f1,4,2,3,a,b)     ! delta_mag+?.tmp=?.als
 
        call j_1(f1,a,b,n)           ! 4 tmp --> all.als  coo_star j_xy4ad
        call lst_4to1(f1,a,b,n)      ! 4 lst --> all.lst  psf_star j_xy6ad
                                     ! according above als_mag, put in
c for fits
	aa=f1(1:k)//'.fit'//char(0)
	call gethead(aa,head)
	call fits_4to1(f1,head,a1,b1,n1,n2)
	write(*,*)'ok'

	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

	subroutine fits_4to1(f1,head,a,b,n1,n2)
	character*69 f1,f2
	character*80 head(72)
        real a(n1,n2),b(n1,n2)
        logical logi
        integer mx(4),my(4)
        data mx/0,2048,0,2048/
        data my/0,0,2016,2016/

	k=lnblnk(f1)
	do 11 m=1,4
	f2=f1(1:k)//char(m+48)//'.fits'//char(0)
        inquire(file=f2,exist=logi)
	if(.not.logi)return
        m1=n1
	m2=n2
        call readfits(f2,b,m1,m2)    ! n1,n2 as para, call c routine, must var.
	do 10 i=mx(m)+1,mx(m)+n1/2
	do 10 j=my(m)+1,my(m)+n2/2
10	a(i,j)=b(i,j)
11      write(*,*)f2
	f2=f1(1:k)//'.fits'//char(0)
	write(*,*)f2
	call swap4(a,n1*n2*4)

	call writefits(f2,head,80*72,1)
	call writefits(f2,a,n1*n2*4,0)
	end

	subroutine gethead(f1,head)
	character*80 f1,head(72)
	open(52,file=f1,status='old',access='direct',recl=5760)
	read(52,rec=1)head
	close (52)
	end

	subroutine lst_4to1(f1,a,b,n)  ! deal with *.lst, psf star
	character*69 f1,aa
	integer a(50000)
	real b(50000,8)
	real d1(750),d2(750),d3(750),d4(750),d5(750),d6(750)
	integer di(750)
	call system('cat '//f1(1:lnblnk(f1))//'?.lst >1.tmp')
	call system('sort 1.tmp >2.tmp')
	open(2,file='2.tmp',status='old')
	j=0
	k=0
10	read(2,'(a)',end=11)aa
	read(aa(1:),*)i
	if(i.eq.j)goto 10
	k=k+1
	di(k)=i
	read(aa(40:),*)d6(k)
	do 13 j=1,n
	if(i.eq.a(j))then
	d1(k)=b(j,1)
	d2(k)=b(j,2)
	d3(k)=b(j,3)
	d4(k)=b(j,4)
	d5(k)=b(j,5)
	goto 101
	endif
13	continue
	goto 10
101	j=i
	goto 10
11	close(2)
	call sortn(k,di,d1,d2,d3,d4,d5,d6)
	open(1,file=f1(1:lnblnk(f1))//'.lst',status='unknown')
	do 12 i=1,k
12	write(1,'(i6,6f10.3)')di(i),d1(i),d2(i),d3(i),d4(i),d5(i),d6(i)	
	close(1)
        write(*,*)'all psf star:',k
	end

	subroutine sortn(k,di,d1,d2,d3,d4,d5,d6)
	real d1(270),d2(270),d3(270),d4(270),d5(270),d6(270)
	integer di(270)
c	  write(*,*)k
	do 10 i=1,k-1
	do 10 j=i+1,k
	if(d3(i).gt.d3(j))then
	  n=di(i)
	  di(i)=di(j)
	  di(j)=n
	  x=d1(i)
	  d1(i)=d1(j)
	  d1(j)=x
	  x=d2(i)
	  d2(i)=d2(j)
	  d2(j)=x
	  x=d3(i)
	  d3(i)=d3(j)
	  d3(j)=x
	  x=d4(i)
	  d4(i)=d4(j)
	  d4(j)=x
	  x=d5(i)
	  d5(i)=d5(j)
	  d5(j)=x
	  x=d6(i)
	  d6(i)=d6(j)
	  d6(j)=x
	endif
10	continue
	end

ccccccccccccccccccccccccccccccc

	subroutine j_23(f1,mm,m0,a,b)
	character*69 f1,f2,f3,c1,c2

	integer a(50000),num(6000)
	real b(50000,8)
	real a1(6000),a2(6000),z(6000),w(6000)

	k=lnblnk(f1)
	f2=f1(1:k)//char(mm+48)//'.als'
	open(1,file=f2,status='old')
	read(1,'(a)')head
	read(1,'(a)')head
	read(1,'(a)')head
	i=0
10	i=i+1
	read(1,*,err=100,end=100)a(i),(b(i,j),j=1,4)
	goto 10
100	n=i-1
	close(1)
	
	f2=f1(1:k)//char(m0+48)//'.tmp'
	open(2,file=f2,status='old')
	read(2,'(a)')head
	read(2,'(a)')head
	read(2,'(a)')head
	ip=0
20	read(2,*,err=200,end=200)k,x,y,xm,ee2
	do 21 i=1,n
	if(k.eq.a(i))then
	  ip=ip+1
	  num(ip)=k
	  a1(ip)=b(i,3)
	  a2(ip)=xm
	  w(ip)=1./sqrt(b(i,4)*b(i,4)+ee2*ee2)
	  z(ip)=a1(ip)-a2(ip)
	endif
21	continue
        goto 20
200	close(2)
	call sigma(z,w,ax,xx,ip)
	write(*,201)n,ip,ax,xx,mm
201	format(2i6,2f7.3,i4)
c	xx=3*xx                 !! notice
	do 40 i=1,ip
40	if(abs(z(i)-ax).gt.xx) num(i)=-1
	j=0
	do 50 i=1,ip
	j=j+1
	z(j)=z(i)
	a1(j)=a1(i)
	a2(j)=a2(i)
	w(j)=w(i)
	num(j)=num(i)
	if(num(j).eq.-1)j=j-1
50	continue
	ip=j
	call sigma(z,w,ax,xx,ip)
        if(ip.eq.0)ax=0
	k=lnblnk(f1)
	f2=f1(1:k)//char(mm+48)//'.als'
	open(2,file=f2,status='old')
	f3=f1(1:k)//char(mm+48)//'.tmp'
	open(3,file=f3,status='unknown')
	read(2,'(a)')c1
	read(2,'(a)')c2
        read(2,'(a)')head
	write(3,'(a)')c1
	write(3,'(a)')c2
	write(3,'(a)')
60	read(2,*,err=70,end=70)k,x1,x2,x3,x4,x5,x6,x7,x8
	x3=x3-ax
	if(x5.lt.9999.99 .and. x5.gt.-999.99)then
	write(3,2)k,x1,x2,x3,x4,x5,x6,x7,x8
2       FORMAT (I6, 3F9.3, F9.4, F9.3, F9.0, 2F9.3)
        ELSE
	write(3,3)k,x1,x2,x3,x4,x5,x6,x7,x8
3       FORMAT (I6, 3F9.3, F9.4, F9.2, F9.0, 2F9.3)
 	endif 
	goto 60
70	close(3)	
	write(*,201)n,ip,ax,xx
	end

cccccccccccccccccccccccccccccccc
	subroutine j_4(f1,mm,m1,m2,a,b)
	character*69 f1,f2,f3,c1,c2
	integer a(50000),num(6000)
	real b(50000,8),a1(6000),a2(6000),z(6000),w(6000)
	
	k=lnblnk(f1)
	f2=f1(1:k)//char(mm+48)//'.als'
	open(1,file=f2,status='old')
	read(1,'(a)')head
	read(1,'(a)')head
	read(1,'(a)')head
	i=0
10	i=i+1
	read(1,*,err=100,end=100)a(i),(b(i,j),j=1,4)
	goto 10
100	n=i-1
	close(1)

        f2=f1(1:k)//'1.als'
        open(2,file=f2,status='old')
        read(2,'(a)')head
        read(2,'(a)')head
        read(2,'(a)')head
        ip=0
20      read(2,*,err=200,end=200)k,x,y,xm,ee2
        do 21 i=1,n
        if(k.eq.a(i))then
          ip=ip+1
          num(ip)=k
          a1(ip)=b(i,3)
          a2(ip)=xm
          w(ip)=1./sqrt(b(i,4)*b(i,4)+ee2*ee2)
          z(ip)=a1(ip)-a2(ip)
          a(i)=-1
        endif
21      continue
        goto 20
200     close(2)

	k=lnblnk(f1)
	f2=f1(1:k)//char(m1+48)//'.tmp'
	open(2,file=f2,status='old')
	read(2,'(a)')head
	read(2,'(a)')head
	read(2,'(a)')head
320	read(2,*,err=300,end=300)k,x,y,xm,ee2
	do 321 i=1,n
	if(k.eq.a(i))then
	  ip=ip+1
	  num(ip)=k
	  a1(ip)=b(i,3)
	  a2(ip)=xm
	  w(ip)=1./sqrt(b(i,4)*b(i,4)+ee2*ee2)
	  z(ip)=a1(ip)-a2(ip)
	  a(i)=-1
	endif
321	continue
        goto 320
300	close(2)

c	write(*,*)'ip ',ip
	k=lnblnk(f1)
	f2=f1(1:k)//char(m2+48)//'.tmp'
	open(2,file=f2,status='old')
	read(2,'(a)')head
	read(2,'(a)')head
	read(2,'(a)')head
420	read(2,*,err=400,end=400)k,x,y,xm,ee2
	do 421 i=1,n
	if(k.eq.a(i))then
	  ip=ip+1
	  num(ip)=k
	  a1(ip)=b(i,3)
	  a2(ip)=xm
	  w(ip)=1./sqrt(b(i,4)*b(i,4)+ee2*ee2)
	  z(ip)=a1(ip)-a2(ip)
	  a(i)=-1
	endif
421	continue
        goto 420
400	close(2)

c	write(*,*)'ip ',ip

	call sigma(z,w,ax,xx,ip)
	write(*,201)n,ip,ax,xx,mm
201	format(2i6,2f7.3,i4)
	      if(xx.gt.0.3)xx=0.3
	do 40 i=1,ip
40	if(abs(z(i)-ax).gt.xx) num(i)=-1
	j=0
	do 50 i=1,ip
	j=j+1
	z(j)=z(i)
	a1(j)=a1(i)
	a2(j)=a2(i)
	w(j)=w(i)
	num(j)=num(i)
	if(num(j).eq.-1)j=j-1
50	continue
	ip=j
	call sigma(z,w,ax,xx,ip)
        if(ip.eq.0)ax=0
	k=lnblnk(f1)
	f2=f1(1:k)//char(mm+48)//'.als'
	open(2,file=f2,status='old')
	f3=f1(1:k)//char(mm+48)//'.tmp'
	open(3,file=f3,status='unknown')
	read(2,'(a)')c1
	read(2,'(a)')c2
        read(2,'(a)')head
	write(3,'(a)')c1
	write(3,'(a)')c2
	write(3,'(a)')
60	read(2,*,err=70,end=70)k,x1,x2,x3,x4,x5,x6,x7,x8
	x3=x3-ax
	if(x5.lt.9999.99 .and. x5.gt.-999.99)then
	write(3,2)k,x1,x2,x3,x4,x5,x6,x7,x8
2       FORMAT (I6, 3F9.3, F9.4, F9.3, F9.0, 2F9.3)
        ELSE
	write(3,3)k,x1,x2,x3,x4,x5,x6,x7,x8
3       FORMAT (I6, 3F9.3, F9.4, F9.2, F9.0, 2F9.3)
 	endif 
	goto 60
70	close(3)	
	write(*,201)n,ip,ax,xx
	end

	subroutine sigma(z,w,ax,xx,ip)
	real z(1),w(1)
	ax=0
	xx=0
	do 40 i=1,ip
	ax=ax+z(i)*w(i)
40	xx=xx+w(i)
	ax=ax/xx
	xx=0
	do 50 i=1,ip
50	xx=xx+(z(i)-ax)**2
	xx=sqrt(xx/ip)
	end

ccccccccccccccccccccccccccccccccc
	
	subroutine j_1(f1,a,b,n)

	character*69 f1,f2,c1,c2
	real b(50000,8)
	integer a(50000)

        integer mx(4),my(4)
        data  mx/0,2048,0,2048/
        data  my/0,0,2016,2016/
	
	i=0
        do 5 iz=1,4
          xx1=mx(iz)
          xx2=xx1+mx(4)
          yy1=my(iz)
          yy2=yy1+my(4)
	k=lnblnk(f1)
	f2=f1(1:k)//char(48+iz)//'.tmp'
	open(1,file=f2,status='old')
	read(1,'(a)')c1
	read(1,'(a)')c2
	read(1,'(a)')j
10	i=i+1
	read(1,*,end=100)a(i),(b(i,j),j=1,8)
         x=b(i,1)
         y=b(i,2)
        if(x.le.xx1.or.x.gt.xx2.or.y.le.yy1.or.y.gt.yy2)i=i-1
	goto 10
100	i=i-1
	close(1)
5	write(*,*)f2,i
        n=i

	call sortm(a,b,n)

	f2=f1(1:k)//'.als'
	open(3,file=f2,status='unknown')
        write(3,'(a)')c1
        write(3,'(a)')c2
        write(3,'(a)')
	do 70 i=1,n
        if(b(i,5).lt.9999.99 .and. x5.gt.-999.99)then
        write(3,2)a(i),(b(i,j),j=1,8)
2       FORMAT (I6, 3F9.3, F9.4, F9.3, F9.0, 2F9.3)
        ELSE
        write(3,3)a(i),(b(i,j),j=1,8)
3       FORMAT (I6, 3F9.3, F9.4, F9.2, F9.0, 2F9.3)
        endif
70	continue
        close(3)
	end

      subroutine sortm(a,b,n)
      integer a(1)
      real b(50000,8),u(8)
      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(a(i).le.a(j))goto 70
      it=a(j)
	do k=1,8
	u(k)=b(j,k)
	enddo
   40 a(j)=a(i)
        do k=1,8
	b(j,k)=b(i,k)
	enddo	
      j=i
      i=i-int
      if(i.le.0)goto 60
      if(a(i).gt.it)goto 40
   60 a(j)=it
	do k=1,8
	b(j,k)=u(k)
	enddo
   70 continue
      if(int.gt.1)goto 20
      return
      end