c jiang zhaoji libary

	subroutine center9(i1,j1,a,kt,ix,iy,iv,n1)
	integer t1(200)
	integer*2 c1(200),a(2048,2048)
	iv=0
	n1=0
	i=i1
	j=j1
10	k=0
	iu=j
	id=j
	if(a(i,j).lt.kt)goto 100	! find end
	k=k+a(i,j)
	a(i,j)=-a(i,j)

20	if(a(i,iu-1).lt.kt)goto 30	! up find end
	iu=iu-1
	  if(iu.lt.1)goto 30
	k=k+a(i,iu)	
	a(i,iu)=-a(i,iu)
	goto 20
30	if(a(i,id+1).lt.kt)goto 40	! down find end
	id=id+1
	  if(id.gt.2048)goto 40
	k=k+a(i,id)
	a(i,id)=-a(i,id)
	goto 30
40	iv=iv+k
	if(n1.lt.200)n1=n1+1
	t1(n1)=k
        c1(n1)=iu+id	
	i=i+1
	j=c1(n1)*.5
	if(a(i,j).lt.a(i,j+1))j=j+1
	goto 10

100     iu=t1(1)
	id=1
	do 110 k=2,n1
	if(t1(k).gt.iu)then
	  iu=t1(k)
	  id=k
	endif
110	continue
	ix=i1+id-1
	iu=0
	do 120 k=1,n1
120	iu=iu+c1(k)
	iy=iu/2/n1	
	end

	subroutine unlink(f1)
	character*(*) f1
	character*80 f2
	f2='rm '//f1
	call system(f2)
	end

	subroutine turnx2(a,n1,n2,it)
	integer*2 a(n1*n2)
	integer*2 b(2048*2048)
	n=n1*n2
	do 110 i=1,n
110	b(i)=a(i)
	k=it/90
	if(k.eq.0.or.k.gt.3)goto 100
	if(k.lt.0)k=4-k
	if(k.eq.8)k=4
	goto (10,20,30,40,50,60,70),k
100	return
10	do 11 j=1,n1
	jj=(j-1)*n2
	kk=1-j
	do 11 i=1,n2
11	a(jj+i)=b(kk+i*n1)
	goto 200
20	kk=n+1
	do 21 i=1,n
21	a(i)=b(kk-i)
	return
30	do 31 j=1,n1
	jj=(j-1)*n2
	kk=n1*n2+j
	do 31 i=1,n2
31	a(jj+i)=b(kk-i*n1)
	goto 200
40	do 41 j=1,n2
	jj=(j-1)*n1
	kk=j*n1+1
	do 41 i=1,n1
41	a(jj+i)=b(kk-i)
	return
50	do 51 j=1,n1
	jj=(j-1)*n2
	kk=j-n1
	do 51 i=1,n2
51	a(jj+i)=b(kk+i*n1)
	goto 200
60	do 61 j=1,n2
	kk=(n2-j)*n1
	jj=(j-1)*n1
	do 61 i=1,n1
61	a(jj+i)=b(kk+i)
	return
70	do 71 j=1,n1
	kk=n1-j+1+n1*n2
	jj=(j-1)*n2
	do 71 i=1,n2
71	a(jj+i)=b(kk-i*n1)
200	n=n1
	n1=n2
	n2=n
	end
 
	subroutine readxfile(file,ihead,map,nbyte)
	byte map(nbyte)
	character*(*) file
	logical qopen
	do 5 lun=99,10,-1
	inquire(unit=lun,opened=qopen)
	if(.not.qopen)goto 6
5	continue
6	open(lun,file=file,status='old',access='direct',recl=512)
	k=1
	iend=nbyte/512
	do 10 i=ihead+1,iend+ihead
	read(lun,rec=i,err=100)(map(j),j=k,k+511)
10	k=k+512
100	close (lun)
	end

	subroutine whitexblack(map,n1,n2,peak,sigma)
	integer*2 map(n1,n2)
	call whitexblackx1(map,n1,n2,peak,sigma)
	end

	subroutine whitexblackx1(map,n1,n2,xpeak,sigma)
	integer*2 map(*)
	integer*4 ihist(32867)
