program modap c Program reads in combined emission-catalogue and aperture corrections and c decides whether object needs flux from a larger aperture to reduce c the size of the original aperture correction. c c c c Input files (emission-line samples) ------------------------------------ c c scaleCoeff (10,3) c slope intercept c 1 2 c c elgLnFlx(3000,(n+3)) (numFitELG) c ID X Y F1 F2 F3 ... Fn c 1 2 3 4 5 6 (n+3) c c elgLnUnc(3000,(n+3)) (numFitELG) c ID X Y dF1 dF2 dF3 ... dFn c 1 2 3 4 5 6 (n+3) c c elgsigdv(3000,(n+3)) **NB: normalised by domin. scatter (numFitELG) c newID X Y sigdv1 sigdv2 sigdv3 ... sigdvn c 1 2 3 4 5 6 (n+3) c c elgLnID(3000,12) (numFitELG) c X Y newID C1 M1 Nfit1 C3 M3 Nfit3 rms/dF devSl1 devSl2 c 1 2 3 4 5 6 7 8 9 10 11 12 c c c c twoMod (10000,20) --> contains Cl=5.0 (2-detect) and 0.0 (1-detect) c ID X Y ap pkWve el FWH stgl APCOR pkSl NoOb avFl avdFl S/N 0.0 0.0 Cl c 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 c c contMod(3000,20) --> contains Cl=2.0 and 3.0 c ID X Y ap ... etc ... avFlx avdFl S/N sigdv Fline dFline Cl waveCor c 1 2 3 4 12 13 14 15 16 17 18 19 c c c c c Output (emission-line samples) ----------------------------------------- c c twoLine(3000,18) --> contains Cl=5.0 (2-detect) and 0.0 (1-detect) c (NB: cols 12-13 and 16-17 identical in this case since no contin) c ID X Y ap pkWv el FWH stgl APCOR pkSl NoOb avF avdF RMS sigdv Fli dFli Cl c 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 c c contLine(3000,18) --> contains Cl=2.0 and 3.0 c ID X Y ap pkWv el FWH stgl APCOR pkSl NoOb avF avdF RMS sigdv Fli dFli Cl c 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 c c - - - - - - - - - - - - - - - - - - - - - - - - - - - c c c contLnFlx(3000,(n+3)) (numFitELG) c ID X Y F1 F2 F3 ... Fn c 1 2 3 4 5 6 (n+3) c c contLnUnc(3000,(n+3)) (numFitELG) c ID X Y dF1 dF2 dF3 ... dFn c 1 2 3 4 5 6 (n+3) c c contsigdv(3000,(n+3)) **NB: normalised by domin. scatter (numFitELG) c newID X Y sigdv1 sigdv2 sigdv3 ... sigdvn c 1 2 3 4 5 6 (n+3) c c contLnID(3000,12) (numFitELG) c X Y newID C1 M1 Nfit1 C3 M3 Nfit3 rms/dF devSl1 devSl2 c 1 2 3 4 5 6 7 8 9 10 11 12 c c - - - - - - - - - - - - - - - - - - - - - - - - - - - c c rejCRid(1000,5) c X Y newID frame where c 1 2 3 4 5 c c DHJ 25/05/99 c ------------------------------------------------------------------------ implicit none character*50 currFile, extn, inWaveFile, inCRfile, inScaleFile character*50 ApmodFile1, ApmodFile2, inApertFile, inRawFile character*50 oldextn, currFile1, currFile2, currFile3, currFile4 integer col, row, Nsl, rowCR, frm, pk, j, Nmin, rej, Nfit integer devsl(2), Nfit(3), numInAp(4), twoInAp(4), intAp(4) integer numIntrErr, numMeasErr, numCRrej, numCRoff, numTotCR integer numHiSl, numLoSl, consecHi, consecLo, prevHi, prevLo integer hiSigSl, loSigSl, iter, found, verbosity real outMean(5000,16), coeff((12-1),2), centWaveL(12) real outFlux(12), outUncert(12), outsigdv((12+3)) real elgLnFlx(5000,(12+3)), elgLnUnc(5000,(12+3)) real elgsigdv(5000,(12+3)), elgLnID(5000,12), rejCRid(3000,5) real*4 t(12), F(12), delF(12), corrWaveL(12), M(3), C(3) real hiFlx(3), loFlx(3), hiSl(3), loSl(3), hiWave(3), loWave(3) real hidF(3), lodF(3), rmsObj(3), rmsErrRat(3) real checkObs, factor real*4 S, Sx, Sy, Stt, Ffit, wave, sigdvcut, fitTol, freq real pkFlx1, pkFlx2, diff, diff2, k, sumd2, n real domScatt, hiSig, loSig, hiAvgWave, loAvgWave real hiSumLine, loSumLine, hiLine, loLine, dist, x1, y1, xc, yc real hiErr, loErr, fOverlap, fTrough, sgn, sigdvnpt real exclEnd, exclNonAdj, nonadj, brtLimit integer numTwo, numCont, colmax, threshCol, chosenAp integer refRow, objRow, objCol, numTotRaw, numObsv, nodev integer numRMSrej, numCorrTwo, numUncorrTwo, numTotTwo integer numCorrCont, numUncorrCont, numTotCont, twoRMS, twoCR real two(10000,20), cont(10000,20), sig(13), apDiam(5) real fluxTh(5,13), corr(5), outtwo(10000,20), outcont(10000,20) real outRaw(100000,20), objbuffer(12,20), CRposns(100000,3) real contLnFlx(3000,(12+3)), contLnUnc(3000,(12+3)) real contsigdv(3000,(12+3)), contLnID(3000,(12+3)) real SNthresh, meanSig, p, corrLim, xx, yy, fluxConst real newFlux, newdF, newAp, oldAp, sigdev, cl real*8 pie pie = 4*atan(1.0) c -------------------------------------------------------------------- c Reading input files: c -------------------------------------------------------------------- freq = 75.0 c --- frequency of input lines to report on progress through input c --- catalogue (freq=500.0 means progress reported every 500th line) write (6,*) ' ' write (6,*) 'Reading input data...' c ----------------------------------- c Reading modap.in: open(unit=20, file='modap.in', status='old') read(20,*,end=14,err=1) oldextn read(20,*,end=14,err=1) extn read(20,*,end=14,err=2) inScaleFile read(20,*,end=14,err=2) ApmodFile1 read(20,*,end=14,err=2) ApmodFile2 read(20,*,end=14,err=2) inWaveFile read(20,*,end=14,err=2) inApertFile read(20,*,end=14,err=2) inRawFile read(20,*,end=14,err=3) xx read(20,*,end=14,err=3) yy read(20,*,end=14,err=3) fluxConst read(20,*,end=14,err=4) threshCol read(20,*,end=14,err=4) sigdvcut read(20,*,end=14,err=5) corrLim read(20,*,end=14,err=5) fitTol read(20,*,end=14,err=6) Nmin read(20,*,end=14,err=6) rej read(20,*,end=14,err=7) fOverlap read(20,*,end=14,err=8) fTrough read(20,*,end=14,err=9) sgn read(20,*,end=14,err=10) inCRfile read(20,*,end=14,err=11) verbosity read(20,*,end=14,err=12) exclEnd read(20,*,end=14,err=12) exclNonAdj read(20,*,end=14,err=12) brtLimit goto 14 1 write (6,*) ' *** Error Reading modap.in: extn ***' goto 13 2 write (6,*) ' *** Error Reading modap.in: filenames ***' goto 13 3 write (6,*) ' *** Error Reading modap.in: xx/yy/fluxConst ***' goto 13 4 write (6,*) ' *** Error Reading modap.in: threshCol/sigdvcut ***' goto 13 5 write (6,*) ' *** Error Reading modap.in: corrLim/fitTol ***' goto 13 6 write (6,*) ' *** Error Reading modap.in: Nmin/rej ***' goto 13 7 write (6,*) ' *** Error Reading modap.in: fOverlap ***' goto 13 8 write (6,*) ' *** Error Reading modap.in: fTrough ***' goto 13 9 write (6,*) ' *** Error Reading modap.in: sgn ***' goto 13 10 write (6,*) ' *** Error Reading modap.in: inCRfile ***' goto 13 11 write (6,*) ' *** Error modap.in: verbosity or exclEnd ***' goto 13 12 write (6,*) ' *** Error modap.in: exclNonAdj or brtLimit ***' goto 13 13 write (6,*) ' *******************************************' 14 continue if (Nmin.lt.3) then write (6,*) ' *******************************************' write (6,*) ' ERROR: Nmin < 3!!' write (6,*) ' --> setting Nmin=3' write (6,*) ' *******************************************' Nmin=3 endif if (rej.gt.2) then write (6,*) ' *******************************************' write (6,*) ' ERROR: max. no. rej pts > 2!!' write (6,*) ' --> setting rej=2' write (6,*) ' *******************************************' rej=2 endif if (abs(sgn).ne.1.0) then write (6,*) ' *******************************************' write (6,*) ' ERROR: sgn is not = +/-1.0 !!' write (6,*) ' --> setting to +1.0' write (6,*) ' *******************************************' sgn=1.0 endif if ((verbosity.ne.1).and.(verbosity.ne.2)) then write (6,*) ' *******************************************' write (6,*) ' ERROR: verbosity is not = 1 or 2' write (6,*) ' --> setting to 2 and continuing' write (6,*) ' *******************************************' verbosity=2 endif if ((exclEnd.ne.1).and.(exclEnd.ne.0)) then write (6,*) ' *******************************************' write (6,*) ' ERROR: exclEnd is not = 0 or 1' write (6,*) ' --> setting to 1 and continuing' write (6,*) ' *******************************************' exclEnd=1 endif if ((exclEnd.ne.1).and.(exclEnd.ne.0)) then write (6,*) ' *******************************************' write (6,*) ' ERROR: exclNonAdj is not = 0 or 1' write (6,*) ' --> setting to 1 and continuing' write (6,*) ' *******************************************' exclNonAdj=1 endif if (verbosity.eq.2) then write (6,*) ' ' write (6,*) 'Output Filename extension: ',extn write (6,*) 'Scale coefficients file: ',inScaleFile write (6,*) '1 or 2-detect EL-cands: ',ApmodFile1 write (6,*) 'Contin-detect EL-cands: ',ApmodFile2 write (6,*) 'Central wavelengths file: ',inWaveFile write (6,*) 'Input aperture sizes/flux thresh file: ',inApertFile write (6,*) 'Input raw fluxes file: ',inRawFile write (6,*) 'X-coord of optical axis centre: ',xx write (6,*) 'Y-coord of optical axis centre: ',yy write (6,*) 'Flux calibration const : ',fluxConst write (6,*) 'Flux thresh column (>=2 => mean values): ',threshCol write (6,*) 'Sigma deviation cut: ',sigdvcut write (6,*) 'Upper limit on correction factors: ',corrLim write (6,*) 'Initial fit precision ' write (6,*) '(as a % of relative flux uncertainty): ',fitTol write (6,*) 'Minimum number of fit pts (all objects): ',Nmin write (6,*) 'Maximum number of rejected pts: ',rej write (6,*) 'Flux scale factor (overlap regions): ',fOverlap write (6,*) 'Flux scale factor (centre trough region): ',fTrough write (6,*) 'Side of fit to reject points from (+/-1.0): ',sgn write (6,*) 'Cosmic ray/ghost file: ',inCRfile write (6,*) 'Output verbosity: ',verbosity write (6,*) 'Exclude end-point devns (1/0=yes/no): ',exclEnd write (6,*) 'Exclude non-adj double-devns (1/0=y/n): ',exclNonAdj write (6,*) 'Reject beyond bright-end (<=0 no-rej): ',brtLimit write (6,*) ' ' endif close(unit=20) c --------------------------------------- c Reading scaleCoeff: c (1st col intercept; 2nd col slope) currFile=inScaleFile row=0 open(unit=20, file=currFile, status='old') 15 row=row+1 read(20,*,end=20,err=20) (coeff(row,col),col=1,2) if (coeff(row,1).eq.1.000) then write (6,*) ' --> Slice ',row,' is the normaliser ...' endif goto 15 20 continue row=row-1 Nsl=row c --- assign Nsl here write (6,*) '(There are ',Nsl,' slices.)' close(unit=20) c -------------------------------------------- c Reading two and cont files: currFile=ApmodFile1 row=0 open(unit=20, file=currFile, status='old') 25 row=row+1 read(20,*,end=30,err=30) (two(row,col),col=1,17) goto 25 30 continue row=row-1 numTwo=row write (6,*) '(There are ',numTwo,' 1 or 2-detect EL-cands.)' close(unit=20) currFile=ApmodFile2 row=0 open(unit=20, file=currFile, status='old') 35 row=row+1 read(20,*,end=40,err=40) (cont(row,col),col=1,19) goto 35 40 continue row=row-1 numCont=row write (6,*) '(There are ',numCont,' continuum-detect EL-cands.)' close(unit=20) c ------------------------------------------------------------ c Reading elgLnFlx, elgLnUnc, elgsigdv and elgLnID files: currFile1='elgLnFlx.'//oldextn currFile2='elgLnUnc.'//oldextn currFile3='elgsigdv.'//oldextn currFile4='elgLnID.'//oldextn open(unit=20, file=currFile1, status='old') open(unit=21, file=currFile2, status='old') open(unit=22, file=currFile3, status='old') open(unit=23, file=currFile4, status='old') row=0 45 row=row+1 read(20,*,end=50,err=50) (elgLnFlx(row,col),col=1,(Nsl+3)) read(21,*,end=50,err=50) (elgLnUnc(row,col),col=1,(Nsl+3)) read(22,*,end=50,err=50) (elgsigdv(row,col),col=1,(Nsl+3)) read(23,*,end=50,err=50) (elgLnID(row,col),col=1,12) goto 45 50 continue row=row-1 if (row.ne.numCont) then write (6,*) '*****************************************' write (6,*) 'ERROR: Mismatch in elg* and numCont no.s!!' write (6,*) '*****************************************' goto 1000 endif close(unit=20) close(unit=21) close(unit=22) close(unit=23) c ------------------------------------- c Reading central wavelength data: j=0 open(unit=20, file=inWaveFile, status='old') 85 j=j+1 read(20,*,end=90,err=90) centwaveL(j) goto 85 90 continue j=j-1 if (j.ne.Nsl) then write (6,*) '*****************************************' write (6,*) 'ERROR: Mismatch in coeff and slice no.s!!' write (6,*) '*****************************************' goto 1000 endif close(unit=20) c ------------------------------------- c Reading CR/ghost positions: j=0 open(unit=20, file=inCRfile, status='old') 95 j=j+1 read(20,*,end=100,err=100) (CRposns(j,col),col=1,3) goto 95 100 continue numTotCR=j-1 write (6,*) '(There are ',numTotCR,' CR/ghost locations.)' close(unit=20) 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 105 row=row+1 read(20,*,end=110,err=110) apDiam(row),(fluxTh(row,col),col=1,Nsl) goto 105 110 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 115 row=row+1 read(20,*,end=120,err=120) apDiam(row),(fluxTh(row,col), $ col=1,colmax) goto 115 120 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 outRaw flux-file lines: j=0 open(unit=20, file=inRawfile, status='old') 135 j=j+1 read(20,*,end=140,err=140) (outRaw(j,col),col=1,16) goto 135 140 continue numTotRaw=j-1 write (6,*) '(There are ',numTotRaw,' raw catalogue lines.)' close(unit=20) c -------------------------------------------------------------------- c 1 or 2-detect EL-cands: c -------------------------------------------------------------------- numCRrej=0 numRMSrej=0 numCorrTwo=0 numUncorrTwo=0 numTotTwo=0 do row=1,4 numInAp(row)=0 enddo write (6,*) ' ' write (6,*) ' ' write (6,*) ' ' write (6,*) '1 or 2-detect EL-cands - - - - - - - - - - - -' write (6,*) ' ' do row=1,numTwo if (two(row,9).le.corrLim) then c ~~~ i.e. <= corrLim, no re-measurement necessary ~ ~ ~ ~ ~ ~ ~ c ~~~ (corrLim = maximum allowable apert corrn) ~~~ c --- 1. Checking for CRs (flag FOUND): -------- xc=two(row,2) yc=two(row,3) pk=int(two(row,10)) found=0 dist=0.0 c --- found=0: no CR; found=1: x,y but not frame (flagged) c --- found=2: CR match in both x,y and frame number rowCR=1 do while ((rowCR.le.numTotCR).and.(found.lt.2)) x1=CRposns(rowCR,1) y1=CRposns(rowCR,2) frm=int(CRposns(rowCR,3)) if ((frm.le.0).or.(frm.gt.Nsl)) then write (6,*) ' *** ERROR: CR file has frm = ',frm frm=0 endif dist=sqrt(((xc-x1)**2)+((yc-y1)**2)) if (dist.le.(0.5*two(row,4))) then c --- i.e. CR pixel falls within aperture radius numCRoff=numCRoff+1 c write (6,*) ' *** Obj',row,' hits CR pix',x1,y1,', but not frame' c found=1 found=2 c --- for 1-detections (cl=0.0): c if ((frm.eq.pk).and.(two(row,17).eq.0.0)) then c found=2 c endif c --- for 2-detections (cl=5.0): c if ((frm.eq.pk).or.(frm.eq.(pk+1))) then c if (two(row,17).eq.5.0) then c found=2 c endif c endif endif rowCR=rowCR+1 enddo if (found.eq.2) then write (6,*) ' Obj ',two(row,1),' *** CR in OLD apert' numCRrej=numCRrej+1 rejCRid(numCRrej,1)=two(row,2) rejCRid(numCRrej,2)=two(row,3) rejCRid(numCRrej,3)=two(row,1) rejCRid(numCRrej,4)=float(frm) rejCRid(numCRrej,5)=1.0 endif if (found.ne.2) then numUncorrTwo=numUncorrTwo+1 numTotTwo=numTotTwo+1 do col=1,13 outtwo(numTotTwo,col)=two(row,col) enddo j=int(two(row,10)) outtwo(numTotTwo,14)=sig(j) outtwo(numTotTwo,15)=two(row,14) outtwo(numTotTwo,16)=two(row,12) outtwo(numTotTwo,17)=two(row,13) outtwo(numTotTwo,18)=two(row,17) c --- i.e. same 18-col output format as the "elgLine" files c --- Includes calibration to physical units and apert corrn c --- i.e. copy direct to output chosenAp=0 col=1 do while ((chosenAp.eq.0).and.(col.le.4)) if (apDiam(col).eq.two(row,4)) then chosenAp=col endif col=col+1 enddo numInAp(chosenAp)=numInAp(chosenAp)+1 c --- i.e. determining aperture number and adding to total endif else c ~~~ i.e. > corrLim, new measurement req'd ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ chosenAp = 0 do col=1,4 p = log(2.0) * ((apDiam(col)/two(row,7))**2) corr(col)=1/(1-exp(-1*p)) if ((chosenAp.eq.0.0).and.(corr(col).lt.corrLim)) then chosenAp = col endif enddo if (chosenAp.eq.0.0) then chosenAp = 4 endif c --- computing and choosing aperture with suitably small corrn c write (6,*) two(row,1),': ',corr(1),corr(2),corr(3),corr(4), c $ ' choose: ',chosenAp c write (6,*) two(row,1),': ',two(row,4),'-->',apDiam(chosenAp) c --- 0. Reassigning apert to *new* value: -------- c --- Also updates with new apert correction value c %%%%%% TEST SECTION TO RUN NEW CODE ON OLD APERTURE FLUXES %%%%% c j=1 c chosenAp=0 c do while (two(row,4).ne.apDiam(j)) c j=j+1 c enddo c chosenAp=j c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% oldAp=two(row,4) newAp=apDiam(chosenAp) two(row,4)=newAp two(row,9)=corr(chosenAp) numInAp(chosenAp)=numInAp(chosenAp)+1 c --- 1. Checking for CRs (flag FOUND): -------- xc=two(row,2) yc=two(row,3) pk=int(two(row,10)) found=0 dist=0.0 c --- found=0: no CR; found=1: x,y but not frame (flagged) c --- found=2: CR match in both x,y and frame number rowCR=1 do while ((rowCR.le.numTotCR).and.(found.lt.2)) x1=CRposns(rowCR,1) y1=CRposns(rowCR,2) frm=int(CRposns(rowCR,3)) if ((frm.le.0).or.(frm.gt.Nsl)) then write (6,*) ' *** ERROR: CR file has frm = ',frm frm=0 endif dist=sqrt(((xc-x1)**2)+((yc-y1)**2)) if (dist.le.(0.5*two(row,4))) then c --- i.e. CR pixel falls within aperture radius numCRoff=numCRoff+1 c write (6,*) ' *** Obj',row,' hits CR pix',x1,y1,', but not frame' c found=1 found=2 c --- for 1-detections (cl=0.0): c if ((frm.eq.pk).and.(two(row,17).eq.0.0)) then c found=2 c endif c --- for 2-detections (cl=5.0): c if ((frm.eq.pk).or.(frm.eq.(pk+1))) then c if (two(row,17).eq.5.0) then c found=2 c endif c endif endif rowCR=rowCR+1 enddo if (found.eq.2) then write (6,*) ' Obj ',two(row,1),' *** CR in new apert' numCRrej=numCRrej+1 rejCRid(numCRrej,1)=two(row,2) rejCRid(numCRrej,2)=two(row,3) rejCRid(numCRrej,3)=two(row,1) rejCRid(numCRrej,4)=float(frm) rejCRid(numCRrej,5)=1.0 endif c --- 2. Extract new flux and compute uncertainty: -------- c --- Initialise object buffer array: do objRow=1,Nsl do objCol=1,16 objbuffer(objRow,objCol)=0.0 enddo enddo refRow=1 numObsv=0 xc=two(row,2) yc=two(row,3) do while ((refRow.le.numTotRaw).and.(numObsv.lt.Nsl)) x1=outRaw(refRow,3) y1=outRaw(refRow,4) dist=sqrt((x1-xc)**2 + (y1-yc)**2) if (dist.le.5.0) then j=int(outRaw(refRow,1)) c --- i.e. slice number if (objbuffer(j,1).ne.0.0) then write (6,*) ' Obj ',two(row,1),'** WARNING: objbuffer overwrite!!' endif do objCol=1,16 objbuffer(j,objCol)=outRaw(refRow,objCol) enddo c --- i.e. copying line to objbuffer numObsv=numObsv+1 endif refRow=refRow+1 enddo c --- Scale by scaleCoeff values: x1=two(row,2) y1=two(row,3) dist=sqrt((x1-xx)**2 + (y1-yy)**2) c --- i.e. dist of object from optical centre do objRow=1,Nsl factor=coeff(objRow,1)+(dist*coeff(objRow,2)) do j=1,4 objbuffer(objRow,2*j+3)=objbuffer(objRow,2*j+3)/factor objbuffer(objRow,2*j+4)=objbuffer(objRow,2*j+4)/factor enddo c --- i.e. correcting fluxes and uncerts in each aperture enddo c --- Computing average values for new flux and uncert: newFlux=0 newdF=0 if (numObsv.eq.0) then write (6,*) '** WARNING: No match found in outRaw **' else do j=1,Nsl newFlux=newFlux+objbuffer(j,2*chosenAp+3) newdF=newdF+objbuffer(j,2*chosenAp+4) enddo newFlux=newFlux/float(numObsv) newdF=newdF/float(numObsv) endif c ---> Do not need to re-compute uncertainties from c --- SMOOTHED SExt data since already done for outRaw by extspect c --- 3. Compute RMS deviation and check (flag NODEV): -------- sigdev=newFlux/newdF c --- i.e. in the case of no continuum, sigdev=S/N ratio nodev=0 if (sigdev.ge.sigdvcut) then nodev=1 else write (6,*) ' *** Object ',two(row,1),' ',objbuffer(j,1),objbuffer(j,3),objbuffer(j,4), c $ objbuffer(j,5),objbuffer(j,6),objbuffer(j,7) c enddo c write (6,*) '.......... new avgFlux, dF = ',newFlux,newdF c write (6,*) '....... radius = ',dist,' sigdev=',sigdev c write (6,*) ' ' c endif endif c ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ enddo write (6,*) ' ' write (6,*) ' ' write (6,*) '1 or 2-detect stats:' write (6,*) ' ' write (6,*) 'Num 1 or 2-detect INPUT: ',numTwo write (6,*) ' --------------' write (6,*) 'Num objects with RMS < sigdvcut: ',numRMSrej write (6,*) 'Num CR rejected: ',numCRrej write (6,*) ' ' write (6,*) 'Num kept same apert: ',numUncorrTwo write (6,*) 'Num with new apert: ',numCorrTwo write (6,*) ' --------------' write (6,*) 'Num 1 or 2-detect OUTPUT: ',numTotTwo write (6,*) ' ' write (6,*) ' ' write (6,*) ' ' c -------------------------------------------------------------------- c Contin-detect EL-cands: c -------------------------------------------------------------------- twoRMS=numRMSrej twoCR=numCRrej c --- i.e store the numbers of these from 1 and 2-detect analysis numCorrCont=0 numUncorrCont=0 numTotCont=0 do row=1,4 twoInAp(row)=numInAp(row) numInAp(row)=0 intAp(row)=int(apDiam(row)) enddo numIntrErr=0 numMeasErr=0 write (6,*) ' ' write (6,*) 'Contin-detect EL-cands - - - - - - - - - - - -' write (6,*) ' ' do row=1,numCont if (cont(row,9).le.corrLim) then c ~~~ i.e. <= corrLim, no re-measurement necessary ~ ~ ~ ~ ~ ~ ~ c ~~~ (corrLim = maximum allowable apert corrn) ~~~ numUncorrCont=numUncorrCont+1 numTotCont=numTotCont+1 do col=1,18 outcont(numTotCont,col)=cont(row,col) enddo j=int(cont(row,10)) outcont(numTotCont,14)=sig(j) c --- i.e. same 18-col output format as the "elgLine" files c --- Includes calibration to physical units and apert corrn do col=1,(Nsl+3) contLnFlx(numTotCont,col)=elgLnFlx(row,col) contLnUnc(numTotCont,col)=elgLnUnc(row,col) contsigdv(numTotCont,col)=elgsigdv(row,col) enddo do col=1,12 contLnID(numTotCont,col)=elgLnID(row,col) enddo c --- i.e. copy direct to output (all files) chosenAp=0 col=1 do while ((chosenAp.eq.0).and.(col.le.4)) if (apDiam(col).eq.cont(row,4)) then chosenAp=col endif col=col+1 enddo numInAp(chosenAp)=numInAp(chosenAp)+1 c --- i.e. determining aperture number and adding to total else c ~~~ i.e. > corrLim, new measurement req'd ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ chosenAp = 0 do col=1,4 p = log(2.0) * ((apDiam(col)/cont(row,7))**2) corr(col)=1/(1-exp(-1*p)) if ((chosenAp.eq.0.0).and.(corr(col).lt.corrLim)) then chosenAp = col endif enddo if (chosenAp.eq.0.0) then chosenAp = 4 endif c --- computing and choosing aperture with suitably small corrn c write (6,*) cont(row,1),': ',corr(1),corr(2),corr(3),corr(4), c $ ' choose: ',chosenAp write (6,*) cont(row,1),': ',cont(row,4),'-->',apDiam(chosenAp) c --- 0. Reassigning apert to *new* value: -------- c --- Also updates with new apert correction value c %%%%%% TEST SECTION TO RUN NEW CODE ON OLD APERTURE FLUXES %%%%% c j=1 c chosenAp=0 c do while (cont(row,4).ne.apDiam(j)) c j=j+1 c enddo c chosenAp=j c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% oldAp=cont(row,4) newAp=apDiam(chosenAp) cont(row,4)=newAp cont(row,9)=corr(chosenAp) numInAp(chosenAp)=numInAp(chosenAp)+1 c --- 1. Checking for CRs (flag FOUND): -------- xc=cont(row,2) yc=cont(row,3) pk=int(cont(row,10)) found=0 dist=0.0 c --- found=0: no CR; found=1: x,y but not frame (flagged) c --- found=2: CR match in both x,y and frame number rowCR=1 do while ((rowCR.le.numTotCR).and.(found.lt.2)) x1=CRposns(rowCR,1) y1=CRposns(rowCR,2) frm=int(CRposns(rowCR,3)) if ((frm.le.0).or.(frm.gt.Nsl)) then write (6,*) ' *** ERROR: CR file has frm = ',frm frm=0 endif dist=sqrt(((xc-x1)**2)+((yc-y1)**2)) if (dist.le.(0.5*cont(row,4))) then c --- i.e. CR pixel falls within aperture radius numCRoff=numCRoff+1 c write (6,*) ' *** Obj',row,' hits CR pix',x1,y1,', but not frame' found=1 c --- for 1 peak (cl=2.0): if ((frm.eq.pk).and.(cont(row,17).eq.2.0)) then found=2 endif c --- for 2 peaks (cl=3.0): if ((frm.eq.pk).or.(frm.eq.(pk+1))) then if (cont(row,17).eq.3.0) then found=2 endif endif endif rowCR=rowCR+1 enddo if (found.eq.2) then write (6,*) ' *** CR found for new apert' numCRrej=numCRrej+1 rejCRid(numCRrej,1)=cont(row,2) rejCRid(numCRrej,2)=cont(row,3) rejCRid(numCRrej,3)=cont(row,1) rejCRid(numCRrej,4)=float(frm) rejCRid(numCRrej,5)=1.0 endif c --- 2. Extract new flux and compute uncertainty: -------- c --- Initialise object buffer array: do objRow=1,Nsl do objCol=1,16 objbuffer(objRow,objCol)=0.0 enddo enddo refRow=1 numObsv=0 xc=cont(row,2) yc=cont(row,3) do while ((refRow.le.numTotRaw).and.(numObsv.lt.Nsl)) x1=outRaw(refRow,3) y1=outRaw(refRow,4) dist=sqrt((x1-xc)**2 + (y1-yc)**2) if (dist.le.5.0) then j=int(outRaw(refRow,1)) c --- i.e. slice number if (objbuffer(j,1).ne.0.0) then write (6,*) '** WARNING: objbuffer over-write !!' endif do objCol=1,16 objbuffer(j,objCol)=outRaw(refRow,objCol) enddo c --- i.e. copying line to objbuffer numObsv=numObsv+1 endif refRow=refRow+1 enddo c --- Scale by scaleCoeff values: x1=cont(row,2) y1=cont(row,3) dist=sqrt((x1-xx)**2 + (y1-yy)**2) c --- i.e. dist of object from optical centre do objRow=1,Nsl factor=coeff(objRow,1)+(dist*coeff(objRow,2)) do j=1,4 objbuffer(objRow,2*j+3)=objbuffer(objRow,2*j+3)/factor objbuffer(objRow,2*j+4)=objbuffer(objRow,2*j+4)/factor enddo c --- i.e. correcting fluxes and uncerts in each aperture enddo c --- Computing average values for new flux and uncert: newFlux=0 newdF=0 if (numObsv.eq.0) then write (6,*) '** WARNING: No match found in outRaw **' else do j=1,Nsl newFlux=newFlux+objbuffer(j,2*chosenAp+3) newdF=newdF+objbuffer(j,2*chosenAp+4) outFlux(j)=objbuffer(j,2*chosenAp+3) outUncert(j)=objbuffer(j,2*chosenAp+4) enddo newFlux=newFlux/float(numObsv) newdF=newdF/float(numObsv) endif c ---> Do not need to re-compute uncertainties from c --- SMOOTHED SExt data since already done for outRaw by extspect c ----------------------------------------------------------------------- c --- Continuum-Fitting Routine: corrected set of wavelength values used c --- (zero flux columns excluded from the fits) c write (6,*) row,'- - - - - - - - - - - - - - - - - - - - - -' devsl(1)=0.0 devsl(2)=0.0 do col=1,3 M(col)=0.0d0 C(col)=0.0d0 rmsObj(col)=0.0 Nfit(col)=0 rmsErrRat(col)=0.0 enddo c ------------------------------------------------------------ c --- A. Initial least-squares fitting: (no rejection) S=0.0 Sx=0.0 Sy=0.0 Stt=0.0 checkObs=0.0 c --- Loop 1: ------------------------------------------ do col=1,Nsl t(col) = 0.0 F(col)=0.0 delF(col)=0.0 corrWaveL(col)=0.0 if (outFlux(col).gt.0.0) then c --- i.e. checking that there is a flux detection c write (6,*), 'col= ',col checkObs=checkObs+1.0 Nfit(1)=Nfit(1)+1 corrWaveL(col) = (1-cont(row,19))*centWaveL(col) F(col) = outFlux(col) delF(col) = outUncert(col) S = S+( 1/(delF(col)*delF(col)) ) Sx = Sx+( corrWaveL(col)/(delF(col)*delF(col)) ) Sy = Sy+( F(col)/(delF(col)*delF(col)) ) endif enddo if (checkObs.ne.cont(row,11)) then c --- i.e. checking that no. of fluxes matches number of obsvns write (6,*) 'ERROR: Mismatch in no. of obsvns!! (at A)' write (6,*) 'checkObs,outMean =',checkObs,cont(row,11) endif c --- Loop 2: -------------------------------------------- checkObs=0.0 do col=1,Nsl if (outFlux(col).gt.0.0) then checkObs=checkObs+1.0 t(col) = (corrWaveL(col)-(Sx/S))/delF(col) Stt = Stt+(t(col)*t(col)) M(1) = M(1)+( (t(col)*F(col))/delF(col) ) endif enddo if (checkObs.ne.cont(row,11)) then c --- i.e. checking that number of obsvns again write (6,*) 'ERROR: Mismatch in no. of obsvns!! (at B)' endif M(1) = M(1)/Stt C(1) = (Sy-Sx*M(1))/S c ---------------------------------------------------- c --- Iterative least-squares fitting with rejection: c ---------------------------------------------------- do iter=1,2 c write (6,*) ' ',iter,':' c write (6,*) ' M, C:',M(iter),C(iter),Nfit(iter) c ------------------------------------------------------------ c --- 2. Measuring RMS from non-zero flux pts (excl deviants) c --- (used to decide whether weight/unweight next LS-fit) sumd2=0.0 n=0.0 do col=1,Nsl if (outFlux(col).gt.0.0) then c --- i.e. checking for a flux detection k=float(col) if ((k.ne.devsl(1)).and.(k.ne.devsl(2))) then c write (6,*), 'k= ',k c --- i.e. checking not a deviant point wave = (1-cont(row,19))*centWaveL(col) Ffit = C(iter)+(wave*M(iter)) c --- i.e. y = C+Mx diff2 = (outFlux(col)-Ffit)**2 sumd2=sumd2+diff2 n=n+1.0 endif endif enddo rmsObj(iter)=sqrt(sumd2/(n-1.0)) rmsErrRat(iter)=rmsObj(iter)/outMean(row,13) c write (6,*) ' rms:',rmsObj(iter),rmsErrRat(iter) c -------------------------------------------------------- c --- 3. Finding slice no.s of 1 or 2 most deviant pt(s) devsl(1)=0.0 devsl(2)=0.0 if (cont(row,11).ge.float(Nmin+1)) then c --- i.e. checking at least 1 observation for rejection pkFlx1 = 0.0 pkFlx2 = 0.0 sigdvnpt=0.0 c --- Calculating dominant source of scatter: if (rmsErrRat(iter).gt.1.0) then domScatt=rmsObj(iter) c --- i.e. measurement errors dominate else domScatt=newdF c --- i.e. intrinsic errors dominate endif do col=1,Nsl if (outFlux(col).gt.0.0) then c --- i.e. checking for a flux detection wave = (1-cont(row,19))*centWaveL(col) Ffit = C(iter)+(wave*M(iter)) c --- i.e. y = C+Mx diff = sgn*(outFlux(col)-Ffit) c --- i.e. sgn=+1.0 biases fit towards +VE devns; c --- sgn=-1.0 biases fit towards -VE devns sigdvnpt = diff/domScatt c --- deviation of the point (in units of sigma) if (sigdvnpt.gt.1.0) then c --- i.e. only if devn is significant wrt scatt if (diff.gt.pkFlx1) then devsl(2)=devsl(1) pkFlx2=pkFlx1 devsl(1)=col pkFlx1=diff else if ((diff.le.pkFlx1).and.(diff.gt.pkFlx2)) then devsl(2)=col pkFlx2=diff endif endif endif enddo if (iter.eq.1) then devsl(2) = 0.0 c --- i.e. only reject 1 point if first time through endif if ((cont(row,11).eq.float(Nmin+1)).or.(rej.eq.1)) then devsl(2)=0.0 c --- i.e. devsl(2)=devsl(1) if only 1 slice to reject endif c write (6,*) ' --> devsl() = ',(devsl(col),col=1,2) c ---------------------------------------------------- c --- 4. Least-squares fitting: (including rejection) S=0.0 Sx=0.0 Sy=0.0 Stt=0.0 checkObs=0.0 c --- Loop 1: --------------------------------------- do col=1,Nsl t(col) = 0.0 F(col)=0.0 delF(col)=0.0 corrWaveL(col)=0.0 if (outFlux(col).gt.0.0) then c --- i.e. checking that there is a flux detection checkObs=checkObs+1.0 if((col.ne.devsl(1)).and.(col.ne.devsl(2)))then c --- checking that slice is not a deviant one Nfit(iter+1)=Nfit(iter+1)+1 corrWaveL(col) = (1-cont(row,19))*centWaveL(col) F(col) = outFlux(col) delF(col) = outUncert(col) if (rmsErrRat(iter).gt.1.0) then delF(col)=1.0 c write (6,*) '**unweighted**' c --- i.e. no weighting if domin by intr err endif S = S+( 1/(delF(col)*delF(col)) ) Sx=Sx+(corrWaveL(col)/(delF(col)*delF(col))) Sy=Sy+(F(col)/(delF(col)*delF(col))) endif endif enddo if (checkObs.ne.cont(row,11)) then write (6,*) 'ERROR: Mismatch in no. of obsvns!!' write (6,*) 'checkObs,outMean =',checkObs,cont(row,11) endif c --- Loop 2: --------------------------------- checkObs=0.0 do col=1,Nsl if (outFlux(col).gt.0.0) then checkObs=checkObs+1.0 if((col.ne.devsl(1)).and.(col.ne.devsl(2)))then c write (6,*) ' col= ',col t(col) = (corrWaveL(col)-(Sx/S))/delF(col) Stt = Stt+(t(col)*t(col)) M(iter+1)= M(iter+1)+( (t(col)*F(col))/delF(col)) endif endif enddo if (checkObs.ne.cont(row,11)) then c --- i.e. checking that number of obsvns again write (6,*) 'ERROR: Mismatch in no. of obsvns!!' endif M(iter+1) = M(iter+1)/Stt C(iter+1) = (Sy-Sx*M(iter+1))/S c ---------------------------------------------------------------------- else if (cont(row,11).eq.float(Nmin)) then c --- i.e. number of obsvn equals the min no. allowed, c --- so no rejection (copy fit parameters from 1st fit) M(iter+1) = M(1) C(iter+1) = C(1) Nfit(iter+1) = Nfit(1) rmsObj(iter+1) = rmsObj(1) rmsErrRat(iter+1) = rmsErrRat(1) endif enddo c --- end of fitting iteration loop c ------------------------------------------------------------ c --- 5. Measuring RMS from final fit (excl deviants) sumd2=0.0 n=0.0 do col=1,Nsl if (outFlux(col).gt.0.0) then c --- i.e. checking for a flux detection k=float(col) if ((k.ne.devsl(1)).and.(k.ne.devsl(2))) then c write (6,*) 'k=',k c --- i.e. checking not a deviant point wave = (1-cont(row,19))*centWaveL(col) Ffit = C(3)+(wave*M(3)) c --- i.e. y = C+Mx diff2 = (outFlux(col)-Ffit)**2 sumd2=sumd2+diff2 n=n+1.0 endif endif enddo rmsObj(3)=sqrt(sumd2/(n-1.0)) rmsErrRat(3)=rmsObj(3)/newdF c --- Because a fit has been made we can determine what aspect c --- dominates the measurement errors if (rmsErrRat(3).gt.1.0) then numIntrErr=numIntrErr+1 else if (rmsErrRat(3).gt.0.0) then numMeasErr=numMeasErr+1 endif c ----------------------------------------------------------------------- c Calculating sigma devn at each wavelength and writing to outsigdv() outsigdv(1)=cont(row,1) outsigdv(2)=cont(row,2) outsigdv(3)=cont(row,3) do col=1,Nsl wave = (1-cont(row,19))*centWaveL(col) Ffit = C(3)+(wave*M(3)) c --- i.e. y = a+bx if (rmsErrRat(3).gt.1.0) then domScatt=rmsObj(3) c --- i.e. measurement errors dominate the continuum scatter else domScatt=newdF c --- i.e. intrinsic errors dominate the continuum scatter endif if (outFlux(col).gt.0.0) then outsigdv(col+3)=(outFlux(col)-Ffit)/domScatt else outsigdv(col+3)=0.0 endif c --- i.e. outsigdv=0.0 if no flux measurement at that point enddo c ------------------------------------------------------------------------ c Checking peak sigma deviation: numHiSl=0 hiSig=0.0 hiSigSl=0 numLoSl=0 loSig=0.0 loSigSl=0 do col=1,3 hiFlx(col)=0.0 hidF(col)=0.0 hiSl(col)=0 hiWave(col)=0.0 loFlx(col)=0.0 lodF(col)=0.0 loSl(col)=0 loWave(col)=0.0 enddo consecHi=0 consecLo=0 prevHi=0 prevLo=0 c --- initialising values hiSumLine=0.0 hiLine=0.0 hiErr=0.0 hiAvgWave=0.0 loSumLine=0.0 loLine=0.0 loErr=0.0 loAvgWave=0.0 c --- to hold weighted-mean values for line flux and wavelength c --- (in both cases of emission and absorption) c --- Filling 3-space high/low deviant arrays; in order of col: do col=1,Nsl if (outsigdv(col+3).ge.sigdvcut) then numHiSl=numHiSl+1 if (numHiSl.le.3) then c --- i.e. only assign latest values if there is space hiSl(numHiSl)=col wave=(1-cont(row,19))*centWaveL(col) hiWave(numHiSl)=wave Ffit=C(3)+(M(3)*wave) hiFlx(numHiSl)=outFlux(col)-Ffit hidF(numHiSl)=outUncert(col) if (outsigdv(col+3).gt.hiSig) then c --- i.e. checking if maximum sigma deviation hiSig=outsigdv(col+3) hiSigSl=col hiLine=hiFlx(numHiSl) hiErr=hidF(numHiSl) c --- i.e. line flux and uncertainty from peak value endif if (numHiSl.eq.1) then consecHi=1 endif if (numHiSl.ge.2) then if ((col-prevHi).eq.1) then consecHi=consecHi+1 endif endif prevHi=col endif else if (outsigdv(col+3).le.(-1*sigdvcut)) then numLoSl=numLoSl+1 if (numLoSl.le.3) then loSl(numLoSl)=col wave=(1-cont(row,19))*centWaveL(col) loWave(numLoSl)=wave Ffit=C(3)+(M(3)*wave) loFlx(numLoSl)=outFlux(col)-Ffit lodF(numLoSl)=outUncert(col) if (outsigdv(col+3).lt.loSig) then loSig=outsigdv(col+3) loSigSl=col loLine=loFlx(numLoSl) loErr=lodF(numLoSl) endif if (numLoSl.eq.1) then consecLo=1 endif if (numLoSl.ge.2) then if ((col-prevLo).eq.1) then consecLo=consecLo+1 endif endif prevLo=col endif endif enddo c --- Finding mean wavelength (weighted by line flux): if (numHiSl.gt.0) then do col=1,numHiSl hiSumLine=hiSumLine+hiFlx(col) c --- unweighted sum of line fluxes hiAvgWave=hiAvgWave+(hiWave(col)*hiFlx(col)) c --- weighted sum of wavelengths enddo hiAvgWave=hiAvgWave/hiSumLine c --- (weighted wavelength sum)/(line flux sum) endif if (numLoSl.gt.0) then do col=1,numLoSl loSumLine=loSumLine+loFlx(col) loAvgWave=loAvgWave+(loWave(col)*loFlx(col)) enddo loAvgWave=loAvgWave/loSumLine endif nonadj=abs(devSl(2)-devSl(1)) c --- i.e. non-adj of double-devn occurs if nonadj.ne.1.0 c --- Classifying nature of high/low deviant sets: if ((((exclEnd.eq.1.0).and.(hiSigSl.gt.1).and. $ (hiSigSl.lt.Nsl))).or.(exclEnd.eq.0.0)) then c --- i.e. checking that devns are not on end-pts if important if (((exclNonAdj.eq.1.0).and.(nonadj.eq.1.0)).or. $ (exclNonAdj.eq.0.0)) then c --- i.e. checking that devns are adj if there are two of them if (((brtLimit.gt.0.0).and.(newFlux.le.brtLimit)).or. $ (brtLimit.le.0.0)) then c --- i.e. checking that EL-cand are .le. bright-end limit if c --- limit is being applied (i.e. limit is set to something +ve) if (numHiSl.eq.1) then c --- class +1 devn: cl=2.0 else if (numHiSl.eq.2) then if (consecHi.eq.2) then c --- class +2 devn: cl=3.0 else c --- class +7 devn: (excluded) cl=8.0 endif else if (numHiSl.eq.3) then if (consecHi.eq.3) then if (hiSigSl.eq.hiSl(2)) then c --- class +3 devn: cl=4.0 else c --- class +9 devn: (excluded) cl=10.0 endif else c --- class +7 devn again: (excluded) cl=8.0 endif else if (numHiSl.eq.0) then c --- class 0 devn: (excluded) cl=1.0 else c --- class +8 devn: (excluded) cl=9.0 endif endif c --- end of brtLimit check if endif c --- end of exclNonAdj check if endif c --- end of exclEndpt check if c ------------------------------------------------------ c --- Deciding whether object still makes the cut as c --- an EL-cand with the fluxes from the new aperture: nodev=0 if (hiSig.ge.sigdvcut) then nodev=1 else write (6,*) ' *** Object ',cont(row,1), $ ' now EXCLUDED from sample' nodev=0 endif c ------------------------------------------------------------------------ c Calibrating and writing to file if object passes if ((found.ne.2).and.(nodev.ge.1)) then numCorrCont=numCorrCont+1 numTotCont=numTotCont+1 if ((cl.eq.2.0).or.(cl.eq.4.0)) then factor=fTrough endif if (cl.eq.3.0) then factor=fOverlap endif c --- outcont file: do col=1,18 outcont(numTotCont,col)=cont(row,col) enddo outcont(numTotCont,4)=newAp*-1.0 outcont(numTotCont,5)=hiAvgWave outcont(numTotCont,10)=hiSigSl outcont(numTotCont,12)=newFlux*fluxConst*cont(row,9) outcont(numTotCont,13)=newdF*fluxConst*cont(row,9) j=int(devSl(1)) c if (j.ne.0) then if ((cl.eq.2.0).or.(cl.eq.4)) then outcont(numTotCont,14)=sig(j) endif if (cl.eq.3.0) then outcont(numTotCont,14)=(sig(j)+sig(int(devSl(2))))/2.0 endif c else c outcont(numTotCont,14)=0.0 c endif outcont(numTotCont,15)=hiSig outcont(numTotCont,16)=hiLine*fluxConst*cont(row,9)/factor outcont(numTotCont,17)=hiErr*fluxConst*cont(row,9)/factor outcont(numTotCont,18)=cl if (cl.ne.cont(row,18)) then write (6,*) ' *** emiss-class changed: ',cont(row,18), $ ' ---> ',cl endif c --- cont(row,9) stores the NEW apert correction c --- ** N.B: Line fluxes modified by factor but not contin flx c --- contLnFlx, contLnUnc, contsigdv, contLnID files: do col=1,3 contLnFlx(numTotCont,col) = cont(row,col) contLnUnc(numTotCont,col) = cont(row,col) contsigdv(numTotCont,col) = cont(row,col) enddo do col=1,Nsl contLnFlx(numTotCont,col+3)=outFlux(col) contLnUnc(numTotCont,col+3)=outUncert(col) contsigdv(numTotCont,col+3)=outsigdv(col+3) enddo contLnID(numTotCont,1) = cont(row,2) contLnID(numTotCont,2) = cont(row,3) contLnID(numTotCont,3) = cont(row,1) contLnID(numTotCont,4) = C(1) contLnID(numTotCont,5) = M(1) contLnID(numTotCont,6) = float(Nfit(1)) contLnID(numTotCont,7) = C(3) contLnID(numTotCont,8) = M(3) contLnID(numTotCont,9) = float(Nfit(3)) contLnID(numTotCont,10) = rmsErrRat(3) contLnID(numTotCont,11) = devsl(1) if ((cl.eq.3.0).or.(cl.eq.4.0)) then contLnID(numTotCont,12) = devsl(2) else contLnID(numTotCont,12) = 0.0 endif endif c ++++++ TEST SECTION +++++++++++++++++++++++++++++++++++++++ c do j=1,Nsl c write (6,*) '=> ',objbuffer(j,1),objbuffer(j,3),objbuffer(j,4), c $ objbuffer(j,2*chosenAp+3),objbuffer(j,2*chosenAp+4) c enddo c write (6,*) '.......... chosen Ap = ',chosenAp c write (6,*) '.......... new avgFlux, dF = ',newFlux,newdF c write (6,*) '....... radius = ',dist,' sigdev=',sigdev c write (6,*) ' ' c +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ endif c ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ enddo write (6,*) ' ' write (6,*) ' ' write (6,*) 'Contin-detect stats:' write (6,*) ' ' write (6,*) 'Num Contin-detect INPUT: ',numCont write (6,*) ' --------------' write (6,*) 'Num objects with RMS < sigdvcut: ',(numRMSrej-twoRMS) write (6,*) 'Num CR rejected: ',(numCRrej-twoCR) write (6,*) ' ' write (6,*) 'Num kept same apert: ',numUncorrCont write (6,*) 'Num with new apert: ',numCorrCont write (6,*) ' --------------' write (6,*) 'Num Contin-detect OUTPUT: ',numTotCont write (6,*) ' ' write (6,*) ' ' write (6,*) 'Num dom by Intr Errors: ',numIntrErr write (6,*) 'Num dom by Meas Errors: ',numMeasErr write (6,*) ' ' write (6,*) ' ' write (6,200) intAp(1),intAp(2),intAp(3),intAp(4) write (6,*) '- - - - - - - - - - - - - - - - - - - - - - - - - -' write (6,200) twoInAp(1),twoInAp(2),twoInAp(3),twoInAp(4), $ numTotTwo write (6,200) numInAp(1),numInAp(2),numInAp(3),numInAp(4), $ numTotCont write (6,*) ' ' write (6,200) (numInAp(1)+twoInAp(1)),(numInAp(2)+twoInAp(2)), $ (numInAp(3)+twoInAp(3)),(numInAp(4)+twoInAp(4)), $ (numTotTwo+numTotCont) write (6,*) '- - - - - - - - - - - - - - - - - - - - - - - - - -' write (6,*) ' ' write (6,*) ' ' 200 format (4(i8), i14) c ---------------------------------------------- c Writing twoLine files: write (6,*) ' ' currFile='twoLine.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numTotTwo,' )' endif open(unit=11, file=currFile, status='unknown') do row=1,numTotTwo outtwo(row,12) = outtwo(row,12)/1e-16 outtwo(row,13) = outtwo(row,13)/1e-16 outtwo(row,16) = outtwo(row,16)/1e-16 outtwo(row,17) = outtwo(row,17)/1e-16 c --- i.e. removing 1e-16 from calibrated fluxes to ease formatting write (11,260) (outtwo(row,col),col=1,18) enddo close(unit=11) c ---------------------------------------------- c Writing contLine files: write (6,*) ' ' currFile='contLine.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numTotCont,' )' endif open(unit=11, file=currFile, status='unknown') do row=1,numTotCont outcont(row,12) = outcont(row,12)/1e-16 outcont(row,13) = outcont(row,13)/1e-16 outcont(row,16) = outcont(row,16)/1e-16 outcont(row,17) = outcont(row,17)/1e-16 c --- i.e. removing 1e-16 from calibrated fluxes to ease formatting write (11,260) (outcont(row,col),col=1,18) enddo close(unit=11) c 220 format (f7.2, f8.2, f8.2, f6.1, f9.2, f7.3, f7.2, f6.2, f7.2, c $ f5.1, f5.0, f8.4, f8.4, f8.2, f8.2, f8.3, f5.1, f5.0) 260 format (f7.2, f8.2, f8.2, f6.1, f9.2, f7.3, f7.2, f6.2, f9.4, $ f7.1, f4.0, f8.3, f7.3, f7.2, f8.2, f8.3, f7.3, f5.0) c ------------------------------------------------------------------ c Writing contLnFlx, contLnUnc, contsigdv and contLnID files: currFile='contLnFlx.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numTotCont,' )' endif open(unit=12, file=currFile, status='unknown') currFile='contLnUnc.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numTotCont,' )' endif open(unit=13, file=currFile, status='unknown') currFile='contsigdv.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numTotCont,' )' endif open(unit=14, file=currFile, status='unknown') currFile='contLnID.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numTotCont,' )' endif open(unit=15, file=currFile, status='unknown') do row=1,numTotCont write (12,440) (contLnFlx(row,col),col=1,(Nsl+3)) write (13,460) (contLnUnc(row,col),col=1,(Nsl+3)) write (14,480) (contsigdv(row,col),col=1,(Nsl+3)) write (15,380) (contLnID(row,col),col=1,12) enddo close(unit=12) close(unit=13) close(unit=14) close(unit=15) 380 format (f8.2, f8.2, f6.0, f14.2, f12.6, f5.0, f14.2, f12.6, $ f5.0, f10.3, f4.0, f4.0) 440 format (f5.0, f8.2, f8.2, 12(f12.3)) 460 format (f5.0, f8.2, f8.2, 12(f12.3)) 480 format (f5.0, f8.2, f8.2, 12(f12.3)) c ---------------------------------------------------------- c Writing rejCRid file: currFile='rejCRid.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numCRrej,' )' endif open(unit=11, file=currFile, status='unknown') do row=1,numCRrej write (11,685) (rejCRid(row,col),col=1,5) enddo close(unit=11) 685 format (f8.2, f8.2, f8.0, f5.0, f5.0) 1000 write (6,*) ' ' write (6,*) ' ' write (6,*) 'Done.' write (6,*) ' ' stop end c ===========================================================================