Skip to content

Commit

Permalink
Major modifications
Browse files Browse the repository at this point in the history
1. Major error in jude_add_frames which underestimated the time by about 10%.
2. Factor of 10 increase in jude_add_frames.
3. Rewrote jude_coadd.pro to treat partial pixels.
4. Add jude_light_curve
  • Loading branch information
jaymurthy committed Dec 9, 2017
1 parent 7c3d09c commit c791a0e
Show file tree
Hide file tree
Showing 3 changed files with 255 additions and 45 deletions.
37 changes: 25 additions & 12 deletions jude_add_frames.pro
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,8 @@
;JM: Aug. 21, 2017: Made Dtime floating instead of double.
;JM: Aug. 27, 2017: Reset ref_frame if required
;JM: Nov. 24, 2017: If par not defined, now reads from jude_params()
;JM: Dec. 08, 2017: Was not taking the reference frame when calculating times.
;JM: Dec. 09, 2017: Major speed increase.
;Copyright 2016 Jayant Murthy
;
; Licensed under the Apache License, Version 2.0 (the "License");
Expand Down Expand Up @@ -113,10 +115,11 @@ if (min_frame eq 0)then min_frame = $
min(where((data.dqi eq 0) and (abs(xoff) lt 1000) and $
(abs(yoff) lt 1000)))

if (n_elements(ref_frame) eq 0)then ref_frame = min_frame
if (ref_frame eq -1)then ref_frame = min_frame
while (((abs(xoff[ref_frame]) gt 1000) or (abs(yoff[ref_frame]) gt 1000)) and $
(ref_frame lt (n_elements(data) - 2))) do ref_frame = ref_frame + 1
if (n_elements(ref_frame) eq 0)then begin
ref_frame = min_frame
while (((abs(xoff[ref_frame]) gt 1000) or (abs(yoff[ref_frame]) gt 1000)) and $
(ref_frame lt (n_elements(data) - 2))) do ref_frame = ref_frame + 1
endif
if (max_frame eq 0) then max_frame = n_elements(data)-1

;If the offsets are not defined, they are set to 0
Expand Down Expand Up @@ -182,6 +185,10 @@ dst = (gx - (256*resolution))^2 + $
q = where(dst lt g_rad^2, nq)
times[q] = 1

old_xindex = 0
old_yindex = 0
shft_times = times
ishft = 0
tstart = systime(1)
for ielem = min_frame,max_frame do begin
if (not(keyword_set(notime)) and ((ielem mod 100) eq 0)) then begin
Expand Down Expand Up @@ -223,14 +230,20 @@ endif
;pixels (modified by the resolution) of the centre. This may have to be refined.
;Note that I don't do this time consuming step if notime is set.
if (not(keyword_set(notime))) then begin
tmp =shift(times, round(xoff[ielem]), round(yoff[ielem]))
xindex = round(xoff[ielem])
yindex = round(yoff[ielem])
if (xindex gt 0)then tmp[0:xindex-1,*] = 0
if (xindex lt 0)then tmp[gxsize + xindex:gxsize - 1, *] = 0
if (yindex gt 0)then tmp[*, 0:yindex-1] = 0
if (yindex lt 0)then tmp[*, gysize + yindex:gxsize - 1] = 0
pixel_time = pixel_time + tmp
xindex = round(xoff[ielem]-xoff_start)
yindex = round(yoff[ielem]-yoff_start)

if ((xindex ne old_xindex) or (yindex ne old_yindex))then begin
pixel_time = pixel_time + ishft*shft_times
shft_times =shift(times, xindex, yindex)
if (xindex gt 0)then shft_times[0:xindex-1,*] = 0
if (xindex lt 0)then shft_times[gxsize + xindex:gxsize - 1, *] = 0
if (yindex gt 0)then shft_times[*, 0:yindex-1] = 0
if (yindex lt 0)then shft_times[*, gysize + yindex:gxsize - 1] = 0
old_xindex = xindex
old_yindex = yindex
ishft = 0
endif else ishft = ishft + 1
endif
endif
endfor
Expand Down
45 changes: 12 additions & 33 deletions jude_coadd.pro
Original file line number Diff line number Diff line change
Expand Up @@ -164,44 +164,23 @@ pro jude_coadd, input_dir, output_file, ra_cent, dec_cent, fov,$
;Use total counts rather than counts/s
time_val = dtime[index]
data_val = data[index]*time_val