c *** clear histogram array
	do 10 i=1,32867
   10	ihist(i)=0
c *** find default levels
	imax=0
	do 60 i=1,n1*n2
 	ii=map(i)
	ii=ii+101
	if(ii.le.0.or.ii.gt.32867)goto 60
	ihist(ii)=ihist(ii)+1
	imax=max(ii,imax)
   60	continue
	call histat(ihist,mode,maxh,xmean,xpeak,sigma,gsigma,imax)
   	xpeak=xpeak-100.0
	sigma=amin1(sigma,gsigma)
	sigma=max(1.0,sigma)
	end

	subroutine histat(ihist,mode,maxh,xmean,xpeak,fsigma,gsigma,n)
	dimension ihist(n),xpar(101),buf(201)
c *** computes statistics of histogram
	rttpi=sqrt(8.0*atan(1.0))
	xcount=0.0
	maxh=0
	mode=0
c *** first get no. of data values in histogram
	do 50 i=1,n
	ii=n+1-i
	if(ihist(ii).eq.0)goto 50
	goto 60
   50 continue
   60 nn=ii
	ihtot=0
	do 70 i=1,nn
   70 ihtot=ihtot+ihist(i)
c *** trap for saturation or holes
c *** only look for mode in 10% to 90% range
	ihsum=0
	do 100 i=1,nn
	ihsum=ihsum+ihist(i)
	if(ihsum.lt.0.1*ihtot.or.ihsum.gt.0.9*ihtot)goto 100
	if(ihist(i).le.maxh)goto 100
	maxh=ihist(i)
	mode=i
  100 continue
	xmean=0.0
	sigmsq=0.0
	do 200 i=1,nn
	temp=float(ihist(i))
	xcount=xcount+temp
	xmean=xmean+float(i)*temp
  200 sigmsq=sigmsq+float(i)**2*temp
	if(xcount.lt.0.5.or.maxh.eq.0)then
	mode=-9999
	maxh=-9999
	xmean=-9999.0
	xpeak=-9999.0
	sigma=0.0
	gsigma=0.0
		return
	endif
c *** median
	xmed=0.0
	do 220 i=1,nn
	xmed=xmed+ihist(i)
	if(xmed.gt.0.5*xcount)then
	xmed=i
	goto 230
	endif
  220 continue
  230 continue
	xmean=xmean/xcount
ccccc	sigma=sqrt(sigmsq/xcount-xmean**2)
c *** parabolic interpolation to find maximum
	gsigma=xcount/(float(maxh)*rttpi)
	if(gsigma.gt.0.1*float(nn))gsigma=0.1*float(nn)
	isigm=max0(3,int(gsigma+0.5))
	isigm=min0(50,isigm)
	nfilt=isigm*4+1
	noff=(nfilt+1)/2
c *** form filter array
	conv=0.5/float(isigm**2)
	do 250 i=1,noff
	xx=i-noff
	buf(i)=exp(-conv*xx**2)
  250 buf(nfilt+1-i)=buf(i)
	xpeak=float(mode)
	do 310 k=1,2
c *** make sure ixpeak bounded by data range
	ixpeak=min(nn,int(xpeak+0.5))
	ixpeak=max(1,ixpeak)
	ilow=max(1,ixpeak-isigm)
	ihih=min(nn,ixpeak+isigm)
	m=ihih-ilow+1
	ii=0
	cogd=0.0
	cogn=0.0
	do 300 i=ilow,ihih
	cogd=cogd+ihist(i)
	cogn=cogn+ihist(i)*float(i)
	temp=0.0
	xnumb=0.0
	do 275 j=1,nfilt
	jj=j-noff
	if(i+jj.lt.1.or.i+jj.gt.nn)goto 275
	temp=temp+ihist(i+jj)*buf(j)
	xnumb=xnumb+buf(j)
  275 continue
	if(xnumb.gt.0.0)then
	temp=temp/xnumb
	else
	temp=0.0
	endif
	ii=ii+1
  300 xpar(ii)=alog(amax1(1.0,temp))
	call parbol(xpar,m,xpeak,sd)
	fsigma=sqrt(amax1(0.0,sd**2-isigm**2))
	xpeak=xpeak+float(ilow-1)
	if(cogd.gt.0.5)then
	  cog=cogn/cogd
	else
	  cog=xpeak
	endif
	if(abs(xpeak-mode).gt.gsigma)xpeak=cog
  310 continue
	mode=mode-1
	xpeak=xpeak-1.0
	xmean=xmean-1.0
	xmed=xmed-1.0
