C Program to take aperture results from individual frames and average C them into single files for each color parameter (maxfilt = 4, maximages=30) real data(15,maximages), tot(maxfilt), sig(maxfilt), data2(15,maximages) real tot2(maxfilt), sig2(maxfilt) character file*64,line*132 character*5 filt(maxfilt) real apcor(maxfilt), expos(maximages) integer i0(maxfilt), i1(maxfilt) character*24 field, suffix real y(maximages,maxfilt), e(maximages,maxfilt) real y2(maximages,maxfilt), e2(maximages,maxfilt) integer ndat(maximages), ndat2(maximages) data apcor/4*0./ C Sgr 1 WFPC2 c Data apcor/0.15,0.245,0.285,0.275, C & 0.20,0.35,0.405,0.380/ C LMC WFPC2 C data apcor/0.15,0.32,0.37,0.34, C & 0.18,0.40,0.42,0.41/ field = 'lmc' print '(1x,''Enter field name: ''$)' read '(a)', field lf = index(field,' ') - 1 print '(1x,''Enter number of filters: ''$)' read *, nfilt do i=1,nfilt print 100, i 100 format(1x,'Enter filter name for filter: ',i3) read '(a)', filt(i) end do print '(1x,''Enter nimages, starting image in each filter: ''$)' read *, nimages, (i0(i), i=1,nfilt), (i1(i),i=1,nfilt) print '(1x,''Enter exp/exp0 for each image: ''$)' read *, (expos(i),i=1,nimages) print '(1x,''Enter fake number: ''$)' read *, ifake print '(1x,''Enter 1 for CTE/NOCTE frames: ''$)' read *, icte if (icte .eq. 1) then ntype = 4 else ntype = 2 end if do itype = 1,ntype if (itype .eq. 1) then suffix = 'ap' else if (itype .eq. 2) then suffix = 'appsf' else if (itype .eq. 3) then suffix = 'apnocte' else if (itype .eq. 4) then suffix = 'appsfnocte' end if ls = index(suffix,' ') - 1 C Open up the output files for the average aperture results do ifilt=1,nfilt l = numchar(filt(ifilt)) write(file,201) field(1:lf),filt(ifilt)(1:l),suffix(1:ls) 201 format(a,a,a,'.nst') open(6+ifilt,file=file,status='unknown') end do C Loop over each of the individual images do i=1,nimages C Open up the aperture results for this file write(file,101) field(1:lf), i, suffix(1:ls) 101 format(a,'_',i2.2,a,'.nst') open(i+10,file=file,status='old') C Read header lines and write them out if this is the first file for C a new filter read(i+10,'(a)') line do ifilt=1,nfilt if (i .eq. i0(ifilt)) write(6+ifilt,'(a)') line end do read(i+10,'(a)') line do ifilt=1,nfilt if (i .eq. i0(ifilt)) write(6+ifilt,'(a)') line end do read(i+10,'(a)') line do ifilt=1,nfilt if (i .eq. i0(ifilt)) write(6+ifilt,'(a)') line end do C Do the same stuff for the ap0 file, which will hold results for the C larger photometric calibration aperture if (itype .lt. 3) then write(file,103) field(1:lf), i 103 format(a,'_',i2.2,'ap0.nst') else write(file,104) field(1:lf), i 104 format(a,'_',i2.2,'ap0nocte.nst') end if open(i+30,file=file,status='old') read(i+30,'(a)') line read(i+30,'(a)') line read(i+30,'(a)') line end do 1 continue C Initialize accumulators for each filter do i=1,nfilt tot(i) = 0. sig(i) = 0. tot2(i) = 0. sig2(i) = 0. ndat(i) = 0 ndat2(i) = 0 end do do j=1,nimages read(j+10,*,end=99) (data(i,j),i=1,9) read(j+30,*,end=99) (data2(i,j),i=1,9) data(5,j) = max(0.0005,data(5,j)) ii=0 do jfilt=1,nfilt do jj=i0(jfilt),i1(jfilt) ii=ii+1 if (ii .eq. j) ifilt = jfilt end do C if (j .ge. i0(jfilt) .and. j .le. i1(jfilt)) ifilt=jfilt end do C Load up arrays with data values C Add in the result for this filter, this is for no rejection if (data(4,j) .lt. 90) then tot(ifilt) = tot(ifilt) + data(4,j)/data(5,j)**2 sig(ifilt) = sig(ifilt) + 1/data(5,j)**2 ndat(ifilt) = ndat(ifilt) + 1 y(ndat(ifilt),ifilt) = data(4,j) e(ndat(ifilt),ifilt) = data(5,j) C weight by exposure times, not observed errors e(ndat(ifilt),ifilt) = 1./sqrt(expos(j)) end if if (data2(4,j) .lt. 90) then tot2(ifilt) = tot2(ifilt) + data2(4,j)/data2(5,j)**2 sig2(ifilt) = sig2(ifilt) + 1/data2(5,j)**2 ndat2(ifilt) = ndat2(ifilt) + 1 y2(ndat2(ifilt),ifilt) = data2(4,j) e2(ndat2(ifilt),ifilt) = data2(5,j) C weight by exposure times, not observed errors e2(ndat2(ifilt),ifilt) = 1./sqrt(expos(j)) end if end do C Get the weighted means for each filter for both aperture sizes C with no rejection. This is used just as a check and also to compute C estimated errors in the means. do ifilt=1,nfilt if (sig(ifilt) .gt. 0) then data(14,i0(ifilt)) = tot(ifilt)/sig(ifilt) data(15,i0(ifilt)) = sqrt(1./sig(ifilt)) else C print *, data(1,i0(ifilt)),tot(ifilt),sig(ifilt) end if if (sig2(ifilt) .gt. 0) then data2(14,i0(ifilt)) = tot2(ifilt)/sig2(ifilt) data2(15,i0(ifilt)) = sqrt(1./sig2(ifilt)) else C print *, data2(1,i0(ifilt)),tot2(ifilt),sig2(ifilt) end if end do C Now do it properly with rejection do ifilt = 1, nfilt call getmean(y(1,ifilt),e(1,ifilt),ndat(ifilt), & data(4,i0(ifilt)),data(5,i0(ifilt))) call getmean(y2(1,ifilt),e2(1,ifilt),ndat2(ifilt), & data2(4,i0(ifilt)),data2(5,i0(ifilt))) C since we weight by exposure times, the errors are incorrect, so C use the errors as computed by the error-weighted computation data(5,i0(ifilt)) = data(15,i0(ifilt)) data2(5,i0(ifilt)) = data2(15,i0(ifilt)) if (abs(data(4,i0(ifilt)) - data(14,i0(ifilt))) .gt. 0.1) then print 707, nint(data(1,i0(ifilt))), ifilt, & data(4,i0(ifilt)), data(15,i0(ifilt)), & data(14,i0(ifilt)), data(15,i0(ifilt)), & data(4,i0(ifilt)) - data(14,i0(ifilt)), & data(5,i0(ifilt)) - data(15,i0(ifilt)) print 708, (y(i,ifilt),e(i,ifilt),i=1,ndat(ifilt)) end if C if (abs(data2(4,i0(ifilt)) - data2(4,i0(ifilt)+1)) .gt. 0.1) then C print 707, nint(data(1,i0(ifilt))), ifilt, C & data2(4,i0(ifilt)), data2(5,i0(ifilt)), C & data2(4,i0(ifilt)+1), data2(5,i0(ifilt)+1), C & data2(4,i0(ifilt)) - data2(4,i0(ifilt)+1), C & data2(5,i0(ifilt)) - data2(5,i0(ifilt)+1) C print 708, (y2(i,ifilt),e2(i,ifilt),i=1,ndat2(ifilt)) C end if 707 format(i5,i3,6f10.3) 708 format(8f10.3) end do do ifilt=1,nfilt C Compute the aperture correction here for the aperture photometry and C put it into column 8 data(8,i0(ifilt)) = data(4,i0(ifilt)) - data2(4,i0(ifilt)) C Compute the aperture correction here for the profile photometry and C put it into column 9 data(9,i0(ifilt)) = data(9,i0(ifilt)) - data2(4,i0(ifilt)) C Add an aperture correction if we have one data(4,i0(ifilt)) = data(4,i0(ifilt)) - apcor(ifilt) C Write out the results j = i0(ifilt) write(6+ifilt,102) nint(data(1,j)),(data(i,j),i=2,9) 102 format(i6,2f9.2,3f9.3,f9.0,f9.2,f9.3) end do goto 1 99 continue do i=1,nfilt close(i+6) end do do i=1,nimages close(i+10) close(i+30) end do end do stop end C subroutine to get the mean value and estimated error from a set of C measurement with individual errors, rejecting outliers in the process subroutine getmean(y,e,n,mean,sig) real mean, fmed, diff(100), y(n), e(n) integer ind(100) parameter (eps=1.e-4) if (n .gt. 100) then pause 'too many points for internal array in getmean' end if if (n .lt. 1) then mean = 99.999 sig = 9.999 return end if C accumulate the sums without any rejection tot = 0. sig = 0. do i=1,n tot = tot + y(i)/e(i)**2 sig = sig + 1./e(i)**2 end do C compute the median value and sort the points with the largest deviant C first call median(y,n,fmed) do i=1,n diff(i) = abs(fmed - y(i)) end do call quick(diff,n,ind) C Now compare value for each with mean and errors from others starting C with the largest deviant do j=n,1,-1 i = ind(j) C mean will contain the mean from all the other points which haven't been C rejected s = sig - 1./e(i)**2 if (s .gt. eps) then mean = (tot - y(i)/e(i)**2) / s if (abs( y(i)-mean ) / e(i) .gt. 5.) then C this is a bad point. Remove it from the sums tot = tot - y(i)/e(i)**2 sig = sig - 1./e(i)**2 end if end if end do C now compute the means after bad values have been removed from the sum if (sig .gt. eps) then mean = tot/sig sig = sqrt(1./sig) else mean = 99.999 sig = 9.999 end if return end SUBROUTINE MEDIAN(X,NPTS,FMEDIA) C This routine finds the median in an array of numbers. C It uses a quicksort algorithm from the text : C 'Algorithms + Data Structures = Programs' by Wirth (1976) page 84. C In the worst case it completely sorts the input array. C (Time is proportional to n * log n.) C Sorting is done on the indices of the array, so as not to C disturb its original ordering. C Input: X The array of points to be examined C NPTS The length of the array C KEY The keys to be partially sorted C Output: FMEDIA The median value C Authors:Mike Fich, Berkeley, April 1981. C B.F. REAL*4 X(NPTS) INTEGER KEY(1000) C Initialize the keys. DO 8701 I = 1, NPTS KEY(I) = I 8701 CONTINUE ILEFT = 1 IRITE = NPTS MIDDLE = (NPTS+1)/2 C Find the middle element, sort everything by exchanging everything C bigger on the left with everything smaller on the right end of the C array. 8702 IF (ILEFT .LT. IRITE) THEN XMID = X(KEY(MIDDLE)) C INDEX1 moves from left to right, INDEX2 moves from right to left C INDEX1 points to something bigger than the key being sorted on. C INDEX2 points to something smaller. INDEX1 = ILEFT INDEX2 = IRITE 8703 IF (INDEX1 .LE. INDEX2) THEN IF (X(KEY(INDEX1)) .LT. XMID) THEN 8704 IF (X(KEY(INDEX1)) .LT. XMID) THEN INDEX1 = INDEX1+1 GO TO 8704 END IF END IF IF (X(KEY(INDEX2)) .GT. XMID) THEN 8705 IF (X(KEY(INDEX2)) .GT. XMID) THEN INDEX2 = INDEX2-1 GO TO 8705 END IF END IF C If INDEX1 is still to the left of INDEX2 then exchange the C keys found out of place and increment pointers. IF (INDEX1 .LE. INDEX2) THEN KTEMP = KEY(INDEX1) KEY(INDEX1) = KEY(INDEX2) KEY(INDEX2) = KTEMP INDEX1 = INDEX1+1 INDEX2 = INDEX2-1 END IF C If not sorted for initial key go back. GO TO 8703 END IF C Choose new limits for sorting remaining block. IF (INDEX2 .LT. MIDDLE) ILEFT = INDEX1 IF (MIDDLE .LT. INDEX1) IRITE = INDEX2 GO TO 8702 END IF FMEDIA = X(KEY(MIDDLE)) RETURN END C======================================================================= SUBROUTINE QUICK (DATUM, N, INDEX) C======================================================================= C C A quick-sorting algorithm suggested by the discussion on pages 114-119 C of THE ART OF COMPUTER PROGRAMMING, Vol. 3, SORTING AND SEARCHING, by C D.E. Knuth, which was referenced in Don Wells' subroutine QUIK. This C is my own attempt at encoding a quicksort-- PBS. C C Arguments C C DATUM (INPUT/OUTPUT) is a vector of dimension N containing randomly C ordered real data upon input. Upon output the elements of C DATUM will be in order of increasing value. C C C INDEX (OUTPUT) is an integer vector of dimension N. Upon return to C the calling program the i-th element of INDEX will tell where C the i-th element of the sorted vector DATUM had been BEFORE C DATUM was sorted. C C======================================================================= C PARAMETER (MAXSTAK=14) C C Parameter C C MAXSTAK is the maximum number of entries the stack can contain. C A limiting stack length of 14 restricts this quicksort C subroutine to vectors of maximum length of order 32,768 C (= 2**15). REAL*4 DATUM(N) INTEGER*4 INDEX(N), STKLO(MAXSTAK), STKHI(MAXSTAK), HI IF (N .LE. 0) THEN PRINT *, 'NO ELEMENTS IN ARRAY TO BE SORTED ' RETURN ENDIF C C Initialize INDEX. C DO 50 I=1,N 50 INDEX(I)=I C C Initialize the pointers. C NSTAK=0 LIMLO=1 LIMHI=N C 100 DKEY=DATUM(LIMLO) IKEY=INDEX(LIMLO) C PRINT *, LIMLO, LIMHI, dkey, ikey C C Compare all elements in the sub-vector between LIMLO and LIMHI with C the current key datum. C LO=LIMLO HI=LIMHI 101 CONTINUE C IF (LO .EQ. HI)GO TO 200 C IF (DATUM(HI) .LE. DKEY) GO TO 109 HI=HI-1 C C The pointer HI is to be left pointing at a datum SMALLER than the C key, which is intended to be overwritten. C GO TO 101 C 109 DATUM(LO)=DATUM(HI) INDEX(LO)=INDEX(HI) LO=LO+1 110 CONTINUE C IF (LO .EQ. HI) GO TO 200 C IF (DATUM(LO) .GE. DKEY) GO TO 119 C LO=LO+1 GO TO 110 C 119 DATUM(HI)=DATUM(LO) INDEX(HI)=INDEX(LO) HI=HI-1 C C The pointer LO is to be left pointing at a datum LARGER than the C key, which is intended to be overwritten. C GO TO 101 C 200 CONTINUE C C LO and HI are equal, and point at a value which is intended to C be overwritten. Since all values below this point are less than C the key and all values above this point are greater than the key, C this is where we stick the key back into the vector. C DATUM(LO)=DKEY INDEX(LO)=IKEY CD DO 1666 I=LIMLO,LO-1 CD1666 PRINT *, DATUM(I) CD PRINT *, DATUM(L0), ' KEY' CD DO 2666 I=LO+1,LIMHI CD2666 PRINT *, DATUM(I) C C At this point in the subroutine, all data between LIMLO and LO-1, C inclusive, are less than DATUM(LO), and all data between LO+1 and C LIMHI are larger than DATUM(LO). C C If both subarrays contain no more than one element, then take the most C recent interval from the stack (if the stack is empty, we're done). C If the larger of the two subarrays contains more than one element, and C if the shorter subarray contains one or no elements, then forget the C shorter one and reduce the other subarray. If the shorter subarray C contains two or more elements, then place the larger subarray on the C stack and process the subarray. C IF (LIMHI-LO .GT. LO-LIMLO) GO TO 300 C C Case 1: the lower subarray is longer. If it contains one or no C elements then take the most recent interval from the stack and go C back and operate on it. C IF (LO-LIMLO .LE. 1) GO TO 400 C C If the upper (shorter) subinterval contains one or no elements, then C process the lower (longer) one, but if the upper subinterval contains C more than one element, then place the lower (longer) subinterval on C the stack and process the upper one. C IF (LIMHI-LO .GE. 2) GO TO 250 C C Case 1a: the upper (shorter) subinterval contains no or one elements, C so we go back and operate on the lower (longer) subinterval. C LIMHI=LO-1 GO TO 100 C 250 CONTINUE C C Case 1b: the upper (shorter) subinterval contains at least two C elements, so we place the lower (longer) subinterval on the stack and C then go back and operate on the upper subinterval. C NSTAK=NSTAK+1 STKLO(NSTAK)=LIMLO STKHI(NSTAK)=LO-1 LIMLO=LO+1 CD DO 3666 I=1,NSTAK CD3666 PRINT *, 'STACK: ', STKLO(I), STKHI(I) GO TO 100 C 300 CONTINUE C C Case 2: the upper subarray is longer. If it contains one or no C elements then take the most recent interval from the stack and C operate on it. C IF (LIMHI-LO .LE. 1) GO TO 400 C C If the lower (shorter) subinterval contains one or no elements, then C process the upper (longer) one, but if the lower subinterval contains C more than one element, then place the upper (longer) subinterval on C the stack and process the lower one. C IF (LO-LIMLO .GE. 2) GO TO 350 C C Case 2a: the lower (shorter) subinterval contains no or one elements, C so we go back and operate on the upper (longer) subinterval. C LIMLO=LO+1 GO TO 100 C 350 CONTINUE C C Case 2b: the lower (shorter) subinterval contains at least two C elements, so we place the upper (longer) subinterval on the stack and C then go back and operate on the lower subinterval. C NSTAK=NSTAK+1 STKLO(NSTAK)=LO+1 STKHI(NSTAK)=LIMHI LIMHI=LO-1 CD DO 4666 I=1,NSTAK CD4666 PRINT *, 'STACK: ', STKLO(I), STKHI(I) GO TO 100 C 400 CONTINUE C C Take the most recent interval from the stack. If the stack happens C to be empty, we are done. C IF (NSTAK .LE. 0) RETURN C ! Normal return LIMLO=STKLO(NSTAK) LIMHI=STKHI(NSTAK) NSTAK=NSTAK-1 GO TO 100 C END