;Some fraction of the counts should be in the central pixel
multx = dmultx[index]
multy = dmulty[index]
grid[ ix, iy, ifile] = grid[ ix, iy, ifile] + $
data_val*multx*multy
times[ix, iy, ifile] = times[ix, iy, ifile] + $
time_val*multx*multy
gadd[ ix, iy, ifile] = gadd[ ix, iy, ifile] + $
multx*multy


;Now the off-center pixels
tx = fix(xf[index] gt 0.5) - fix(xf[index] le 0.5)
ty = fix(yf[index] gt 0.5) - fix(yf[index] le 0.5)
for itx = -1, 1 do for ity = -1, 1 do begin
case 1 of
itx eq tx: multx = 1 - dmultx[index]
itx eq 0: multx = dmultx[index]
else: multx = 0d
endcase
case 1 of
ity eq ty: multy = 1 - dmulty[index]
ity eq 0: multy = dmulty[index]
else: multy = 0d
endcase
if ((itx eq 0) and (ity eq 0))then begin
multx = 0d
multy = 0d
endif
grid[ ix + itx, iy + ity, ifile] = grid[ix + itx, iy + ity, ifile] + $
data_val*multx*multy
times[ix + itx, iy + ity, ifile] = times[ix + itx, iy + ity, ifile] + $
time_val*multx*multy
gadd[ ix + itx, iy + ity, ifile] = gadd[ix + itx, iy + ity, ifile] + $
multx*multy

ity = [-1, -1, -1, 0, 0, 0, 1, 1, 1]
itx = [-1, 0, 1, -1, 0, 1, -1, 0, 1]
multx = (itx eq tx)*(1 - dmultx[index]) + (itx eq 0)*dmultx[index]
multy = (ity eq ty)*(1 - dmulty[index]) + (ity eq 0)*dmulty[index]
grid[ ix-1:ix+1, iy-1:iy+1, ifile] = grid[ix-1:ix+1, iy-1:iy+1, ifile] + $
data_val*multx*multy
times[ix-1:ix+1, iy-1:iy+1, ifile] = times[ix-1:ix+1, iy-1:iy+1, ifile] + $
time_val*multx*multy
gadd[ ix-1:ix+1, iy-1:iy+1, ifile] = gadd[ix-1:ix+1, iy-1:iy+1, ifile] + $
multx*multy
;str = string(nqy) + string(xf[index]) + string(yf[index]) + string(itx) + string(ity) + string(multx) + string(multy)
;print,strcompress(str)
endfor
endfor; iqy
endfor ;iy
endif ;nqx
Expand Down
218 changes: 218 additions & 0 deletions jude_light_curve.pro
Original file line number Diff line number Diff line change
@@ -0,0 +1,218 @@
;+
; NAME: JUDE_CENTROID
; PURPOSE: Improves registration by centroiding on a single star
; CALLING SEQUENCE:
; jude_centroid, events_file, grid2, params, xstar, ystar, $
; xoff = xoff, yoff = yoff ,$
; boxsize = boxsize,$
; init_size = init_size, medsiz = medsiz, $
; test = test, cent_file = cent_file, display = display, $
; nosave = nosave, defaults = defaults, new_star = new_star
; INPUTS:
; Events_file: Name of the input photon list (Level 2) file.
;
; OPTIONAL INPUTS: (If not defined beforehand, defined in the program)
; Params: parameter list
; Xstar: Known position of star. If blank, then position used
; ycent: Ystar: in program is returned.
; OUTPUT:
; Grid2: Array containing final image.
; KEYWORDS:
; Xoff: If the spacecraft offsets are known, they are applied to
; Yoff: the data; if not, they are calculated and passed back
; Boxsize: The search box size for the star. If there is significant
; spacecraft motion, boxsize may have to be larger.
; Init_size: The initisal search box for a centroid in case the selection
; is poor.
; Medsize: I filter out noise using a median filter of Medsize.
; Test: Stops every N frame to check the centroiding
; Cent_file: Writes a file containing centroids and fluxes.
; Defaults: If set, goes through without asking any questions
; Display: If set, each frame is displayed.
; Nosave: The default is to save the resulting events and images. If
; this is set, they are not saved.
; new_star: Forces selection of a star each time instead of reading from
; the default file.
; RESTRICTIONS:
; none
; NOTES:
;
;Modification history
;JM: Apr. 07, 2017
;JM: Apr. 10, 2017: Generalized for FUV or NUV
;JM: Apr. 11, 2017: Was not passing parameters.
;JM: Apr. 14, 2017: Was not taking starting frame into account
;JM: Apr. 16, 2017: Added write to parameters
;JM: May 13, 2017: Improvements in centroiding by using known offsets.
;JM: May 23, 2017: Version 3.1
;JM: Jun 10, 2017: Check to see if centroid stars are defined in header
;JM: Jun 10, 2017: Code cleanup and fixes
;JM: Jun 23, 2017: Corrected edge effect where the subarray was too small.
;JM: Jul 24, 2017: Added option to not display array.
;JM: Aug. 02, 2017: Explicitly print number of frames.
;JM: Aug. 03, 2017: Added option to quit.
;JM: Aug. 11, 2017: Nbin is redundant (params.fine_bin) so removed the option
;JM: Aug. 21, 2017: Fixed an inconsistency in passing offsets
;JM: Sep. 14, 2017: Fixed problem if the offsets were not defined.
;JM: Nov. 7, 2017: Cosmetic changes.
;JM: Nov. 25, 2017: Correct if a1 was not finite.
;Copyright 2016 Jayant Murthy
;
; Licensed under the Apache License, Version 2.0 (the "License");
; you may not use this file except in compliance with the License.
; You may obtain a copy of the License at
;
; http://www.apache.org/licenses/LICENSE-2.0
;
; Unless required by applicable law or agreed to in writing, software
; distributed under the License is distributed on an "AS IS" BASIS,
; WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
; See the License for the specific language governing permissions and
; limitations under the License.
;-

