c c this differs from _i4 because it c properly accounts for zero c #define _LINUX_ .true. #define _PXDIMX_ 4096 #define _PXDIMY_ 4096 subroutine WFCFLTREAD(FILE,pix) implicit none real pix(4096,4096) character*80 FILE call readfits_WFC3(FILE,pix(0001,0001),1) call readfits_WFC3(FILE,pix(0001,2049),4) return end subroutine readfits_WFC3(FILE,pix,nimg) implicit none character*80 FILE real pix(4096,2051) character*70 INFO(10) common / fitsinfo / INFO integer nimg integer naxes integer laxis(3) character*8 field character*20 stream integer nbyte0 integer nbyteE integer nbyte1 integer nbyte2 integer nbper integer i,ios, k character*2880 buffc byte buffb(2880) equivalence (buffc,buffb) integer ifirst, i1, i2 integer np1, np2, npt integer nextend integer nread real bscale, bzero integer bitpix logical LINUX data LINUX/.true./ character*70 HDR(25) common/HDR/HDR open(10,file=FILE,status='old', . err=900,recl=2880,form='UNFORMATTED', . access='DIRECT') bscale = 1.0 bzero = 0.0 naxes = -1 laxis(1) = 1 laxis(2) = 1 laxis(3) = 1 nextend = 0 do i = 1, 10 INFO(i) = ' ' enddo do i = 1, 25 HDR(i) = ' ' enddo i = 0 nread = 0 100 continue i = i + 1 read(10,rec=i,iostat=ios) buffc do k = 0, 35, 1 field = buffc(k*80+01:k*80+08) stream = buffc(k*80+11:k*80+31) if (field.eq.'NAXIS ') read(stream,*) naxes if (field.eq.'NAXIS1 ') read(stream,*) laxis(1) if (field.eq.'NAXIS2 ') read(stream,*) laxis(2) if (field.eq.'NAXIS3 ') read(stream,*) laxis(3) if (field.eq.'NEXTEND ') read(stream,*) nextend if (field.eq.'BITPIX ') read(stream,*) bitpix if (field.eq.'BSCALE ') read(stream,*) bscale if (field.eq.'BZERO ') read(stream,*) bzero if (field.eq.'EXPTIME ') INFO(01) = stream if (field.eq.'FILTER ') INFO(02) = stream if (field.eq.'FILENAME') INFO(03) = stream if (field.eq.'DATE-OBS') INFO(04) = stream if (field.eq.'TIME-OBS') INFO(05) = stream if (field.eq.'DEC_TARG') INFO(06) = stream if (field.eq.'RA_TARG ') INFO(07) = stream if (field.eq.'PA_V3 ') INFO(08) = stream if (field.eq.'PROPOSID') INFO(09) = stream if (field.eq.'CCDGAIN ') INFO(10) = stream if (field.eq.'CRPIX1 ') HDR(01) = stream if (field.eq.'CRPIX2 ') HDR(02) = stream if (field.eq.'CRVAL1 ') HDR(03) = stream if (field.eq.'CRVAL2 ') HDR(04) = stream if (field.eq.'CTYPE1 ') HDR(05) = stream if (field.eq.'CTYPE2 ') HDR(06) = stream if (field.eq.'CD1_1 ') HDR(07) = stream if (field.eq.'CD1_2 ') HDR(08) = stream if (field.eq.'CD2_1 ') HDR(09) = stream if (field.eq.'CD2_2 ') HDR(10) = stream if (field.eq.'ORIENTAT') HDR(11) = stream if (field.eq.'PA_APER ') HDR(12) = stream if (field.eq.'PA_V3 ') HDR(13) = stream if (field.eq.'DATE-OBS') HDR(14) = stream if (field.eq.'TIME-OBS') HDR(15) = stream if (field.eq.'EXPTIME ') HDR(16) = stream if (field.eq.'ROOTNAME') HDR(17) = stream if (field.eq.'TARGNAME') HDR(18) = stream if (field.eq.'RA_TARG ') HDR(19) = stream if (field.eq.'DEC_TARG') HDR(20) = stream if (field.eq.'PROPOSID') HDR(21) = stream if (field.eq.'FILTER ') HDR(22) = stream if (field.eq.'FILTER ') HDR(23) = stream if (field.eq.'VAFACTOR') HDR(24) = stream if (field.eq.'CCDGAIN ') HDR(25) = stream if (field.eq.'END ') goto 101 enddo goto 100 101 continue nread = nread + 1 ifirst = i+1 i1 = i i2 = i if (naxes.eq.0) then ! maybe multiple images stored as extensions if (nextend.eq.0) then print*,'THIS IS A NULL IMAGE: ' print*,'NAXES: ',NAXES print*,'NEXND: ',NEXTEND stop endif endif if (nread.ne.nimg+1) then nbper = abs(BITPIX/8)*laxis(1)*laxis(2) if (NAXES.eq.0) nbper = 0 i = i + 1.0*nbper/2880 + 0.9999 goto 100 endif if (laxis(1).ne.4096.or. . laxis(2).ne.2051) then print*,' laxis1: ',laxis(1) print*,' laxis2: ',laxis(2) print*,' 4096: ',4096 print*,' 2051: ',2051 stop endif if (naxes.eq.2) then ! nimg is irrelevant; ignore nbper = 4*laxis(1)*laxis(2) npt = laxis(1)*laxis(2) nbyte1 = 1 nbyte2 = nbper i1 = i+1 + nbyte1/2880 i2 = i+1 + nbyte2/2880 endif if (naxes.eq.3) then nbper = 4*laxis(1)*laxis(2) nbyte1 = 1 + nbper*(nimg-1) nbyte2 = nbper*(nimg ) i1 = i+1 + nbyte1/2880 i2 = i+1 + nbyte2/2880 npt = laxis(1)*laxis(2) endif do i = i1, i2, 1 read(10,rec=i,iostat=ios) buffc nbyte0 = (i-ifirst)*2880+ 1 nbyteE = (i-ifirst)*2880+2880 np1 = (nbyte0-nbyte1)/4 + 1 np2 = (nbyteE-nbyte1)/4 + 1 call buff2pix_r4(buffb,pix,np1,npt) enddo return 900 continue print*,'readfits_WFC: READFITS ERROR' stop end subroutine writfits_i2_32K(FILE,pix,PXDIMX,PXDIMY) implicit none integer PXDIMX,PXDIMY character*80 FILE integer*2 pix(PXDIMX,PXDIMY) integer nbyte0 integer nbyteE integer nbyte1 integer nbyte2 integer nbper integer i,ios character*2880 buffc byte buffb(2880) equivalence (buffb,buffc) integer ifirst, i1, i2 integer np1, np2, npt character*70 HDR(25) common/HDR/HDR character*80 FILEU FILEU = FILE do i = 75,2,-1 if (FILE(i:i+4).eq.'.fits') FILEU = FILE(1:i+4) enddo write(*,'(''ENTER writfits_i2_32K: '',80a)') FILEU open(10,file=FILEU, . status='unknown', . err =900, . recl =2880, . form ='UNFORMATTED', . access='DIRECT') i = 1 write(buffc( 0*80+1: 1*80),'(''SIMPLE = T '')') write(buffc( 1*80+1: 2*80),'(''BITPIX = 16 '')') write(buffc( 2*80+1: 3*80),'(''NAXIS = '',i12)') 2 write(buffc( 3*80+1: 4*80),'(''NAXIS1 = '',i12)') PXDIMX write(buffc( 4*80+1: 5*80),'(''NAXIS2 = '',i12)') PXDIMY write(buffc( 5*80+1: 6*80),'(''DATATYPE= '',9a)') . ' ''INTEGER*2 ''' write(buffc( 6*80+1: 7*80),'(''DATE = '',11a)') . ' ''28/01/00''' write(buffc( 7*80+1: 8*80),'(''BSCALE = '',i12)') 00001 write(buffc( 8*80+1: 9*80),'(''BZERO = '',i12)') 32000 write(buffc(09*80+1:10*80),'(''CRPIX1 = '',a20)') HDR(01) write(buffc(10*80+1:11*80),'(''CRPIX2 = '',a20)') HDR(02) write(buffc(11*80+1:12*80),'(''CRVAL1 = '',a20)') HDR(03) write(buffc(12*80+1:13*80),'(''CRVAL2 = '',a20)') HDR(04) write(buffc(13*80+1:14*80),'(''CTYPE1 = '',a20)') HDR(05) write(buffc(14*80+1:15*80),'(''CTYPE2 = '',a20)') HDR(06) write(buffc(15*80+1:16*80),'(''CD1_1 = '',a20)') HDR(07) write(buffc(16*80+1:17*80),'(''CD1_2 = '',a20)') HDR(08) write(buffc(17*80+1:18*80),'(''CD2_1 = '',a20)') HDR(09) write(buffc(18*80+1:19*80),'(''CD2_2 = '',a20)') HDR(10) write(buffc(19*80+1:20*80),'(''ORIENTAT= '',a20)') HDR(11) write(buffc(20*80+1:21*80),'(''PA_APER = '',a20)') HDR(12) write(buffc(21*80+1:22*80),'(''PA_V3 = '',a20)') HDR(13) write(buffc(22*80+1:23*80),'(''DATE-OBS= '',a20)') HDR(14) write(buffc(23*80+1:24*80),'(''TIME-OBS= '',a20)') HDR(15) write(buffc(24*80+1:25*80),'(''EXPTIME = '',a20)') HDR(16) write(buffc(25*80+1:26*80),'(''ROOTNAME= '',a20)') HDR(17) write(buffc(26*80+1:27*80),'(''TARGNAME= '',a20)') HDR(18) write(buffc(27*80+1:28*80),'(''RA_TARG = '',a20)') HDR(19) write(buffc(28*80+1:29*80),'(''DEC_TARG= '',a20)') HDR(20) write(buffc(29*80+1:30*80),'(''PROPOSID= '',a20)') HDR(21) write(buffc(30*80+1:31*80),'(''FILTER = '',a20)') HDR(22) write(buffc(31*80+1:32*80),'(''FILTER = '',a20)') HDR(23) write(buffc(32*80+1:33*80),'(''VAFACTOR= '',a20)') HDR(24) write(buffc(33*80+1:34*80),'(''COMMENT '')') write(buffc(34*80+1:35*80),'(''COMMENT '')') write(buffc(35*80+1:36*80),'(''END '')') write(10,rec=i,iostat=ios) buffc ifirst = i+1 i1 = i i2 = i nbper = 2*PXDIMX*PXDIMY npt = PXDIMX*PXDIMY nbyte1 = 1 nbyte2 = nbper i1 = i+1 + nbyte1/2880 i2 = i+1 + nbyte2/2880 do i = i1, i2, 1 nbyte0 = (i-ifirst)*2880+ 1 nbyteE = (i-ifirst)*2880+2880 np1 = (nbyte0-nbyte1)/4 + 1 np2 = (nbyteE-nbyte1)/4 + 1 np1 = (nbyte0-nbyte1)/2 + 1 np2 = (nbyteE-nbyte1)/2 + 1 !if (i.lt.0010) print*,'i: ',i,i1,i2,np1,np2,nbyte0 !if (i.gt.2900) print*,'i: ',i,i1,i2,np1,np2,nbyte0 call pix2buff_i2(buffb,pix,np1,npt) write(10,rec=i,iostat=ios) buffc enddo close(10) return 900 continue print*,'WRITFITS.f ERROR' stop end subroutine buff2pix_r4(buff,pix,n1,nt) implicit none byte buff(2880) real pix(*) integer n1,nt byte b(4) real r equivalence(r,b) integer i, npu, nbu do i = 1, 720 npu = n1+i-1 nbu = (i-1)*4 if (.not.(_LINUX_)) then b(1) = buff(nbu+1) b(2) = buff(nbu+2) b(3) = buff(nbu+3) b(4) = buff(nbu+4) endif if ((_LINUX_)) then b(4) = buff(nbu+1) b(3) = buff(nbu+2) b(2) = buff(nbu+3) b(1) = buff(nbu+4) endif if (npu.ge.1.and.npu.le.nt) pix(npu) = r enddo return end subroutine pix2buff_r4(buff,pix,n1,nt) implicit none byte buff(2880) real*4 pix(*) integer n1,nt byte b(4) real*4 r equivalence(r,b) integer i, npu, nbu do i = 1, 720 npu = n1+i-1 nbu = (i-1)*4 if (npu.ge.1.and.npu.le.nt) r = pix(npu) if (.not.(_LINUX_)) then buff(nbu+1) = b(1) buff(nbu+2) = b(2) buff(nbu+3) = b(3) buff(nbu+4) = b(4) endif if ((_LINUX_)) then buff(nbu+1) = b(4) buff(nbu+2) = b(3) buff(nbu+3) = b(2) buff(nbu+4) = b(1) endif enddo return end subroutine buff2pix_i4(buff,pix,n1,nt) implicit none byte buff(2880) integer*4 pix(*) integer n1,nt byte b(4) integer ii equivalence(ii,b) integer i, npu, nbu do i = 1, 720 npu = n1+i-1 nbu = (i-1)*4 if (.not.(_LINUX_)) then b(1) = buff(nbu+1) b(2) = buff(nbu+2) b(3) = buff(nbu+3) b(4) = buff(nbu+4) endif if ((_LINUX_)) then b(4) = buff(nbu+1) b(3) = buff(nbu+2) b(2) = buff(nbu+3) b(1) = buff(nbu+4) endif if (npu.ge.1.and.npu.le.nt) pix(npu) = ii enddo return end subroutine buff2pix_i2(buff,pix,n1,nt) implicit none byte buff(2880) integer*2 pix(*) integer n1,nt byte b(2) integer ii equivalence(ii,b) integer i, npu, nbu do i = 1, 1440 npu = n1+i-1 nbu = (i-1)*2 if (.not.(_LINUX_)) then b(1) = buff(nbu+1) b(2) = buff(nbu+2) endif if ((_LINUX_)) then b(2) = buff(nbu+1) b(1) = buff(nbu+2) endif if (npu.ge.1.and.npu.le.nt) pix(npu) = ii enddo return end subroutine pix2buff_i2(buff,pix,n1,nt) implicit none byte buff(2880) integer*2 pix(*) integer n1,nt byte b(2) integer*2 ii equivalence(ii,b) integer i, npu, nbu do i = 1, 1440 npu = n1+i-1 nbu = (i-1)*2 if (npu.ge.1.and.npu.le.nt) ii = pix(npu) if (.not.(_LINUX_)) then buff(nbu+1) = b(1) buff(nbu+2) = b(2) endif if ((_LINUX_)) then buff(nbu+1) = b(2) buff(nbu+2) = b(1) endif enddo return end subroutine pix2buff_i4(buff,pix,n1,nt) implicit none byte buff(2880) integer*4 pix(*) integer n1,nt byte b(4) integer ii equivalence(ii,b) integer i, npu, nbu do i = 1, 720 npu = n1+i-1 nbu = (i-1)*4 if (npu.ge.1.and.npu.le.nt) ii = pix(npu) if (.not.(_LINUX_)) then buff(nbu+1) = b(1) buff(nbu+2) = b(2) buff(nbu+3) = b(3) buff(nbu+4) = b(4) endif if ((_LINUX_)) then buff(nbu+1) = b(4) buff(nbu+2) = b(3) buff(nbu+3) = b(2) buff(nbu+4) = b(1) endif enddo return end program convert_WFC implicit none real pixr(4096,4096) real pixh(4096,4096) integer*2 pixi(4096,4096) integer*4 pixii integer i integer j integer NARGs, NARG, iargc character*80 FILEI character*80 FILEO character*70 HDR(25) common/HDR/HDR data HDR / . ' ', . ' ', . ' ', . ' ', . ' ', . ' ', . ' ', . ' ', . ' ', . ' ', . ' ', . ' ', . ' ', . ' ', . ' ', . ' ', . ' ', . ' ', . ' ', . ' ', . ' ', . ' ', . ' ', . ' ', . ' '/ NARGs = iargc() 144 format(1x,i4,1x,i4,1x,i7,1x,i7) c open(44,file='convert_WFC.SATN',status='unknown') c open(45,file='convert_WFC.USAT',status='unknown') do NARG = 1, NARGs call getarg(NARG,FILEI) write(*,'(''ENTER WFCFLTREAD: '',a18)') FILEI call WFCFLTREAD(FILEI,pixr) FILEO = FILEI(01:09) // '_WJ2.fits' do i = 0001, 4096 pixr(i,2046) = -750.0 pixr(i,2047) = -750.0 pixr(i,2048) = -750.0 pixr(i,2049) = -750.0 pixr(i,2050) = -750.0 pixr(i,2051) = -750.0 enddo c do i = 0001, 4096 c do j = 0001, 4096 c if (pixr(i,j).gt.40000) write(99,'(i4.4,1x,i4.4,1x,f8.1)') c . i,j,pixr(i,j) c enddo c enddo c c if (.true.) stop call sat_fix(pixr) do i = 1, 4096 do j = 1, 4096 pixii = int(1000+pixr(i,j)+0.5) pixii = pixii - 1000 if (pixii.gt.55000) then pixii = 55000 + (pixii-55000+2)/5 endif pixii = pixii - 32000 if (pixii.lt.-32750) pixii = -32750 if (pixii.gt. 32750) pixii = 32750 pixi(i,j) = pixii enddo enddo write(*,'(''ENTER writfits_i2_32K: '',a18)') FILEO call writfits_i2_32K(FILEO,pixi,4096,4096) print*,' ' write(*,110) (i,i=0888-5,0888+5) do j = 0963+10,0963-10,-1 write(*,111) j,(pixi(i,j)+32000,i=0888-5,0888+5) enddo print*,' ' write(*,110) (i,i=0888-5,0888+5) print*,' ' do j = 0963+10,0963-10,-1 write(*,111) j,(int(pixr(i,j)+0.5),i=0888-5,0888+5) enddo print*,' ' write(*,110) (i,i=0888-5,0888+5) print*,' ' print*,' ' 110 format(4x,1x,11i6) 111 format(i4,1x,11i6) enddo stop end subroutine sat_fix(pixr) implicit none real pixr(4096,4096) real pixs(4096,4096) real pixx(4096,4096) integer*2 pixi(4096,4096) integer*4 pixii byte pixu(4096,4096) integer i integer j integer NARGs, NARG, iargc integer L, LL, Ls integer iil(9999) integer jjl(9999) real frl(9999) real rrl(9999) integer nrl(9999) real srl(9999) logical trip character*80 FILEI character*80 FILEO real UPR_LIM integer iu real pixe integer il(99999), jl(99999) integer NSATs integer NCENs real sr, mbar_sky real fr, fru, fr0 real apphot_sat real apphot_satX integer nr, nr0 integer ir real reff integer ireff integer ii, jj integer ncirc NARGs = iargc() 144 format(1x,i4,1x,i4,1x,i7,1x,i7) print*,' ' print*,'ASSUME GAIN OF TWO.' print*,' ' UPR_LIM = 62500.0 do i = 0001, 4096 pixr(i,2048) = -750 pixr(i,2049) = -750 enddo do i = 0001, 4096 do j = 0001, 4096 pixu(i,j) = 0 pixs(i,j) = pixr(i,j) pixx(i,j) = pixr(i,j) enddo enddo NSATs = 0 NCENs = 0 do i = 1, 4096 if (i.eq.i/128*128) write(*,'('' i: '',i4.4,1x,2i9)') . i,NSATs,NCENs do j = 1, 4096 if (pixr(i,j).gt.UPR_LIM.and.pixu(i,j).eq.0) then call find_sat(i,j,pixr,pixu,UPR_LIM,il,jl,Ls) endif if (pixr(i,j).ge.99998) NSATs = NSATs + 1 if (pixr(i,j).eq.99999) NCENs = NCENs + 1 enddo enddo do i = 0001, 4096 do j = 0001, 4096 if (pixr(i,j).ge.99998) pixx(i,j) = -750 enddo enddo Ls = 0 do i = 0001+10, 4096-10 do j = 0001+10, 4096-10 if (pixr(i,j).eq.99999) then sr = mbar_sky(i,j,15,20,pixs) fr = apphot_sat(i,j,pixs,sr,nr) Ls = Ls + 1 iil(Ls) = i jjl(Ls) = j frl(Ls) = fr nrl(Ls) = nr srl(Ls) = sr reff = sqrt(fr/90000/3.14159) rrl(Ls) = reff endif enddo enddo call cosort_riiirr(frl,iil,jjl,nrl,srl,rrl,Ls) do L = 1, Ls write( *,888) iil(L),jjl(L), . -2.5*log10(max(frl(L),1)), . srl(L),frl(L),nrl(L),L c write(98,888) iil(L),jjl(L), c . -2.5*log10(max(frl(L),1)), c . srl(L),frl(L),nrl(L),L 888 format(i4.4,1x,i4.4,1x,f10.4,1x,f8.2, . 1x,f12.1,1x,i8,1x,i5.5) enddo do L = 1, Ls i = iil(L) j = jjl(L) fr = frl(L) reff = sqrt(fr/90000/3.14159) ireff = reff + 0.90 if (ireff.lt.1) ireff = 1 ncirc = 0 do ii = -ireff, +ireff do jj = -ireff, +ireff if (ii**2+jj**2.le.(ireff+0.5)**2) . ncirc = ncirc + 1 enddo enddo fru = fr-100000 do ii = -ireff-1, +ireff+1 do jj = -ireff-1, +ireff+1 if (ii**2+jj**2.le.(ireff+1.5)**2) . pixx(i+ii,j+jj) = -750 enddo enddo do LL = L+1, Ls if ((iil(L)-iil(LL))**2+ . (jjl(L)-jjl(LL))**2.lt. . (rrl(L)+rrl(LL)+3)**2) goto 777 enddo if (abs(j-2048.5).lt.reff+5) goto 777 if (j+0001.lt.reff+5) goto 777 if (4095-j.lt.reff+5) goto 777 do ii = -ireff, +ireff do jj = -ireff, +ireff if (ii**2+jj**2.le.(ireff+0.5)**2) . pixx(i+ii,j+jj) = fru/(ncirc-1) enddo enddo pixx(i,j) = 100000 777 continue enddo if (.false.) then open(55,file='LOG.peaksat_CEN.reg',status='unknown') write(55,'(''# Region file format: DS9 version 3.0'')') write(55,'(''global color=yellow'')') open(56,file='LOG.peaksat_ALL.reg',status='unknown') write(56,'(''# Region file format: DS9 version 3.0'')') write(56,'(''global color=yellow'')') do i = 1, 4096 do j = 1, 4096 if (pixr(i,j).eq.99999) write(55,112) i,j,1.0,'red' if (pixr(i,j).ge.99998) write(56,112) i,j,0.5,'yellow' enddo enddo close(55) close(56) 112 format('image;circle( ',i4,',',i4,',',f6.3,') # color=',a7) endif do i = 1, 4096 do j = 1, 4096 pixr(i,j) = pixx(i,j) enddo enddo c if (.true.) then c print*,'---> re-anal: ' c do i = 0001+10, 4096-10 c do j = 0001+10, 4096-10 c if (pixx(i,j).gt.95000) then c ir = 1 c 3 ir = ir + 1 c fr = 0 c trip = .false. c do ii = -ir, ir c do jj = -ir, ir c if (ii**2+jj**2.le.(ir+0.5)**2) then c if (pixx(i+ii,j+jj).gt.0) c . fr = fr + pixx(i+ii,j+jj) c if (pixx(i+ii,j+jj).lt.0) trip = .true. c endif c enddo c enddo c print*,i,j,ir,fr,trip c if (.not.trip) goto 3 c write(99,120) i,j,-2.5*log10(max(fr,1)) c 120 format(i4.4,1x,i4.4,1x,f10.4) c endif c enddo c enddo c endif return end c---------------------------------------------------- c c c c subroutine procsatn(pixarr,UPR_LIM) implicit none real pixarr(_PXDIMX_,_PXDIMY_) real UPR_LIM integer i, j real rpsf_photijk real dx, dy real ri real pixval,pixtot real psfval,psftot real zlist(99999) integer nl, nt integer ii, jj integer ir real reff integer iru real fru integer nsat c c this is the radial profile of the WFC F675W PSF... c c from r = 0 to r = 50 _flt pixels... c real rprof(51) data rprof / . 0.27000000,0.07940824,0.02189085,0.00431118,0.00219986, . 0.00184045,0.00074173,0.00039569,0.00032109,0.00023005, . 0.00017415,0.00016298,0.00013543,0.00010287,0.00008278, . 0.00007402,0.00005836,0.00004262,0.00003057,0.00002460, . 0.00002157,0.00001737,0.00001614,0.00001478,0.00001383, . 0.00001180,0.00001164,0.00001120,0.00000997,0.00000976, . 0.00000882,0.00000834,0.00000775,0.00000671,0.00000619, . 0.00000547,0.00000536,0.00000503,0.00000473,0.00000469, . 0.00000431,0.00000448,0.00000425,0.00000403,0.00000381, . 0.00000336,0.00000298,0.00000237,0.00000193,0.00000173, . 0.00000040 / write(*,'(''# '')') write(*,'(''# ENTER PROCSATN: '')') write(*,'(''# UPR_LIM: '',f8.2)') UPR_LIM write(*,'(''# '')') nsat = 0 do j = _PXDIMY_-2, 3, -1 if (j.eq.j/064*064) .write(*,'(''# '',i4.4,'' out of 4096. NSAT: '',i6)') j,nsat do i = 3, _PXDIMX_-2 if (pixarr(i,j).gt.UPR_LIM) then nsat = nsat + 1 ir = 1 5 continue if (ir.lt.10) ir = ir + 1 nl = 0 nt = 0 psftot = 0. pixtot = 0. do ii = i-ir, i+ir do jj = j-ir, j+ir ri = sqrt(1.*(ii-i)**2+(jj-j)**2) if (ii.ge.1+3.and.ii.le.(_PXDIMX_-3).and. . jj.ge.1+3.and.jj.le.(_PXDIMY_-3).and. . ri.lt.ir) then nt = nt + 1 if (pixarr(ii,jj).gt.UPR_LIM) goto 6 if (pixarr(ii+1,jj).gt.UPR_LIM.and. . pixarr(ii+2,jj).gt.UPR_LIM.and. . pixarr(ii+3,jj).gt.UPR_LIM) goto 6 if (pixarr(ii-1,jj).gt.UPR_LIM.and. . pixarr(ii-2,jj).gt.UPR_LIM.and. . pixarr(ii-3,jj).gt.UPR_LIM) goto 6 if (pixarr(ii,jj+1).gt.UPR_LIM.and. . pixarr(ii,jj+2).gt.UPR_LIM.and. . pixarr(ii,jj+3).gt.UPR_LIM) goto 6 if (pixarr(ii,jj-1).gt.UPR_LIM.and. . pixarr(ii,jj-2).gt.UPR_LIM.and. . pixarr(ii,jj-3).gt.UPR_LIM) goto 6 nl = nl + 1 if (nl.gt.99999) then print*,' nl wants to be >99 999...' stop endif pixval = pixarr(ii,jj) reff = ri if (reff.gt.3.5) reff = 3.5 + (reff-3.5)*0.5 iru = int(ri) fru = ri-int(ri) psfval = rprof(iru+1)+ . (fru)*(rprof(iru+2)-rprof(iru+1)) pixtot = pixtot + pixval psftot = psftot + psfval zlist(nl) = pixval/psfval 6 continue endif enddo enddo if (nl.lt.nt/2.and.ir.lt.10) goto 5 pixarr(i,j) = pixtot/psftot*rprof(1) if (psftot.le.0) pixarr(i,j) = UPR_LIM if (.not.(pixarr(i,j).gt.UPR_LIM)) pixarr(i,j) = UPR_LIM endif enddo!i enddo!j return end subroutine find_sat(i0,j0,pixr,pixu,UPR_LIM,il,jl,Ls) implicit none integer i0 integer j0 real pixr(4096,4096) byte pixu(4096,4096) real UPR_LIM integer il(99999) integer jl(99999) integer ml(99999) integer Ls integer Lo, l integer jmin, jmax integer ii,jj integer iii,jjj integer oi, oj, m integer imax, mmax Ls = 0 if (pixr(i0,j0).lt.UPR_LIM) return Ls = 1 il(1) = i0 jl(1) = j0 pixu(i0,j0) = 1 jmin = 0001 jmax = 2048 if (j0.ge.2049) then jmin = 2049 jmax = 4096 endif 1 continue Lo = Ls do l = 1, Lo do ii = max(il(L)-1,0001), min(il(L)+1,4096) do jj = max(jl(L)-1,jmin), min(jl(L)+1,jmax) if (pixu(ii,jj).eq.0) then if (pixr(ii,jj).gt.55000) then pixu(ii,jj) = 1 Ls = Ls + 1 if (Ls.gt.99999) then print*,' Ls wants to be > 999999... ' endif il(Ls) = ii jl(Ls) = jj m = 0 3 continue m = m + 1 iii = ii + oi(m) jjj = jj + oj(m) if (iii.ge.0001.and.iii.le.4096.and. . jjj.ge.jmin.and.jjj.le.jmax) then if (pixr(iii,jjj).gt.55000) goto 3 endif ml(Ls) = m endif endif enddo enddo enddo if (Ls.gt.Lo) goto 1 imax = il(1) jmax = jl(1) mmax = ml(1) do L = 1, Ls if (ml(L).gt.mmax) then imax = il(L) jmax = jl(L) mmax = ml(L) endif enddo do L = 1, Ls pixr(il(L),jl(L)) = 99998 enddo pixr(imax,jmax) = 99999 return end c-------------------------------------------------- c c This routine orders the pixels by successive distance c from each other. If you step through the array c like this: c c c do m = 1, 99 c i = ixc + oi(m) c j = iyc + oj(m) c ... c enddo c c c Then, the pixels (i,j) will be in order of increasing c distance from the central pixel (ixc,iyc). c c This can be useful for many applications. c c integer function oi(m) implicit none integer m integer qm, fm integer im integer mm integer oim, ojm integer nm integer oil(91) data oil / 0,1,1,2,2,1,2,3,3,1,3,2,4,4, . 1,3,4,2,3,5,4,5,1,5,2,4, . 3,5,6,1,6,6,2,5,4,3,6,7,5,1, . 7,6,4,2,7,7,3,6,5,8,1,4, . 8,7,8,2,6,8,3,7,5,4,8,9,1,9, . 9,7,2,6,5,8,3,9,4,9,7, . 9,5,8,6,8,7,8,9,6,9,7,9,8,9 / integer ojl(91) data ojl / 0,0,1,0,1,2,2,0,1,3,2,3,0,1,4, . 3,2,4,4,0,3,1,5,2,5,4, . 5,3,0,6,1,2,6,4,5,6,3,0,5,7,1, . 4,6,7,2,3,7,5,6,0,8,7, . 1,4,2,8,6,3,8,5,7,8,4,0,9,1,2, . 6,9,7,8,5,9,3,9,4,7, . 5,9,6,8,7,8,8,6,9,7,9,8,9,9 / qm = m+2 - 4*int((m+2)/4) fm = int((m+6)/4) if (fm.le.91) then if (qm.eq.0) oi = oil(fm) if (qm.eq.1) oi = -oil(fm) if (qm.eq.2) oi = -ojl(fm) if (qm.eq.3) oi = ojl(fm) return endif im = (sqrt(1.*m-1)-1)/2 mm = (2*im+1)**2+1 nm = (m-mm)/4 oim = nm-im ojm = im if (qm.eq.0) oi = oim if (qm.eq.1) oi = -oim if (qm.eq.2) oi = -ojm if (qm.eq.3) oi = ojm return end c---------------------------------------------------------------- c c the y-analog to oi(m) c integer function oj(m) implicit none integer m integer qm, fm integer im, mm, nm, oim, ojm integer oil(91) data oil / 0,1,1,2,2,1,2,3,3,1,3,2,4,4, . 1,3,4,2,3,5,4,5,1,5,2,4, . 3,5,6,1,6,6,2,5,4,3,6,7,5,1, . 7,6,4,2,7,7,3,6,5,8,1,4, . 8,7,8,2,6,8,3,7,5,4,8,9,1,9, . 9,7,2,6,5,8,3,9,4,9,7, . 9,5,8,6,8,7,8,9,6,9,7,9,8,9 / integer ojl(91) data ojl / 0,0,1,0,1,2,2,0,1,3,2,3,0,1,4, . 3,2,4,4,0,3,1,5,2,5,4, . 5,3,0,6,1,2,6,4,5,6,3,0,5,7,1, . 4,6,7,2,3,7,5,6,0,8,7, . 1,4,2,8,6,3,8,5,7,8,4,0,9,1,2, . 6,9,7,8,5,9,3,9,4,7, . 5,9,6,8,7,8,8,6,9,7,9,8,9,9 / qm = m+2 - 4*int((m+2)/4) fm = int((m+6)/4) if (fm.le.91) then if (qm.eq.0) oj = ojl(fm) if (qm.eq.1) oj = -ojl(fm) if (qm.eq.2) oj = oil(fm) if (qm.eq.3) oj = -oil(fm) return endif im = (sqrt(1.*m-1)-1)/2 mm = (2*im+1)**2+1 nm = (m-mm)/4 oim = nm-im ojm = im if (qm.eq.0) oj = ojm if (qm.eq.1) oj = -ojm if (qm.eq.2) oj = oim if (qm.eq.3) oj = -oim return end c---------------------------------------------------------------- c c should be obvious c real function orm(m) implicit none integer m integer oi, oj orm = sqrt(1.*(oi(m)**2 + oj(m)**2)) return end c--------------------------------------------------- c c This routine will find an iteratively clipped mean sky c value from an annulis about a point. The point is (ixc,iyc). c It will take the pixels between r=inr and r=ior and will c analyze them for a robust average. c real function mbar_sky(ixc,iyc,inr,ior,pixarr) implicit none integer ixc, iyc integer inr, ior real pixarr(_PXDIMX_,_PXDIMY_) integer i , j integer ixh, iyh real rij integer nuse logical inimage real sklist(99999) integer nsk real bar, sig integer HIFLAG common / HIFLAG_ / HIFLAG integer LOFLAG common / LOFLAG_ / LOFLAG nsk = 0 do i = -ior+1, ior-1 do j = -ior+1, ior-1 rij = i**2+j**2 if (rij.ge.ior**2) goto 333 if (rij.lt.inr**2) goto 333 ixh = ixc + i iyh = iyc + j if (.not.inimage(ixh,iyh)) goto 333 if (pixarr(ixh,iyh).le.LOFLAG) goto 333 if (pixarr(ixh,iyh).ge.HIFLAG) goto 333 nsk = nsk + 1 if (nsk.gt.99999) goto 333 sklist(nsk) = pixarr(ixh,iyh) 333 continue enddo enddo call rbarsigs(sklist,nsk,bar,sig,nuse,1.75) mbar_sky = bar return end c---------------------------------------- c c Returns "true" if a given pixels is within c the confines of the image. c logical function inimage(i,j) implicit none integer i, j inimage = .true. if (i.lt. 1) inimage = .false. if (i.gt._PXDIMX_) inimage = .false. if (j.lt. 1) inimage = .false. if (j.gt._PXDIMY_) inimage = .false. return end subroutine rbarsigs(xlist,NTOT,bar,sig,NUSE,SIGCLIP) implicit none integer NTOT real xlist(NTOT) real bar real sig integer NUSE real SIGCLIP integer n real*8 bsum, ssum integer nsum, nsumo integer NIT nsum = 0 bar = 0.e0 sig = 9e9 do NIT = 1, 20 nsumo = nsum bsum = 0. ssum = 0. nsum = 0. do n = 1, NTOT if (abs(xlist(n)-bar).le.SIGCLIP*sig) then bsum = bsum + xlist(n) ssum = ssum + abs(xlist(n)-bar) nsum = nsum + 1 endif enddo if (nsum.gt.0) bar = bsum / nsum if (nsum.gt.1) sig = ssum/(nsum-1) if (nsum.lt.0.35*NTOT.and.NIT.ge.3) return if (nsum.eq.nsumo.and.NIT.ge.3) goto 1 enddo 1 continue NUSE = nsum if (nsum.le.1) sig = 0.999 return end subroutine zero_fixq(pix) implicit none real pix(4096,4096) integer hist1(25), hist2(25), hist3(25) integer i, j integer PM2, PM1, P00, PP1, PP2 integer EM2, EM1, E00, EP1, EP2 integer ERR real f real A, B, C real fmin, emin real rand print*,' ' print*,' ' print*,'ZERO_FIX: (for antiquated way of storing images)' print*,' ' print*,' ' call find_hist(pix,hist1) do i = 0001, 4096 do j = 0001, 4096 !write(*,'(2i5,1x,f10.3)') i,j,pix(i,j) if (pix(i,j).lt.-0.9) pix(i,j) = pix(i,j)-1 enddo enddo call find_hist(pix,hist2) c print*,' ' c do h = 01, 25 c hz = h-13 c write(*,'(i3,1x,i3,1x,2i10)') h,hz,hist1(h),hist2(h) c enddo c print*,' ' fmin = 0 emin = 9e9 c print*,' ' do f = 0.00, 1.00, 0.01 PM2 = hist2(11) PM1 = hist2(13)*( f ) P00 = hist2(13)*(1-f) PP1 = hist2(14) PP2 = hist2(15) A = P00 B = ((PP1-PM1)/2.0 + (PP2-PM2)/4.0)/2.0 C = (PP2+PM2-2*A)/8.0 c----------------------------------------------------- c A = -0.08571*PM2 c . +0.34286*PM1 c . +0.48572*P00 c . +0.34286*PP1 c . -0.08571*PP2 c B = -0.20000*PM2 c . -0.10000*PM1 c . +0.00000*P00 c . +0.10000*PP1 c . +0.20000*PP2 c C = +0.71430*PM2 c . -0.35710*PM1 c . -0.71430*P00 c . -0.35710*PP1 c . +0.71430*PP2 c print*,' ' c write(*,'(30x,3f15.2)') A,B,C c A = P00 c B = ((PP1-PM1)/2.0 + (PP2-PM2)/4.0)/2.0 c C = (PP2+PM2-2*A)/8.0 c write(*,'(30x,3f15.2)') A,B,C c print*,' ' c----------------------------------------------------- EM2 = A - 2*B + 4*C EM1 = A - B + C E00 = A EP1 = A + B + C EP2 = A + 2*B + 4*C ERR = abs(PM2-EM2) . + abs(PM1-EM1) . + abs(P00-E00) . + abs(PP1-EP1) . + abs(PP2-EP2) c write(*,'(f8.4,10x,5i8,10x,5i8,10x,i12)') f, c . PM2,PM1,P00,PP1,PP2, c . EM2,EM1,E00,EP1,EP2, c . ERR if (ERR.lt.emin) then emin = ERR fmin = f endif enddo c print*,' ' c print*,'fmin: ',fmin c print*,' ' do i = 0001, 4096 do j = 0001, 4096 if (int(pix(i,j)).eq.0) then if (rand().lt.fmin) pix(i,j) = -1 endif enddo enddo call find_hist(pix,hist3) c print*,' ' c do h = 01, 25 c hz = h-13 c write(*,'(i3,1x,i3,1x,3i10)') h,hz,hist1(h),hist2(h),hist3(h) c enddo c print*,' ' return end c----------------------------------- c c subroutine find_hist(pix,hist) implicit none real pix(4096,4096) integer hist(25) integer h integer i, j do h = 01, 25 hist(h) = 0 enddo do i = 0001, 4096 do j = 0001, 4096 h = 13 + pix(i,j) if (h.ge.01.and.h.le.25) hist(h) = hist(h)+1 enddo enddo return end c---------------------------------------------- c c this routine will search for a field in a fits header c subroutine query_hdr(filename,FIELDX,streamx) implicit none character filename*80 character*8 field character*20 stream character*8 fieldx character*20 streamx integer i integer ios, k character*2880 buff integer nread c----------------------------------------------- open(10,file=FILENAME,status='old', . err=900,recl=2880,form='UNFORMATTED', . access='DIRECT') streamx = '0' i = 0 nread = 0 100 continue i = i + 1 read(10,rec=i,iostat=ios) buff if (ios.lt.0) goto 900 do k = 0, 35, 1 field = buff(k*80+01:k*80+08) stream = buff(k*80+11:k*80+31) if (field.eq.fieldx) streamx = stream(1:20) if (field.eq.'END ') goto 101 109 continue enddo goto 100 101 continue close(10) return 900 continue print*,' ' print*,'imginfo.e ERROR EXIT. ' print*,' ' print*,'ONE OF THE IMAGES WAS NOT IN STANDARD' print*,'HST FITS FORMAT.' print*,' ' write(*,'(''PROBLEM FILE: '',a80)') FILENAME print*,' ' stop end c real function apphot_satX(ii,jj,pixc,sr,nsat) c implicit none c integer ii,jj c real pixc(4096,4096) c real sr c integer nsat c c integer u0(101,101) c integer u1(101,101) c integer u2(101,101) c integer nadd, naddo c integer iu,ju c integer iii,jjj c integer nadd2 c real ftot, f5 c integer i,j c integer ir c c nsat = 0 c ftot = 0 c do i = -5, 5 c do j = -5, 5 c if (i**2+j**2.lt.5.5**2) then c ftot = ftot + (pixc(ii+i,jj+j)-sr) c if (pixc(ii+i,jj+j).gt.70000) then c nsat = nsat + 1 c endif c endif c enddo c enddo c apphot_satX = ftot c if (nsat.lt.15) return c c c do i = 001, 101 c do j = 001, 101 c u1(i,j) = 0 c u2(i,j) = 0 c enddo c enddo c c do i = -50, 50 c do j = -50, 50 c u0(51+i,51+j) = 0 c if (i+ii.gt.0004.and.i+ii.lt.4093.and. c . j+jj.gt.0004.and.j+jj.lt.4093) then c if (pixc(ii+i,jj+j).gt.50000) u0(51+i,51+j) = 1 c endif c enddo c enddo c c do i = -2, 2 c do j = -2, 2 c u1(51+i,51+j) = 1 c enddo c enddo c c c ir = 2 c 3 continue c ir = ir + 1 c nadd = 0 c do i = -ir, ir c do j = -ir, ir c if (i**2+j**2.le.(ir+0.5)**2.and. c . u0(51+i,51+j).eq.1.and. c . u1(51+i,51+j).eq.0) then c if (u1(51+i-1,51+j )+ c . u1(51+i+1,51+j )+ c . u1(51+i ,51+j+1)+ c . u1(51+i ,51+j-1)+ c . u1(51+i-1,51+j-1)+ c . u1(51+i+1,51+j-1)+ c . u1(51+i-1,51+j+1)+ c . u1(51+i+1,51+j+1).ge.1) then c u1(51+i,51+j) = 1 c nadd = nadd + 1 c endif c endif c enddo c enddo c if (nadd.gt.0.and.ir.le.48) goto 3 c c nadd2 = 0 c do i = 002, 100 c do j = 002, 100 c u2(i,j) = 0 c if (u1(i ,j ).eq.1.or. c . u1(i+1,j ).eq.1.or. c . u1(i-1,j ).eq.1.or. c . u1(i ,j+1).eq.1.or. c . u1(i ,j-1).eq.1.or. c . u1(i+1,j+1).eq.1.or. c . u1(i-1,j+1).eq.1.or. c . u1(i+1,j-1).eq.1.or. c . u1(i-1,j-1).eq.1) u2(i,j) = 1 c if ((i-51)**2+(j-51)**2.le.5.5**2) u2(i,j) = 1 c nadd2 = nadd2 + u2(i,j) c enddo c enddo c c f5 = 0 c ftot = 0 c nsat = 0 c do i = 001, 101 c do j = 001, 101 c if (u2(i,j).eq.1.and. c . i+ii.gt.0001.and.i+ii.lt.4096.and. c . j+jj.gt.0001.and.j+jj.lt.4096) then c if (pixc(ii+(i-51),jj+(j-51))-sr.gt.0) c . ftot = ftot + (pixc(ii+(i-51),jj+(j-51))-sr) c if (pixc(ii+(i-51),jj+(j-51)).gt.70000) nsat = nsat + 1 c endif c enddo c enddo c c apphot_satX = ftot c c return c end c c c c-------------------------------------------------------- c c sorts in ASCENDING ORDER!!! c subroutine cosort_riiirr(r1,i2,i3,i4,r5,r6,NTOT) implicit none real r1(1) integer i2(1), i3(1), i4(1) real r5(1), r6(1) integer NTOT integer n logical change 777 continue change = .false. do n = 1, NTOT-1 if (r1(n).gt.r1(n+1)) then call swap(r1(n),r1(n+1)) call swip(i2(n),i2(n+1)) call swip(i3(n),i3(n+1)) call swip(i4(n),i4(n+1)) call swap(r5(n),r5(n+1)) call swap(r6(n),r6(n+1)) change = .true. endif enddo if (change) goto 777 end subroutine swap(x,y) implicit none real x, y, t t = x x = y y = t return end subroutine swip(x,y) implicit none integer x, y, t t = x x = y y = t return end real function apphot_sat(ii,jj,pixc,sr,nsat) implicit none integer ii,jj real pixc(4096,4096) real sr integer nsat integer u1(101,101) integer u2(101,101) integer nadd, naddo integer iu,ju integer iii,jjj integer nadd2 real ftot, f5 integer i,j nsat = 0 ftot = 0 do i = -5, 5 do j = -5, 5 if (i**2+j**2.lt.5.5**2) then ftot = ftot + (pixc(ii+i,jj+j)-sr) if (pixc(ii+i,jj+j).gt.70000) nsat = nsat + 1 endif enddo enddo apphot_sat = ftot if (nsat.le.3) return do i = 001, 101 do j = 001, 101 u1(i,j) = 0 u2(i,j) = 0 enddo enddo do i = -2, 2 do j = -2, 2 if (pixc(ii+i,jj+j).gt.50000) u1(51+i,51+j) = 1 enddo enddo naddo = 0 nadd = 0 1 continue naddo = nadd do i = -49, 49 do j = -49, 49 if (i+ii.le.0002) goto 3 if (j+jj.le.0002) goto 3 if (i+ii.ge.4095) goto 3 if (j+jj.ge.4095) goto 3 iu = max(min(i+ii,4095),2) ju = max(min(j+jj,4095),2) if (u1(51+i,51+j).eq.1) then do iii = -1,+1 do jjj = -1,+1 if (pixc(iu+iii,ju+jjj).gt.50000) u1(51+i+iii,51+j+jjj)=1 enddo enddo endif 3 continue enddo enddo nadd = 0 do i = 001, 101 do j = 001, 101 nadd = nadd + u1(i,j) enddo enddo if (nadd.gt.naddo) goto 1 nadd2 = 0 do i = 002, 100 do j = 002, 100 u2(i,j) = 0 if (u1(i ,j ).eq.1.or. . u1(i+1,j ).eq.1.or. . u1(i-1,j ).eq.1.or. . u1(i ,j+1).eq.1.or. . u1(i ,j-1).eq.1.or. . u1(i+1,j+1).eq.1.or. . u1(i-1,j+1).eq.1.or. . u1(i+1,j-1).eq.1.or. . u1(i-1,j-1).eq.1) u2(i,j) = 1 if ((i-51)**2+(j-51)**2.le.5.5**2) u2(i,j) = 1 nadd2 = nadd2 + u2(i,j) enddo enddo f5 = 0 ftot = 0 nsat = 0 do i = 001, 101 do j = 001, 101 if ((i-51)**2+(j-51)**2.lt.5.5**2) then f5 = f5 + (pixc(ii+(i-51),jj+(j-51))-sr) endif if (u2(i,j).eq.1.and. . ii+(i-51).gt.0002.and. . jj+(j-51).gt.0002.and. . ii+(i-51).lt.4095.and. . jj+(j-51).lt.4095) then ftot = ftot + (pixc(ii+(i-51),jj+(j-51))-sr) if (pixc(ii+(i-51),jj+(j-51)).gt.70000) nsat = nsat+1 endif enddo enddo apphot_sat = ftot nsat = nadd2 return end