diff --git a/README.md b/README.md deleted file mode 100644 index 6420074..0000000 --- a/README.md +++ /dev/null @@ -1,4 +0,0 @@ -# JUDE -Alternative pipeline for UVIT -Includes a number of GDL procedures to obtain Level 2 scientific products from the Level 1 UVIT data. -Reference is at https://ui.adsabs.harvard.edu/#abs/2016ascl.soft07007M/abstract diff --git a/jude_add_frames.pro b/jude_add_frames.pro index 874ada1..49b18f5 100644 --- a/jude_add_frames.pro +++ b/jude_add_frames.pro @@ -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. @@ -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"); @@ -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 diff --git a/jude_apply_astrometry.pro b/jude_apply_astrometry.pro new file mode 100644 index 0000000..9820981 --- /dev/null +++ b/jude_apply_astrometry.pro @@ -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 \ No newline at end of file diff --git a/jude_apply_time.pro b/jude_apply_time.pro new file mode 100644 index 0000000..3dcdf8a --- /dev/null +++ b/jude_apply_time.pro @@ -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 \ No newline at end of file diff --git a/jude_call_astrometry.pro b/jude_call_astrometry.pro new file mode 100644 index 0000000..eab3d91 --- /dev/null +++ b/jude_call_astrometry.pro @@ -0,0 +1,91 @@ +;+ +; NAME: JUDE_CALL_ASTROMETRY +; PURPOSE: Adds frames into grid for UVIT data +; CALLING SEQUENCE: +; jude_call_astrometry, inp_dir, out_dir, min_exp_time = min_exp_time, $ +; noupdate = noupdate +; INPUTS: +; inp_dir: Directory containing image files to be corrected +; OPTIONAL INPUTS: +; out_dir: Directory to put files from astrometry.net +; OPTIONAL KEYWORDS: +; min_exp_time: Ignore those with lower exp times +; noupdate: Don't update the files. +; OUTPUTS: +; Updates original image file with astrometry information. +; RESTRICTIONS: +; Astrometry.net must be installed. Some image files may not +; be recognized by astrometry.net +; NOTES: This program is only a shell which calls the +; local astrometry.net installation. +; +;Modification history +;JM: July 19, 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_call_astrometry, inp_dir, out_dir, min_exp_time = min_exp_time, $ + noupdate = noupdate, fuv = fuv, new = new + +;Check all the files + files=file_search(inp_dir, "*.fits*",count=nfiles) + if (n_elements(min_exp_time) eq 0)then min_exp_time = 10 + if (n_elements(out_dir) eq 0)then spawn, "mkdir ",out_dir + + + for ifile = 0, nfiles - 1 do begin + im = mrdfits(files[ifile], 0, im_hdr, /silent) + exp_time = sxpar(im_hdr, "EXP_TIME") + ra_pnt = sxpar(im_hdr, "RA_PNT") + dec_pnt = sxpar(im_hdr, "DEC_PNT") + if (keyword_set(new))then astr_done = "FALSE" else $ + astr_done = strcompress(sxpar(im_hdr, "ASTRDONE"),/rem) + + if ((exp_time gt min_exp_time) and (astr_done ne "TRUE"))then begin + str = "solve-field" + if (keyword_set(fuv))then $ + str = "solve-field --backend-config /Volumes/UVIT_Data/astrometric/astrometry.cfg" + str = str + " --downsample 4 --scale-units" + str = str + " + str = str + " degwidth --scale-low .4 --scale-high .6 " + str = str + " --no-plots --continue" + str = str + " --dir " + out_dir + if ((ra_pnt ne 0) and (dec_pnt ne 0))then begin + str = str + " --ra " + string(ra_pnt) + " --dec " + str = str + string(dec_pnt) + " --radius 1 " + endif + str = strcompress(str + " " + files[ifile] + " > solve.txt") + print,"Beginning solve of file no ", ifile, " at ",systime(0) + spawn,str + print,"Ending solve at ",systime(0) + + if (not(keyword_set(noupdate)))then begin + fname = file_basename(files[ifile]) + fname = strmid(fname,0,strlen(fname)-8) + new_file = out_dir + "/" + fname + ".new" + tst = file_test(new_file) + + if (tst gt 0)then begin + out_file = inp_dir + fname + ".fits" + times = mrdfits(files[ifile], 1, thdr, /silent) + asthdr = headfits(new_file, /silent) + sxaddpar, asthdr, "ASTRDONE", "TRUE", "Is astrometry done" + mwrfits, im, out_file, asthdr, /create + mwrfits, times, out_file, thdr + spawn, "gzip -fv " + out_file + endif + endif + endif + endfor +end diff --git a/jude_centroid.pro b/jude_centroid.pro index bded8fa..3bf7de9 100644 --- a/jude_centroid.pro +++ b/jude_centroid.pro @@ -51,6 +51,7 @@ ;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. ;Copyright 2016 Jayant Murthy ; ; Licensed under the Apache License, Version 2.0 (the "License"); @@ -161,6 +162,7 @@ pro jude_centroid, events_file, grid2, params, xstar, ystar, $ ystar = sxpar(data_hdr0, "YCENT", count = nystar) xstar = xstar*params.resolution ystar = ystar*params.resolution + if ((nxstar eq 0) or (nystar eq 0))then begin x = 0 thresh = .0005 @@ -172,6 +174,7 @@ pro jude_centroid, events_file, grid2, params, xstar, ystar, $ s = reverse(sort(flux)) xstar = x[s[0]] ystar = y[s[0]] + endif ans = 'n' while (ans eq 'n') do begin @@ -186,8 +189,12 @@ pro jude_centroid, events_file, grid2, params, xstar, ystar, $ tv,bytscl(rebin(h1,siz[1]*5, siz[2]*5), 0, max_im_value), 512, 0 print,"Width is ",a1[2], a1[3] ans = 'y' - if (not(keyword_set(defaults)))then read,"Is this star ok? ", ans - if (ans ne 'n')then ans = 'y' + if (not(keyword_set(defaults)))then $ + read,"Is this star ok (s for single frame)? ", ans + if (ans eq 'd')then begin + display = 0 + ans = 'y' + endif if (ans eq 'n')then begin print,"Select star" cursor,a,b,/dev @@ -196,6 +203,31 @@ pro jude_centroid, events_file, grid2, params, xstar, ystar, $ print,"boxsize is now: ", boxsize read,"Enter new boxsize: ",boxsize endif + if (ans eq 's')then begin + par = params + endif + while (ans eq 's')do begin + g2 = grid2*0 + while (max(g2) eq 0)do begin + par.min_frame = par.min_frame + par.fine_bin + par.max_frame = par.min_frame + par.fine_bin + print,par.min_frame,string(13b),format="(i6,a,$)" + nframes = jude_add_frames(data_l2, g2, pixel_time, par, $ + xoff*params.resolution, $ + yoff*params.resolution, /notime) + endwhile + tv,bytscl(rebin(g2,512,512),0,max_im_value) + read, "Is this frame ok? ", ans + if (ans eq 'y')then begin + print,"Select star" + cursor,a,b,/dev + xstar = a*params.resolution + ystar = b*params.resolution + print,"boxsize is now: ", boxsize + read,"Enter new boxsize: ",boxsize + ans = 'n' + endif else ans = 's' + endwhile endwhile q = where(finite(a1) eq 0, nq) if ((xstar eq 0) or (nq gt 0)) then begin @@ -298,7 +330,8 @@ pro jude_centroid, events_file, grid2, params, xstar, ystar, $ params.max_frame = ndata_l2 - 1 nframes = jude_add_frames(data_l2, grid2, pixel_time, params, $ xoff, yoff, /notime) - if (display eq 1)then tv,bytscl(rebin(grid2,512,512),0,max_im_value) + if ((display eq 1) and (max(grid2) gt 0))then $ + tv,bytscl(rebin(grid2,512,512),0,max_im_value) h1 = set_limits(grid2, xstar_first, ystar_first, boxsize, params.resolution) r1 = mpfit2dpeak(h1, a1) siz = size(h1, /dimensions) @@ -314,7 +347,7 @@ pro jude_centroid, events_file, grid2, params, xstar, ystar, $ sxaddpar,data_hdr0, "XCENT", xstar_first/params.resolution, "XPOS of centroid star" sxaddpar,data_hdr0, "YCENT", ystar_first/params.resolution, "YPOS of centroid star" temp_file = strmid(events_file,0,strlen(events_file)-3) - print,"writing image file to ", temp_file + print,"writing events file to ", temp_file mwrfits,data_l2,temp_file,data_hdr0,/create,/no_comment mwrfits,offsets,temp_file,off_hdr,/no_comment diff --git a/jude_driver_uv.pro b/jude_driver_uv.pro index 5ae6d58..8faf354 100644 --- a/jude_driver_uv.pro +++ b/jude_driver_uv.pro @@ -80,6 +80,7 @@ ; JM: Dec. 11, 2016 : Cleaning up ; JM: May 23, 2017 : Version 3.1 ; JM: Jun 27, 2017 : Added overwrite option +; JM: Jul 27, 2017 : Add reference frame to files for astrometry ;Copyright 2016 Jayant Murthy ; ; Licensed under the Apache License, Version 2.0 (the "License"); @@ -372,10 +373,10 @@ endif ;The notimes is much faster if (keyword_set(notime))then begin nframes = JUDE_ADD_FRAMES(data_l2, grid, pixel_time, par, $ - xoff, yoff, /notime) + xoff, yoff, ref_frame = ref_frame, /notime) endif else begin nframes = JUDE_ADD_FRAMES(data_l2, grid, pixel_time, par, $ - xoff, yoff) + xoff, yoff, ref_frame = ref_frame) endelse ;Write PNG file @@ -405,6 +406,10 @@ endif ;Starting and ending frames for image production sxaddpar, out_hdr,"MINFRAME", params.min_frame,"Starting frame" sxaddpar, out_hdr,"MAXFRAME", params.max_frame,"Ending frame" + sxaddpar, out_hdr, "REFFRAME", ref_frame, "Reference frame." + sxaddpar, bout_hdr,"MINFRAME", params.min_frame,"Starting frame" + sxaddpar, bout_hdr,"MAXFRAME", params.max_frame,"Ending frame" + sxaddpar, bout_hdr, "REFFRAME", ref_frame, "Reference frame." ;Calibration factor cal_factor = jude_apply_cal(detector, nom_filter) diff --git a/jude_get_xy.pro b/jude_get_xy.pro index 30c7ed4..2c7f7de 100644 --- a/jude_get_xy.pro +++ b/jude_get_xy.pro @@ -37,6 +37,7 @@ ; JM: Aug 10, 2016: Integer overflow. ; JM: Aug 21, 2016: Add filter information to files. ; JM: May 23, 2017: Version 3.1 +; JM: Jul 21, 2017: Added RA and DEC to Level 2 data format. ; COPYRIGHT: ;Copyright 2016 Jayant Murthy ; @@ -91,6 +92,8 @@ function jude_get_xy,data_l1, data_l1a, data_l2, out_hdr nevents:0, $ ;Number of events in the frame x:fltarr(nevents), $ ;X position y:fltarr(nevents), $ ;Y position + ra:fltarr(nevents), $ ;RA of every event + dec:fltarr(nevents),$ ;Dec of every event mc:intarr(nevents), $ ;minimum value in centroid dm:intarr(nevents), $ ;maximum value in centroid time: 0d, $ ;Time from L1 file diff --git a/jude_interactive.pro b/jude_interactive.pro index 278574f..a7b02c0 100644 --- a/jude_interactive.pro +++ b/jude_interactive.pro @@ -34,6 +34,8 @@ ; JM: Jun 10, 2017 : Save star centroids, if defined ; JM: Jun 11, 2017 : Changes in centroid calls. ; JM: Jun 23, 2017 : Handled case where centroids were not defined. +; JM: Jul 21, 2017 : Was writing centroids wrong. +; JM: Jul 27, 2017 : Add reference frame to files for astrometry ;Copyright 2016 Jayant Murthy ; ; Licensed under the Apache License, Version 2.0 (the "License"); @@ -216,6 +218,12 @@ pro plot_diagnostics, data_l2, offsets, data_hdr0, im_hdr, fname, grid, $ if (nq gt 1)then oplot,offsets[q].time - data_l2[0].time, yoff_vis[q], $ col=65535,psym=3,symsize=3 +;What is the time per frame? + ndata = n_elements(data) + dtimes = (data_l2[1:ndata-1].time - data_l2[0:ndata-2].time) + h=histogram(dtimes,min=0,bin=.00001,max=.1) + dtime = where(h eq max(h))*.00001 + ;Print diagnostics l2_tstart = data_l2[0].time ndata_l2 = n_elements(data_l2) @@ -223,10 +231,11 @@ pro plot_diagnostics, data_l2, offsets, data_hdr0, im_hdr, fname, grid, $ diff = l2_tend - l2_tstart if (n_elements(im_hdr) gt 0) then exp_time = sxpar(im_hdr, "exp_time") $ else exp_time = 0 + good = where(data_l2.dqi eq 0, ngood) + str = fname - str = str + " " + string(long(l2_tstart)) - str = str + " " + string(long(l2_tend)) str = str + " " + string(long(diff)) + str = str + " " + string(ngood * dtime) str = str + " " + string(long(exp_time)) str = str + string(mode) str = str + " " + string(fix(total(grid))) @@ -298,27 +307,34 @@ pro jude_interactive, data_file, uv_base_dir, data_l2, grid, offsets, params = p events_dir = uv_base_dir + params.events_dir png_dir = uv_base_dir + params.png_dir image_file = image_dir + fname + ".fits.gz" - - + ;Read the VIS offsets if they exist. If they don't exist, set up a ;dummy array and header. - offsets = mrdfits(data_file, 2, off_hdr) + offsets = mrdfits(data_file, 2, off_hdr, /silent) if (n_elements(offsets) le 1)then begin offsets = replicate({offsets, time:0d, xoff:0., yoff:0., att:0}, ndata_l2) offsets.time = data_l2.time offsets.att = 1 fxbhmake,off_hdr,ndata_l2,/initialize sxaddhist,"No offsets from visible",off_hdr, /comment + print,"No visible offsets available." endif detector = strcompress(sxpar(data_hdr0,"detector"), /remove) calc_uv_offsets, offsets, xoff_vis, yoff_vis, detector xoff_uv = data_l2.xoff yoff_uv = data_l2.yoff -;Create image - nframes = jude_add_frames(data_l2, grid, pixel_time, params, $ +;Does the image exist + if (file_test(image_file) gt 0)then begin + print,"Reading image from file." + grid = mrdfits(image_file, 0, im_hdr) + if (strcompress(sxpar(im_hdr, "ASTRDONE"), /rem) eq "TRUE") then $ + print,"Astrometry done." + endif else begin + nframes = jude_add_frames(data_l2, grid, pixel_time, params, $ xoff_uv*params.resolution, yoff_uv*params.resolution,$ /notime) + endelse check_diag: ;Plot the image plus useful diagnostics @@ -427,6 +443,7 @@ if (ans ne "d")then begin ;****************************** CENTROIDING *************************** if (run_centroid eq 'y')then begin +print,"Starting centroid" nbin = params.fine_bin if (strupcase(detector) eq "FUV")then $ thg1 = params.ps_threshold_fuv else $ @@ -476,7 +493,8 @@ if (ans ne "d")then begin ;Until getting the time issues sorted out I will leave this as no time calculation. data_l2.dqi = dqi nframes = jude_add_frames(data_l2, grid, pixel_time, params, $ - xoff_sc*params.resolution, yoff_sc*params.resolution, /notime) + xoff_sc*params.resolution, yoff_sc*params.resolution, /notime,$ + ref_frame = ref_frame) ;File definitions fname = file_basename(data_file) @@ -518,6 +536,7 @@ if (ans ne "d")then begin 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" + sxaddpar, out_hdr, "REFFRAME", ref_frame, "Reference frame." sxaddhist,"Times are in Extension 1", out_hdr, /comment if (run_centroid eq 'y')then $ sxaddhist, "Centroiding run on image.", out_hdr @@ -540,14 +559,15 @@ if (ans ne "d")then begin "Recommended starting frame" sxaddpar, data_hdr0,"MAXFRAME", params.max_frame,$ "Recommended ending frame" + sxaddpar, data_hdr0, "REFFRAME", ref_frame, "Reference frame." sxaddpar, data_hdr0, "CALF", cal_factor, $ "Ergs cm-2 s-1 A-1 (cps)-1" if ((run_centroid eq 'y') and $ (n_elements(xcent) gt 0) and (n_elements(ycent) gt 0))then begin sxaddhist, "Centroiding run for s/c motion correction.", data_hdr0 - sxaddpar,data_hdr0, "XCENT", xcent/params.resolution, "XPOS of centroid star" - sxaddpar,data_hdr0, "YCENT", ycent/params.resolution, "YPOS of centroid star" + sxaddpar,data_hdr0, "XCENT", xcent, "XPOS of centroid star" + sxaddpar,data_hdr0, "YCENT", ycent, "YPOS of centroid star" endif mwrfits,temp,t,data_hdr0,/create,/no_comment diff --git a/jude_mosaic.pro b/jude_mosaic.pro new file mode 100644 index 0000000..006d09f --- /dev/null +++ b/jude_mosaic.pro @@ -0,0 +1,132 @@ ++ +; NAME: JUDE_MOSAIC +; PURPOSE: Creates mosaic of desired region from multiple UVIT data sets. +; CALLING SEQUENCE: +; jude_mosaic, events_dir, out_file, ra_cent, dec_cent, fov, pixel_size +; INPUTS: +; events_dir: Directory containing multiple event files. The +; files must have astrometric information. +; out_file: FITS image file +; OPTIONAL INPUTS: +; ra_cent: Central RA in degrees +; dec_Cent: Central Dec in degrees +; Fov: Field of view in degrees +; pixel_size: Size of each pixel in degrees +; OPTIONAL KEYWORDS: +; OUTPUTS: +; RESTRICTIONS: +; NOTES: +; +;Modification history +;JM: July 19, 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_mosaic, events_dir, out_file, ra_cent, dec_cent, fov, pixel_size + +;Read parameters if they are not set + if (n_elements(ra_cent) eq 0) then read, $ + "Please enter field center: ", ra_cent, dec_cent + if (n_elements(fov) eq 0)then read, $ + "Please enter field width: ", fov + if (n_elements(pixel_size) eq 0)then read, $ + "Please enter pixel size: ", pixel_size + +;The UVIT field is 28 arcminutes in diameter. This is used in calculating the +;exposure time + instr_size = 14d/60./pixel_size + +;Set up WCS parameters + naxis = round(2*fov/pixel_size) + crpix = double([naxis/2, naxis/2]) + crval = double([ra_cent, dec_cent]) + delt = double([-pixel_size, pixel_size]) + make_astr, astr, delta = delt, crpix = crpix, crval = crval + +;Output variables + grid = fltarr(naxis, naxis) + times = fltarr(naxis, naxis) + xaxis = lindgen(naxis, naxis) mod naxis + yaxis = lindgen(naxis, naxis) / naxis + xref_save = 0 + yref_save = 0 + +;List of event files + files = file_search(events_dir, "*.fits.gz", count = nfiles) + for ifile = 0, nfiles - 1 do begin + data_l2 = mrdfits(files[ifile], 1, d2_hdr, /silent) + ndata = n_elements(data_l2) + +;Select good data + q = where((data_l2.dqi eq 0) and (data_l2.nevents gt 0), ndata) + if (ndata gt 0)then begin + data_l2 = data_l2[q] + +;What is the time per frame? + dtimes = (data_l2[1:ndata-1].time - data_l2[0:ndata-2].time) + h=histogram(dtimes,min=0,bin=.00001,max=.1) + dtime = where(h eq max(h))*.00001 + print,"Using ", files[ifile]," with timing: ",dtime + endif else print,"Not using ",files[ifile] + +;Run through each frame + for idata = 0, ndata - 1 do begin + print,idata,string(13b),format="(i5, a, $)" + ra = data_l2[idata].ra + dec = data_l2[idata].dec + q = where((ra ne 0) and (dec ne 0), nq) + +;If there is good data + if (nq gt 0)then begin + ra = ra[q] + dec = dec[q] + ad2xy,ra,dec,astr,x,y + x = round(x) + y = round(y) + q = where((x ge 0) and (x lt naxis) and (y ge 0) and (y lt naxis), nq) + +;If the data are within the grid + if (nq gt 0)then begin + x = x[q] + y = y[q] +;Add events into the grid + for i = 0l, nq - 1 do grid[x[i], y[i]] = grid[x[i], y[i]] + 1 +;Boresight position + ra_ref = data_l2[idata].roll_ra + dec_ref = data_l2[idata].roll_dec + ad2xy, ra_ref, dec_ref, astr, xref, yref + +;Find where the distance is within the UVIT field + if ((round(xref) ne xref_save) or (round(yref) ne yref_save))then begin +;This saves time in case the boresight doesn't move + xref_save = round(xref) + yref_save = round(yref) + dst = (xaxis - xref_save)^2 + (yaxis - yref_save)^2 + qdst = where(dst lt instr_size^2, ndst) + endif +;Set exposure time + if (ndst gt 0) then times(qdst) = times[qdst] + dtime + endif + endif + endfor + endfor + +;Write out file + q = where(times gt 0, nq) + if (nw gt 0)then grid = grid/times + mkhdr, out_hdr, grid + putast, out_hdr, astr + mwrfits, grid, out_file, out_hdr, /create + mwrfits, times, out_file + +end \ No newline at end of file diff --git a/jude_obs_log.pro b/jude_obs_log.pro index be3c6a4..28471bf 100644 --- a/jude_obs_log.pro +++ b/jude_obs_log.pro @@ -36,7 +36,7 @@ pro jude_obs_log, output_file, uv_base_dir, params ;Which files are available files = file_search(uv_base_dir + params.events_dir, "*", count = nfiles) ;Header titles - printf,obs_lun,"Level2 Source Filter OBSDATE TSTART TEND RA DEC Exp_Time IM_Exp" + printf,obs_lun,"Level2 Source Filter OBSDATE TSTART TEND RA DEC Exp_Time Good_time IM_Exp" for ifile = 0, nfiles - 1 do begin print,ifile," of ",nfiles-1,string(13b),format="(i5,a,i5, a, $)" @@ -48,14 +48,22 @@ pro jude_obs_log, output_file, uv_base_dir, params im_name = uv_base_dir + params.image_dir + im_name + ".fits.gz" q = where(data_l2.dqi eq 0, nq) if (nq gt 0)then begin - tstart = min(data_l2.time) - tend = max(data_l2.time) + tstart = min(data_l2[q].time) + tend = max(data_l2[q].time) +;What is the time per frame? + ndata = n_elements(data_l2) + dtimes = (data_l2[1:ndata-1].time - data_l2[0:ndata-2].time) + h=histogram(dtimes,min=0,bin=.00001,max=.1) + dtime = where(h eq max(h))*.00001 + if (n_elements(dtime) gt 1)then dtime=reform(dtime[1],1) endif else begin tstart = -1 tend = -1 + dtime = 0 printf,rm_lun,"rm " + l2_file printf,rm_lun,"rm " + im_name endelse + str = l2_file str = str + " " + strcompress(string(sxpar(hdr, "sourceid")),/remove) str = str + " " + string(sxpar(hdr, "filter")) @@ -64,13 +72,15 @@ pro jude_obs_log, output_file, uv_base_dir, params str = str + " " + string(tend, form = "(f15.4)") str = str + " " + string(sxpar(hdr, "ra_pnt")) str = str + " " + string(sxpar(hdr, "dec_pnt")) - str = str + " " + string(tend - tstart, form = "(f15.4)") + str = str + " " + string(tend - tstart, form = "(f15.2)") + str = str + " " + string(nq*dtime, form = "(f15.2)") if (file_test(im_name))then begin im = mrdfits(im_name, 0, im_hdr,/silent) str = str + " " + string(sxpar(im_hdr, "exp_time")) endif else str = str + " 0" str = strcompress(str) printf,obs_lun,str + endfor free_lun,obs_lun free_lun,rm_lun diff --git a/write_astrometry_dot_net.pro b/write_astrometry_dot_net.pro deleted file mode 100644 index 131d2e9..0000000 --- a/write_astrometry_dot_net.pro +++ /dev/null @@ -1,26 +0,0 @@ -;Assuming that Astrometry.net is locally installed, this will write a shell file -;to process all the image files in a directory -;We have to downsample by some factor which, empirically, we choose to be 4 -;(1024x1024 arrays) - -f=file_search("*.fits.gz",count=nf) -openw,1,"solve.sh" -for i = 0, nf-1 do begin - im = mrdfits(f[i], 0, hdr) - exp_time = sxpar(hdr, "EXP_TIME") - ra_pnt = sxpar(hdr, "RA_PNT") - dec_pnt = sxpar(hdr, "DEC_PNT") - if (exp_time gt 10)then begin - str = "solve-field --downsample 2 --scale-units" - str = str + " degwidth --scale-low .2 --scale-high .75 " - str = str + " --no-plots --continue" - if ((ra_pnt ne 0) and (dec_pnt ne 0))then begin - str = str + " --ra " + string(ra_pnt) + " --dec " - str = str + string(dec_pnt) + " --radius 10 " - endif - str = str + " " + f[i] - printf,1,strcompress(str) - endif -endfor -close,1 -end