c *** another trap for saturation or holes
	if(xpeak.gt.nn-gsigma)xpeak=xmed
	if(xpeak.lt.gsigma)xpeak=xmed
c *** find sigma by 68% cutoff level
	ixpeak=int(xpeak+1.5)
	iwid=0
	icount=ihist(ixpeak)
	ilim=int(0.68*xcount+0.5)
	do 400 i=1,nn
	ii=ixpeak+i
	if(ii.gt.nn)goto 350
	iwid=iwid+1
	icount=icount+ihist(ii)
	if(icount.gt.ilim)goto 500
  350 ii=ixpeak-i
	if(ii.lt.1)goto 375
	iwid=iwid+1
	icount=icount+ihist(ii)
  375 if(icount.gt.ilim)goto 500
  400 continue
  500 gsigma=0.5*iwid
	return
	end

	subroutine solve(a,b,m)
	dimension a(25,25),b(25)
c *** gauss elimination to solve ax=b
	iu=m-1
	do 200 i=1,iu
c *** find largest remaining term in ith column for pivot
	big=0.0
	do 50 k=i,m
	rmax=abs(a(k,i))
	if(rmax-big)50,50,30
   30 big=rmax
	l=k
   50 continue
c *** check for non-zero term
	if(big.eq.0.0)then
	do ib=1,m
	b(ib)=0.0
	enddo
    	write(*,*) ' Zero determinant from SOLVE'
	return
	endif
	if(i.eq.l)goto 120
c *** switch rows
	do 100 j=1,m
	temp=a(i,j)
	a(i,j)=a(l,j)
  100 a(l,j)=temp
	temp=b(i)
	b(i)=b(l)
	b(l)=temp
c *** pivotal reduction
  120 pivot=a(i,i)
	jl=i+1
	do 200 j=jl,m
	temp=a(j,i)/pivot
	b(j)=b(j)-temp*b(i)
	do 200 k=i,m
  200 a(j,k)=a(j,k)-temp*a(i,k)
c *** back substitution for solution
	do 500 i=1,m
	ir=m+1-i
	if(a(ir,ir).eq.0.0)goto 490
	temp=b(ir)
	if(ir.eq.m)goto 480
	do 450 j=2,i
	k=m+2-j
  450 temp=temp-a(ir,k)*b(k)
  480 b(ir)=temp/a(ir,ir)
	goto 500
  490 b(ir)=0.0
  500 continue
	return
	end

	subroutine parbol(a,n,xpeak,sd)
	dimension a(n),cmatx(25,25),bvect(25)
c *** least-squares parabolic curve fitting
	do 100 i=1,25
	bvect(i)=0.0
	do 100 j=1,25
  100 cmatx(i,j)=0.0
	cmatx(1,1)=float(n)
	suma=0.0
	sumay=0.0
	sumayy=0.0
	sumy1=0.0
	sumy2=0.0
	sumy3=0.0
	sumy4=0.0
	do 200 i=1,n
	aa=a(i)
	fi=float(i)
	suma=suma+aa
	sumay=sumay+aa*fi
	sumayy=sumayy+aa*fi**2
	sumy1=sumy1+fi
	sumy2=sumy2+fi**2
	sumy3=sumy3+fi**3
  200 sumy4=sumy4+fi**4
	cmatx(1,2)=sumy1
	cmatx(1,3)=sumy2
	cmatx(2,1)=sumy1
	cmatx(2,2)=sumy2
	cmatx(2,3)=sumy3
	cmatx(3,1)=sumy2
	cmatx(3,2)=sumy3
	cmatx(3,3)=sumy4
	bvect(1)=suma
	bvect(2)=sumay
	bvect(3)=sumayy
	call solve(cmatx,bvect,3)
	if(bvect(3).eq.0.0)then
	xpeak=0.0
	sd=0.0
	return
	endif
	xpeak=-bvect(2)/(2.0*bvect(3))
	sdsq=1.0/(2.0*abs(bvect(3)))
	if(sdsq.lt.0.0)sdsq=0.0
	sd=sqrt(sdsq)
	return
	end

	subroutine toms(aa,c,k)
