c combine 2048*2048 pic. parameter (nfile=24) parameter (nb=5259) integer*2 sx(nfile),sy(nfile) integer*2 xx1(nfile),xx2(nfile),yy1(nfile),yy2(nfile) character*80 f0,f1,f2,head(72) real*8 a8(8),adcoef(2,3) real a(2048*2048),wx(nfile),wy(nfile) real b(2048*32,nfile) real ab(2048*2432) real ba(2048*2048) equivalence (ab,b),(ab(2048*384+1),a),(ab,ba) if(iargc().lt.2)then write(*,*) write(*,*)' ******** Combine 2048*2048 Pic. ****2000.7' write(*,*) write(*,*) write(*,*)' sh-- mean fragment operate' write(*,*)' cob origion' write(*,*)' cob50 if > 29000, put 500000' write(*,*)' cob51 if >29000, process as bad point' write(*,*)' cobad add frames direct' write(*,*) write(*,*)' cob batch_file out_file' write(*,*) write(*,*)' 1996,12,19' 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.3)stop 'too less files, must>=2' close(51) m=i-1 if(m.gt.nfile)m=nfile open(51,file=f0,status='old') write(*,2) 2 format(35x,'shift(x,y) UT. exp. sky'/78('-')) nt=0 np=0 do 30 i=1,m xx1(i)=1 xx2(i)=2048 yy1(i)=1 yy2(i)=2048 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) wx(i)=x2-sx(i) wy(i)=y2-sy(i) ipos=indexpos(head,'NAXIS ') read(head(ipos)(65:69),'(i5)')n12 if(n12.ne.0)then ipos=indexpos(head,'NAXIS1 ') read(head(ipos)(65:74),'(2i5)')xx1(i),xx2(i) ipos=indexpos(head,'NAXIS2 ') read(head(ipos)(65:74),'(2i5)')yy1(i),yy2(i) else n12=1 endif np=np+n12 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 write(*,1)i,f2(1:22),x2,y2,sx(i),sy(i), * head(ipos1)(16:27),head(ipos2)(24:30),white 1 format(i2,': ',a,2f6.1,2x,2i4,2a,f8.2) 30 continue close(51) ******* get each shift end ****** open(51,file=f0,status='old') do 120 i=1,m read(51,'(a)')f2 call readata(a,f2) call shfile(a,ba,wx(i),wy(i),sx(i),sy(i),i) 120 write(*,*)f2(1:22) close(51) ********************************** do 215 i=1,m 215 open(50+i,file=char(64+i)//'.tmp', * status='old',access='direct',recl=2048*32*4) k=0 do 220 l=1,64 call dispdot() do 230 i=1,m 230 call readrec(50+i,l,b(1,i)) do 240 i=1,2048*32 xad=0. do 235 j=1,m 235 xad=xad+b(i,j) k=k+1 240 a(k)=xad 220 continue do 250 i=1,m 250 close(50+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,iy2i,xx1,xx2,yy1,yy2) 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:69),"('frame',i5)")np 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=char(i+64)//'.tmp' 200 call unlink(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,xx1,xx2,yy1,yy2) real a(2048,2048) integer*2 sx(1),sy(1) integer*2 xx1(1),xx2(1),yy1(1),yy2(1) do 5 i=1,m xx1(i)=xx1(1)-sx(i) xx2(i)=xx2(i)-sx(i) yy1(i)=yy1(i)-sy(i) 5 yy2(i)=yy2(i)-sy(i) ix1=1 ix2=2048 iy1=1 iy2=2048 do 10 i=1,m if(xx1(i).gt.ix1)ix1=xx1(i) if(xx2(i).lt.ix2)ix2=xx2(i) if(yy1(i).gt.iy1)iy1=yy1(i) 10 if(yy2(i).lt.iy2)iy2=yy2(i) 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 subroutine readrec(k,l,b) real b(2048*32) read(k,rec=l)b end subroutine shfile(a,b,wx,wy,ix,iy,k) real a(2048,2048),b(2048,2048) integer*2 ix,iy i=iy*2048+ix open(52,file=char(k+64)//'.tmp',status='unknown', * access='direct',recl=2048*2048*4) if(k.ne.1)then call shift(a,i) call p_shift(a,b,wx,wy) write(52,rec=1)b else write(52,rec=1)a endif 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 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 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 p_shift(a,b,wx,wy) real a(2048,2048),b(2048,2048) i1=1 if(wx.lt.0.)i1=-1 j1=1 if(wy.lt.0.)j1=-1 if(wx.lt.0.)wx=-wx if(wy.lt.0.)wy=-wy zx=1.-wx zy=1.-wy x1=zx*zy x2=wx*zy x3=zx*wy x4=wx*wy do 10 j=1,2048 jj1=j+j1 if(jj1.lt.1)jj1=1 if(jj1.gt.2048)jj1=2048 do 10 i=1,2048 ii1=i+i1 if(ii1.lt.1)ii1=1 if(ii1.gt.2048)ii1=2048 10 b(i,j)=a(i,j)*x1+a(ii1,j)*x2+a(i,jj1)*x3+a(ii1,jj1)*x4 return end