function set_limits, grid2, xstar, ystar, boxsize, resolution,$
xmin = xmin, ymin = ymin, display = display
siz = size(grid2)
ndim = siz[1]
xmin = xstar - boxsize*resolution
xmax = xstar + boxsize*resolution
ymin = ystar - boxsize*resolution
ymax = ystar + boxsize*resolution

xmin = 0 > xmin < (ndim - 1)
xmax = (xmin + 1) > xmax < (ndim - 1)
ymin = 0 > ymin < (ndim - 1)
ymax = (ymin + 1) > ymax < (ndim - 1)

if (keyword_set(display))then begin
plots,/dev, [xmin, xmin, xmax, xmax, xmin]/resolution,$
[ymin, ymax, ymax, ymin, ymin]/resolution,$
col=255,thick=2
endif
if (((xmax - xmin) lt 5) or ((ymax - ymin) lt 5)) then begin
h1 = fltarr(2*boxsize*resolution, 2*boxsize*resolution)
endif else h1 = grid2[xmin:xmax, ymin:ymax]
return,h1
end

pro jude_light_curve, events_file, out_file, $
xoff = xoff, yoff = yoff , params = params, $
nbin = nbin,$
max_im_value = max_im_value, mask = mask

;*********************** INITIALIZATION **************************
if (n_elements(nbin) eq 0)then $
read,"How many frames to bin over? ",nbin
;MAX_IM_VALUE for different images
if (n_elements(max_im_value) eq 0)then max_im_value = 0.00001 ;For display
device,window_state = window_state
if (window_state[0] eq 0)then $
window, 0, xs = 1024, ys = 512, xp = 10, yp = 500
if (n_elements(params) eq 0)then params = jude_params()
naxis = 512 * params.resolution
;************************* END INITIALIZATION ***************************

openw,write_lun,out_file,/get
;Read data
data_l2 = mrdfits(events_file,1,data_hdr0, /silent)
ndata_l2 = n_elements(data_l2)

