program extspect c Program is a refined version of the program spectext which itself was a c more sophisticated version of bandspectra which originally grew out of c bandsort. The program here: c c - groups individual observations of an object together c - sorts object observations in order of increasing wavelength c - discards objects with dud (negative) fluxes c - combines the SExtractor flags for each observation of an object c into a single number c - searches for and flags calibration stars in the field and uses c their raw fluxes to produce mean normalisation factors at each c wavelength c - assigns fresh ID numbers to each object c - calculates wavelength correction based on distance of object c from optical axis c - chooses fixed aperture size that optimizes S/N for each object c c c Input format is as follows (4 files in total): -------------------------- c c Input file name is 'extspect.in' c c Input format requires the following 16 columns: c slice no., SExt ID no., X, Y, flux(ap1), err(ap1), flux (ap2), err(ap2) c flux(ap3), err(ap3), flux(ap4), err(ap4), ellipticity, FWHM, c star/gal classification, SExt flag c Assumes objects to be matched are sorted into order of increasing X in c X,Y being the first two entries in each column. c c c c Output format is as follows (12 files in total): ------------------------- c c outRaw(60000,16) **N.B: starRaw() contains a larger sample of stars than c the other star files (since it is only selected on stargal). c sl ID X Y F1 dF1 F2 dF2 F3 dF3 F4 dF4 ell FWHM stgl Flg c 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 c c outMean(10000,13) c newID X Y ap wavCor ell FWHM strgl Flags rad numOb avFlx avdFl S/N c 1 2 3 4 5 6 7 8 9 10 11 12 13 14 c c outFlux(10000,(n+3)) c newID X Y F1 F2 F3 ... Fn c 1 2 3 4 5 6 (n+3) c c outUncert(10000,(n+3)) c newID X Y dF1 dF2 dF3 ... dFn c 1 2 3 4 5 6 (n+3) c c outSN(10000,(n+3)) c newID X Y F1/dF1 F2/dF2 F3/dF3 ... Fn/dFn c 1 2 3 4 5 6 (n+3) c c starNorm() and starErrNrm() c newID X Y rad S/N F1/F1 F2/F1 F3/F1 ... Fn/F1 c 1 2 3 4 5 6 7 8 (n+5) c c outID(10000,3) c X Y numObsv c 1 2 3 c c c Identical subset files with only the stars (filename prefix 'star'). c DHJ 26/10/98 c ------------------------------------------------------------------------ c Max number of observations = 100000 (must be > N) c Max number of objects = 20000 c Max number of slices = 12 (must be > Nsl) c Number of observation fields (input) = 16 implicit none character*50 inFile, inWaveFile, inApertFile character*50 currFile, outFilePref, inMatchFile integer col, row, inpRow, N, numObj, numObsv, outRow, objRow integer inpStart, bestCol, Nmatches, nextmatch, Nsl, starRow integer refRow, rowDepth, distflag1, distflag2, maxRowDepth, k integer threshCol, colmax, numStarRaw, numStar, strgalObs integer numLostSN, strApNum, numSpill, obsvReqd, normSlice real input(100000,17), objbuffer(13,17), meanRow(17) real apDiam(5), bestSN, currSN real mat(5000,3), centwaveL(13), fluxTh(5,13) real outRaw(100000,17), outMean(20000,15) real outFlux(20000,(13+3)), outUncert(20000,(13+3)) real starIDnos(4000) real Xprev(10*13), Yprev(10*13), sig(13) real reorder(40000,5*2), rat(13+3), err(13+3), err2(13+3) real dist, maxDist, Xref, Yref, X, Y, Xcen, Ycen, realRow real Dpix, beta, Rpix2, placeholder, lowSN, matchflag real gn, SNthresh, meansig real distXY, distXYrefs, freq, strgalCut real flux, Area, sig2, pie, Nstr, a, b, c, d, err1 real strMaxFlx, strMinFlx, strMaxFWHM, strMinFWHM real fracFluxLost, minLost, maxLost, currLost, S real*8 Rap, Robj c maxRowDepth = 400 c --- maximum search distance allowed (in rows) from reference obsvn c --- (=2*2 pixel distance *6 frames *20 objects in that narrow strip) c k = 3 c --- maximum comparison distance to go when checking current obsvn c --- against those previous (in units of Nsl; i.e. k=3 => dist=3*Nsl) freq = 250.0 c freq = 10.0 c --- frequency of input lines to report on progress through input c --- catalogue (freq=500.0 means progress reported every 500th line) pie = 4*atan(1.0) open(unit=20, file='extspect.in') read(20,*,end=10,err=10) inFile read(20,*,end=10,err=10) inMatchFile read(20,*,end=10,err=10) inWaveFile read(20,*,end=10,err=10) inApertFile read(20,*,end=10,err=10) threshCol read(20,*,end=10,err=10) normSlice read(20,*,end=10,err=10) maxDist read(20,*,end=10,err=10) S read(20,*,end=10,err=10) strgalCut read(20,*,end=10,err=10) strgalObs read(20,*,end=10,err=10) strMaxFlx read(20,*,end=10,err=10) strMinFlx read(20,*,end=10,err=10) strMaxFWHM read(20,*,end=10,err=10) strMinFWHM read(20,*,end=10,err=10) strApNum read(20,*,end=10,err=10) gn read(20,*,end=10,err=10) Dpix read(20,*,end=10,err=10) beta read(20,*,end=10,err=10) Xcen read(20,*,end=10,err=10) Ycen read(20,*,end=10,err=10) outFilePref read(20,*,end=10,err=10) k read(20,*,end=10,err=10) maxRowDepth read(20,*,end=10,err=10) obsvReqd 10 continue close(unit=20) if (k.gt.10) then write (6,*) ' ***********************************************' write (6,*) ' ERROR: k > 10 !!' write (6,*) ' --> setting k=3' write (6,*) ' ***********************************************' k=10 endif if (maxRowDepth.lt.30) then write (6,*) ' ***********************************************' write (6,*) ' ERROR: maxRowDepth < 30 !!' write (6,*) ' --> setting maxRowDepth=200' write (6,*) ' ***********************************************' maxRowDepth=200 endif write (6,*) ' ' write (6,*) 'Input (near-sorted) catalogue: ',inFile write (6,*) 'Optional (sorted) object match file: ',inMatchFile write (6,*) 'Input (sorted) central wavelengths file: ',inWaveFile write (6,*) 'Input aperture sizes/flux thresh file: ',inApertFile write (6,*) 'Flux thresh column (>=2 => mean values): ',threshCol write (6,*) 'Slice number for stellar normalisations: ',normSlice write (6,*) 'Maximum pix-dist for an object match: ',maxDist write (6,*) 'Sigma threshold for variations in seeing: ',S write (6,*) 'Star/gal. cut-off for stellar sample: ',strgalCut write (6,*) 'Maximum no. of missed stellar observns: ',strgalObs write (6,*) 'Maximum stellar flux: ',strMaxFlx write (6,*) 'Minimum stellar flux: ',strMinFlx write (6,*) 'Maximum stellar FWHM: ',strMaxFWHM write (6,*) 'Minimum stellar FWHM: ',strMinFWHM write (6,*) 'Stellar aperture to use (1-4): ',strApNum write (6,*) 'Gain of the CCD (e-/ADU): ',gn write (6,*) 'Side length (pix) of CCD: ',Dpix write (6,*) 'Angle subtended by CCD: ',beta write (6,*) 'X-coord of optical axis centre: ',Xcen write (6,*) 'Y-coord of optical axis centre: ',Ycen write (6,*) 'Output catalogue extension: ',outFilePref write (6,*) ' ' write (6,*) 'k: ',k write (6,*) 'maxRowDepth: ',maxRowDepth write (6,*) 'Min. number of required obsvns: ',obsvReqd c -------------------------------------------------------------------- c Reading input file and filling input array: c -------------------------------------------------------------------- write (6,*) ' ' write (6,*) 'Reading input data...' N=0 open(unit=20, file=inFile, status='old') 15 N=N+1 read(20,*,end=20,err=20) (input(N,col),col=1,16) do col=1,8 reorder(N,col)=input(N,(col+4)) enddo input(N,5)=reorder(N,1) input(N,6)=reorder(N,5) input(N,7)=reorder(N,2) input(N,8)=reorder(N,6) input(N,9)=reorder(N,3) input(N,10)=reorder(N,7) input(N,11)=reorder(N,4) input(N,12)=reorder(N,8) c --- work-around for SExtractor re-ordering of fluxes and uncertainties goto 15 20 continue N=N-1 write (6,*) '(There are ',N,' lines of observations.)' close(unit=20) c -------------------------------------------------------------------- c Reading central wavelength data: c -------------------------------------------------------------------- Nsl=0 open(unit=20, file=inWaveFile, status='old') 25 Nsl=Nsl+1 read(20,*,end=30,err=30) centwaveL(Nsl) goto 25 30 continue Nsl=Nsl-1 write (6,*) '(There are ',Nsl,' slices.)' c -------------------------------------------------------------------- c Reading aperture sizes/noise data: c -------------------------------------------------------------------- row=0 open(unit=20, file=inApertFile, status='old') read(20,*) SNthresh,(sig(col),col=1,Nsl),meansig c --- i.e. reading first line which contains sigma rather than sky values if (threshCol.le.1) then 35 row=row+1 read(20,*,end=40,err=40) apDiam(row),(fluxTh(row,col),col=1,Nsl) goto 35 40 continue write (6,*) 'Flux threshold calculated for individual slices.' end if c --- i.e. if threshold column = 0,1 etc, read threshold values for ALL sl if (threshCol.gt.1) then colmax = threshCol-1 45 row=row+1 read(20,*,end=50,err=50) apDiam(row),(fluxTh(row,col), $ col=1,colmax) goto 45 50 continue write (6,*) 'Mean flux threshold taken for all slices.' end if c --- i.e. if threshold column >= 2, read only single column of threshol c --- (ordinarily some sort of mean threshold varying only with aperture) c --- c --- Mean fluxes in fluxTh(row,(threshCol-1)) colmax=row-1 write (6,*) '(Detection threshold, S/N = ',SNthresh,'.)' write (6,*) '(Mean sky RMS = ',meansig,' ADU rms.)' write (6,*) '(There are ',colmax,' apertures.)' close(unit=20) c -------------------------------------------------------------------- c Reading match coordinates: c -------------------------------------------------------------------- Nmatches=0 open(unit=20, file=inMatchFile, status='old') 55 Nmatches=Nmatches+1 read(20,*,end=60,err=60) mat(Nmatches,1),mat(Nmatches,2) goto 55 60 continue Nmatches=Nmatches-1 write (6,*) '(There are ',Nmatches,' objects to match.)' write (6,*) ' ' close(unit=20) c -------------------------------------------------------------------- c *** Re-calculating photometric uncertainty c from UNSMOOTHED rms data (inApertFile): *** c -------------------------------------------------------------------- do row=1,N do col=1,colmax flux = input(row,2*col+3) Area = pie*((apDiam(col)/2.0)**2.0) sig2 = (sig( int(input(row,1)) )**2.0) input(row,2*col+4) = sqrt(abs((Area*sig2)+(flux/gn))) end do end do c --- N.B: take abs() inside the sqrt() for cases of negative flux c ===================================================================== c Object Loop: inpRow = 1 inpStart = 1 numObj = 0 outRow = 1 numObsv = 0 nextmatch = 1 matchflag = 0 numLostSN = 0 numSpill = 0 write (6,*) ' ------------------------------------------------' c --- Copy X and Y positions and input row number for reference; also c --- initialise previous values of (Xref,Yref). do row=1,Nsl Xprev(row)=0.0 Yprev(row)=0.0 enddo Xref = input(inpRow,3) Yref = input(inpRow,4) refRow = inpRow c --- New object: initialise object buffer array: do row=1,Nsl objbuffer(row,1)=row do col=2,16 objbuffer(row,col)=0.0 enddo enddo c --- Copy first entry to object buffer: objRow = int(input(inpRow,1)) do col=1,16 objbuffer(objRow,col)=input(inpRow,col) enddo inpRow=inpRow+1 numObsv = numObsv+1 c ------------------------------------------------------------- do while (inpRow.le.(N+1)) c cccccccccc write (6,*) 'A: ',inpRow, numObj, inpStart X = input(inpRow,3) Y = input(inpRow,4) rowDepth = inpRow - refRow dist = sqrt((X-Xref)**2+(Y-Yref)**2) c --- i.e. calculating distances between current observation and c --- current (Xref,Yref), as well the distance from the reference row. c --- Checking current observation against previous ones: distflag1 = 0 distflag2 = 0 c --- means that neither the current (X,Y) or (Xref,Yref) are nearby c --- any of the Nsl previous (Xref,Yref) points. row=1 do while (row.le.(k*Nsl)) distXY = sqrt((X-Xprev(row))**2+(Y-Yprev(row))**2) distXYrefs = sqrt((Xref-Xprev(row))**2+(Yref-Yprev(row))**2) if (distXY.lt.maxDist) then distflag1 = 1 endif c --- i.e. current (X,Y) is no good for adding to object buffer since it c --- is nearby one of the Nsl previous (Xref,Yref). if (distXYrefs.lt.maxDist) then distflag2 = 1 endif c --- i.e. current (Xref,Yref) is no good to write to either output file c --- since it is nearby one of the Nsl previous (Xref,Yref). if ((distflag1.eq.1).and.(distflag2.eq.1)) then row = (k*Nsl)+1 endif c --- i.e. break out of loop. row=row+1 enddo c ccccccccc write (6,*) 'B:' if (inpRow.eq.(N+1)) then rowDepth = maxRowDepth+1 end if c --- i.e. special detour straight to writing out object buffer if c --- last time through (final input line) if (distflag1.eq.1) then dist = 1.1*maxDist end if c --- i.e. also special detour straight past distance comparison c --- if current matches previous two obsvns. c write (6,*) inpRow, Xref, Yref, X, Y, int(input(inpRow,1)), dist if (dist.le.maxDist) then c write (6,*) ' ---> obsvn copied across !!' objRow = int(input(inpRow,1)) do col=1,16 objbuffer(objRow,col)=input(inpRow,col) enddo numObsv = numObsv+1 c --- i.e. starting position for next round moves down one since c --- object was copied out and across to the object buffer. c --- ***N.B: This is the only place where inpStart+1 occurs. endif c ===================================================================== c === Only enters writing-out phase if obsvns thought to be complete: if ((rowDepth.gt.maxRowDepth).or.(numObsv.eq.Nsl)) then c cccccccccccc write (6,*) 'C: ' inpStart=inpStart+1 realRow = inpStart*1.0 if ( (inpStart/freq).eq.(int(inpStart/freq)) ) then write (6,*) ' Passing input line: ',inpStart, $ ' Objects found: ',numObj endif c --- i.e. reporting on progress made on working through input catalogue if (distflag2.eq.0) then c --- i.e. only permit newly encountered (Xref,Yref) obsvns past this c --- point for writing of object buffer to output files. c write (6,*) ' ** in the assessing/output stage ...' c ccccccccc write (6,*) 'D: ' c --- Copy contents of object buffer to outRaw array: -- do row=1,Nsl do col=1,16 outRaw(outRow,col)=objbuffer(row,col) enddo outRow = outRow+1 enddo c cccccccccc write (6,*) 'D1: ' c --- Calculating mean values for outMean array: ---------------- do col=1,16 meanRow(col) = 0.0 enddo do row=1,Nsl do col=1,16 meanRow(col) = meanRow(col)+objbuffer(row,col) enddo enddo do col=1,16 meanRow(col) = meanRow(col)/(float(numObsv)) enddo c --- calculating mean value of each column in object buffer c --- (automatically excludes non-detections from the average) c cccccccc write (6,*) 'D2: ' bestCol = 0 bestSN = 0.0 do col=1,4 currSN = meanRow(2*col+3)/meanRow(2*col+4) c --- Calculating difference in fractional flux lost c --- (outside of the aperture presently examined): Rap = apDiam(col)/2.0d0 minLost = 0.0 maxLost = 0.0 do row=1,Nsl Robj = objbuffer(row,14)/2.0d0 if (Robj.gt.0.0) then currLost=exp(-0.693147d0*(Rap/Robj)*(Rap/Robj)) c --- above gives fraction flux lost outside aperture if((minLost.eq.0.0).and.(maxLost.eq.0.0))then c --- i.e. initialising max/min. 1st time thru minLost=currLost maxLost=currLost endif if (currLost.lt.minLost) then minLost=currLost else if (currLost.gt.maxLost) then maxLost=currLost endif endif enddo fracFluxLost = maxLost-minLost c write (6,*) 'col = ',col,' bestCol =',bestCol c write (6,*) ' currSN=',currSN, ' bestSN=',bestSN c --- Testing if aperture is the best found so far: if (currSN.gt.bestSN) then if (fracFluxLost.lt.(S/currSN)) then c --- i.e only include if both criteria satisfied bestSN = currSN bestCol = col endif endif enddo c --- selecting aperture that maximises S/N c cccccccc write (6,*) 'D3: ' if (bestCol.eq.0) then c --- i.e. none of the aperture sizes are suitable for this obj; c --- set to largest aperture and flag with -ve star/gal param bestCol=4 bestSN=currSN meanRow(15)=-1.0*meanRow(15) endif c ccccccccccc write (6,*) 'D4: ' a=meanRow(15) b=numObsv c=meanRow(2*strApNum+3) d=meanRow(14) c --- a=strgal; b=numObsv; c=avgFlx; d=fwhm: to select stars if ((a.gt.strgalCut).and.(b.ge.Nstr)) then if ((c.le.strMaxFlx).and.(c.ge.strMinFlx)) then if ((d.le.strMaxFWHM).and.(d.ge.strMinFWHM)) then bestCol=strApNum bestSN=meanRow(2*bestCol+3)/meanRow(2*bestCol+4) endif endif endif c --- modifying best aperture if object fits stellar criteria c --- (since use a constant large aperture across all stars) if (nextmatch.le.Nmatches) then dist =(meanRow(3)-mat(nextmatch,1))**2 dist = dist+(meanRow(4)-mat(nextmatch,2))**2 dist = sqrt(dist) if (dist.le.maxDist) then write (6,*) ' Obj ',(numObj+1),' is match No. ' $ ,nextmatch,'.' matchflag = 9*(10**Nsl) nextmatch = nextmatch+1 endif endif c --- Checking if current object matches one of the objects in the c --- input match-up list; flagging it (with a "9") if it does. lowSN = 0 row = 1 do while (row.le.Nsl) flux = objbuffer(row,2*bestCol+3) if (flux.lt.0) then lowSN = 1 row = Nsl endif row = row+1 enddo c --- checking for neg flux; flagging this (0=no,1=yes) c ccccccccccc write (6,*) 'D5: ' c --- Deciding if current object is legit: -------------- c --- Criteria: At least 2 or 3 obsvns; SN +ve; SExt flag<8 (neighbs only) if ((numObsv.ge.obsvReqd).and.(lowSN.eq.0).and. $ (meanRow(16).lt.8)) then c ccccccccccc write (6,*) 'E: ' c write (6,*) ' ==> real object; written to catalogue !!' numObj=numObj + 1 outMean(numObj,4) = apDiam(bestCol) do row=1,Nsl outFlux(numObj,row+3) = objbuffer(row,2*bestCol+3) outUncert(numObj,row+3) = objbuffer(row,2*bestCol+4) if (outFlux(numObj,row+3).le.0.0) then currSN = 0.0 else currSN = outFlux(numObj,row+3)/outUncert(numObj,row+3) endif c --- checking for zero fluxes/flux uncertainties if (currSN.le.SNthresh) then outFlux(numObj,row+3)=0.0 c if (threshCol.le.1) then c outFlux(numObj,row+3)=fluxTh(bestCol,row) c endif c if (threshCol.gt.1) then c outFlux(numObj,row+3)=fluxTh(bestCol,threshCol-1) c endif c outUncert(numObj,row+3)=outFlux(numObj,row+3)/SNthresh outUncert(numObj,row+3)=1.0 c currSN = SNthresh currSN = 0.0 endif c NEW METHOD: c --- Assigns flux=0.0 for all measurements less than the threshold. c --- Assigns fluxUncertainty=1.0 for the same objects. c c OLD METHOD: c --- i.e. assigns the threshold flux if that measured by SExtractor c --- is less than the threshold sigma level c --- (this includes zero flux measurements recorded by SExtractor) end do c --- Copying best aperture data from object buffer (to outFlux) c --- along with corresponding uncertainty (to outUncert). outMean(numObj,11)=0.0 do col=1,Nsl if (outFlux(numObj,col+3).gt.0.0) then outMean(numObj,11)=outMean(numObj,11)+1.0 endif enddo c if (outMean(numObj,11).lt.numObsv) then c numObsv=outMean(numObj,11) c outMean(numObj,11)=-1.0*outMean(numObj,11) c endif numObsv=outMean(numObj,11) c --- recounting number of observations; c --- OLDER METHOD: noting with a -ve if some c --- have been lost due to falling below the SNthresh threshold. c --- CURRENT METHOD: c --- numObsv is set to the new value if some obsvns have been discarded. c -------------------------------------------------------------------- c Writing to main output arrays (not the raw output): c -------------------------------------------------------------------- outMean(numObj,1) = numObj outFlux(numObj,1) = numObj outUncert(numObj,1) = numObj do col=2,3 outMean(numObj,col) = meanRow(col+1) outFlux(numObj,col) = meanRow(col+1) outUncert(numObj,col) = meanRow(col+1) enddo c --- i.e. writing ID no., X, Y, star/gal class, ellipticity, FWHM c --- to outMean (in that order of the program lines above). outMean(numObj,8) = meanRow(15) if (meanRow(15).lt.0.0) then c --- i.e. checking if object leaks flux from largest ap numSpill = numSpill+1 endif outMean(numObj,6) = meanRow(13) outMean(numObj,7) = meanRow(14) Rpix2 = (meanRow(3)-Xcen)**2.0+(meanRow(4)-Ycen)**2.0 outMean(numObj,5) = ((beta/Dpix)**2)*Rpix2*0.5 c --- copying wavelength correction to outMean row outMean(numObj,10)=sqrt(Rpix2) c outMean(numObj,11)=numObsv outMean(numObj,12)=meanRow(2*bestCol+3) outMean(numObj,13)=meanRow(2*bestCol+4) outMean(numObj,14)=outMean(numObj,12)/outMean(numObj,13) c --- Also writing radius, mean flux, mean flux uncert and mean S/N c --- ID no., X, and Y also written to outFlux and outUncert. c --- c --- All columns of current row of outMean have been written except c --- for col 9 (SExtractor flag; calculated below). outMean(numObj,9) = 0.0 do row=1,Nsl placeholder = 10**(row-1) outMean(numObj,9) = outMean(numObj,9) + $ (placeholder*objbuffer(row,16)) outMean(numObj,9) = outMean(numObj,9)+matchflag matchflag = 0 enddo c --- Combining SExtractor flags into a single number; a 9 in leftmost c --- position means that object is one of the match-up objects if ((abs(numObsv)).lt.obsvReqd) then c i.e. min numObsv for *any* object is obsvReqd numObj=numObj-1 numLostSN=numLostSN+1 endif c --- i.e. discard object just written if it has only < obsvReqd obsvns c --- (after fluxes are checked against S/N threshold). c --- numLostSN = no. of objects lost because S/N cut reduced to c --- less than the minimum number of required observations c --- (abs value to account for -ve numObsv number (if used to signal c --- below-SNthresh slice removal -- not used currently in program.) end if c ccccccccc write (6,*) 'F' end if c ccccccccc write (6,*) 'G' c --- Set-up for Next Object: --------------------------- inpRow=inpStart do row=0,(k*Nsl-2) c write (6,*) 'row =',row,' k= ',k,' ',(k*Nsl-2) Xprev(k*Nsl-row) = Xprev(k*Nsl-row-1) Yprev(k*Nsl-row) = Yprev(k*Nsl-row-1) enddo Xprev(1) = Xref Yprev(1) = Yref Xref = input(inpRow,3) Yref = input(inpRow,4) refRow = inpRow c --- copy X and Y positions and input row number for reference c ccccccccc write (6,*) 'G1' do row=1,Nsl objbuffer(row,1)=row do col=2,16 objbuffer(row,col)=0.0 c write (6,*) row,col,objbuffer(row,col) enddo enddo c --- initialise object buffer array c ccccccccc write (6,*) 'G2' objRow = int(input(inpRow,1)) if (objRow.lt.1) then objRow=1 c --- i.e. have a zero inpRow the final time through; c --- to avoid problems this causes the objbuffer() array endif do col=1,16 objbuffer(objRow,col)=input(inpRow,col) enddo c --- copy first entry to object buffer c ccccccccc write (6,*) 'G3' numObsv = 1 endif c ccccccccc write (6,*) 'H' inpRow = inpRow+1 c --- ***N.B: This is the only place where inpRow+1 occurs. enddo numObj = numObj-1 write (6,*) ' ------------------------------------------------' c ===================================================================== if (nextmatch.le.Nmatches) then write (6,*) ' ' write (6,*) ' *** NOT ALL MATCH-UP OBJECTS DETECTED ***' write (6,*) ' --> ONLY ',(nextmatch-1),'/',Nmatches,' FOUND.' write (6,*) ' ******************************************' endif c -------------------------------------------------------------------- c Writing output files (12 of them): c -------------------------------------------------------------------- Nstr = Nsl - strgalObs c --- Nstr = minimum number of stellar observations required per object write (6,*) ' ' write (6,*) 'Writing to files...' currFile='outRaw.'//outFilePref write (6,*) ' ',currFile open(unit=11, file=currFile, status='unknown') currFile='starRaw.'//outFilePref write (6,*) ' ',currFile open(unit=12, file=currFile, status='unknown') numStarRaw=0 do row=1,outRow write (11,100) (outRaw(row,col),col=1,16) if (outRaw(row,15).gt.strgalCut) then numStarRaw=numStarRaw+1 write (12,100) (outRaw(row,col),col=1,16) endif enddo close(unit=11) close(unit=12) c --- i.e. long version of all obsvns, vertical sort style 100 format (f4.0, f7.0, f8.2, f8.2, 4(f12.3, f12.4), f7.3, f6.2, $f6.2, f4.0) currFile='outMean.'//outFilePref write (6,*) ' ',currFile open(unit=11, file=currFile, status='unknown') currFile='starMean.'//outFilePref write (6,*) ' ',currFile open(unit=12, file=currFile, status='unknown') numStar=0 do row=1,numObj write (11,120) (outMean(row,col),col=1,14) a=outMean(row,8) b=outMean(row,11) c=outMean(row,12) d=outMean(row,7) c --- a=strgal; b=numObsv; c=avgFlx; d=fwhm: for stellar selection if ((a.gt.strgalCut).and.(b.ge.Nstr)) then c write (6,*) ' --> a=',a,' (',strgalCut,') b=',b if ((c.le.strMaxFlx).and.(c.ge.strMinFlx)) then if ((d.le.strMaxFWHM).and.(d.ge.strMinFWHM)) then numStar=numStar+1 write (12,120) (outMean(row,col),col=1,14) starIDnos(numStar)=outMean(row,1) endif endif endif enddo close(unit=11) close(unit=12) c --- mean value catalogue (all parameters except those flux-related) 120 format (f5.0, f8.2, f8.2, f6.1, f12.8, f7.3, f7.2, f6.2, f14.0, $ f8.2, f5.0, f12.2, f8.2, f8.2) currFile='outFlux.'//outFilePref write (6,*) ' ',currFile open(unit=11, file=currFile, status='unknown') currFile='starFlux.'//outFilePref write (6,*) ' ',currFile open(unit=12, file=currFile, status='unknown') starRow=1 do row=1,numObj write (11,140) (outFlux(row,col),col=1,(Nsl+3)) if (outMean(row,1).eq.starIDnos(starRow)) then write (12,140) (outFlux(row,col),col=1,(Nsl+3)) starRow=starRow+1 endif c --- i.e. checking if object is also a star enddo if (starRow.ne.(numStar+1)) then write (6,*) ' *** ERROR: Not all stars found for starFlux()' write (6,*) ' *********************************************' write (6,*) starRow,(numStar+1) endif close(unit=11) close(unit=12) c --- flux catalogue 140 format (f5.0, f8.2, f8.2, 12(f12.3)) currFile='outUncert.'//outFilePref write (6,*) ' ',currFile open(unit=11, file=currFile, status='unknown') currFile='starUncert.'//outFilePref write (6,*) ' ',currFile open(unit=12, file=currFile, status='unknown') starRow=1 do row=1,numObj write (11,160) (outUncert(row,col),col=1,(Nsl+3)) if (outMean(row,1).eq.starIDnos(starRow)) then write (12,160) (outUncert(row,col),col=1,(Nsl+3)) starRow=starRow+1 endif enddo if (starRow.ne.(numStar+1)) then write (6,*) ' *** ERROR: Not all stars found for starUncert()' write (6,*) ' ***********************************************' write (6,*) starRow,(numStar+1) endif close(unit=11) close(unit=12) c --- flux uncertainties catalogue 160 format (f5.0, f8.2, f8.2, 12(f12.3)) currFile='outSN.'//outFilePref write (6,*) ' ',currFile open(unit=11, file=currFile, status='unknown') currFile='starSN.'//outFilePref write (6,*) ' ',currFile open(unit=12, file=currFile, status='unknown') starRow=1 do row=1,numObj write (11,180) (outFlux(row,col),col=1,3), $ ((outFlux(row,col)/outUncert(row,col)),col=4,(Nsl+3)) if (outMean(row,1).eq.starIDnos(starRow)) then write (12,180) (outFlux(row,col),col=1,3), $ ((outFlux(row,col)/outUncert(row,col)),col=4,(Nsl+3)) starRow=starRow+1 endif enddo if (starRow.ne.(numStar+1)) then write (6,*) ' *** ERROR: Not all stars found for starSN()' write (6,*) ' *******************************************' write (6,*) starRow,(numStar+1) endif close(unit=11) close(unit=12) c --- signal-to-noise ratio catalogue 180 format (f5.0, f8.2, f8.2, 12(f12.3)) currFile='starNorm.'//outFilePref write (6,*) ' ',currFile open(unit=11, file=currFile, status='unknown') currFile='starErrNrm.'//outFilePref write (6,*) ' ',currFile open(unit=12, file=currFile, status='unknown') starRow=1 do row=1,numObj if (outMean(row,1).eq.starIDnos(starRow)) then c --- i.e. if object is a star, then: if (outFlux(row,normSlice+3).eq.0.0) then err1=0.0 else err1=outUncert(row,normSlice+3)/outFlux(row,normSlice+3) endif c --- i.e. relative error on the slice used for flux normalisation do col=4,(Nsl+3) if ((outFlux(row,col).eq.0.0).or. $ (outFlux(row,normSlice+3).eq.0.0)) then rat(col)=0.0 err(col)=0.0 c --- i.e. zero flux entry case (in cols 4 on) else rat(col) = outFlux(row,col)/outFlux(row,normSlice+3) err2(col) = outUncert(row,col)/outFlux(row,col) err(col) = rat(col)*(sqrt((err1)**2 + (err2(col))**2)) endif end do write (11,200) (outFlux(row,col),col=1,3), $ outMean(row,10), outMean(row,14), (rat(col),col=4,(Nsl+3)) write (12,200) (outFlux(row,col),col=1,3), $ outMean(row,10), outMean(row,14), (err(col),col=4,(Nsl+3)) starRow=starRow+1 endif enddo if (starRow.ne.(numStar+1)) then write (6,*) ' *** ERROR: Not all stars found for starNorm()' write (6,*) ' *********************************************' write (6,*) starRow,(numStar+1) endif close(unit=11) close(unit=12) c --- normalised stellar fluxes and normalised uncertainties catalogues 200 format (f5.0, f8.2, f8.2, f12.2, f8.2, 12(f7.3)) currFile='outID.'//outFilePref write (6,*) ' ',currFile open(unit=11, file=currFile, status='unknown') currFile='starID.'//outFilePref write (6,*) ' ',currFile open(unit=12, file=currFile, status='unknown') starRow=1 do row=1,numObj write (11,220) (outMean(row,col),col=2,3),outMean(row,11) c --- record numObsv for each object in the ID version of the full catalogue if (outMean(row,1).eq.starIDnos(starRow)) then write (12,220) (outMean(row,col),col=2,3),outMean(row,1) starRow=starRow+1 endif enddo if (starRow.ne.(numStar+1)) then write (6,*) ' *** ERROR: Not all stars found for starID()' write (6,*) ' *********************************************' write (6,*) starRow,(numStar+1) endif close(unit=11) close(unit=12) c --- ID catalogue (for use with TVMARK) 220 format (f8.2, f8.2, f9.0) c write (6,*) ' ' write (6,*) ' outRaw entries: ',outRow write (6,*) ' starRaw entries: ',numStarRaw write (6,*) ' outMean entries: ',numObj write (6,*) ' starMean entries: ',numStar write (6,*) ' Total objects lost by S/N cut: ',numLostSN write (6,*) ' Total objects with signif flux spill: ',numSpill write (6,*) ' ' write (6,*) ' Check contents of starMean.' write (6,*) ' (k = ',k,'; maxRowDepth = ',maxRowDepth,')' write (6,*) ' ' c write (6,*) ' ' write (6,*) 'Done.' write (6,*) ' ' stop end c ===========================================================================