program scaleflux c Program reads in scale factors fit by TSCALE and output from extspect.f c and applies scale factors to all flux-related data. c The selects ELG candidates on the basis of these scaled fluxes. c c Input/Output format ('Raw' files left untouched) ----------------------- c c outMean(5000,16) (numObj/numStar) c ID X Y ap wvCor ell FWHM strgl Flgs rad NoOb avFlx avdFl S/N rms/dF Cl c 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 c c outFlux(5000,(n+3)) (numObj/numStar) c newID X Y F1 F2 F3 ... Fn c 1 2 3 4 5 6 (n+3) c c outUncert(5000,(n+3)) (numObj/numStar) c newID X Y dF1 dF2 dF3 ... dFn c 1 2 3 4 5 6 (n+3) c c outSN(5000,(n+3)) (numObj/numStar) c newID X Y F1/dF1 F2/dF2 F3/dF3 ... Fn/dFn c 1 2 3 4 5 6 (n+3) c c outsigdv(5000,(n+3)) **NB: normalised by domin. scatter (numFit/numStar) c newID X Y sigdv1 sigdv2 sigdv3 ... sigdvn c 1 2 3 4 5 6 (n+3) c c outID(10000,12) (** note the extra columns) (numObj/numStar) 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 Identical subset files with only the stars (filename prefix 'star'). c c c Additional Output (emission-line samples) ------------------------------- c c elgRed(3000,17) --> contains Flux Devn classes +1 -> +5 (numTotELG) c (NB: cols 5, 10 and 16-17 differ from outMean) c ID X Y ap pkWve el FWH stgl Flgs pkSl NoOb avFl avdFl S/N rms/dF sigdv Cl c 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 c c elgID(3000,3) (subset of elgRed() catalogue) (numTotELG) c X Y sigdv c 1 2 3 c c - - - - - - - - - - - - - - - - - - - - - - - - - - - c c elgLine(3000,18) --> contains Flux Devn classes +1 -> +3,+5 (numFitELG) c (NB: cols 15-18 differ from outMean) c ID X Y ap ... etc ... avFlx avdFl S/N sigdv Fline dFline Cl c 1 2 3 4 12 13 14 15 16 17 18 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 rejCRid(1000,5) c X Y newID frame where c 1 2 3 4 5 c c DHJ 10/1/99 c ------------------------------------------------------------------------ implicit none character*50 currFile, extn, inWaveFile, inCRfile integer col, row, numObj, numStar, Nsl, numFit, rowCR integer adjFlag, numFitELG, numTotELG, starRow, ID, frm, pk integer j, Nmin, rej, Nfit, consecStats(12), flxdnStats(13) integer devsl(2), fullStats(11), redStats(6), Nfit(3) integer numIntrErr, numMeasErr, numCRrej, numCRoff, numTotCR integer numHiSl, numLoSl, consecHi, consecLo, prevHi, prevLo integer hiSigSl, loSigSl, iter, found, verbosity, numCRfit real outMean(5000,16), coeff((12-1),2), centWaveL(12) real outFlux(5000,(12+3)), outUncert(5000,(12+3)) real outsigdv(5000,(12+3)), outSN(5000,(12+3)), outID(5000,12) real starMean(500,16), CRposns(100000,3) real starFlux(500,(12+3)), starUncert(500,(12+3)) real starsigdv(500,(12+3)), starSN(500,(12+3)), starID(500,12) real elgRed(5000,17), elgID(3000,3), elgLine(5000,18) 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 flux, uncert, checkObs, factor, pkSl, pkWave real*4 S, Sx, Sy, Stt, Ffit, wave, sigdvcut, fitTol, freq real pkFlx1, pkFlx2, diff, diff2, k, sumd2, n real domScatt, hiSig, loSig, fl, hiAvgWave, loAvgWave real hiSumLine, loSumLine, hiLine, loLine, dist, x1, y1, xc, yc real hiErr, loErr, fOverlap, fTrough, sgn, sigdvnpt, w1, w2 real exclEnd, exclNonAdj, nonadj, brtLimit 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 scaleflux.in: open(unit=20, file='scaleflux.in', status='old') read(20,*,end=14,err=1) extn read(20,*,end=14,err=2) inWaveFile read(20,*,end=14,err=3) sigdvcut read(20,*,end=14,err=4) fitTol read(20,*,end=14,err=5) 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=11) exclEnd read(20,*,end=14,err=12) exclNonAdj read(20,*,end=14,err=12) brtLimit goto 14 1 write (6,*) ' *** Error Reading scaleflux.in: extn ***' goto 13 2 write (6,*) ' *** Error Reading scaleflux.in: inWaveFile ***' goto 13 3 write (6,*) ' *** Error Reading scaleflux.in: sigdvcut ***' goto 13 4 write (6,*) ' *** Error Reading scaleflux.in: fitTol ***' goto 13 5 write (6,*) ' *** Error Reading scaleflux.in: Nmin ***' goto 13 6 write (6,*) ' *** Error Reading scaleflux.in: rej ***' goto 13 7 write (6,*) ' *** Error Reading scaleflux.in: fOverlap ***' goto 13 8 write (6,*) ' *** Error Reading scaleflux.in: fTrough ***' goto 13 9 write (6,*) ' *** Error Reading scaleflux.in: sgn ***' goto 13 10 write (6,*) ' *** Error Reading scaleflux.in: inCRfile ***' goto 13 11 write (6,*) ' *** Error scaleflux.in: verbosity or exclEnd ***' goto 13 12 write (6,*) ' *** Error scaleflux.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,*) 'Filename extension: ',extn write (6,*) 'Central wavelengths file: ',inWaveFile write (6,*) 'Sigma deviation cut: ',sigdvcut 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='scaleCoeff.'//extn 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 outMean, outFlux and outUncert: currFile='outMean.'//extn row=0 open(unit=20, file=currFile, status='old') 25 row=row+1 read(20,*,end=30,err=30) (outMean(row,col),col=1,14) goto 25 30 continue row=row-1 numObj=row write (6,*) '(There are ',numObj,' objects.)' close(unit=20) currFile='outFlux.'//extn open(unit=20, file=currFile, status='old') currFile='outUncert.'//extn open(unit=21, file=currFile, status='old') do row=1,numObj read(20,*,end=40,err=40) (outFlux(row,col),col=1,(Nsl+3)) read(21,*,end=40,err=40) (outUncert(row,col),col=1,(Nsl+3)) enddo 40 continue close(unit=20) close(unit=21) c ---------------------------------------------- c Reading starMean, starFlux and starUncert: currFile='starMean.'//extn row=0 open(unit=20, file=currFile, status='old') 65 row=row+1 read(20,*,end=70,err=70) (starMean(row,col),col=1,14) goto 65 70 continue row=row-1 numStar=row write (6,*) '(There are ',numStar,' stars.)' close(unit=20) currFile='starFlux.'//extn open(unit=20, file=currFile, status='old') currFile='starUncert.'//extn open(unit=21, file=currFile, status='old') do row=1,numStar read(20,*,end=80,err=80) (starFlux(row,col),col=1,(Nsl+3)) read(21,*,end=80,err=80) (starUncert(row,col),col=1,(Nsl+3)) enddo 80 continue close(unit=20) close(unit=21) 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,*) 'ERROR: Mismatch in coeff and slice no.s!!' 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 Perfoming scaling on each slice: c -------------------------------------------------------------------- write (6,*) ' ' write (6,*) 'Performing scaling... ' c --------------------------------------------------------------------- c Scaling all objects and revising mean flux, uncertainty and S/N: do row=1,numObj do col=1,Nsl c do col=2,Nsl c factor=coeff(col-1,1)+(outMean(row,10)*coeff(col-1,2)) factor=coeff(col,1)+(outMean(row,10)*coeff(col,2)) outFlux(row,col+3)=outFlux(row,col+3)/factor if (outUncert(row,col+3).le.1.0) then outUncert(row,col+3)=0.0 endif outUncert(row,col+3)=outUncert(row,col+3)/factor if (outUncert(row,col+3).eq.0.0) then outUncert(row,col+3)=1.0 endif enddo c --- i.e. col starts at 2 since fluxes normalised to first slice. c --- Special section of code to retain 1.0 for unobserved fluxUncertainties c --- Calculating mean values: flux=0.0 uncert=0.0 checkObs=0.0 do col=1,Nsl flux=flux+outFlux(row,col+3) if (outUncert(row,col+3).gt.1.0) then c --- i.e. need to check for outUncert=1.0 for dud entries uncert=uncert+outUncert(row,col+3) checkObs=checkObs+1.0 endif enddo outMean(row,12)=flux/outMean(row,11) outMean(row,13)=uncert/outMean(row,11) if (checkObs.ne.abs(outMean(row,11))) then write (6,*) 'Obj',outMean(row,1),': numObs do not match!!' write (6,*) ' --> (numObs,checkObs)=(',outMean(row,11), $ checkObs,')' endif outID(row,1)=outMean(row,2) outID(row,2)=outMean(row,3) outID(row,3)=outMean(row,1) c --- copying to ID output do col=1,3 outSN(row,col)=outFlux(row,col) enddo do col=4,(Nsl+3) if (outUncert(row,col).gt.0.0) then outSN(row,col)=outFlux(row,col)/outUncert(row,col) else outSN(row,col)=0.0 endif enddo c --- copying S-N measurements enddo c --------------------------------------------------------------- c Scaling stars and revising mean flux, uncertainty and S/N: do row=1,numStar do col=1,Nsl c do col=2,Nsl c factor=coeff(col-1,1)+(starMean(row,10)*coeff(col-1,2)) factor=coeff(col,1)+(outMean(row,10)*coeff(col,2)) starFlux(row,col+3)=starFlux(row,col+3)/factor starUncert(row,col+3)=starUncert(row,col+3)/factor enddo c --- i.e. col starts at 2 since fluxes normalised to first slice. c --- Special section of code to retain 1.0 for unobserved fluxUncertainties c --- Calculating mean values: flux=0.0 uncert=0.0 checkObs=0.0 do col=1,Nsl flux=flux+starFlux(row,col+3) if (starUncert(row,col+3).gt.1.0) then c --- i.e. need to check for starUncert=1.0 for dud entries uncert=uncert+starUncert(row,col+3) checkObs=checkObs+1 endif enddo starMean(row,12)=flux/starMean(row,11) starMean(row,13)=uncert/starMean(row,11) if (checkObs.ne.abs(starMean(row,11))) then write (6,*) 'Star',starMean(row,1),': numObs do not match!!' write (6,*) ' --> (numObs,checkObs)=(',starMean(row,11), $ checkObs,')' endif do col=1,3 starSN(row,col)=starFlux(row,col) enddo do col=4,(Nsl+3) if (starUncert(row,col).gt.0.0) then starSN(row,col)=starFlux(row,col)/starUncert(row,col) else outSN(row,col)=0.0 endif enddo c --- copying S-N measurements enddo c -------------------------------------------------------------------- c Selecting emission-line candidates: c -------------------------------------------------------------------- if (verbosity.eq.2) then write (6,*) ' ' endif write (6,*) 'Selecting emission-line candidates: ' numFit=0 numFitELG=0 numTotELG=0 numCRrej=0 numCRoff=0 numCRfit=0 do col=1,11 fullStats(col)=0 if ((col.ge.1).and.(col.le.6)) then redStats(col)=0 endif enddo do col=1,Nsl consecStats(col)=0 flxdnStats(col)=0 enddo flxdnStats(Nsl+1)=0 c --- i.e. clearing stats arrays (of differing length) c ------------------------------------------------------------------- c 1st pass -- selecting all the 3 adjacent-slice ELG detections: if (verbosity.eq.2) then write (6,*) ' --> the 3 adjacent-slice candidates... ' endif do row=1,numObj if (outMean(row,11).eq.3) then c --- i.e. if object has only 3 detections (anywhere) adjFlag=0 pkSl=0.0 pkWave=0.0 c --- adjFlag=0 for no observations; 1 for three adjacent col=2 do while ((col.lt.Nsl).and.(adjFlag.eq.0)) factor=outFlux(row,col+4)*outFlux(row,col+3)*outFlux(row,col+2) if (factor.gt.0.0) then adjFlag=1 c pkSl=(float(col)+float(col-1))/2.0; 2-adj case pkSl=float(col) c pkWave=(centWaveL(col)+centWaveL(col-1))/2.0; 2-adj case pkWave=centWaveL(col) pkWave=(1-outMean(row,5))*pkWave c --- i.e. applying waveL correction factor to mean waveL endif col=col+1 enddo c --- i.e. tests adjacency by multiplying adjacent fluxes c --- Checking against CR list if object meets adjacency criterion: if (adjFlag.eq.1) then xc=outMean(row,2) yc=outMean(row,3) pk=int(pkSl) 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*outMean(row,4))) then c --- i.e. CR pixel falls within aperture radius c write (6,*) ' *** Obj',row,' hits CR pix',x1,y1,', but not frame' numCRoff=numCRoff+1 found=1 if ((frm.eq.pk).or.(frm.eq.(pk+1)).or.(frm.eq.(pk-1))) then found=2 endif endif rowCR=rowCR+1 enddo if (found.eq.2) then numCRrej=numCRrej+1 rejCRid(numCRrej,1)=outMean(row,2) rejCRid(numCRrej,2)=outMean(row,3) rejCRid(numCRrej,3)=outMean(row,1) rejCRid(numCRrej,4)=float(frm) rejCRid(numCRrej,5)=1.0 else c --- genuine EL object numTotELG=numTotELG+1 outMean(row,16)=4.0 elgRed(numTotELG,17)=4.0 fullStats(5)=fullStats(5)+1 redStats(5)=redStats(5)+1 consecStats(int(pkSl))=consecStats(int(pkSl))+1 do col=1,14 elgRed(numTotELG,col)=outMean(row,col) enddo c --- Writing peak location (slicenumber and wavelength): elgRed(numTotELG,10)=pkSl elgRed(numTotELG,5)=pkWave elgRed(numTotELG,15)=0.0 elgID(numTotELG,1)=outMean(row,2) elgID(numTotELG,2)=outMean(row,3) elgID(numTotELG,3)=0.0 endif else c --- non-consec 2-detection outMean(row,16)=10.0 fullStats(11)=fullStats(11)+1 c &&&& endif endif c consecStats(int(pkSl))=consecStats(int(pkSl))+1 endif enddo c ----------------------------------------------------------------------- c --- Continuum-Fitting Routine: corrected set of wavelength values used c --- (zero flux columns excluded from the fits) if (verbosity.eq.2) then write (6,*) ' --> weighted least-squares continuum fit...' write (6,*) ' --> (rejecting up to ',rej,' pts maximum)' endif numIntrErr=0 numMeasErr=0 do row=1,numObj 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 if ((outMean(row,11).ge.float(Nmin)).and. $ (outMean(row,8).ge.0.0)) then c --- i.e. checking that object does not flux-spill (strgal<0) c ------------------------------------------------------------ c --- 1. 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(row,col+3).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-outMean(row,5))*centWaveL(col) F(col) = outFlux(row,col+3) delF(col) = outUncert(row,col+3) 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.outMean(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,outMean(row,11) endif c --- Loop 2: -------------------------------------------- checkObs=0.0 do col=1,Nsl if (outFlux(row,col+3).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.outMean(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(row,col+3).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-outMean(row,5))*centWaveL(col) Ffit = C(iter)+(wave*M(iter)) c --- i.e. y = C+Mx diff2 = (outFlux(row,col+3)-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 (outMean(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=outMean(row,13) c --- i.e. intrinsic errors dominate endif do col=1,Nsl if (outFlux(row,col+3).gt.0.0) then c --- i.e. checking for a flux detection wave = (1-outMean(row,5))*centWaveL(col) Ffit = C(iter)+(wave*M(iter)) c --- i.e. y = C+Mx diff = sgn*(outFlux(row,col+3)-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 ((outMean(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(row,col+3).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-outMean(row,5))*centWaveL(col) F(col) = outFlux(row,col+3) delF(col) = outUncert(row,col+3) 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.outMean(row,11)) then write (6,*) 'ERROR: Mismatch in no. of obsvns!!' write (6,*) 'checkObs,outMean =',checkObs,outMean(row,11) endif c --- Loop 2: --------------------------------- checkObs=0.0 do col=1,Nsl if (outFlux(row,col+3).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.outMean(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 (outMean(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(row,col+3).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-outMean(row,5))*centWaveL(col) Ffit = C(3)+(wave*M(3)) c --- i.e. y = C+Mx diff2 = (outFlux(row,col+3)-Ffit)**2 sumd2=sumd2+diff2 n=n+1.0 endif endif enddo rmsObj(3)=sqrt(sumd2/(n-1.0)) rmsErrRat(3)=rmsObj(3)/outMean(row,13) c write (6,*) ' 3:' c write (6,*) ' M, C:',M(3),C(3),Nfit(3) c write (6,*) ' rms: ',rmsObj(3),rmsErrRat(3) c write (6,*) ' ' endif c --- end of min no. of pts (for *any* fit) criterion c ---------------------------------------------------- c --- 6. Copying fit results to output files: c --- (objects with too few fit pts already taken care of) outID(row,4) = C(1) outID(row,5) = M(1) outID(row,6) = float(Nfit(1)) outID(row,7) = C(3) outID(row,8) = M(3) outID(row,9) = float(Nfit(3)) outID(row,10) = rmsErrRat(3) outID(row,11) = devsl(1) outID(row,12) = devsl(2) outMean(row,15) = rmsErrRat(3) if (outID(row,4).ne.0.0) then c --- i.e. object has been fit if non-zero fit value if (rmsErrRat(3).gt.1.0) then numIntrErr=numIntrErr+1 else if (rmsErrRat(3).gt.0.0) then numMeasErr=numMeasErr+1 endif endif enddo c --- end of main object loop c ----------------------------------------------------------------------- c Calculating sigma devn at each wavelength and writing to outsigdv() if (verbosity.eq.2) then write (6,*) ' --> Sigma deviation at each wavelength... ' endif j=0 c --- i.e. row counter through outsigdv (=numFit at end) starRow=1 do row=1,numObj if (outID(row,4).ne.0.0) then c --- i.e. object has been fit if non-zero fit value j=j+1 outsigdv(j,1)=outMean(row,1) outsigdv(j,2)=outMean(row,2) outsigdv(j,3)=outMean(row,3) do col=1,Nsl wave = (1-outMean(row,5))*centWaveL(col) Ffit = outID(row,7)+(wave*outID(row,8)) c --- i.e. y = a+bx rmsErrRat(3)=outMean(row,15) rmsObj(3)=rmsErrRat(3)*outMean(row,13) c --- i.e. recalculating rmsObj from its ratio and dF if (rmsErrRat(3).gt.1.0) then domScatt=rmsObj(3) c --- i.e. measurement errors dominate the continuum scatter else domScatt=outMean(row,13) c --- i.e. intrinsic errors dominate the continuum scatter endif if (outFlux(row,col+3).gt.0.0) then outsigdv(j,col+3)=(outFlux(row,col+3)-Ffit)/domScatt else outsigdv(j,col+3)=0.0 endif c --- i.e. outsigdv=0.0 if no flux measurement at that point enddo if (starMean(starRow,1).eq.outsigdv(j,1)) then do col=1,(Nsl+3) starsigdv(starRow,col) = outsigdv(j,col) enddo starRow=starRow+1 endif c --- i.e. copying current entry to starsigdv if one of the stars; endif enddo numFit=j if (starRow.ne.(numStar+1)) then write (6,*) ' ** ERROR: Not all stars relocated !!' write (6,*) ' starRow,(numStar+1) =',starRow,(numStar+1) endif c ------------------------------------------------------------------- c 2nd pass -- selecting all the fit-deviating ELG candidates: if (verbosity.eq.2) then write (6,*) ' --> objects deviating from the fit...' endif k=0 c --- counter through the fit ELG sample (=numFitELG at end) do row=1,numFit ID=int(outsigdv(row,1)) 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(row,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-outMean(ID,5))*centWaveL(col) hiWave(numHiSl)=wave Ffit=outID(ID,7)+(outID(ID,8)*wave) hiFlx(numHiSl)=outFlux(ID,col+3)-Ffit hidF(numHiSl)=outUncert(ID,col+3) if (outsigdv(row,col+3).gt.hiSig) then c --- i.e. checking if maximum sigma deviation hiSig=outsigdv(row,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(row,col+3).le.(-1*sigdvcut)) then numLoSl=numLoSl+1 if (numLoSl.le.3) then loSl(numLoSl)=col wave=(1-outMean(ID,5))*centWaveL(col) loWave(numLoSl)=wave Ffit=outID(ID,7)+(outID(ID,8)*wave) loFlx(numLoSl)=outFlux(ID,col+3)-Ffit lodF(numLoSl)=outUncert(ID,col+3) if (outsigdv(row,col+3).lt.loSig) then loSig=outsigdv(row,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 c if ((numHiSl.ne.0).or.(numLoSl.ne.0)) then c write (6,*) ID,' (',row,') - - - - - - - - - - - - - - - - - ' c write (6,*) ' ',numHiSl,':',(hiFlx(col),col=1,3) c write (6,*) ' ',hiSig,':',(hiWave(col),col=1,3) c write (6,*) ' ',(hiFlx(col),col=1,3) c write (6,*) ' ',(hidF(col),col=1,3) c write (6,*) ' -----> ',hiSumLine, hiAvgWave c write (6,*) ' =====> ',hiLine, hiErr,' (',hiSigSl,')' c write (6,*) ' ' c write (6,*) ' ',numLoSl,':',(loFlx(col),col=1,3) c write (6,*) ' ',loSig,':',(loWave(col),col=1,3) c write (6,*) ' -----> ',loSumLine, loAvgWave c write (6,*) ' =====> ',loLine, loErr,' (',loSigSl,')' c write (6,*) '(consecLo, consecHi): ',consecLo, consecHi c endif if (outID(ID,12).eq.0.0) then nonadj=1.0 c --- i.e. non-adj of double-devn irrelevant when only 1 devn; c --- set to 1.0 anyway. else nonadj=abs(outID(ID,11)-outID(ID,12)) c --- i.e. non-adj of double-devn occurs if nonadj.ne.1.0 endif c --- Classifying nature of high/low deviant sets: c ---- Emission Cases: if ((hiSig.gt.0.0).and.(loSig.eq.0.0)) then 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.(outMean(ID,12).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: numTotELG=numTotELG+1 k=k+1 outMean(ID,16)=1.0 elgRed(numTotELG,17)=1.0 elgLine(k,18)=1.0 fullStats(2)=fullStats(2)+1 redStats(2)=redStats(2)+1 elgRed(numTotELG,5)=hiAvgWave elgLine(k,16)=hiLine/fTrough elgLine(k,17)=hiErr/fTrough else if (numHiSl.eq.2) then if (consecHi.eq.2) then c --- class +2 devn: numTotELG=numTotELG+1 k=k+1 outMean(ID,16)=2.0 elgRed(numTotELG,17)=2.0 elgLine(k,18)=2.0 fullStats(3)=fullStats(3)+1 redStats(3)=redStats(3)+1 elgRed(numTotELG,5)=hiAvgWave elgLine(k,16)=hiLine/fOverlap elgLine(k,17)=hiErr/fOverlap else c --- class +7 devn: (excluded) outMean(ID,16)=7.0 fullStats(8)=fullStats(8)+1 endif else if (numHiSl.eq.3) then if (consecHi.eq.3) then if (hiSigSl.eq.hiSl(2)) then c --- class +3 devn: numTotELG=numTotELG+1 k=k+1 outMean(ID,16)=3.0 elgRed(numTotELG,17)=3.0 elgLine(k,18)=3.0 fullStats(4)=fullStats(4)+1 redStats(4)=redStats(4)+1 elgRed(numTotELG,5)=hiAvgWave elgLine(k,16)=hiLine/fTrough elgLine(k,17)=hiErr/fTrough else c --- class +9 devn: (excluded) outMean(ID,16)=9.0 fullStats(10)=fullStats(10)+1 endif else c --- class +7 devn again: (excluded) outMean(ID,16)=7.0 fullStats(8)=fullStats(8)+1 endif else if (numHiSl.eq.0) then c --- class 0 devn: (excluded) outMean(ID,16)=0.0 fullStats(1)=fullStats(1)+1 else c --- class +8 devn: (excluded) outMean(ID,16)=8.0 fullStats(9)=fullStats(9)+1 endif endif c --- end of brtLimit check if endif c --- end of exclNonAdj check if endif c --- end of exclEndpt check if c ---- Absorption Cases: ---------------------------------------- else if ((hiSig.eq.0.0).and.(loSig.lt.0.0)) then if ((((exclEnd.eq.1.0).and.(loSigSl.gt.1).and. $ (loSigSl.lt.Nsl))).or.(exclEnd.eq.0.0)) then 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.(outMean(ID,12).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 (numLoSl.eq.1) then c --- class -1 devn: outMean(ID,16)=-1.0 fullStats(2)=fullStats(2)+1 c absRed(numTotELG,5)=loAvgWave c absLine(k,16)=loLine/fTrough c absLine(k,17)=loErr/fTrough c &&&& NO FLUX/WAVELENGTH CALCULATION FOR ABSORPTION CASE YET &&&& else if (numLoSl.eq.2) then if (consecLo.eq.2) then c --- class -2 devn: outMean(ID,16)=-2.0 fullStats(3)=fullStats(3)+1 c absRed(numTotELG,5)=loAvgWave c absLine(k,16)=loLine/fOverlap c absLine(k,17)=loErr/fOverlap c &&&& NO FLUX/WAVELENGTH CALCULATION FOR ABSORPTION CASE YET &&&& else c --- class -7 devn: (excluded) outMean(ID,16)=-7.0 fullStats(8)=fullStats(8)+1 endif else if (numLoSl.eq.3) then if (consecLo.eq.3) then if (loSigSl.eq.LoSl(2)) then c --- class -3 devn: outMean(ID,16)=-3.0 fullStats(4)=fullStats(4)+1 c absRed(numTotELG,5)=loAvgWave c absLine(k,16)=loAvgLine/fTrough c absLine(k,17)=loAvgLine/fTrough c &&&& NO FLUX/WAVELENGTH CALCULATION FOR ABSORPTION CASE YET &&&& else c --- class -9 devn: (excluded) outMean(ID,16)=-9.0 fullStats(10)=fullStats(10)+1 endif else c --- class -7 devn again: (excluded) outMean(ID,16)=-7.0 fullStats(8)=fullStats(8)+1 endif else if (numLoSl.eq.0) then c --- class 0 devn: (excluded) outMean(ID,16)=0.0 fullStats(1)=fullStats(1)+1 else c --- class -8 devn: (excluded) outMean(ID,16)=-8.0 fullStats(9)=fullStats(9)+1 endif c --- Simultaneous emission and absorption: else if ((hiSig.gt.0.0).and.(loSig.lt.0.0)) then if (hiSig.ge.(abs(loSig))) then c --- class 5 devn: numTotELG=numTotELG+1 k=k+1 outMean(ID,16)=5.0 elgRed(numTotELG,17)=5.0 elgLine(k,18)=5.0 fullStats(6)=fullStats(6)+1 redStats(6)=redStats(6)+1 elgRed(numTotELG,5)=hiAvgWave if ((numHiSl.eq.1).or.(numHiSl.eq.3)) then elgLine(k,16)=hiLine/fTrough elgLine(k,17)=hiErr/fTrough else if (numHiSl.eq.2) then elgLine(k,16)=hiLine/fOverlap elgLine(k,17)=hiErr/fOverlap else elgLine(k,16)=0.0 elgLine(k,17)=0.0 endif else c --- class 6 devn (excluded) outMean(ID,16)=6.0 fullStats(7)=fullStats(7)+1 c absRed(numTotELG,5)=loAvgWave c if ((numLoSl.eq.1).or.(numLoSl.eq.3)) then c absLine(k,16)=loLine/fTrough c absLine(k,17)=loErr/fTrough c else if (numHiSl.eq.2) then c absLine(k,16)=loLine/fOverlap c absLine(k,17)=loErr/fOverlap c else c absLine(k,16)=0.0 c absLine(k,17)=0.0 c endif c &&&& NO FLUX/WAVELENGTH CALCULATION FOR ABSORPTION CASE YET &&&& endif endif c --- end of brtLimit check if endif c --- end of exclNonAdj check if endif c --- end of exclEndpt check if c --- No deviation: else c --- class 0 (no devn): outMean(ID,16)=0.0 fullStats(1)=fullStats(1)+1 redStats(1)=redStats(1)+1 endif c --- Checking against CR list if an emission-line suspect: fl=outMean(ID,16) if ( ((fl.ge.1).and.(fl.le.3)).or.(fl.eq.5) ) then c --- this excludes 2-detection ELGs since already written c --- Checking against CR list: xc=outMean(ID,2) yc=outMean(ID,3) 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*outMean(row,4))) then c --- i.e. CR pixel falls within aperture radius c write (6,*) ' *** Obj',row,' hits CR pix',x1,y1,', but not frame' numCRoff=numCRoff+1 found=1 do col=1,numHiSl if (frm.eq.(hiSl(col))) then found=2 endif enddo endif rowCR=rowCR+1 enddo if (found.eq.2) then numTotELG=numTotELG-1 k=k-1 c --- i.e. removing current EL cand from the count c --- (and hence, previously written data from ELG file) numCRrej=numCRrej+1 numCRfit=numCRfit+1 rejCRid(numCRrej,1)=outMean(ID,2) rejCRid(numCRrej,2)=outMean(ID,3) rejCRid(numCRrej,3)=outMean(ID,1) rejCRid(numCRrej,4)=float(frm) rejCRid(numCRrej,5)=2.0 else c --- i.e. genuine EL object c --- Writing EL cand information to file: do col=1,14 if (col.ne.5) then elgRed(numTotELG,col)=outMean(ID,col) endif c --- i.e. do not copy column 5 as this will overwrite c --- what has been written about the line (immed. above) elgLine(k,col)=outMean(ID,col) if (col.le.12) then elgLnID(k,col)=outID(ID,col) endif enddo elgRed(numTotELG,15)=outMean(ID,15) elgRed(numTotELG,10)=hiSigSl elgRed(numTotELG,16)=hiSig elgID(numTotELG,1)=outMean(ID,2) elgID(numTotELG,2)=outMean(ID,3) elgID(numTotELG,3)=hiSig elgLine(k,15)=hiSig do col=1,(Nsl+3) elgLnFlx(k,col)=outFlux(ID,col) elgLnUnc(k,col)=outUncert(ID,col) elgsigdv(k,col)=outsigdv(row,col) enddo endif endif c --- i.e. end of EL-suspect routine enddo c --- i.e end of main object loop numFitELG=k c ---------------------------------------------------------------- c Compiling histogram of emission wavelengths flxdnStats(); c (only from the sample of fit ELG candidates): k=numTotELG-numFitELG do row=1,numFitELG col=0 found=0 c --- flag: found=0 if wavelength bin yet to be located; =1 if has do while ((col.le.Nsl).and.(found.eq.0)) if (col.eq.0) then c -- first slice is a special EXTRA bin towards the blue, c -- to catch emission phase shifted blueward of centre wave. c -- (NOTE: no lower limit to bin.) w1=0 w2=centWaveL(2)-(1.5*(centWaveL(2)-centWaveL(1))) else if ((col.gt.0).and.(col.lt.Nsl)) then w1=centWaveL(col)-(0.5*(centWaveL(col+1)-centWaveL(col))) w2=centWaveL(col)+(0.5*(centWaveL(col+1)-centWaveL(col))) else if (col.eq.Nsl) then w1=centWaveL(col)-(0.5*(centWaveL(col)-centWaveL(col-1))) w2=centWaveL(col)+(0.5*(centWaveL(col)-centWaveL(col-1))) endif if((elgRed(k+row,5).gt.w1).and.(elgRed(k+row,5).le.w2))then found=1 flxdnStats(col+1)=flxdnStats(col+1)+1 else col=col+1 endif enddo if ((col.ge.(Nsl+1)).and.(found.eq.0)) then c -- checking for objects not falling within any of the bins write (6,*) '*** Obj ',elgRed(k+row,1),' misses histogram' endif enddo c ---------------------------------------------------------------- c Copying outID() contents to starID() for star entries; c (also writing values to starMean(row,15)(row,16) columns): starRow=1 do row=1,numObj if (starMean(starRow,1).eq.outMean(row,1)) then do col=1,12 starID(starRow,col)=outID(row,col) enddo starMean(starRow,15)=outMean(row,15) starMean(starRow,16)=outMean(row,16) starRow=starRow+1 endif enddo if (starRow.ne.(numStar+1)) then write (6,*) ' ** ERROR: Not all stars relocated !!' write (6,*) ' starRow,(numStar+1) =',starRow,(numStar+1) endif c -------------------------------------------------------------------- c Writing output files (20 of them): c -------------------------------------------------------------------- if (verbosity.eq.2) then write (6,*) ' ' endif write (6,*) 'Writing to files...' c ------------------------------------------------------------- c Writing outMean2, outFlux2, outUncert2, outID2 and outSN2: currFile='outMean2.'//extn if (verbosity.eq.2) then write (6,*) ' ' write (6,*) ' ',currFile,'(',numObj,' )' endif open(unit=11, file=currFile, status='unknown') currFile='outFlux2.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numObj,' )' endif open(unit=12, file=currFile, status='unknown') currFile='outUncert2.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numObj,' )' endif open(unit=13, file=currFile, status='unknown') currFile='outID2.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numObj,' )' endif open(unit=14, file=currFile, status='unknown') currFile='outSN2.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numObj,' )' endif open(unit=15, file=currFile, status='unknown') do row=1,numObj outMean(row,16)=outMean(row,16)+1.0 c --- since these class no.s have increased by 1 since first writing write (11,120) (outMean(row,col),col=1,16) write (12,140) (outFlux(row,col),col=1,(Nsl+3)) write (13,160) (outUncert(row,col),col=1,(Nsl+3)) write (14,170) (outID(row,col),col=1,12) write (15,180) (outSN(row,col),col=1,(Nsl+3)) enddo close(unit=11) close(unit=12) close(unit=13) close(unit=14) close(unit=15) c ------------------------------------------------ c Writing outsigdv (of different length): currFile='outsigdv.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numFit,' )' endif open(unit=11, file=currFile, status='unknown') do row=1,numFit write (11,180) (outsigdv(row,col),col=1,(Nsl+3)) enddo close(unit=11) c ----------------------------------------------------------------- c Writing starMean2, starFlux2, starUncert2, starID2, starSN2 c and starsigdv: write (6,*) ' ' currFile='starMean2.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numStar,' )' endif open(unit=11, file=currFile, status='unknown') currFile='starFlux2.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numStar,' )' endif open(unit=12, file=currFile, status='unknown') currFile='starUncert2.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numStar,' )' endif open(unit=13, file=currFile, status='unknown') currFile='starID2.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numStar,' )' endif open(unit=14, file=currFile, status='unknown') currFile='starSN2.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numStar,' )' endif open(unit=15, file=currFile, status='unknown') currFile='starsigdv.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numStar,' )' endif open(unit=16, file=currFile, status='unknown') do row=1,numStar starMean(row,16)=starMean(row,16)+1.0 c --- since these class no.s have increased by 1 since first writing write (11,120) (starMean(row,col),col=1,16) write (12,140) (starFlux(row,col),col=1,(Nsl+3)) write (13,160) (starUncert(row,col),col=1,(Nsl+3)) write (14,170) (starID(row,col),col=1,12) write (15,180) (starSN(row,col),col=1,(Nsl+3)) write (16,180) (starsigdv(row,col),col=1,(Nsl+3)) enddo close(unit=11) close(unit=12) close(unit=13) close(unit=14) close(unit=15) close(unit=16) c --- i.e. new star catalogues 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, f8.3, f5.0) 140 format (f5.0, f8.2, f8.2, 12(f12.3)) 160 format (f5.0, f8.2, f8.2, 12(f12.3)) 170 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) 180 format (f5.0, f8.2, f8.2, 12(f12.3)) c 200 format (f5.0, f8.2, f8.2, f8.4, f8.2, 12(f8.3)) c ---------------------------------------------- c Writing elgRed and elgID files: write (6,*) ' ' currFile='elgRed.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numTotELG,' )' endif open(unit=11, file=currFile, status='unknown') currFile='elgID.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numTotELG,' )' endif open(unit=12, file=currFile, status='unknown') do row=1,numTotELG elgRed(row,17)=elgRed(row,17)+1.0 c --- since these class no.s have increased by 1 since first writing write (11,220) (elgRed(row,col),col=1,17) write (12,240) (elgID(row,col),col=1,3) enddo close(unit=11) close(unit=12) 220 format (f5.0, f8.2, f8.2, f6.1, f9.2, f7.3, f7.2, f6.2, f14.0, $ f5.1, f5.0, f12.2, f8.2, f8.2, f8.2, f8.3, f5.1, f5.0) 240 format (f8.2, f8.2, f10.3) c ---------------------------------------------------------- c Writing elgLine, elgLnFlx, elgLnUnc, elgsigdv and c elgLnID files: currFile='elgLine.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numFitELG,' )' endif open(unit=11, file=currFile, status='unknown') currFile='elgLnFlx.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numFitELG,' )' endif open(unit=12, file=currFile, status='unknown') currFile='elgLnUnc.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numFitELG,' )' endif open(unit=13, file=currFile, status='unknown') currFile='elgsigdv.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numFitELG,' )' endif open(unit=14, file=currFile, status='unknown') currFile='elgLnID.'//extn if (verbosity.eq.2) then write (6,*) ' ',currFile,'(',numFitELG,' )' endif open(unit=15, file=currFile, status='unknown') do row=1,numFitELG elgLine(row,18)=elgLine(row,18)+1.0 c --- since these class no.s have increased by 1 since first writing write (11,260) (elgLine(row,col),col=1,18) write (12,140) (elgLnFlx(row,col),col=1,(Nsl+3)) write (13,160) (elgLnUnc(row,col),col=1,(Nsl+3)) write (14,180) (elgsigdv(row,col),col=1,(Nsl+3)) write (15,280) (elgLnID(row,col),col=1,12) enddo close(unit=11) close(unit=12) close(unit=13) close(unit=14) close(unit=15) 260 format (f5.0, f8.2, f8.2, f6.1, f12.8, f7.3, f7.2, f6.2, f14.0, $ f6.0, f4.0, f12.2, f8.2, f8.2, f7.2, f12.2, f8.2, f4.0) 280 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) 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,285) (rejCRid(row,col),col=1,5) enddo close(unit=11) 285 format (f8.2, f8.2, f8.0, f5.0, f5.0) c ---------------------------------------------------------- c Writing scaleflux.out (break down of numbers): currFile='scaleflux.out' if (verbosity.eq.2) then write (6,*) ' ' write (6,*) ' ',currFile endif open(unit=11, file=currFile, status='unknown') write (11,*) 'Filename extension: ',extn write (11,*) 'Central wavelengths file: ',inWaveFile write (11,*) 'Sigma deviation cut: ',sigdvcut write (11,*) 'Initial fit precision ' write (11,*) '(as a % of relative flux uncertainty): ',fitTol write (11,*) 'Minimum number of fit pts (all objects): ',Nmin write (11,*) 'Maximum number of rejected pts: ',rej write (11,*) 'Flux scale factor (overlap regions): ',fOverlap write (11,*) 'Flux scale factor (centre trough region): ',fTrough write (11,*) 'Side of fit to reject points from (+/-1.0): ',sgn write (11,*) 'Cosmic ray/ghost file: ',inCRfile write (11,*) 'Output verbosity: ',verbosity write (11,*) 'Exclude end-point devns (1/0=yes/no): ',exclEnd write (11,*) ' ' write (11,*) ' ' write (11,*) 'Total Objects: ',numObj,' (',numStar,' stars)' j=int(100*(float(numFit)/float(numObj))) write (11,*) 'Number fit: ',numFit,' (',j,'% total)' j=int(100*(float(numIntrErr)/float(numFit))) write (11,*) 'Dom intr/meas: ',numIntrErr,',',numMeasErr, $ ' (',j,'%,',(100-j),'% those fit)' j=int(100*(float(numTotELG)/float(numObj))) write (11,*) 'Total ELGs: ',numTotELG,' (',j,'% total)' j=int(100*(float(numFitELG)/float(numObj))) write (11,*) 'Fit ELGs: ',numFitELG,' (',j,'% total)' write (11,*) ' ' write (11,*) ' Tot EL reject: ',numCRrej j=int(100*(float(numCRfit)/float(numCRrej))) write (11,*) ' Fit EL reject: ',numCRfit,' (',j,'% tot CR)' write (11,*) ' CR mismatches: ',numCRoff write (11,*) ' ' write (11,*) '- - - - - - - - - - - - - - - - - - - - - - - - - -' write (11,300) (col, col=1,11) write (11,*) '- - - - - - - - - - - - - - - - - - - - - - - - - -' write (11,300) (fullStats(col), col=1,11) write (11,300) (redStats(col), col=1,6) write (11,300) ((fullStats(col)-redStats(col)), col=1,6) write (11,*) ' ' write (11,*) ' 1: no flux deviations at all from a fit' write (11,*) ' 2: 1 deviant band from a fit' write (11,*) ' 3: 2 consecutive deviants from a fit' write (11,*) ' 4: 3 consecutive deviants from a fit' write (11,*) ' 5: 3 consecutive detections and no others (no fit)' write (11,*) ' 6: both + & -ve devns in which the +ve is stronger' write (11,*) ' 7: both + & -ve devns in which the -ve is stronger' write (11,*) ' 8: 2 or more non-consecutive deviants from a fit' write (11,*) ' 9: 4 or more deviants (non-/consecutive) from fit' write (11,*) '10: 3 consecutive fit deviants with non-central max' write (11,*) '11: 3 non-consecutive detections (no fit)' write (11,*) ' All columns above count CRs except 5' write (11,*) ' All columns below exclude CRs' write (11,*) ' ' write (11,*) '- - - - - - - - - - - - - - - - - - - - - - - - - -' write (11,290) ((col+0.5),col=1,(Nsl-1)) write (11,295) (consecStats(col),col=1,(Nsl-1)) write (11,*) ' ==> ',fullStats(5) c --- this is the range of col for 2-consec case write (11,*) ' ' w1=centWaveL(1)-(centWaveL(2)-centWaveL(1)) write (11,290) (w1/10.0), (centWaveL(col)/10.0,col=1,Nsl) write (11,295) (flxdnStats(col),col=1,(Nsl+1)) write (11,*) ' ==> ',numFitELG write (11,*) '- - - - - - - - - - - - - - - - - - - - - - - - - -' close(unit=11) c --- Writing to terminal: if (verbosity.eq.2) then write (6,*) ' ' write (6,*) ' ' endif write (6,*) ' Total Objects: ',numObj,' (',numStar,' stars)' j=int(100*(float(numFit)/float(numObj))) write (6,*) ' Number fit: ',numFit,' (',j,'% total)' j=int(100*(float(numIntrErr)/float(numFit))) write (6,*) ' Dom intr/meas: ',numIntrErr,',',numMeasErr, $ ' (',j,'%,',(100-j),'% those fit)' j=int(100*(float(numTotELG)/float(numObj))) write (6,*) ' Total ELGs: ',numTotELG,' (',j,'% total)' j=int(100*(float(numFitELG)/float(numObj))) write (6,*) ' Fit ELGs: ',numFitELG,' (',j,'% total)' write (6,*) ' ' write (6,*) ' Tot EL reject: ',numCRrej if (numCRrej.gt.0.0) then j=int(100*(float(numCRfit)/float(numCRrej))) else j = 0 endif write (6,*) ' Fit EL reject: ',numCRfit,' (',j,'% tot CR)' write (6,*) ' CR mismatches: ',numCRoff write (6,*) ' ' write (6,*) ' Minimum fit pts (all objects): ',Nmin write (6,*) ' Maximum no. rejected pts: ',rej write (6,*) ' ' write (6,*) '- - - - - - - - - - - - - - - - - - - - - - - - - -' write (6,300) (col, col=1,11) write (6,*) '- - - - - - - - - - - - - - - - - - - - - - - - - -' write (6,300) (fullStats(col), col=1,11) write (6,300) (redStats(col), col=1,6) write (6,300) ((fullStats(col)-redStats(col)), col=1,6) write (6,*) ' ' if (verbosity.eq.2) then write (6,*) ' 1: no flux deviations at all from a fit' write (6,*) ' 2: 1 deviant band from a fit' write (6,*) ' 3: 2 consecutive deviants from a fit' write (6,*) ' 4: 3 consecutive deviants from a fit' write (6,*) ' 5: 3 consecutive detections and no others (no fit)' write (6,*) ' 6: both + & -ve devns in which the +ve is stronger' write (6,*) ' 7: both + & -ve devns in which the -ve is stronger' write (6,*) ' 8: 2 or more non-consecutive deviants from a fit' write (6,*) ' 9: 4 or more deviants (non-/consecutive) from fit' write (6,*) '10: 3 consecutive fit deviants with non-central max' write (6,*) '11: 3 non-consecutive detections (no fit)' write (6,*) ' All columns above count CRs except 5' write (6,*) ' All columns below exclude CRs' endif write (6,*) ' ' write (6,*) '- - - - - - - - - - - - - - - - - - - - - - - - - -' write (6,290) ((col+0.5),col=1,(Nsl-1)) write (6,295) (consecStats(col),col=1,(Nsl-1)) write (6,*) ' ==> ',fullStats(5) c --- this is the range of col for 2-consec case write (6,*) ' ' w1=centWaveL(1)-(centWaveL(2)-centWaveL(1)) write (6,290) (w1/10.0), (centWaveL(col)/10.0,col=1,Nsl) write (6,295) (flxdnStats(col),col=1,(Nsl+1)) write (6,*) ' ==> ',numFitELG write (6,*) '- - - - - - - - - - - - - - - - - - - - - - - - - -' 290 format (11(f7.1)) 295 format (11(i7)) 300 format (11(i5)) write (6,*) ' ' write (6,*) ' ' write (6,*) 'Done.' write (6,*) ' ' stop end c ===========================================================================