;+ ; PURPOSE: ; Code to accompany script ghostholes_ind ; ; DESCRIPTION: ; Takes and source removed event file, and a source list ; Fills the gaps in the event list, by sampling and copy ; events around each circular gap, randomising the positions ; in DETX and DETY. The SAS task attcalc must be used ; afterwards (in the accompanying script ghostholes_ind) ; to readjust the positions in X and Y (sky coordinates). ; ; INPUT: ; Environment variables as set by above code ; ; OUTPUT: ; Event file, FITS, with removed circular ; holes filled ; ; AUTHOR: ; J. A. Carter ; Department of Physics and Astronomy ; University of Leicester ; Leicester, LE1 7RH ; U.K. ; jac48@star.le.ac.uk ; ;- FUNCTION randomgen, xc, yc, rc xminc = DOUBLE(xc - rc) xmaxc = DOUBLE(xc + rc) yminc = DOUBLE(yc - rc) ymaxc = DOUBLE(yc + rc) wxc = xmaxc - xminc wyc = ymaxc - yminc x = DOUBLE(RANDOMU(a, /DOUBLE)) y = DOUBLE(RANDOMU(b, /DOUBLE)) xe = xminc + (wxc * x) ye = yminc + (wyc * y) RETURN, [xe, ye] END ; -------------------------------- ; -------------------------------- FUNCTION DRAWBOX, boxx, boxy, boxr, steps Stepx = (2*boxr)/(steps-1) Stepy = Stepx fvalx = boxx-boxr fvaly = boxy-boxr vxs = (FINDGEN(steps)*Stepx) + fvalx vys = (FINDGEN(steps)*Stepy) + fvaly coords = [[vxs], [vys]] RETURN, coords END ;------------------------------- ;------------------------------- FUNCTION FILLBOX, fillxc, fillyc, fillr, n, values, fbarray, valuefill, offPlot=offPlot, holePlot=holePlot, type=type, replace=replace IF KEYWORD_SET(type) THEN typeV = type ELSE typeV = 'A' boxXs = values[*, 0] boxYs = values[*, 1] boxVs = fbarray ;IF KEYWORD_SET(replace) THEN PRINT, 'Replace set ' IF KEYWORD_SET(offPlot) THEN BEGIN ; Diagnostic area ENDIF ELSE BEGIN FOR i=0, n-1 DO BEGIN testDistAll = ((boxXs[i]-fillxc)^2+(boxYs-fillyc)^2) ; Is this correct?? inindex = WHERE(testDistAll LE (fillr^2), ci) rowBvs = fbarray[i, *] IF ci NE 0 THEN BEGIN IF KEYWORD_SET(replace) THEN BEGIN rowBvs[inindex] = valuefill ENDIF ELSE BEGIN rowBvs[inindex] = rowBvs[inindex] + valuefill ENDELSE ENDIF boxVs[i, *] = rowBvs ENDFOR ENDELSE RETURN, boxVs END ;------------------------------- ;------------------------------- FUNCTION calhrarea, subx, suby, subr, subs, relX=relX, relY=relY, relR=relR, relType=relType, type=type, fy=fy, rl=rl, printInfo=printInfo, plothole=plothole, n=n IF KEYWORD_SET(n) THEN BEGIN n = n ENDIF ELSE BEGIN n=99999 ENDELSE ; Calculate hole and ring areas basicunit = 0.20 dim = 100 ;dim = FIX(basicunit * 2 * subs) unitsize = (dim/(2*subs))^2 boxarray = MAKE_ARRAY(dim, dim, value=0) coords = DRAWBOX(subx, suby, subs, dim) ringTag = 17 holeTag = 33 relsTag = 1 offPlot = 500 goneintorecalc = 'No' ; Fill boxarray thisFill = FILLBOX(subx, subY, subs, dim, coords, boxarray, ringTag) boxarray = thisFill thisFill = FILLBOX(subx, suby, subr, dim, coords, boxarray, holeTag, /holePlot) boxarray = thisFill ;thisFill = FILLBOX(subx, suby, subr, dim, coords, boxarray, offPlot, /offPlot) ;boxarray = thisFill pureHole = WHERE(thisFill EQ (holeTag + ringTag)) IF SIZE(pureHole, /N_DIMENSIONS) NE 0 THEN pureHole = N_ELEMENTS(pureHole) ELSE pureHole = 0 pureRing = WHERE(thisFill EQ ringTag) IF SIZE(pureRing, /N_DIMENSIONS) NE 0 THEN pureRing = N_ELEMENTS(pureRing) ELSE pureRing = 0 ; Count areas IF type EQ 'S' THEN BEGIN ; Simple hole pureHole = pureHole ENDIF extremeH = 7000L extreme = 1000L holeTagSpecial = 578 ringTagSpecial = 633 IF type EQ 'F' OR type EQ 'X' THEN BEGIN ; First fill - only adapt ring, not hole area pureHole = pureHole IF KEYWORD_SET(relX) THEN BEGIN FOR i = 0, N_ELEMENTS(relX)-1 DO BEGIN IF type EQ 'F' THEN thisFill = FILLBOX(relX[i], relY[i], relR[i], dim, coords, boxarray, relsTag, type='F') IF type EQ 'X' THEN thisFill = FILLBOX(relX[i], relY[i], relR[i], dim, coords, boxarray, relsTag, type='X', /holePlot) boxarray = thisFill ENDFOR ENDIF pureRing = WHERE(thisFill EQ ringTag) IF SIZE(pureRing, /N_DIMENSIONS) NE 0 THEN pureRing = N_ELEMENTS(pureRing) ELSE pureRing = 0 IF type EQ 'X' AND KEYWORD_SET(fy) AND KEYWORD_SET(rl) AND KEYWORD_SET(relType) THEN BEGIN ; PRINT, 'Calculating type X hole area incorporating unfilled crossovers' ; Need to re-fill do the calculation to include in the hole ; area where others holes have not been filled toTake = 999999 positions = 999999 FOR g = 0, N_ELEMENTS(rl)-1 DO BEGIN ti = WHERE(rl[g] NE fy, notdone) value = N_ELEMENTS(fy) IF notdone EQ value AND relType[g] EQ 2 THEN BEGIN toTake = [toTake, g] positions= [positions, rl[g]] ENDIF ENDFOR IF N_ELEMENTS(toTake) GE 2 THEN BEGIN ;PRINT, 'Recalculating hole size' goneintorecalc = 'Yes' toTake = toTake[1:*] mask = REPLICATE(0, N_ELEMENTS(rl)) mask[toTake] = 1 others = WHERE(mask NE 0, count) avoids = WHERE(mask NE 1, counta) IF count GT 0 THEN BEGIN newX = relX[others] newY = relY[others] newR = relR[others] newF = relType[others] ; Take x,y, r and fill box with extreme value for not-yet ; overlapped hole FOR i=0, N_ELEMENTS(newX)-1 DO BEGIN IF type EQ 'X' THEN thisFill = FILLBOX(newX[i], newY[i], newR[i], dim, coords, thisFill, extreme, type='X', /replace) ENDFOR thisFill = FILLBOX(subx, suby, subr, dim, coords, thisFill, holeTagSpecial) ; Need to do other overlaps again IF counta GT 0 THEN BEGIN avoidX = relX[avoids] avoidY = relY[avoids] avoidR = relR[avoids] FOR i=0, counta-1 DO BEGIN IF type EQ 'X' THEN thisFill = FILLBOX(avoidX[i], avoidY[i], avoidR[i], dim, coords, thisFill, relsTag, type='X') ENDFOR ENDIF total1 = holeTagSpecial + holeTag + ringTag total2 = extreme + holeTagSpecial pureHole = WHERE((thisFill EQ total1) OR (thisFill EQ total2)) ENDIF ELSE BEGIN PRINT, 'No overlap with case not filled X tested' pureHole = WHERE(thisFill EQ (holeTag + ringTag)) ENDELSE ENDIF ELSE BEGIN ;PRINT, 'No overlap with case not filled X at all' pureHole = WHERE(thisFill EQ (holeTag + ringTag)) ENDELSE IF SIZE(pureHole, /N_DIMENSIONS) NE 0 THEN pureHole = N_ELEMENTS(pureHole) ELSE pureHole = 0 ENDIF ELSE BEGIN IF type EQ 'X' AND goneintorecalc EQ 'No' THEN pureHole = WHERE(thisFill EQ (holeTag + ringTag)) ENDELSE ENDIF IF KEYWORD_SET(plothole) AND (n EQ 6937) AND KEYWORD_SET(rl) THEN BEGIN window,xsize=640,ysize=640 minX = subX - subS maxX = subX + subS minY = subY - subS maxY = subY + subS uS = dim/(2 * subs) axisX = (FINDGEN(dim+1)*((2*subs)/FLOAT(dim))) + minX PRINT, minX, maxX, minY, maxY SET_PLOT, 'PS' thisDevice = !D.Name ; Color table loadct, 13 ; Rainbow Device, /Color, FILENAME='Lookatspecifichole.ps' PLOT, axisX, thisFill[0, *], XRANGE=[minX, maxX], YRANGE=[minY, maxY], /NODATA, PSYM=2, xstyle=1, ystyle=1 TVCIRCLE, subs, subx, suby, COLOR=255, /data, thick=2 TVCIRCLE, subr, subx, suby, /data, thick=2 FOR i = 0, N_ELEMENTS(newX)-1 DO BEGIN dcX = subX - newX[i] dcY = subY - newY[i] coord1 = dim/2 - (dcX * uS) coord2 = dim/2 - (dcY * uS) PRINT, 'Coordinate 1, 2: ', newX[i], newY[i], coord1, coord2 ;TVCIRCLE, (subr * uS), coord1, coord2, COLOR=255, /DATA TVCIRCLE, subr, newX[i], newY[i], /data, color=255 ; STOP ENDFOR XYOUTS, subx, suby, 'Centre' FOR i=0, dim-1 DO BEGIN index1 = WHERE(thisFill[*, i] EQ total1) index2 = WHERE(thisFill[*, i] EQ total2) ;index3 = WHERE(thisFill[*, i] EQ ringTag + holeTag) ;index4 = WHERE(thisFill[*, i] EQ ringTag + holeTag + holeTagSpecial + relsTag) ; index4 = WHERE((thisFill[*, i] EQ (holeTag + ringTag)) OR (thisFill GE (extreme + holeTag + ringTag))) IF SIZE(index1, /N_DIMENSIONS) GT 0 THEN BEGIN index1 = axisX[index1] OPLOT, index1, REPLICATE((i/uS)+minY, N_ELEMENTS(index1)), psym=2 ENDIF IF SIZE(index2, /N_DIMENSIONS) GT 0 THEN BEGIN index2 = axisX[index2] OPLOT, index2, REPLICATE((i/uS)+minY, N_ELEMENTS(index2)), psym=2, color=255 ENDIF IF SIZE(index3, /N_DIMENSIONS) GT 0 THEN BEGIN index3 = axisX[index3] OPLOT, index3, REPLICATE((i/uS)+minY, N_ELEMENTS(index3)), psym=4 ENDIF IF SIZE(index4, /N_DIMENSIONS) GT 0 THEN BEGIN index4 = axisX[index4] OPLOT, index4, REPLICATE((i/uS)+minY, N_ELEMENTS(index4)), psym=4, color=255 ENDIF ENDFOR DEVICE, /CLOSE_FILE SET_PLOT, 'X' ENDIF calhrareaval = PTR_NEW([pureRing/unitSize, pureHole/unitSize]) RETURN, calhrareaval END ;------------------------------- ;------------------------------- FUNCTION relations, xc, yc, r, s, regionInfo, num ; Return new relationships newCrossX = 99999 newCrossY = 99999 newCrossR = 99999 newFlags = 99999 newRels = 99999 caseSwitch=0 caseType = 0L ; Comparison hole FOR i=0, N_ELEMENTS(regionInfo)-1 DO BEGIN ; caseSwitch = 0 ; Somewhat controversial - check!! caseType = 0L cX = regionInfo[i].DETX cY = regionInfo[i].DETY cR = regionInfo[i].R cX = cX[0] cY = cY[0] cR = cR[0] ; Calculate distance between holes dh = SQRT((cX-xc)^2+(cY-yc)^2) IF i EQ num THEN BEGIN ; PRINT, 'Same hole' ENDIF ELSE BEGIN ; PRINT, 'Comparing hole ', num ,' with ', i ; CASE 1: simple hole IF ((dh GT (r+cR)) AND (dh GE (s + cR))) THEN BEGIN ENDIF ; CASE 2: hole-hole IF (dh LT (r+cR)) THEN BEGIN newCrossX = [newCrossX, cX] newCrossY = [newCrossY, cY] newCrossR = [newCrossR, cR] caseType = 2L caseSwitch = 1 newFlags = [newFlags, caseType] newRels = [newRels, i] ENDIF ; CASE 3: ring-hole overlap IF ((dh LT (s + cR)) AND caseType NE 2L) THEN BEGIN newCrossX = [newCrossX, cX] newCrossY = [newCrossY, cY] newCrossR = [newCrossR, cR] caseType = 3L caseSwitch = 1 newFlags = [newFlags, caseType] newRels = [newRels, i] ENDIF ENDELSE ENDFOR IF caseSwitch EQ 1 THEN BEGIN newCrossX = newCrossX[1:*] newCrossY = newCrossY[1:*] newCrossR = newCrossR[1:*] newFlags = newFlags[1:*] newRels = newRels[1:*] ENDIF ELSE BEGIN ENDELSE relations = {RELATES, holeNo: num, newCrossX:PTR_NEW(newCrossX), newCrossY:PTR_NEW(newCrossY), newCrossR:PTR_NEW(newCrossR), newFlags:PTR_NEW(newFlags), newRelates:PTR_NEW(newRels), caseS: caseSwitch} RETURN, relations END ;================== Filling function ============================= FUNCTION filling, subx, suby, subr, subs, theseevents, xts=xts, yts=yts, rts=rts, type=type, n=n newevent = 0 ringindex = WHERE(((theseevents.DETX - subx)^2+(theseevents.DETY - suby)^2 LE (subs^2)) AND (theseevents.flag EQ 0)) ringindex = WHERE(((theseevents.DETX-subx)^2)+((theseevents.DETY-(suby))^2) LT (subs^2)) ; PRINT, 'Filling: ', N_ELEMENTS(ringIndex) IF KEYWORD_SET(n) THEN n=n ELSE n=999999 IF (SIZE(ringindex, /DIMENSIONS) NE 0) THEN BEGIN ring = theseevents[ringindex] ; Back to taking all the events newevent = {TIME:DOUBLE(0.0), RAWX:0L, RAWY:0L, DETX:0.0, DETY:0.0, X:0L, Y:0L, PHA:0L, PI:0L, FLAG:0L, PATTERN:0L, CCDNR:0L} newevent = REPLICATE(newevent, N_ELEMENTS(ring)) neX = 0 neY = 0 loopc = 0 ; Find new coordinates WHILE (N_ELEMENTS(neX)-1 LT N_ELEMENTS(ring)) DO BEGIN loopc = loopc + 1 c = 0 newcoords = RANDOMGEN(subx, suby, subr) xe = newcoords[0] ye = newcoords[1] IF (((xe-(subx))^2+(ye-(suby))^2) LT (subr^2)) THEN BEGIN IF type EQ 'S' OR type EQ 'F' THEN BEGIN neX = [neX, xe] neY = [neY, ye] ENDIF ELSE BEGIN IF type EQ 'X' THEN BEGIN FOR f=0, N_ELEMENTS(xts)-1 DO BEGIN IF (((xe - xts[f])^2+(ye - yts[f])^2 GT (rts[f]^2))) THEN BEGIN ;If the new x, y values fall outside all the overlapping holes c = c + 1 ENDIF ENDFOR IF (c EQ N_ELEMENTS(xts)) THEN BEGIN neX = [neX, xe] neY = [neY, ye] ENDIF ENDIF ENDELSE ENDIF ENDWHILE IF (N_ELEMENTS(neX) GT 1) THEN BEGIN neX = neX[1:*] neY = neY[1:*] nebadindex = WHERE(((neX-subx)^2+(neY-suby)^2 GE (subr^2))) ; Assign these values to a real event newevent.TIME = ring.TIME newevent.RAWX = ring.RAWX newevent.RAWY = ring.RAWY newevent.DETX = neX newevent.DETY = neY newevent.X = ring.X newevent.Y = ring.Y newevent.PHA = ring.PHA newevent.PI = ring.PI newevent.FLAG = ring.FLAG newevent.PATTERN = ring.PATTERN newevent.CCDNR = ring.CCDNR ENDIF RETURN, newevent ENDIF ELSE BEGIN RETURN, 0 ENDELSE END ;============================================================= FUNCTION cleanregsoutbounds, origevents, inregs, instr, cleanB ; Apply boundaries to holes depending on instrument xmin = FLOAT(MIN(origevents.DETX)) xmax = FLOAT(MAX(origevents.DETX)) ymin = FLOAT(MIN(origevents.DETY)) ymax = FLOAT(MAX(origevents.DETY)) IF cleanB EQ 'YES' THEN BEGIN PRINT, 'Cleaning boundaries' IF (instr EQ 'PN') THEN BEGIN PRINT, 'PN boundaries set' exX1 = -18360L ; extreme X, 1 exX2 = 14000L; extreme X, 2 exY1 = -16800L; extreme Y, 1 exY2 = 14600L; extreme Y, 2 midX1 = exX1 ; mid X, 1 midX2 = exX2 ; mid X, 2 midY1 = exY1 ; mid Y, 1 midY2 = exY2 ; mid Y, 2 ENDIF ELSE BEGIN exX1 = -19720L ; extreme X, 1 exX2 = 19880L ; extreme X, 2 exY1 = -19960L ; extreme Y, 1 exY2 = 19800L ; extreme Y, 2 midX1 = -13110. ; mid X, 1 midX2 = 13161. ; mid X, 2 midY1 = -6959. ; mid Y, 1 midY2 = 6640. ; mid Y, 2 ENDELSE exX1 = -20399L ; extreme X, 1 exX2 = xmax ; extreme X, 2 exY1 = ymin ; extreme Y, 1 exY2 = ymax ; extreme Y, 2 xcd = ABS((exX2 - exX1)/2) xcLine = xcd + exX1 ycd = ABS((exY2 - (exY1))/2) ycLine = ycd + exY1 index = WHERE((inregs.DETX GE exX1) AND $ (inregs.DETX LE exX2) AND $ (inregs.DETY GE exY1) AND $ (inregs.DETY LE exY2)) index = WHERE(((ABS(ABS(inregs.DETX - xcLine) - inregs.R)) LT xcd) AND $ ((ABS(ABS(inregs.DETY - ycLine) - inregs.R)) LT ycd)) inregs = inregs[index] midxcd = (midX2 - midX1)/2 midxcLine = midxcd + midX1 midycd = (midY2 - midY1)/2 midycLine = midycd + midY1 crossindex1 = WHERE(((ABS(ABS(inregs.DETX - midxcLine) - inregs.R)) LT midxcd)) crossindex2 = WHERE(((ABS(ABS(inregs.DETY - midycLine) - inregs.R)) LT midycd)) crossindex = [crossindex1, crossindex2] crossindex = crossindex[SORT(crossindex)] crossindex = crossindex[UNIQ(crossindex)] outevents = inregs[crossindex] ENDIF ELSE BEGIN outevents = inregs ENDELSE RETURN, outevents END ;================== Program body ============================= PRO simple_ghostit_ind PRINT, 'Program called is simple_ghostit_ind' time = SYSTIME(1) start_mem = MEMORY(/CURRENT) instr = (GETENV('IDL_ARG1')) event = (GETENV('IDL_ARG2')) region = (GETENV('IDL_ARG3')) resultfile = (GETENV('IDL_ARG4')) plotit = (GETENV('IDL_ARG5')) FITS_INFO, event, /SILENT, N_ext = next eventfile = event regionsfile = region event = MRDFITS(event, 1, headerevents, /SILENT) head = HEADFITS(eventfile, /SILENT) regions = MRDFITS(region, 1, /SILENT) originalevent = event takeevent = event finalresults = originalevent[0] ; New region structure and fill new structure newregionsbase = {SHAPE: STRING('C'), DETX: 0.0, DETY: 0.0, R: 0.0, ROTANG: 0.0, COMPONENT: 0L} newregions = REPLICATE(newregionsbase, N_ELEMENTS(regions)) xmin = FLOAT(MIN(originalevent.DETX)) xmax = FLOAT(MAX(originalevent.DETX)) ymin = FLOAT(MIN(originalevent.DETY)) ymax = FLOAT(MAX(originalevent.DETY)) FOR r=0, N_ELEMENTS(regions)-1 DO BEGIN rDX = regions[r].DETX rDX = rDX[0] rDY = regions[r].DETY rDY = rDY[0] rr = regions[r].R rr = rr[0] ro = regions[r].ROTANG ro = ro[0] ; Fill new structure newregions[r].SHAPE = regions[r].SHAPE newregions[r].DETX = rDX newregions[r].DETY = rDY newregions[r].R = rr newregions[r].ROTANG = ro newregions[r].COMPONENT = regions[r].COMPONENT ENDFOR ;instm1 = cleanregsoutbounds(originalevent, newregions, instr, 'NO') instm1 = cleanregsoutbounds(originalevent, newregions, instr, 'YES') ;---------------- Main processing ---------------------------- ;sortHoles = SORT(instm1.DETX) ;PRINT, 'N_ELEMENTS(instm1) pre sort: ', N_ELEMENTS(instm1) ;instm1 = instm1[sortHoles] ;PRINT, 'N_ELEMENTS(instm1) post sort: ', N_ELEMENTS(instm1) ;STOP fillList = 7898757 ringList = 99999 xList = 99999 yList = 99999 incr = 4.0 oincr = incr big = 1.08 small = 0.70 vsmall = 0.55 limit = 2000 addingMyNewCount = 0 zeroCount = 0 alarmList = 0 mydeficits = 0 beenfilled = 'No' defArray = 0 bigCount = 0 PRINT, 'Total number of holes to ghost: ', N_ELEMENTS(instm1) start = 0 ;fin = 50 fin = N_ELEMENTS(instm1)-1 IF start NE 0 THEN fillList = FIX(FINDGEN((start))) ;print, 'Fill list: ', fillList FOR n=start, fin DO BEGIN xc = instm1[n].DETX yc = instm1[n].DETY r = instm1[n].R s = SQRT(DOUBLE(2.)) * DOUBLE(r) xList = [xList, xc] yList = [yList, yc] ; Find relationships to this hole and act accordingly firstRels = relations(xc, yc, r, s, regions, n) holeClass = firstRels.caseS PRINT, 'Hole : ', n, ' of:', N_ELEMENTS(instm1)-1,' xc:', xc, ' yc: ', yc, ' class:', holeClass ;PRINT, 'Radius of hole, ring, relations: ', r, s, N_ELEMENTS(*firstRels.newcrossX) ;PRINT, 'Xs: ', *firstRels.newcrossX ;PRINT, 'Ys: ', *firstRels.newcrossY ;PRINT, 'Rels: ', *firstRels.newRelates ;PRINT, 'Flags: ', *firstRels.newFlags IF holeClass NE 0 THEN BEGIN nowcrossXt = *firstRels.newcrossX nowcrossYt = *firstRels.newcrossY nowcrossRt = *firstRels.newcrossR nowFlagsFt = *firstRels.newFlags nowRelates = *firstRels.newRelates initRels = firstRels ninitRels = N_ELEMENTS(*firstRels.newcrossX) bitsFilled = 'No' FOR g = 0, N_ELEMENTS(nowRelates)-1 DO BEGIN IF SIZE(WHERE(nowRelates[g] EQ fillList), /DIMENSIONS) GT 0 AND nowFlagsFT[g] EQ 2 THEN BEGIN bitsFilled = 'Yes' ENDIF ENDFOR IF bitsFilled EQ 'No' THEN BEGIN ; Need to increase ring radius as area reduced holeArea = *(calhrarea(xc, yc, r, s, relX=nowcrossXt, relY=nowcrossYt, relR=nowcrossRt, type='F')) currHole = holeArea[1] currRing = holeArea[0] initHole = currHole currRingRadius = s incr = oincr WHILE currRing LT currHole AND (ABS(currRing-currHole) GT limit) DO BEGIN currRingRadius = currRingRadius + incr newRels = relations(xc, yc, r, currRingRadius, regions, n) initRels = newRels ringcrossX = *newRels.newcrossX ringcrossY = *newRels.newcrossY ringcrossR = *newRels.newcrossR holeArea = *(calhrarea(xc, yc, r, currRingRadius, relX=ringcrossX, relY=ringcrossY, relR=ringcrossR, type='F')) currRing = holeArea[0] currHole = holeArea[1] currHole = initHole ;PRINT, 'F:INCR:Rrad, Ring, Hole, def: ', incr, currRingRadius, currRing, currHole, currRing-currHole upp = 120000 mid = 30000 IF ABS(currHole-currRing) GT upp THEN incr = incr * big IF ABS(currHole-currRing) GT mid AND ABS(currHole-currRing) LT upp THEN incr = incr IF ABS(currHole-currRing) LT mid THEN incr = incr * vsmall IF incr LT 0.5 OR ABS(currHole-currRing) LT 20000 THEN incr = 0.5 IF ABS(currHole-currRing) LT 6500 THEN incr = 0.6 IF incr GT 30 THEN incr = 30.0 ENDWHILE incr = oincr * 0.50 incr = 1.0 WHILE currRing GT currHole AND currRingRadius GT r AND ABS(currRing-currHole) GT limit DO BEGIN preChange = currRingRadius currRingRadius = currRingRadius - incr newRels = relations(xc, yc, r, currRingRadius, regions, n) initRels = newRels ringcrossX = *newRels.newcrossX ringcrossY = *newRels.newcrossY ringcrossR = *newRels.newcrossR ringRelate = *newRels.newRelates ringFlags = *newRels.newFlags holeArea = *(calhrarea(xc, yc, r, currRingRadius, relX=ringcrossX, relY=ringcrossY, relR=ringcrossR, type='F')) currRing = holeArea[0] currHole = holeArea[1] currHole = initHole ;PRINT, 'F:DECR:Rrad, Ring, Hole, def: ', incr, currRingRadius, currRing, currHole, currRing-currHole upp = 120000 mid = 30000 IF ABS(currHole-currRing) GT upp THEN incr = incr * big IF ABS(currHole-currRing) GT mid AND ABS(currHole-currRing) LT upp THEN incr = incr IF ABS(currHole-currRing) LT mid THEN incr = incr * vsmall IF incr LT 1.0 OR ABS(currHole-currRing) LT 20000 THEN incr = 1.0 IF incr GT 30.0 THEN incr = 30.0 ENDWHILE defArray = [defArray, (currRing-currHole)] calRingRadius = SQRT(currRing/!PI) ;PRINT, 'Type F new ring: ', currRing, currHole, currRing-currHole myNew = filling(xc, yc, r, currRingRadius, takeevent, type='F') ringList = [ringList, currRingRadius] ENDIF ELSE BEGIN holeArea = *(calhrarea(xc, yc, r, s, relX=nowcrossXt, relY=nowcrossYt, relR=nowcrossRt, relType=nowFlagsFt, type='X', fy=fillList, rl=nowRelates, n=n, /printInfo, /plothole)) currHole = holeArea[1] currRing = holeArea[0] currRingRadius = s initHole = currHole ; Can use this value as hole area shouldn't change incrCount = 0 beenIncr = 'No' incr = oincr lastT = 'X' WHILE currRing LT currHole AND (ABS(currRing-currHole) GT limit) DO BEGIN beenIncr = 'Yes' currRingRadius = currRingRadius + incr newRels = relations(xc, yc, r, currRingRadius, regions, n) initRels = newRels ringcrossX = *newRels.newcrossX ringcrossY = *newRels.newcrossY ringcrossR = *newRels.newcrossR ringRelate = *newRels.newRelates ringFlags = *newRels.newFlags holeArea = *(calhrarea(xc, yc, r, currRingRadius, relX=ringcrossX, relY=ringcrossY, relR=ringcrossR, type='X', n=n, /plothole)) currRing = holeArea[0] currHole = holeArea[1] currHole = initHole incrCount = incrCount + 1 ;PRINT, 'X:INCR:Rrad, Ring, Hole, def: ', incr, currRingRadius, currRing, currHole, currRing-currHole upp = 120000 mid = 30000 IF ABS(currHole-currRing) GT upp THEN incr = incr * big IF ABS(currHole-currRing) GT mid AND ABS(currHole-currRing) LT upp THEN incr = incr IF ABS(currHole-currRing) LT mid THEN incr = incr * vsmall IF incr LT 1 OR ABS(currHole-currRing) LT 20000 THEN incr = 1.0 IF incr GT 30 THEN incr = 30.0 ENDWHILE preChange = currRingRadius IF beenIncr EQ 'Yes' THEN incr = incr * 0.3 ELSE incr = oincr * 0.50 incr = oincr * 0.50 incr = 1.0 incrCount = 0 lastDef = 0 WHILE currRing GT currHole AND currRingRadius GT r AND ABS(currRing-currHole) GT limit AND currHole NE 0 DO BEGIN preChange = currRingRadius currRingRadius = currRingRadius - incr newRels = relations(xc, yc, r, currRingRadius, regions, n) initRels = newRels ringcrossX = *newRels.newcrossX ringcrossY = *newRels.newcrossY ringcrossR = *newRels.newcrossR ringRelate = *newRels.newRelates ringFlags = *newRels.newFlags thisT = 'X' holeArea = *(calhrarea(xc, yc, r, currRingRadius, relX=ringcrossX, relY=ringcrossY, relR=ringcrossR, type=thisT, n=n, /plothole)) currRing = holeArea[0] currHole = holeArea[1] currHole = initHole ; PRINT, 'X:DECR:Rrad, Ring, Hole, df: ', incr, currRingRadius, currRing, currHole, currRing-currHole upp = 150000 mid = 50000 low = 30000 IF ABS(currRing-currHole) GT upp THEN incr = incr * big IF ABS(currRing-currHole) GT mid AND ABS(currRing-currHole) LT upp THEN incr = incr IF ABS(currRing-currHole) GT low AND ABS(currRing-currHole) LT mid THEN incr = incr*small IF ABS(currHole-currRing) LT low THEN incr = incr * vsmall IF incr LT 1 OR ABS(currHole-currRing) LT 20000 THEN incr = 1.0 IF incr GT 30 THEN incr = 30.0 incrCount = incrCount + 1 lastDef = (currRing-currHole) ENDWHILE ;PRINT, 'Decreasing ring area count: ', incrCount finalRingRadius = currRingRadius ; Sort out holes to avoid filling nowcrossXt = *firstRels.newcrossX nowcrossYt = *firstRels.newcrossY nowcrossRt = *firstRels.newcrossR nowFlagsFt = *firstRels.newFlags nowRelates = *firstRels.newRelates ;nowRelates = ringRelate ;nowFlags = ringFlags relsFilled = 99999 FOR g = 0, N_ELEMENTS(nowRelates)-1 DO BEGIN thisFill = WHERE(nowRelates[g] EQ fillList, countFilled) IF countFilled GT 0 AND nowFlagsFt[g] EQ 2 THEN BEGIN relsFilled = [relsFilled, g] ENDIF ENDFOR IF N_ELEMENTS(relsFilled) GE 2 THEN BEGIN relsFilled = relsFilled[1:*] badcrossX = nowcrossXt[relsFilled] badcrossY = nowcrossYt[relsFilled] badcrossR = nowcrossRt[relsFilled] badcrossF = nowFlagsFt[relsFilled] badcrossRels = nowRelates[relsFilled] ENDIF ELSE BEGIN badcrossX = nowcrossXt badcrossY = nowcrossYt badcrossR = nowcrossRt badcrossF = nowFlagsFt badcrossRels = nowRelates ENDELSE ; Fill holes mem = MEMORY(/CURRENT) IF currHole NE 0 THEN BEGIN defArray = [defArray, (currRing-currHole)] myNew = filling(xc, yc, r, finalRingRadius, takeevent, type='X', xts=badcrossX, yts=badcrossY, rts=badcrossR, n=n) ENDIF ELSE BEGIN defArray = [defArray, 0.0] zeroCount = zeroCount + 1 ENDELSE ringList = [ringList, finalRingRadius] ENDELSE ENDIF ELSE BEGIN ; Simple hole defArray = [defArray, 0.0] myNew = filling(xc, yc, r, s, takeevent, type='S') ringList = [ringList, s] ENDELSE ;PRINT, '**************This hole completed*************' fillList = [fillList, n] IF bigCount EQ 0 AND start EQ 0 THEN fillList = fillList[1:*] IF SIZE(myNew, /DIMENSIONS) NE 0 AND N_ELEMENTS(myNew) GT 1 THEN BEGIN addingMyNewCount = addingMyNewCount + 1 resultstruct = finalresults[0] resultstruct = REPLICATE(resultstruct, N_ELEMENTS(myNew)) resultstruct.TIME = myNew.TIME resultstruct.RAWX = myNew.RAWX resultstruct.RAWY = myNew.RAWY resultstruct.DETX = myNew.DETX resultstruct.DETY = myNew.DETY resultstruct.X = myNew.X resultstruct.Y = myNew.Y resultstruct.PHA = myNew.PHA resultstruct.PI = myNew.PI resultstruct.FLAG = myNew.FLAG resultstruct.PATTERN = myNew.PATTERN resultstruct.CCDNR = myNew.CCDNR finalresults = [finalresults, resultstruct] beenfilled = 'Yes' ENDIF ; This is very important to reset myNew to zero myNew = 0 bigCount = bigCount + 1 PTR_FREE, PTR_VALID() ENDFOR IF beenfilled EQ 'Yes' THEN BEGIN PRINT, 'Holes filled, create final structure' finalresults = finalresults[1:*] finalresults = [originalevent, finalresults] finalresults = finalresults[SORT(finalresults.TIME)] ENDIF ; Cleaning up IF (instr EQ 'PN') THEN BEGIN exX1 = -18360L ; extreme X, 1 exX2 = 14000L; extreme X, 2 exY1 = -16800L; extreme Y, 1 exY2 = 14600L; extreme Y, 2 midX1 = exX1 ; mid X, 1 midX2 = exX2 ; mid X, 2 midY1 = exY1 ; mid Y, 1 midY2 = exY2 ; mid Y, 2 ENDIF ELSE BEGIN exX1 = -19720L ; extreme X, 1 exX2 = 19880L ; extreme X, 2 exY1 = -19960L ; extreme Y, 1 exY2 = 19800L ; extreme Y, 2 midX1 = -13110. ; mid X, 1 midX2 = 13161. ; mid X, 2 midY1 = -6959. ; mid Y, 1 midY2 = 6640. ; mid Y, 2 ENDELSE exX1 = -20399L ; extreme X, 1 exX2 = xmax ; extreme X, 2 exY1 = ymin ; extreme Y, 1 exY2 = ymax ; extreme Y, 2 inside = WHERE((finalresults.DETX GE exX1) AND (finalresults.DETX LE exX2) AND $ (finalresults.DETY GE exY1) AND (finalresults.DETY LE exY2)) finalresults = finalresults[inside] midxcd = (midX2 - midX1)/2 midxcLine = midxcd + midX1 midycd = (midY2 - midY1)/2 midycLine = midycd + midY1 crossindex1 = WHERE(((ABS(ABS(finalresults.DETX - midxcLine))) LT midxcd)) crossindex2 = WHERE(((ABS(ABS(finalresults.DETY - midycLine))) LT midycd)) crossindex = [crossindex1, crossindex2] crossindex = crossindex[SORT(crossindex)] crossindex = crossindex[UNIQ(crossindex)] finalresults = finalresults[crossindex] ; Write to file MWRFITS, 0, resultfile, head, /CREATE, /No_comment SXADDHIST, 'IDL events added to source holes', headerevents MWRFITS, finalresults, resultfile, headerevents, /SILENT, /No_comment FOR j = 2, next DO BEGIN ; PRINT, "Extension to add: ", j ext = MRDFITS(eventfile, j, thisHeader, /SILENT) MWRFITS, ext, resultfile, thisHeader, /SILENT ENDFOR PTR_FREE, PTR_VALID() time2 = SYSTIME(1) PRINT, 'Time started: ', time, ' Differences in time: ', time2 - time, (time2 - time)/60.0 , (time2 - time)/(60.0*60.0) END