Skip to content

Commit

Permalink
Added analysis files
Browse files Browse the repository at this point in the history
1) Changed Level 2 event files to include RA and Dec.
2) Added astrometry through astrometry.net.
3) Minor corrections
  • Loading branch information
jaymurthy committed Aug 2, 2017
1 parent 5a8647d commit 6c6a884
Show file tree
Hide file tree
Showing 12 changed files with 609 additions and 51 deletions.
4 changes: 0 additions & 4 deletions README.md

This file was deleted.

6 changes: 5 additions & 1 deletion jude_add_frames.pro
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,8 @@
; spacecraft.
; REF_FRAME: Because of spacecraft motion, the frame chosen by the program
; as the reference frame may not be optimal. This allows the
; frame to be changed to another frame number.
; frame to be changed to another frame number. Also used in
; marking which frame is used as the reference frame.
; APPLY_DIST: We have derived a distortion matrix which may be applied to
; every frame. I recommend that we apply distortion only to the
; coadded image.
Expand All @@ -57,6 +58,7 @@
;JM: Apr. 10, 2017: Fixed problem with choice for ref frame.
;JM: May 1, 2017: Option to apply distortion coefficients.
;JM: May 22, 2017: Version 3.1
;JM: Jul 27, 2017: Corrected reference frame.
;Copyright 2016 Jayant Murthy
;
; Licensed under the Apache License, Version 2.0 (the "License");
Expand Down Expand Up @@ -110,6 +112,8 @@ if (min_frame eq 0)then min_frame = $
(abs(yoff) lt 1000)))

if (n_elements(ref_frame) eq 0)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))) do ref_frame = ref_frame + 1
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
157 changes: 157 additions & 0 deletions jude_apply_astrometry.pro
Original file line number Diff line number Diff line change
@@ -0,0 +1,157 @@
;+
; NAME: JUDE_APPLY_ASTROMETRY
; PURPOSE: Puts astrometry information from image into photon list
; CALLING SEQUENCE:
; jude_apply_astrometry, l2_file, params, iamge_dir = image_dir, $
; image_file = image_file
; INPUTS:
; l2_file: Events list from Level 2 data
; OPTIONAL INPUTS:
; Params: parameters for pipeline operation
; OPTIONAL KEYWORDS:
; Image_dir: Directory in which to look for image files.
; Image_file: If image file has a different name than expected
; OUTPUTS:
; RESTRICTIONS:
; none
; NOTES:
; The astrometry is taken from the image and fed back into the
; events file using the x and y offsets. This will associate an
; RA and Dec with every photon.
;
;Modification history
;JM: July 21, 2017
;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.
;-
pro jude_apply_astrometry, l2_file, params, data_dir = data_dir, new = new

;Begin by reading from Level 2 files
if (file_test(l2_file))then begin
data_l2 = mrdfits(l2_file, 1, hdr_l2, /silent)
offsets = mrdfits(l2_file, 2, hdr_off, /silent)
endif else begin
print, "Could not find ", l2_file
goto, noproc
endelse

if (keyword_set(new))then astr_done = "FALSE" else $
astr_done = strcompress(sxpar(hdr_l2, "ASTRDONE"),/rem)
if (astr_done eq "TRUE")then begin
print,"Astrometry already done."
goto, noproc
endif
;If the parameters are not set, we use the defaults.
if (n_elements(params) eq 0)then params = jude_params()

;Finding the default location for the image files.
detector = strcompress(sxpar(hdr_l2, "DETECTOR"), /rem)
if (detector eq "NUV")then uv_base_dir = params.def_nuv_dir $
else uv_base_dir = params.def_fuv_dir
if (n_elements(data_dir) gt 0)then uv_base_dir = data_dir + uv_base_dir

;Name definitions
fname = file_basename(l2_file)
f1 = strpos(fname, "level1")
f2 = strpos(fname, "_", f1+8)
fname = strmid(fname, 0, f2)
if (n_elements(image_dir) eq 0)then $
image_dir = uv_base_dir + params.image_dir
if (n_elements(image_file) eq 0)then $
image_file = image_dir + fname + ".fits.gz"