c hour or degree to char_line
c if k=1 (hour case) in **:**:**.** format
c    k=0 (degree        **:**:**.*  format
c    k=2 (degree)  in  ***:**:**.*  format
	character*12 c
	a=aa
	c(1:1)=' '
	if(a.lt.0.)c(1:1)='-'
	if(a.lt.0.)a=-a
	i=a
	b=(a-i)*60.
	j=b
	x=(b-j)*60.
	if(j.eq.60)then
	j=0
	i=i+1
	endif
	if(k.eq.1)write(c(2:),1)i,j,x
1	format(i2.2,':',i2.2,':',f5.2)
	if(k.eq.0)write(c(2:),2)i,j,x
2       format(i2.2,':',i2.2,':',f4.1)
        if(k.eq.2)write(c(1:),3)i,j,x
3       format(i3.2,':',i2.2,':',f4.1)
	if(c(8:8).eq.' ')c(8:8)='0'
	if(c(9:9).eq.' ')c(9:9)='0'
c solve 60 problem, ex 12:09:00 alway display 12:08:60.0
	if(c(8:8).eq.'6')then
	c(8:8)='0'
	j=j+1
	write(c(5:6),'(i2.2)')j
	if(c(5:5).eq.'6')then
	c(5:5)='0'
	i=i+1
	if(k.ne.2)write(c(2:3),'(i2.2)')i
	if(k.eq.2)write(c(1:3),'(i3.2)')i
	endif
	endif
	end

	function itohd(cc,alpha,delta,ep)
c 12:23 -23:00:2.3 1950 in a char_line, give value to 3 avariable
c 12:23 -23:00:2.3       is also ok, but last avariable get zero value
c 12:34:23.5 to 12.578567
	character*(*) cc
	character*20 dd
	byte bell
	bell=7
	itohd=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
	ep=0.
	read(dd(k+1:),*,end=40,err=100)ep	
40	if(ep.eq.0.or.(ep.ge.1900..and.ep.le.2100.))return
100	itohd=0
	write(*,'(1x,a1)')bell
	end

      subroutine astprs (ra1, dec1, ep1, ra2, dec2, ep2)
c come from IRAF in spp language.
c epoch operate
c alpha in hour
c delta in degree
      real*8 epoch1,epoch2,r0(3),r1(3),p(3,3),arc
	epoch1=ep1
	epoch2=ep2
	arc=45d0/atan(1d0)
            ra2 = ra1
            dec2 = dec1
	 if(epoch1.eq.epoch2)return
         ra2 = ra1 *15d0/arc
         dec2 =dec1/arc
         r0(1) = cos (ra2) * cos (dec2)
         r0(2) = sin (ra2) * cos (dec2)
         r0(3) = sin (dec2)
         if (epoch1 .eq. 2000d0) goto 120
            call astrox (epoch1, p)
            r1(1) = p(1, 1) * r0(1) + p(1, 2) * r0(2) + p(1, 3) * r0(3)
            r1(2) = p(2, 1) * r0(1) + p(2, 2) * r0(2) + p(2, 3) * r0(3)
            r1(3) = p(3, 1) * r0(1) + p(3, 2) * r0(2) + p(3, 3) * r0(3)
            r0(1) = r1(1)
            r0(2) = r1(2)
            r0(3) = r1(3)
120     if (epoch2 .eq. 2000d0) goto 130
            call astrox (epoch2, p)
            r1(1) = p(1, 1) * r0(1) + p(2, 1) * r0(2) + p(3, 1) * r0(3)
            r1(2) = p(1, 2) * r0(1) + p(2, 2) * r0(2) + p(3, 2) * r0(3)
            r1(3) = p(1, 3) * r0(1) + p(2, 3) * r0(2) + p(3, 3) * r0(3)
            r0(1) = r1(1)
            r0(2) = r1(2)
            r0(3) = r1(3)
130      ra2 = atan2 (r0(2), r0(1))/15d0*arc
         dec2 = asin (r0(3))*arc
         if (ra2.lt.0d0)ra2=ra2+24.
      end

      subroutine astrox (epoch, p)
c sub of asrprs
      real*8 epoch,p(3,3)
      real*8 t,a,b,c,ca,cb,cc,sa,sb,sc,astjuy,arc
	arc=45d0/atan(1d0)
         t = (astjuy (epoch) - 2451545.0) / 36525
         a = t * (0.6406161 + t * (0.0000839 + t * 0.0000050))
         b = t * (0.6406161 + t * (0.0003041 + t * 0.0000051))
         c = t * (0.5567530 - t * (0.0001185 + t * 0.0000116))
         ca = cos (a/arc)
         sa = sin (a/arc)
         cb = cos (b/arc)
         sb = sin (b/arc)
         cc = cos (c/arc)
         sc = sin (c/arc)
         p(1,1) = ca * cb * cc - sa * sb
         p(2,1) = -sa * cb * cc - ca * sb
         p(3,1) = -cb * sc
         p(1,2) = ca * sb * cc + sa * cb
         p(2,2) = -sa * sb * cc + ca * cb
         p(3,2) = -sb * sc
         p(1,3) = ca * sc
         p(2,3) = -sa * sc
         p(3,3) = cc
      end

      function astjuy (epoch)
c sub of astrox of astprs
      real*8 astjuy,epoch,jd
      integer year,centuy
         year = int (epoch) - 1
         centuy = year / 100
         jd =1721425.5+365*year-centuy+int(year/4)+int(centuy/4)
         astjuy = jd + (epoch - int(epoch)) * 365.25
      end

	subroutine standc(rc,dc,ra,de,xi,xn)
c alpha delta --------> X Y
c rc dc : center of plate
c rc dc ra de in arc
	sd=sin(dc)
	cd=cos(dc)
	td=tan(de)
	co=cos(ra-rc)
	tmp=sd*td+cd*co
	xi=sin(ra-rc)/tmp
	xn=(cd*td-sd*co)/tmp
	end

	subroutine astand(rc,dc,xi,xn,ra,de)
c X Y -------------> alpha delta
	cd=cos(dc)
	td=tan(dc)
	ra=atan(xi/cd/(1.-xn*td))
	de=atan((xn+td)*cos(ra)/(1.-xn*td))
	ra=ra+rc
	if(abs(dc-de).gt.2.)then
	de=-de
	ra=ra+4.*atan(1.)
	endif
	end
	
	subroutine xytoad(x,a)
c 2 set 6 plate coef convert
	real*8 x(2,6),a(2,6),z
	z=x(1,1)*x(2,2)-x(1,2)*x(2,1)
	a(1,1)=x(2,2)/z
	a(1,2)=-x(1,2)/z
	a(1,3)=(x(1,2)*x(2,3)-x(1,3)*x(2,2))/z
	a(2,1)=-x(2,1)/z
	a(2,2)=x(1,1)/z
	a(2,3)=(x(1,3)*x(2,1)-x(1,1)*x(2,3))/z
	end

	subroutine plate(ad,xy,w,n,c,nc)
c come from Mike
	real*8 ad(2,n),xy(2,n),c(2,6),a(6,6),b(6),wk(6)
	byte w(n)
	do 100 ixy=1,2
	do 20 i=1,6
	b(i)=0.d0
	do 20 j=1,6
20	a(i,j)=0.d0
	do 40 k=1,n
	if(w(k).ne.1)goto 40
	wk(1)=xy(1,k)
	wk(2)=xy(2,k)
	wk(3)=1.d0
	z=ad(1,k)
	if(ixy.eq.2)z=ad(2,k)
	do 41 j=1,nc
	b(j)=b(j)+wk(j)*z
	do 41 i=1,nc
41	a(i,j)=a(i,j)+wk(i)*wk(j)
40	continue
	call solve8(a,b,nc)
	do 100 i=1,nc
100	c(ixy,i)=b(i)
	end

      subroutine solve8(a,b,m)
      real*8 a(6,6),b(6)
      real*8 big,rmax,temp,pivot
      iu=m-1
      do 200 i=1,iu
      big=0.d0
      do 50 k=i,m
      rmax=abs(a(k,i))
      if(rmax-big)50,50,30
   30 big=rmax
      l=k
   50 continue
      if(big.eq.0.0)then
	do ib=1,m
	b(ib)=0.0
	enddo
    	write(*,*)'Zero determinant from SOLVE'
	return
      endif
      if(i.eq.l)goto 120
      do 100 j=1,m
      temp=a(i,j)
      a(i,j)=a(l,j)
  100 a(l,j)=temp
      temp=b(i)
      b(i)=b(l)
      b(l)=temp
  120 pivot=a(i,i)
      jl=i+1
      do 200 j=jl,m
      temp=a(j,i)/pivot
      b(j)=b(j)-temp*b(i)
      do 200 k=i,m
  200 a(j,k)=a(j,k)-temp*a(i,k)
      do 500 i=1,m
      ir=m+1-i
      if(a(ir,ir).eq.0.0)goto 490
      temp=b(ir)
      if(ir.eq.m)goto 480
      do 450 j=2,i
      k=m+2-j
  450 temp=temp-a(ir,k)*b(k)
  480 b(ir)=temp/a(ir,ir)
      goto 500
  490 b(ir)=0.0
  500 continue
      end

	subroutine correct216(xi,xn)
        real d(4)
        data d/6.66661e-8,-4.69485e-3,-8.64714e-3,1.40871e2/
        r=sqrt(xi*xi+xn*xn)
        z=1.
        if(r.ne.0.)z=1.+d(1)/r+d(2)+d(3)*r+d(4)*r*r
        xi=xi*z
        xn=xn*z
        end

	function lnblnk(a)
	character*(*) a
	k=len(a)
10	if(a(k:k).le.' ')then
	k=k-1
	goto 10
	endif
	lnblnk=k
	end

	subroutine swap4(a,n)
        byte a(n),c
        do 10 i=1,n,4
        c=a(i)
        a(i)=a(i+3)
        a(i+3)=c
        c=a(i+1)
        a(i+1)=a(i+2)
        a(i+2)=c
10      continue
        end

	subroutine swap2(a,n)
        byte a(n),c
        do 10 i=1,n,2
        c=a(i)
        a(i)=a(i+1)
        a(i+1)=c
10      continue
        end	

	function ior(a,b)
	integer*2 a,b
	integer c,d
	integer*2 e(16)
	c=a
	d=b
	if(c.lt.0)c=c+65536
	if(d.lt.0)d=d+65536
	do 10 i=1,16
	e(i)=mod(d,2)+mod(c,2)
	c=c/2
10	d=d/2
	k=0
	do 20 i=16,1,-1
	k=k*2
20	if(e(i).ne.0)k=k+1	
	ior=k
	end

	function iand(a,b)
	integer*2 a,b
	integer c,d
	integer*2 e(16)
	c=a
	d=b
	if(c.lt.0)c=c+65536
	if(d.lt.0)d=d+65536
	do 10 i=1,16
	e(i)=mod(d,2)+mod(c,2)
	c=c/2
10	d=d/2
	k=0
	do 20 i=16,1,-1
	k=k*2
20	if(e(i).eq.2)k=k+1	
	iand=k
	end

	function inot(a)
	integer*2 a
	inot=65535-a
	end