;Set and check parameters
if (params.max_counts eq 0)then begin
q = where(data_l2.dqi eq 0, nq)
if (nq gt 10)then begin
dave = median(data_l2[q].nevents)
dstd = sqrt(dave)
params.max_counts = dave + dstd*3
endif else params.max_counts = 1000
endif
if (params.max_frame eq 0)then params.max_frame = ndata_l2 - 1
ngrid = (params.max_frame - params.min_frame)/nbin
ref_frame = sxpar(data_hdr0, "REFFRAME")

;If we haven't defined xoff and yoff set it to data_l2.xoff
if (n_elements(xoff) eq 0)then xoff = data_l2.xoff
if (n_elements(yoff) eq 0)then yoff = data_l2.yoff
q = where(abs(xoff) gt 500, nq)
if (nq gt 0)then xoff[q]= 0
q = where(abs(yoff) gt 500, nq)
if (nq gt 0)then yoff[q]= 0

;Select ROI
if (n_elements(mask) eq 0)then begin
nframes = jude_add_frames(data_l2, grid2, pixel_time, params, $
xoff*params.resolution, $
yoff*params.resolution, /notime, $
ref_frame = ref_frame)
ans = 'n'
while (ans eq 'n')do begin
tv,bytscl(rebin(grid2, 512, 512),0, max_im_value)
read,"Is the scaling ok? ",ans
if (ans eq 'n')then read, "New max value: ",max_im_value
endwhile
if (max(grid2) eq 0)then goto, noproc
print,"Please select midpoint of mask."
cursor,xstar,ystar,/dev
xstar = xstar*params.resolution
ystar = ystar*params.resolution
read, "Enter c for a circular mask; else square", ans
if (ans eq 'c')then begin
read,"Enter radius: ",boxsize
xmask = lindgen(naxis, naxis) mod naxis
ymask = lindgen(naxis, naxis)/naxis
dist = sqrt((xmask - xstar)^2 + (ymask - ystar)^2)
q = where(dist le boxsize, nq)
mask = fltarr(naxis, naxis)
if (nq gt 0)then mask[q] = 1
endif else begin
read,"Enter half width of box side: ", boxsize
mask = fltarr(naxis, naxis)
mask[(xstar-boxsize)>0:(xstar+boxsize)<(naxis - 1), $
(ystar-boxsize)>0:(ystar+boxsize)<(naxis - 1)] = 1
endelse
endif
stop
;Initial definitions for centroid region
xstar_first = xstar
ystar_first = ystar

;Loop through data to centroid
start_frame = params.min_frame
end_frame = params.max_frame
time0 = systime(1)
for i=1l, ngrid - 1 do begin
str = string(i*nbin)
str = str + " S/C time: "
str = str + string(data_l2[i*nbin].time - data_l2[0].time)
str = str + " Time Left: "
str = str + string((systime(1) - time0)/float(i)*float(ngrid - i))
print,strcompress(str), string(13b),format="(a,a,$)"
params.min_frame = start_frame + i*nbin
params.max_frame = params.min_frame + nbin - 1
dqi = where(data_l2[params.min_frame:params.max_frame].dqi eq 0,ndqi)
if (ndqi gt 3)then begin
nframes = jude_add_frames(data_l2, grid2, pixel_time, params, $
xoff*params.resolution, $
yoff*params.resolution, /notime, ref_frame = ref_frame)

tv,bytscl(rebin(grid2,512,512),0,max_im_value)

;Define a small array around star
carray = set_limits(grid2, xstar, ystar, boxsize, params.resolution, $
xmin = xmin, ymin = ymin, display = display)
siz = size(carray, /dimensions)
xmax = xmin + siz[0]
ymax = ymin + siz[1]
if (siz[0] lt 2)then stop
tv,bytscl(rebin(carray, siz[0]*(512/siz[0]),siz[1]*(512/siz[0])),0,max_im_value),512,0
plots,/dev,[xmin, xmin, xmax, xmax, xmin]/params.resolution,$
[ymin, ymax, ymax, ymin, ymin]/params.resolution,col=255
tcent = total(grid2*mask)
printf,write_lun, params.min_frame, $
tcent,sqrt(tcent)/sqrt(float(nframes)), $
nframes
endif else if (ndqi gt 3)then nlost = nlost + 1
endfor
noproc:
end

0 comments on commit c791a0e

Please sign in to comment.