;Read image file if it exists; if it doesn't we don't have astrometric info
if (file_test(image_file))then begin
grid = mrdfits(image_file, 0, hdr_im, /silent)
endif else begin
print, "Could not find ", image_file
goto, noproc
endelse

;Check resolution. This is primarily used as a sanity check. I have to know
;the image resolution to convert into the coordinates for the photon list so
;I'm just confirming that it is what I think it is.
resolution = params.resolution
if ((resolution * 512) ne (sxpar(hdr_im, "NAXIS1")))then begin
print,"Mismatch in coordinates"
goto, noproc
endif

;Check to make sure astrometry has been done on the image
astr_done = strcompress(sxpar(hdr_im, "ASTRDONE"),/rem)
if (astr_done ne "TRUE")then begin
print,"no astrometry done on image"
goto, noproc
endif

;We have to use the same parameters as were used in the image production. The
;coordinate transformation depends on the coordinate transformation from x and
;y which is dependent on the image details
min_frame = sxpar(hdr_im, "MINFRAME")
max_frame = sxpar(hdr_im, "MAXFRAME")
xoff = data_l2.xoff * params.resolution
yoff = data_l2.yoff * params.resolution

;This is the code from jude_add_frames so I reproduce it here.
if (min_frame eq 0)then min_frame = $
min(where((data_l2.dqi eq 0) and (abs(xoff) lt 1000) and $
(abs(yoff) lt 1000)))

if (max_frame eq 0)then $
max_frame = params.max_frame < (n_elements(data_l2)-1)
;This is the reference name forr the x and y offsets.
ref_frame = min_frame

;Read astrometric data from image file
extast,hdr_im, astr

;Initialization
xoff_start = xoff[ref_frame]
yoff_start = yoff[ref_frame]

;Run through the photon list to get the coordinates for every photon.
for ielem = 0, n_elements(data_l2) - 1 do begin
if (data_l2[ielem].dqi eq 0)then begin

;Translate each event into its position in the image grid
q = where(data_l2[ielem].x gt 0,nq)
x = data_l2[ielem].x(q)*resolution + xoff(ielem) - xoff_start
y = data_l2[ielem].y(q)*resolution + yoff(ielem) - yoff_start
;I define the boresight as the central pixel in the events list.
xref = 256d*resolution + xoff[ielem] - xoff_start
yref = 256d*resolution + yoff[ielem] - yoff_start

;Get the RA and Dec of each photon and feed it back into the events list.
xy2ad, x, y, astr, r, d
data_l2[ielem].ra[q] = r
data_l2[ielem].dec[q] = d
;RA and Dec of the boresight.
xy2ad, xref, yref, astr, ra_ref, dec_ref
data_l2[ielem].roll_ra = ra_ref
data_l2[ielem].roll_dec = dec_ref
endif
endfor

;Write modified data out. Overwrite the original data.
fname = strmid(l2_file, 0, strlen(l2_file) - 3)
sxaddpar, hdr_l2, "ASTRDONE", "TRUE", "Astrometry done"
sxaddhist, "ROLL_RA,ROLL_DEC contain central ra, dec",hdr_l2,/comment
mwrfits, data_l2, fname, hdr_l2, /create, /no_comment
mwrfits, offsets, fname, hdr_off, /no_comment
spawn, "gzip -fv " + fname

;Go here if the file cannot be processed for any reason.
noproc:

end
133 changes: 133 additions & 0 deletions jude_apply_time.pro
Original file line number Diff line number Diff line change
@@ -0,0 +1,133 @@
;+
; NAME: JUDE_APPLY_TIME
; PURPOSE: Apply times to iamges.
; CALLING SEQUENCE:
; jude_apply_time, data_file, uv_base_dir
; INPUTS:
; Data_file :The photon list to be processed (Level 2 file).
; UV_Base_Dir :The top level UV directory (typically FUV or NUV)
; OUTPUTS:
; OPTIONAL INPUT KEYWORDS:
; NOTES:
; Data files are written out as per the original names
; MODIFICATION HISTORY:
; JM: July 12, 2017
;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.
;
;-

pro jude_apply_time, data_file, uv_base_dir
;Define bookkeeping variables
exit_success = 1
exit_failure = 0
version_date = "July 11, 2017"
print,"Software version: ",version_date

