c combine 2048*2048 pic. parameter (nfile=89) real white(nfile),scale(nfile),d(nfile),s(nfile),g(nfile) integer*2 sx(nfile),sy(nfile),badtype(nfile) integer*2 bx(30000),by(30000) character*80 f0,f1,f2,head(72) real*8 a8(8),adcoef(2,3) real a(2048*2048) real starx(400),stary(400),starv(400) real p0(100) integer*2 px(100),py(100),ac real pz0(100),pz1(100) real b(2048*32,nfile) logical logi ac=iargc() if(ac.lt.2)then write(*,*) write(*,*)' ******** Combine 2048*2048 Pic. ***2003.5' write(*,*) write(*,*)' cob batch_file out_file [!]' write(*,*) write(*,*)' if cob.par is not existed in default directory' write(*,*)' default value will be taken as following' write(*,*)' 3.4 ! gain, before 2003,2,5' write(*,*)' 6.0 ! sigam for star' write(*,*)' 6.0 ! sigma for background' write(*,*)' 500000 ! default value, if saturating' write(*,*)' 8.26 ! gain, after 2003,2,5' write(*,*)' if default value is not suitable, you can edit' write(*,*)' a file named cob.par inclued 5 line in your' write(*,*)' default directory, ' write(*,*)' 1997,11' write(*,*)' after 2003,2,5 fits, at end, data*=8.26/3.4' write(*,*)' 2003,5; 2003,8 modify' write(*,*)' auto selected 3 bad_pixel tables' write(*,*)' zjzx.bad ljl.bad ccd5.bad' stop endif call getarg(1,f0) open(51,file=f0,status='old') i=0 20 i=i+1 read(51,'(a)',end=10)f2 goto 20 10 if(i.lt.4)stop 'too less files, must>=3' close(51) m=i-1 if(m.gt.nfile)m=nfile inquire(file='cob.par',exist=logi) if(logi)then open(51,file='cob.par',status='old') read(51,*)gain read(51,*)s61 read(51,*)s62 read(51,*)s50 read(51,*,end=501)gain2003 502 close (51) write(*,*)'parameters come from default directory' goto 1234 501 gain2003=8.0 goto 502 endif inquire(file='/vega2/rhbin/cob.par',exist=logi) if(logi)then open(51,file='/vega2/rhbin/cob.par',status='old') read(51,*)gain read(51,*)s61 read(51,*)s62 read(51,*)s50 read(51,*)gain2003 close (51) write(*,*)'parameters come from /vega2/rhbin' goto 1234 endif gain=3.4 s61=6. s62=6. s50=500000. gain2003=8.0 write(*,*)'default parameters:' 1234 n49=s50-1000.+0.1 write(*,*)'gain=',gain,' n1(star)=',s61,' n2(background)=',s62 write(*,*)'gain_2003=',gain2003 write(*,*)'saturated data=',s50,' high good data=',n49 open(51,file=f0,status='old') write(*,2) 2 format(35x,'shift(x,y) UT. exp. sky'/78('-')) nt=0 do 30 i=1,m if(i.eq.1)then read(51,'(a)')f1 ! mother file call gethead(f1,head) call geta8(head,a8) call xytoad(a8,adcoef) x=a8(7) y=a8(8) f2=f1 else read(51,'(a)')f2 call gethead(f2,head) endif ipos=indexpos(head,'A87 ') read(head(ipos)(20:),*)alpha read(head(ipos+1)(20:),*)delta call standc(x,y,alpha,delta,xi,xn) x2=adcoef(1,1)*xi+adcoef(1,2)*xn+adcoef(1,3)-1024.5 y2=adcoef(2,1)*xi+adcoef(2,2)*xn+adcoef(2,3)-1024.5 sx(i)=nint(x2) sy(i)=nint(y2) ipos1=indexpos(head,'TIME ') ipos2=indexpos(head,'EXPOSURE') read(head(ipos2)(26:),*)j nt=nt+j ipos3=indexpos(head,'VOLT2 ') if(head(ipos3)(42:46).ne."PIXEL")then write(*,*)"do statsky first !" stop endif read(head(ipos3)(22:),*)white(i) write(*,1)i,f2(1:22),x2,y2,sx(i),sy(i), * head(ipos1)(16:27),head(ipos2)(24:30),white(i) 1 format(i2,': ',a,2f6.1,2x,2i4,2a,f8.2) ipos=indexpos(head,'INSTRUME') if(ipos.gt.70)stop 'not our fits file**' badtype(i)=0 ! zjzx.bad if(head(ipos)(13:13).eq.'1')badtype(i)=1 ! new controller, ljl.bad if(indexpos(head,'MPP ').gt.70)badtype(i)=0 if(head(ipos)(13:13).eq.'5')badtype(i)=5 ! after 2003 if(head(ipos)(17:17).eq.'4')badtype(i)=4 ! thin ccd, no bad pixel if(ac.ge.3)badtype(i)=-1 c 20030808 ipos=indexpos(head,'DATE-OBS') g(i)=gain if(ipos.lt.70)then read(head(ipos)(22:29),'(i2,1x,i2,1x,i2)')id,im,iy if(iy.lt.90)then j=iy*10000+im*100+id if(j.gt.30205)g(i)=gain2003 endif endif 30 continue close(51) if(ac.ge.3)then write(*,*)" no bad_pixel corrected" goto 12345 endif inquire(file='zjzx.bad',exist=logi) if(logi)then write(*,*)'read zjzx.bad from default directory' open(52,file='zjzx.bad',status='old') else write(*,*)'read zjzx.bad from /vega2/rhbin' open(52,file='/vega2/rhbin/zjzx.bad',status='old') endif nb=1 301 read(52,*,end=302)bx(nb),by(nb) if(bx(nb).lt.1.or.bx(nb).gt.2048)goto 301 if(by(nb).lt.1.or.by(nb).gt.2048)goto 301 if(nb.lt.10000)nb=nb+1 goto 301 302 nb=nb-1 write(*,*)"readin badpixel table from file zjzx.bad:",nb close (52) inquire(file='ljl.bad',exist=logi) if(logi)then write(*,*)'read ljl.bad from default directory' open(52,file='ljl.bad',status='old') else write(*,*)'read ljl.bad from /vega2/rhbin' open(52,file='/vega2/rhbin/ljl.bad',status='old') endif lb=10001 3011 read(52,*,end=3021)bx(lb),by(lb) if(bx(lb).lt.1.or.bx(lb).gt.2048)goto 3011 if(by(lb).lt.1.or.by(lb).gt.2048)goto 3011 if(lb.lt.20000)lb=lb+1 goto 3011 3021 lb=lb-1 write(*,*)"readin badpixel table from file ljl.bad:",lb-10000 close (52) inquire(file='ccd5.bad',exist=logi) if(logi)then write(*,*)'read ccd5.bad from default directory' open(52,file='ccd5.bad',status='old') else write(*,*)'read ccd5.bad from /vega2/rhbin' open(52,file='/vega2/rhbin/ccd5.bad',status='old') endif mb=20001 4011 read(52,*,end=4021)bx(mb),by(mb) if(bx(mb).lt.1.or.bx(mb).gt.2048)goto 4011 if(by(mb).lt.1.or.by(mb).gt.2048)goto 4011 if(lb.lt.30000)mb=mb+1 goto 4011 4021 mb=mb-1 write(*,*)"readin badpixel table from file ccd5.bad:",mb-20000 close (52) ******* get each shift end ****** ****** get mother pic star begin ***** 12345 call readata(a,f1) call shrink4(a,a) nstar=400 call getstar(a,a(2048*1024+1), * starx,stary,starv,nstar,white(1),sx,sy,m,g(1),gain) write(*,*)'nstar: ',nstar k1=nstar/3. k2=k1+99 if(k2.gt.nstar)k2=nstar np=0 do 110 i=k1,k2 np=np+1 px(np)=starx(i)+0.5 110 py(np)=stary(i)+0.5 call readata(a,f1) write(*,*)f1(1:22),'platen file' call star100(a,sx(1),sy(1),white(1),px,py,pz0,np) call badfile(a,nb,lb,mb,bx,by,badtype(1),0,0,1) scale(1)=1. open(51,file=f0,status='old') read(51,'(a)')f2 c call pgbegin(0,'/su',2,2) do 120 i=2,m read(51,'(a)')f2 call readata(a,f2) call star100(a,sx(i),sy(i),white(i),px,py,pz1,np) call badfile(a,nb,lb,mb,bx,by,badtype(i),sx(i),sy(i),i) c call pgenv(0.,120.,-10.,10.1,0,0) do 130 j=1,np p0(j)=j 130 pz1(j)=pz1(j)/pz0(j) c call pgline(np,p0,pz1) scale(i)=xmedian(pz1,np) 120 write(*,*)f2(1:22),'scale:',scale(i) c call pgend close(51) ******** find scale end ******* do 215 i=1,m 215 open(10+i,file=char(i/10+64)//char(mod(i,10)+64)//'.tmp', * status='old',access='direct',recl=2048*32*4) do 216 i=1,m c 20030808 d(i)=sqrt(white(i)/g(i)) 216 s(i)=s62*d(i)/scale(i) sigma=xmedian(d,m) k=0 do 220 l=1,64 call dispdot() do 230 i=1,m 230 call readrec(10+i,l,b(1,i)) do 240 i=1,2048*32 do 235 j=1,m 235 d(j)=b(i,j) k=k+1 240 a(k)=allc(d,white,scale,m,sigma,s,s61,s50,g,gain) 220 continue do 250 i=1,m 250 close(10+i) write(*,*) open(51,file=f1,status='old',access='direct',recl=5760) read(51,rec=1)head close(51) ipos=indexpos(head,'END ') if(ipos.lt.72)then f2=head(ipos)(1:80) head(ipos)(1:80)=head(72)(1:80) head(72)(1:80)=f2 endif if(ipos.eq.73)head(72)(1:8)='END ' write(head(ipos2)(26:30),'(i5)')nt head(ipos2)(54:63)='(combined)' ipos=indexpos(head,'TTIME ') if(ipos.ne.73)write(head(ipos)(25:30),"('>',i5)")nt x1=statt(a,b,256,256) x2=statt(a,b,256,1280) x3=statt(a,b,768,768) x4=statt(a,b,1280,256) x5=statt(a,b,1280,1280) sky=amin1(x1,x2,x3,x4,x5) ipos=indexpos(head,'VOLT2 ') write(head(ipos)(22:),"(f9.3,' / SKY ADU/PIXEL')")sky c delete 4 edge call edge(a,sx,sy,m,ix1,ix2,iy1,iy2) write(*,*)' Edge_frame: ccd(i,j)=ccd(i,j)+5000.' write(*,*)' line_number: col~~',ix1,ix2,' raw~~',iy1,iy2 c add frame numbers, ix1,ix2, iy1,iy2 96,12 ipos=indexpos(head,'NAXIS ') if(ipos.ne.73)write(head(ipos)(60:77),"('frame',i5,i8)")m,n49 ipos=indexpos(head,'NAXIS1 ') if(ipos.ne.73)write(head(ipos)(60:74),"('range',2i5)")ix1,ix2 ipos=indexpos(head,'NAXIS2 ') if(ipos.ne.73)write(head(ipos)(60:74),"('range',2i5)")iy1,iy2 call getarg(2,f2) write(*,*)' combined file: '//f2(1:25) k=2048*2048*4+5760 call f_blank(a,ix1,ix2,iy1,iy2) call swap4(a,2048*2048*4) open(50,file=f2,status='unknown',access='direct',recl=k) write(50,rec=1)head,a close(50) do 200 i=1,m f0='rm '//char(i/10+64)//char(mod(i,10)+64)//'.tmp' 200 call system(f0) write(*,*)'re coordinating!' f0='coord '//f2(1:25)//' !' call system(f0) stop end subroutine f_blank(a,ix1,ix2,iy1,iy2) real a(2048,2048) if(ix1.gt.1)then do 10 j=1,2048 do 10 i=1,ix1 10 a(i,j)=0. endif if(ix2.lt.2048)then do 20 j=1,2048 do 20 i=ix2,2048 20 a(i,j)=0. endif if(iy1.gt.1)then do 30 i=1,2048 do 30 j=1,iy1 30 a(i,j)=0. endif if(iy2.lt.2048)then do 40 i=1,2048 do 40 j=iy2,2048 40 a(i,j)=0. endif end subroutine edge(a,sx,sy,m,ix1,ix2,iy1,iy2) real a(2048,2048) integer*2 sx(1),sy(1) ix1=1 ix2=-1 iy1=1 iy2=-1 do 10 i=1,m if(sx(i).lt.ix1)ix1=sx(i) if(sx(i).gt.ix2)ix2=sx(i) if(sy(i).lt.iy1)iy1=sy(i) 10 if(sy(i).gt.iy2)iy2=sy(i) ix1=-ix1 ix2=2048-ix2 iy1=-iy1 iy2=2048-iy2 k=iy1 if(k.lt.1)k=1 do 20 j=k,iy2 if(ix1.ne.0)a(ix1,j)=a(ix1,j)+5000. 20 if(ix2.ne.2048)a(ix2,j)=a(ix2,j)+5000. k=ix1 if(k.lt.1)k=1 do 30 i=k,ix2 if(iy1.ne.0)a(i,iy1)=a(i,iy1)+5000. 30 if(iy2.ne.2048)a(i,iy2)=a(i,iy2)+5000. if(ix1.eq.0)ix1=1 if(iy1.eq.0)iy1=1 end function statt(a,b,m1,m2) real a(2048,2048) integer*2 b(512,512) y=amin1(a(m1,m2),a(m1+512,m2+512),a(m1,m2+512), * a(m1+512,m2),a(m1+256,m2+256)) y=amin1(y,a(m1+256,m2),a(m1+512,m2+256), * a(m1+256,m2+512),a(m1,m2+256)) y=y*0.8 j1=m2 do 10 j=1,512 i1=m1 j1=j1+1 do 10 i=1,512 i1=i1+1 x=a(i1,j1)-y if(x.lt.-32767.)x=-32767. if(x.gt. 32767.)x= 32767. 10 b(i,j)=x call whitexblack(b,512,512,x,s) stat=x+y end function allc(d,white,scale,m,sigma,s,s61,s50,g,gain) real d(1),white(1),scale(1),s(1),g(1) real d1(99) ! normal value real d2(99) ! bad working unit bad=-999. ! s(i)=6*sigma(i)/scale(i) x=bad c 80717 j=0 k=0 do 211 i=1,m if(d(i).lt.-10000.)j=j+1 if(d(i).gt. 10000.)k=k+1 211 continue if(j.gt.0 .and. ((j+k)*3.ge.m*2))then c write(*,*) c write(*,*)(d(i),i=1,m) goto 220 endif c bad pixel n=0 do 5 i=1,m 5 if(d(i).lt.white(i)*.5)d(i)=bad do 10 i=1,m 10 if(d(i).le.bad)n=n+1 ! eq if(n.eq.m)goto 220 do 20 i=1,m 20 if(d(i).gt.bad)d1(i)=(d(i)-white(i))/scale(i) !ne if(n.eq.0)goto 100 n=0 do 30 i=1,m if(d(i).le.bad)goto 30 ! eq n=n+1 d2(n)=d1(i) 30 continue x=xmedian(d2,n) do 40 i=1,m 40 if(d(i).le.bad)d1(i)=x ! eq c cosmos ray 100 do 110 i=1,m 110 d2(i)=d1(i) x=xmedian(d2,m) if(x.gt.sigma)then ! star if(d2(1).lt.s61*x)goto 200 do 120 i=1,m 120 if(d1(i).eq.d2(1))d1(i)=x else ! background do 130 i=1,m 130 if( d1(i).eq.d2(1) .and. d1(i).gt.s(i) )d1(i)=x endif c combine 200 x=0. do 210 i=1,m y= d1(i)*scale(i)+white(i) if(y.ge.29000.0)then x=s50 goto 220 endif c if after 20030205's data, y*=8.0/3.4 c 20030808 y=y*g(i)/gain 210 x=x+y 220 allc=x end subroutine readrec(k,l,b) real b(2048*32) read(k,rec=l)b end subroutine badfile(a,nb,lb,mb,bx,by,badtype,ix,iy,k) real a(2048,2048) integer*2 bx(1),by(1),badtype,ix,iy if(badtype.eq.0)then do 10 i=1,nb 10 a(bx(i),by(i))=-999. endif if(badtype.eq.1)then do 20 i=10001,lb 20 a(bx(i),by(i))=-999. endif if(badtype.eq.5)then do 30 i=20001,mb 30 a(bx(i),by(i))=-999. endif i=iy*2048+ix i=iy*2048+ix call shift(a,i) open(52,file=char(k/10+64)//char(mod(k,10)+64)//'.tmp', * status='unknown',access='direct',recl=2048*2048*4) write(52,rec=1)a close(52) end subroutine shift(a,k) real a(2048*2048) mm=2048*2048 if(k.gt.0)then j=k do 10 i=1,mm-k j=j+1 10 a(i)=a(j) endif if(k.lt.0)then j=mm+k+1 do 20 i=mm,1-k,-1 j=j-1 20 a(i)=a(j) endif end function xmedian(a,n) real a(n) do 10 i=1,n-1 do 10 j=i+1,n if(a(i).gt.a(j))goto 10 x=a(i) a(i)=a(j) a(j)=x 10 continue m=n/2 xmedian=a(m+1) if(m+m.ne.n)return xmedian=(a(m)+a(m+1))/2 end subroutine star100(a,ix,iy,white,px,py,pz,n) integer*2 ix,iy,px(1),py(1) real a(2048,2048),pz(1) do 10 i=1,n z=0. do 20 ii=px(i)+ix-7,px(i)+ix+7 do 20 jj=py(i)+iy-7,py(i)+iy+7 20 z=z+a(ii,jj)-white 10 pz(i)=z end subroutine readata(a,f1) integer c(720),a(1) character*80 f1,ch*2880 equivalence (c,ch) nn=720 open(1,file=f1,status='old',access='direct',recl=nn*4) read(1,rec=1)ch n1=2048 n2=2048 n=1 55 k=index(ch,'END ') if(k.eq.0)then n=n+1 read(1,rec=n)ch goto 55 endif i=0 k1=n+1 k2=n+(n1*n2-1)/nn do 60 k=k1,k2 if(mod(k,200).eq.0)call dispdot() read(1,rec=k)ch do 60 j=1,nn i=i+1 60 a(i)=c(j) close (1) call swap4(a,n1*n2*4) end subroutine shrink4(a,b) real a(1) integer*2 b(1) do 10 i=1,2048*2048 x=a(i) if(x.gt.32767.)x=32767. if(x.lt.0.)x=0. 10 b(i)=x j=2048*2048 do 20 i=1,2048*2048 20 b(j+i)=b(i) end subroutine getstar(a,b,starx,stary,starv,nstar, c white,sx,sy,m,g1,gain) integer*2 a(2048,2048),b(1) real starx(1),stary(1),starv(1) integer*2 sx(1),sy(1) ix1=1 ix2=-1 iy1=1 iy2=-1 do 5 i=1,m if(sx(i).lt.ix1)ix1=sx(i) if(sx(i).gt.ix2)ix2=sx(i) if(sy(i).lt.iy1)iy1=sy(i) 5 if(sy(i).gt.iy2)iy2=sy(i) ix1=-ix1 ix2=2048-ix2 iy1=-iy1 iy2=2048-iy2 ix1=ix1+20 ix2=ix2-20 iy1=iy1+20 iy2=iy2-20 c 20030808 kt=2000*gain/g1 kk=0 ktw=kt+white do 10 i=ix1,ix2 do 10 j=iy1,iy2 if(a(i,j).gt.ktw)then call center9(i,j,a,ktw,ix,iy,iv,n9) if(iv.lt.ktw*2)goto 10 x=ix y=iy call star(b,2048,2048,x,y,4,ierr,white) if(ierr.ne.0)x=-99. z=(ix-x)**2+(iy-y)**2 if(z.gt.5.)goto 10 if(kk.eq.nstar)goto 11 kk=kk+1 starx(kk)=x stary(kk)=y starv(kk)=iv if(mod(kk,10).eq.1)call dispdot() endif 10 continue 11 nstar=kk do 20 i=1,nstar-1 do 20 j=i+1,nstar if(starv(i).ge.starv(j))goto 20 x=starx(i) starx(i)=starx(j) starx(j)=x x=stary(i) stary(i)=stary(j) stary(j)=x x=starv(i) starv(i)=starv(j) starv(j)=x 20 continue if(nstar.gt.300)nstar=300 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 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 xy_rad(x,y,a8,alpha,delta) c 2048*2048 ccd's X Y (in pixel) -------------> alpha delta (in h,d) real x,y,alpha,delta real*8 a8(8) ! from fits_head "A81"--"A88" pi=4.*atan(1.) xi=a8(1)*x+a8(3)*y+a8(5) xn=a8(2)*x+a8(4)*y+a8(6) rc=a8(7) ! center of ccd (in rad.) dc=a8(8) cd=cos(dc) td=tan(dc) alpha=atan(xi/cd/(1.-xn*td)) delta=atan((xn+td)*cos(alpha)/(1.-xn*td)) alpha=alpha+rc if(abs(dc-delta).gt.2.)then delta=-delta alpha=alpha+pi endif alpha=alpha*12./pi delta=delta*180./pi end function indexpos(head,f1) character*80 head(72),f1*8 do 10 indexpos=1,72 10 if(head(indexpos)(1:8).eq.f1)return end subroutine star(map,n1,n2,xx,yy,ir,ierr,sky) c when ir=4, 15*15 for center integer*2 map(n1,n2),maps(-25:25,-25:25) real b(-25:25),xmap(-25:25,-25:25) c sky first n25=25 x=xx y=yy kx=nint(x) ky=nint(y) if(kx.le.n25)kx=n25+1 if(kx.gt.n1-n25)kx=n1-n25 if(ky.le.n25)ky=n25+1 if(ky.gt.n2-n25)ky=n2-n25 do 10 j=-n25,n25 do 10 i=-n25,n25 10 maps(i,j)=map(kx+i,ky+j) n51=n25+n25+1 ierr=1 ! error =1 i5=ir*1.1+2 do 1000 ll=1,7 if(x.lt.1..or.x.gt.n1)return if(y.lt.1..or.y.gt.n2)return ix=nint(x) iy=nint(y) cccccc x center weight_center method do 30 i=ix-i5,ix+i5 z=0. do 40 j=iy-i5,iy+i5 40 z=z+map(i,j)-sky 30 b(i-ix)=z ! for weight_center cal. z=0. s=0. do 50 i=-i5,i5 z=z+b(i) 50 s=s+b(i)*i s=s/z x=ix+s cccccc y center do 130 j=iy-i5,iy+i5 z=0. do 140 i=ix-i5,ix+i5 140 z=z+map(i,j)-sky 130 b(j-iy)=z z=0. s=0. do 150 i=-i5,i5 z=z+b(i) 150 s=s+b(i)*i s=s/z y=iy+s 1000 continue ix=nint(x) iy=nint(y) do 101 j=-n25,n25 do 101 i=-n25,n25 101 xmap(i,j)=map(ix+i,iy+j) x=ix-25+center_1(xmap,sky,51,1) y=iy-25+center_1(xmap,sky,51,0) if(x.lt.1..or.x.gt.n1)return if(y.lt.1..or.y.gt.n2)return ierr=0 xx=x yy=y return end function center_1(a,sky,n,ip) real a(51,51) integer b(51) skyn=sky*n if(ip.eq.1)then do 10 i=1,n x=0. do 20 j=1,n 20 x=x+a(i,j) 10 b(i)=x-skyn else do 110 i=1,n x=0. do 120 j=1,n 120 x=x+a(j,i) 110 b(i)=x-skyn endif call histat(b,mode,maxh,xmean,xpeak,fsigma,gsigma,n) center_1=xpeak end