;**************************INITIALIZATION**************************

;The parameters are read using JUDE_PARAMS
;Assuming the path is set correctly, a personalized file can be in
;the current directory.
if (n_elements(params) eq 0)then $
params = jude_params()
;Error log: will append to existing file.
jude_err_process,"errors.txt",data_file

;************************LEVEL 2 DATA *********************************
data_l2 = mrdfits(data_file,1,data_hdr0,/silent)
ndata_l2 = n_elements(data_l2)

;if the keywords exist, read the starting and ending frame from the
;data file
params.min_frame = sxpar(data_hdr0, "MINFRAME")
params.max_frame = sxpar(data_hdr0, "MAXFRAME")
if ((params.max_frame eq 0) or (params.max_frame gt (ndata_l2 -1)))then $
params.max_frame = ndata_l2 - 1
start_frame = params.min_frame
end_frame = params.max_frame
save_dqi = data_l2.dqi
dqi = data_l2.dqi

;Calculate the median from the data. This is photon counting so sigma
; = sqrt(median). I allow 5 sigma.
q = where((data_l2.dqi eq 0) and (data_l2.nevents gt 0), nq)

if (nq gt 10)then begin
dave = median(data_l2[q].nevents)
dstd = sqrt(dave)
params.max_counts = dave + dstd*5

;Name definitions
fname = file_basename(data_file)
f1 = strpos(fname, "level1")
f2 = strpos(fname, "_", f1+8)
fname = strmid(fname, 0, f2)
image_dir = uv_base_dir + params.image_dir
image_file = image_dir + fname + ".fits.gz"

xoff_uv = data_l2.xoff
yoff_uv = data_l2.yoff

;Create image
nframes = jude_add_frames(data_l2, grid, pixel_time, params, $
xoff_uv*params.resolution, yoff_uv*params.resolution)

print,"Total of ",nframes," frames ",nframes*.035," seconds"

;File definitions
fname = file_basename(data_file)
fname = strmid(fname,0,strlen(fname)-8)
imname = file_basename(image_file)
imname = strmid(imname, 0, strlen(imname) - 8)
if (file_test(image_dir) eq 0)then spawn,"mkdir " + image_dir

;Make the basic header
mkhdr, out_hdr, grid
jude_create_uvit_hdr,data_hdr0,out_hdr

;Write FITS image file
nom_filter = strcompress(sxpar(out_hdr, "filter"),/remove)
sxaddpar,out_hdr,"NFRAMES",nframes,"Number of frames"

;Calibration factor
cal_factor = jude_apply_cal(detector, nom_filter)
sxaddpar, out_hdr, "CALF", cal_factor, "Ergs cm-2 s-1 A-1 (cps)-1"

;Check the exposure time
q = where(data_l2.dqi eq 0, nq)
if (nq gt 0)then begin
avg_time = $
(max(data_l2[q].time) - min(data_l2[q].time))/(max(q) - min(q))
endif else avg_time = 0

sxaddpar,out_hdr,"EXP_TIME",nframes * avg_time, "Exposure Time in seconds"
print,"Total exposure time is ",nframes * avg_time
nom_filter = nom_filter[0]
sxaddpar,out_hdr,"FILTER",nom_filter
sxaddpar,out_hdr,"MINEVENT",params.min_counts,"Counts per frame"
sxaddpar,out_hdr,"MAXEVENT",params.max_counts,"Counts per frame"
sxaddpar, out_hdr,"MINFRAME", params.min_frame,"Starting frame"
sxaddpar, out_hdr,"MAXFRAME", params.max_frame,"Ending frame"
sxaddhist,"Times are in Extension 1", out_hdr, /comment
if (run_centroid eq 'y')then $
sxaddhist, "Centroiding run on image.", out_hdr
sxaddhist,fname,out_hdr
t = uv_base_dir + params.image_dir + imname + ".fits"
print,"writing image file to ",t
mwrfits,grid,t,out_hdr,/create
mwrfits,pixel_time,t
spawn,"gzip -fv " + t
noproc:
end
Loading

0 comments on commit 6c6a884

Please sign in to comment.