;+ ; NAME: ; even_light ; PURPOSE: ; Program used for the certification of the SMEI cameras using ; the light box data take in the 3th floor clean room. ; CATEGORY: ; stuff ; CALLING SEQUENCE: ; INPUTS: ; OPTIONAL INPUT PARAMETERS: ; OUTPUTS: ; OPTIONAL OUTPUT PARAMETERS: ; CALLS: ; InitVar, gridgen, Inside_Wedge, SuperArray, BadValue ; flat_centerofmass, wedge_bounding_box, TimeSystem, ArrayLocation ; bin_read, bin_write, destroyvar, TimeOp, TimeUnit ; PROCEDURE: ; MODIFICATION HISTORY: ; ???-200?, Paul Hick (UCSD/CASS; pphick@ucsd.edu) ;- pro even_light_geometric_fix, zval, zphi, optaxis0, doptaxis @compile_opt.pro ; On error, return to caller zphi = temporary(zphi)-optaxis0[0] zval = temporary(zval)/cos(zphi) zval = temporary(zval)*(1+2*doptaxis[0]*sin(zphi)) return & end pro even_light_view, center, phi, r @compile_opt.pro ; On error, return to caller cphi = cos(phi) sphi = sin(phi) tphi = transpose([[cphi],[sphi]]) cntr = center#[1,1] plots, cntr+[cphi[0],sphi[0]]#[0,r[1]] , color=255 plots, cntr+[cphi[1],sphi[1]]#[0,r[1]] , color=255 plots, cntr+tphi*([1,1]# r ) , color=255 plots, cntr+tphi*([1,1]#reverse(r)) , color=255 return & end pro even_light_firstguess, spot_centroid, sym_spot_nr, good_spot, bins_axis, center0, optaxis0, binwidth, dbinwidth @compile_opt.pro ; On error, return to caller x = reform( spot_centroid[0,good_spot] ) y = reform( spot_centroid[1,good_spot] ) p = lsqCircleFit(x,y) center0 = center0+p[0:1] ; Center coordinates optaxis0 = p[2] ; Radius p = transpose([[x-p[0]], [y-p[1]]]) p = cv_coord(from_rect=p, /to_polar) p = reform(p[0,*]) ; Azimuth for all spots in image p = lsqQuadraticFit(sym_spot_nr[good_spot], p) ;p = poly_fit(sym_spot_nr[good_spot],p,2) optaxis0 = [p[0], optaxis0] binwidth = p[1] dbinwidth = p[2] return & end pro even_light_prorate, nx, ny, center, p_box, b_box, weight, $ prorate=prorate, zphi=zphi, exclude_p_box=exclude_p_box @compile_opt.pro ; On error, return to caller InitVar, prorate, /key ix = b_box[2]-b_box[0]+1 iy = b_box[3]-b_box[1]+1 xyvec = gridgen([ix,iy],origin=center-b_box[0:1]) prvec = cv_coord(from_rect=xyvec, /to_polar) prvec = reform(prvec, 2, ix, iy, /overwrite) zphi = reform(prvec[0,*,*]) in_center = Inside_Wedge(p_box, prvec, exclude_p_box=exclude_p_box) prcor = cv_coord(from_rect=gridgen([ix,iy]+1,origin=center-b_box[0:1]+0.5), /to_polar) prcor = reform(prcor, 2, ix+1, iy+1, /overwrite) in_corner = Inside_Wedge(p_box, prcor, exclude_p_box=exclude_p_box) inside = in_center eq in_corner[0:ix-1,0:iy-1] and $ in_center eq in_corner[1:ix ,0:iy-1] and $ in_center eq in_corner[0:ix-1,1:iy ] and $ in_center eq in_corner[1:ix ,1:iy ] across = where(inside eq 0B, n_across) inside = where(inside and in_center, n_inside) n_across = prorate*n_across case n_across of 0 : weight = intarr(ix,iy) else: weight = dblarr(ix,iy) endcase if n_inside gt 0 then weight[inside] = 1 ; The prorating is done by dividing each pixel on the boundary in a group of 10x10 subpixels if n_across gt 0 then begin mx = 10 my = 10 ; x,y coordinates of subpixels between -0.5,+0.5 r = [mx,my] r_small = gridgen( r, range=(0.5*[-1,1])#((r-1.)/r) ) ; x,y coordinate relative to origin of pixels to be prorated r = xyvec[*,across] ; Set up an mx x my x n_across array for all subpixels in all pixels to be prorated r = SuperArray(r, mx*my, after=1)+SuperArray(r_small,n_across,/trail) r = reform(r,2,mx*my*n_across, /overwrite) ; Convert the subpixel centers to polar coordinates ; Determine which subpixels lie inside the wedge r = Inside_Wedge(p_box, cv_coord(from_rect=r, /to_polar), exclude_p_box=exclude_p_box) r = reform(r, mx*my, n_across, /overwrite) r = total(r,1)/(mx*my) ; Fraction of subpixels lying inside the wedge weight[across] = r endif return & end function even_light_shift, hhvv, i @compile_opt.pro ; On error, return to caller ; Input: ; hhvv string identifying image sequence ; i image number in sequence ; ; Output: ; i_shift[2] # shifts relative to first image y_run = strmid(hhvv, strlen(hhvv)-1,1) eq 'y' i_unit = [1-y_run, y_run] if i le 5 then $ i_shift = i_unit $ else if i le 45 then $ i_shift = ( i-5)*i_unit $ else if i le 47 then $ i_shift = (45-5)*i_unit $ else if i le 50 then $ i_shift = i_unit return, i_shift & end function even_light_bin_phi, optaxis0, j, doptaxis, binwidth, dbinwidth, bins_step ; j=0 is center spot @compile_opt.pro ; On error, return to caller dphi_j = j*(binwidth+j*dbinwidth) ; Include linear correction to binwidth across image phi_j = optaxis0[0]+dphi_j ; Azimuth of spot j in 1st image of the sequence phi_j = phi_j+doptaxis[0]*cos(dphi_j) ; Azimuth of spot j in shifted image if n_elements(bins_step) ne 0 then $ phi_j = 0.5*( [even_light_bin_phi(optaxis0, j-bins_step, doptaxis, binwidth, dbinwidth), $ even_light_bin_phi(optaxis0, j+bins_step, doptaxis, binwidth, dbinwidth)] + phi_j ) return, phi_j & end function even_light_bin_rad, optaxis0, j, doptaxis ; j=0 is center spot @compile_opt.pro ; On error, return to caller rad_j = optaxis0[1] ; Radius of spot j in 1st image of the sequence rad_j = rad_j+doptaxis[1] ; Radius of spot j in shifted image return, rad_j & end function even_light_centroid, nx, ny, nimg, nimg_min, nimg_max, images, hhvv, $ center0, dcenter, optaxis0, doptaxis_unit, binwidth, dbinwidth, $ sym_spot_nr, visible_spot, bins_step, move_grid, dr_fov, $ gain, P1, P2, by_center=by_center, spot_center=spot_center @compile_opt.pro ; On error, return to caller InitVar, by_center, /key box_edge = 4d0*[(8*0.0001d0),1] ; ~ 4 pixels at radius 1200 dr_max = 0.01 nn_max = 10 nvis = n_elements(visible_spot) nspot = n_elements(sym_spot_nr) if nimg gt 0 then begin spot_centroid = replicate(BadValue(0d0),2,nspot,nimg) spot_center = replicate(BadValue(0d0),2,nimg) nimg_beg = nimg_min nimg_end = nimg_max endif else begin spot_centroid = replicate(BadValue(0d0),2,nspot) spot_center = replicate(BadValue(0d0),2) nimg_beg = -nimg nimg_end = -nimg endelse for i=nimg_beg,nimg_end do begin ; Loop over all images print, i i_img = i-nimg_min zfov = even_light_background(images[*,*,i_img], gain[i_img], P1[i_img], P2[i_img]) zmin = min(zfov, max=zmax) view, in= (alog10(zfov-zmin+1) > 0) < alog10(zmax-zmin+1), /stretch, title=hhvv i_units = even_light_shift(hhvv, i) center = center0+dcenter*i_units doptaxis = i_units*doptaxis_unit ; We calculate centroids for all visible bins, including 'disqualified' bins ; (i.e. bins not to be used for fitting circles or normalization). for j=0,nvis-1 do begin ; Loop over all bins spot = visible_spot[j] jj = sym_spot_nr[spot] phi = even_light_bin_phi(optaxis0, jj, move_grid*doptaxis, binwidth, dbinwidth, bins_step) rad = even_light_bin_rad(optaxis0, jj, move_grid*doptaxis) rad = rad+dr_fov even_light_view, center, phi, rad q_box = transpose([[phi], [rad]]) ; outer box p_box = q_box+box_edge#[1,-1] ; inner box centroid = flat_centerofmass(p_box, /polar) ; centroid for inner box phi = even_light_bin_phi(optaxis0, 0, move_grid*doptaxis, binwidth, dbinwidth, bins_step) q_box = transpose([[phi], [rad]]) ; outer box p_box = q_box+box_edge#[1,-1] ; inner box box = flat_centerofmass(p_box, /polar, /shift_origin) nn = 0 repeat begin prev_centroid = centroid ; Centroid in polar coordinates p_box = flat_centerofmass(box, prev_centroid, /polar) ; outer box q_box = p_box+box_edge#[-1,1] ; inner box b_box = wedge_bounding_box(q_box, center=center) b_box = (b_box > 0) < ([nx,ny]-1)#[1,1] even_light_prorate, nx, ny, center, q_box, b_box, dr, $ prorate=prorate, exclude_p_box=p_box, zphi=zphi zval = zfov[b_box[0]:b_box[2], b_box[1]:b_box[3]] if by_center then zphi = prev_centroid[0] even_light_geometric_fix, zval, zphi, optaxis0, doptaxis dr = dr/total(dr) faint = total(dr*zval) even_light_prorate, nx, ny, center, p_box, b_box, dr, prorate=prorate last_centroid = b_box[0:1]+CenterOfMass(dr*(zval-faint)) last_centroid = last_centroid-center centroid = last_centroid dr = centroid-cv_coord(from_polar=prev_centroid, /to_rect) centroid = cv_coord(from_rect=centroid, /to_polar) dr = sqrt( total(dr*dr) ) nn = nn+1 endrep until dr lt dr_max or nn ge nn_max case nimg gt 0 of 0: spot_centroid[*,spot ] = last_centroid 1: spot_centroid[*,spot,i_img] = last_centroid endcase endfor case nimg gt 0 of 0: spot_center = center 1: spot_center [*,i_img] = center endcase endfor return, spot_centroid & end pro even_light, $ ny = ny, $ ; 300 <= ny <= 600 (default: 600) check = check, $ ; check single image fix_grid=fix_grid, $ ; don't let the grid move prorate = prorate, $ ; prorate pixels for centroid calculations hhvv = hhvv, $ ; name of sequence sloped = sloped, $ ; allow pedestal P1 to depend linearly on frame number bmp = bmp, $ ; write bmp files as output paper = paper, $ ; write to printer (hardcopy output) raw = raw, $ ; subtract raw background for gain calculations indir = indir, $ ; input directory is indir+hhvv outdir = outdir, $ ; output directory for all .pph and .bmp files by_center=by_center,$ ; make geometric corrections based on centroid of spot flatfield=flatfield @compile_opt.pro ; On error, return to caller InitVar, sloped , /key InitVar, raw , /key InitVar, by_center, /key InitVar, flatfield, /key InitVar, bmp , /key InitVar, paper, /key if n_elements(hhvv) eq 0 then hhvv = 'flt_1x' ; Output files (*.pph and *.bmp) to $TUB unless otherwise specified if n_elements(indir ) eq 0 then indir = getenv('TUB') if n_elements(outdir) eq 0 then outdir = getenv('TUB') dirnic = filepath(root=indir, sub=hhvv) start_time = TimeSystem() nx = 1280L if n_elements(ny) eq 0 then ny = 600L ny = ny > 300L InitVar, fix_grid, /key move_grid = 1-fix_grid InitVar, prorate, /key fmt = '(I2.2)' if hhvv eq 'dot20' then fmt = '(I3.3)' type = '.nic' outer_spot = 20L ; 20 spots left and right of center spot nspot = 2*outer_spot+1 ; # number of spots ; Spots are numbered left to right across the image. ; There are 41 spots (i=0,40); we introduce a symmetric numbering ; which assigns 0 to the center spot. sym_spot_nr = -outer_spot+lindgen(nspot) ; j=0,40 -> sym_spot_nr[j]=-20,20; center spot has sym_spot_nr[20] = 0 bins_axis = float( where( sym_spot_nr eq 0 ) ) ; Index for optical axis ; 'visible_spot' lists the visible spots among spots 0..40. ; We assume that only the 20 odd-numbered spots are used. bins_step = 2 visible_spot = sym_spot_nr/bins_step*bins_step ne sym_spot_nr ; 'good_spot' identifies the spots to be used in fitting circles to ; the centroids. This will be a subset of 'visible_spot'. ; Currently only spot 15 (which has bad column 502 in it) for the UCSD ; engineering CCD is excluded. ; We don't want to exclude this spot for the Birmingham runs. good_spot = visible_spot if strpos(strlowcase(hhvv), 'em_' ) eq -1 and $ strpos(strlowcase(hhvv), 'bflt' ) eq -1 and $ strpos(strlowcase(hhvv), 'z3flt') eq -1 then $ good_spot[15] = 0 $ else $ message, /info, ' Birmingham run : '+hhvv ; Indices for visible and good spots (can be used as index to sym_spot_nr) visible_spot= where(visible_spot) good_spot = where(good_spot ) nimg_min = 1 nimg_max = 50 nimg = nimg_max-nimg_min+1 img_nr = nimg_min+indgen(nimg) center0 = [632d0, 1296d0] dcenter = [ 0d0, 0d0] optaxis0 = [-89.5d0*(!dpi/180), 1186] doptaxis_unit = [0d0, 0d0] binwidth = 1.5*(!dpi/180) dbinwidth = 0.0d0 moving_img = [ 5,45]-nimg_min dr_fov = [-50 +15, +50 -15] if n_elements(check) ne 0 then begin i = (check[0] < nimg_max) > nimg_min i_img = i-nimg_min files = file_search(filepath(root=dirnic,hhvv+string(format=fmt,i)+type),count=nfiles) if img_read(files[0],img) then begin img = float(img[*,0:ny-1]) spot_centroid = even_light_centroid(nx, ny, 0, 0, 0, img, hhvv, $ center0, dcenter, optaxis0, doptaxis_unit, binwidth, dbinwidth, $ sym_spot_nr, visible_spot, bins_step, move_grid, dr_fov, $ 1, 0, 0, by_center=by_center) even_light_firstguess, spot_centroid, sym_spot_nr, good_spot, bins_axis, center0, optaxis0, binwidth, dbinwidth zmin = min(img, max=zmax) view, in=(alog10(img-zmin+1) > 0) < alog10(zmax-zmin+1), /stretch, title=hhvv i_units = even_light_shift(hhvv, i) doptaxis = i_units*doptaxis_unit center = center0+dcenter*i_units phi = even_light_bin_phi(optaxis0, 0, move_grid*doptaxis, binwidth, dbinwidth) ; Azimuth center spot rad = even_light_bin_rad(optaxis0, 0, move_grid*doptaxis) ; Radius center spot rad = rad+dr_fov plots, center#[1,1]+[cos(phi),sin(phi)]#[0,rad[0]], color=255 a = visible_spot[0] phi1 = even_light_bin_phi(optaxis0, sym_spot_nr[a] , move_grid*doptaxis, binwidth, dbinwidth, bins_step) a = visible_spot[n_elements(visible_spot)-1] phi2 = even_light_bin_phi(optaxis0, sym_spot_nr[a] , move_grid*doptaxis, binwidth, dbinwidth, bins_step) n = 100 a = gridgen(n, range=[phi1[0],phi2[1]]) a = transpose([[cos(a)],[sin(a)]]) cntr = center#replicate(1,n) plots, cntr+rad[0]*a, color=255 plots, cntr+rad[1]*a, color=255 for j=0,n_elements(good_spot)-1 do begin phi = even_light_bin_phi(optaxis0, sym_spot_nr[good_spot[j]] , move_grid*doptaxis, binwidth, dbinwidth, bins_step) even_light_view, center, phi, rad endfor endif return endif ; Read all the images into a floating point array images = uintarr(nx, ny, nimg, /nozero) blnk_raw = fltarr(nimg, /nozero) dark_raw = fltarr(nimg, /nozero) back_raw = fltarr(nimg, /nozero) for i=nimg_min,nimg_max do begin i_img = i-nimg_min files = file_search(filepath(root=dirnic,hhvv+string(format=fmt,i)+type),count=nfiles) stat = img_read(files[0],img) images[*,*,i_img] = img[*,0:ny-1] ; Fix for anomalous images in EM_ runs made at Birmingham, following Chris' suggestion ; that an extra pixel gets inserted at the start of each image. Shift moves the extra ; pixel to the end. if ( hhvv eq 'em_1x' and (i eq 17 or i eq 23 or i eq 35 ) ) or $ ( hhvv eq 'em_1y' and (i eq 17 or i eq 23 or i eq 47 ) ) or $ ( hhvv eq 'em_2x' and (i eq 2 or i eq 45 ) ) or $ ( hhvv eq 'em_2y' and (i eq 10 or i eq 12 or i eq 18 or i eq 32) ) or $ ( hhvv eq 'em_3x' and (i eq 5 or i eq 19 ) ) or $ ( hhvv eq 'em_3y' and (i eq 29 or i eq 36 ) ) or $ ( hhvv eq 'bflt1x' and (i eq 11 ) ) or $ ( hhvv eq 'bflt2x' and (i eq 22 or i eq 41 ) ) or $ ( hhvv eq 'bflt2y' and (i eq 3 or i eq 11 or i eq 35 ) ) or $ ( hhvv eq 'bflt3x' and (i eq 42 ) ) or $ ( hhvv eq 'bflt3y' and (i eq 2 or i eq 49 ) ) then $ images[*,*,i_img] = shift(images[*,*,i_img], -1) if hhvv eq 'ucsdxe' and (37 le i and i le 44 ) then $ images[*,*,i_img] = shift(images[*,*,i_img], 0, 24) if max(images[*,*,i-nimg_min],pos) eq uint(-1) then begin print, ArrayLocation(pos,size=size(images[*,*,i-nimg_min])) message, /info, 'saturated' endif blnk_raw[i_img] = mean(images[1270:1279,20:(570 < (ny-1)),i_img]) dark_raw[i_img] = mean(images[ 17: 20,20:(570 < (ny-1)),i_img]) back_raw[i_img] = mean(images[ 50: 130,30:(100 < (ny-1)),i_img]) endfor status = bin_read( filepath(root=indir, 'flatfield.pph'), flat) if flatfield and status then begin flat = 1+flat[*,0:ny-1] images = float(images) for i=nimg_min,nimg_max do begin i_img = i-nimg_min images[*,*,i_img] = images[*,*,i_img]*flat endfor endif else $ flatfield = 0 boxes = [ [ 1, 21], [ 4, 24], $ ; 0: blank pix at bottom left [ 13, 21], [ 16, 24], $ ; 1: blank pix at bottom left [ 17, 21], [ 20, 24], $ ; 2: dark pix at bottom left [ 66, 21], [ 69, 24], $ ; 3: background at bottom left [ 644, 21], [ 647, 24], $ ; 4: background at bottom center [1216, 21], [1219, 24], $ ; 5: background at bottom right [1265, 21], [1268, 24], $ ; 6: dark pix at bottom right [1271, 21], [1274, 24], $ ; 7: blank pix at bottom right [ 1,551], [ 4,554], $ ; 8: blank pix at top left [ 13,551], [ 16,554], $ ; 9: blank pix at top left [ 17,551], [ 20,554], $ ; 10: dark pix at top left [ 66,551], [ 69,554], $ ; 11: background at top left [ 644,551], [ 647,554], $ ; 12: background at top center [1216,551], [1219,554], $ ; 13: background at top right [1265,551], [1268,554], $ ; 14: dark pix at top right [1271,551], [1274,554], $ ; 15: blank pix at top right [ 644,197], [ 647,200], $ ; 16: center, above fov [ 644,109], [ 647,112] ] ; 17: center, inside fov nbox = n_elements(boxes)/4 ; # boxes boxes = reform(boxes,2,2,nbox,/overwrite) box_stat = replicate(BadValue(0.0), 16, nbox, nimg) for i=0,nbox-1 do $ if boxes[1,0,i] le ny-1 then $ box_stat[*,i,*] = images[boxes[0,0,i]:boxes[0,1,i], boxes[1,0,i]:(boxes[1,1,i] < (ny-1)),*] box_stat = reform(box_stat, 16*nbox, nimg, /overwrite) ; First calculate the centroids of the image in the middle in the sequence to get ; a first correction on center and radius spot_centroid = even_light_centroid(nx, ny, -(nimg_min+nimg_max)/2, nimg_min, nimg_max, images, hhvv, $ center0, dcenter, optaxis0, doptaxis_unit, binwidth, dbinwidth, $ sym_spot_nr, visible_spot, bins_step, move_grid, dr_fov, $ replicate(1,nimg), back_raw, replicate(0,nimg), $ by_center=by_center, spot_center=spot_center) ; At this point center0=spot_center even_light_firstguess, spot_centroid, sym_spot_nr, good_spot, bins_axis, center0, optaxis0, binwidth, dbinwidth ; Note that the first image is not used so image 0 is file with number 01 ; Calculate centroids to refine other geometrical corrections for the grid ; Note also that this is done on the raw images. We don't have the gain yet ; so we set it to 1. We do subtract the background value 'back_raw', but ; this does not affect the centroid since the centroid calculation subtracts ; a local background based on a 4 pixel rim along the edge of each spot. ; At this point doptaxis_unit (indicating the motion across the sequence) is still set to zero), ; so the bins used to calculate the centroids do not follow the motion. spot_centroid = even_light_centroid(nx, ny, nimg, nimg_min, nimg_max, images, hhvv, $ center0, dcenter, optaxis0, doptaxis_unit, binwidth, dbinwidth, $ sym_spot_nr, visible_spot, bins_step, move_grid, dr_fov, $ replicate(1,nimg), back_raw, replicate(0,nimg), $ by_center=by_center, spot_center=spot_center) ; spot_centroid refers to the same centroid for all images: center0 (=spot_center) even_light_corrections, hhvv, spot_centroid, sym_spot_nr, good_spot, bins_axis, $ bins_step, img_nr, center0, dcenter, optaxis0, doptaxis_unit, binwidth, dbinwidth, $ y_run, bmp=bmp, paper=paper, outdir=outdir file = hhvv+'_centroids_rough.pph' bin_write, filepath(root=outdir, file), $ ; Each record has 41 spots [ transpose(img_nr), $ ; image number spot_center, $ reform(spot_centroid, 2*nspot, nimg)] ; (x,y)-pairs for each spot file = hhvv+'_geometry.pph' bin_write, filepath(root=outdir, file), $ [center0, dcenter, optaxis0, doptaxis_unit, binwidth, dbinwidth, y_run] ; Calculate the total intensity and total number of pixels in the SMEI fov. ; This includes only spots listed in 'good_spot' ; We use the intensity from the raw images with the background 'back_raw' subtracted. ; The remaining intensity are corrected for geometrical affects (cosine- and ; 1/r^2 effect) adu_in_fov = fltarr(nimg, /nozero) pix_in_fov = fltarr(nimg, /nozero) zfov = fltarr(nx, ny) wfov = fltarr(nx, ny) ngood = n_elements(good_spot) for i=nimg_min,nimg_max do begin print, i case prorate of 0B: i_fix = (nimg_min+nimg_max)/2 1B: i_fix = i endcase i_units = even_light_shift(hhvv, i_fix) doptaxis = i_units*doptaxis_unit center = center0+dcenter*i_units i_img = i-nimg_min wfov[*] = 0. zfov[*] = 0. for j=0,ngood-1 do begin jj = sym_spot_nr[good_spot[j]] phi = even_light_bin_phi(optaxis0, jj, move_grid*doptaxis, binwidth, dbinwidth, bins_step) rad = even_light_bin_rad(optaxis0, jj, move_grid*doptaxis) r = rad+dr_fov p_box = transpose([[phi], [r]]) b_box = wedge_bounding_box(p_box, center=center) b_box = (b_box > 0) < ([nx,ny]-1)#[1,1] even_light_prorate, nx, ny, center, p_box, b_box, fraction, prorate=prorate, zphi=zphi zval = images[b_box[0]:b_box[2], b_box[1]:b_box[3],i_img]-back_raw[i_img] if by_center then $ zphi = even_light_bin_phi(optaxis0, jj, doptaxis, binwidth, dbinwidth) even_light_geometric_fix, zval, zphi, optaxis0, doptaxis wfov[b_box[0]:b_box[2], b_box[1]:b_box[3]] = wfov[b_box[0]:b_box[2], b_box[1]:b_box[3]] + fraction zfov[b_box[0]:b_box[2], b_box[1]:b_box[3]] = zfov[b_box[0]:b_box[2], b_box[1]:b_box[3]] + fraction*zval endfor view, in=bytscl((wfov ne 0)*alog10(zfov-min(zfov)+1))+(wfov eq 0)*255B, /stretch, title=hhvv pix_in_fov[i_img] = total(wfov) adu_in_fov[i_img] = total(zfov)/total(wfov) endfor destroyvar, zval, fraction, zfov, wfov, zphi, spot_centroid even_light_pedestal, adu_in_fov, back_raw, sloped=sloped, $ gain, P1pP2, P1, P2, P1pP2sigma, P1sigma, P2sigma ; Final normalized background (should be zero on average) back_pix = even_light_background(back_raw, gain, P1, P2) ; Write the normalized images to disk for i=nimg_min,nimg_max do begin i_img = i-nimg_min zfov = even_light_background(images[*,*,i_img], gain[i_img],$ (1-raw)*P1[i_img]+raw*back_raw[i_img], $ (1-raw)*P2[i_img]) zfov = uint( (zfov > 0) < uint(-1) ) trailer = 'Gain='+strcompress(gain[i_img], /rem) + $ ' Pedestal1='+strcompress(P1[i_img], /rem) + $ ' Pedestal2='+strcompress(P2[i_img], /rem) ;bin_write, filepath(root=outdir, hhvv+string(format=fmt,i)+'.pph'), zfov, trailer=trailer endfor ; Calculate intensity and number of pixels for every spot. ; We now use the normalized intensity (corrected for gain). nvis = n_elements(visible_spot) pix_per_bin = replicate(BadValue(1.0),nspot,nimg) adu_per_bin = replicate(BadValue(1.0),nspot,nimg) for i=nimg_min,nimg_max do begin ; Loop over all images print, i i_img = i-nimg_min zfov = even_light_background(images[*,*,i_img], gain[i_img], $ (1-raw)*P1[i_img]+raw*back_raw[i_img], $ (1-raw)*P2[i_img]) zmin = min(zfov, max=zmax) view, in=(alog10(zfov-zmin+1) > 0) < alog10(zmax-zmin+1), /stretch, title=hhvv i_units = even_light_shift(hhvv, i) doptaxis = i_units*doptaxis_unit center = center0+dcenter*i_units for j=0,nvis-1 do begin ; Loop over all visible spots spot = visible_spot[j] jj = sym_spot_nr[spot] phi = even_light_bin_phi(optaxis0, jj, move_grid*doptaxis, binwidth, dbinwidth, bins_step) rad = even_light_bin_rad(optaxis0, jj, move_grid*doptaxis) rad = rad+dr_fov even_light_view, center, phi, rad q_box = transpose([[phi], [rad]]) b_box = wedge_bounding_box(q_box, center=center) b_box = (b_box > 0) < ([nx,ny]-1)#[1,1] even_light_prorate, nx, ny, center, q_box, b_box, fraction, prorate=prorate, zphi=zphi zval = zfov[b_box[0]:b_box[2], b_box[1]:b_box[3]] if by_center then $ zphi = even_light_bin_phi(optaxis0, jj, doptaxis, binwidth, dbinwidth) even_light_geometric_fix, zval, zphi, optaxis0, doptaxis pix_per_bin[spot,i_img] = total(fraction) adu_per_bin[spot,i_img] = total(fraction*zval) endfor endfor destroyvar, zval, zphi, fraction file = hhvv+'_photometry.pph' bin_write, filepath(root=outdir, file), $ ; Each record has 90 numbers: [ transpose( [[img_nr ] , $ ; image number [adu_in_fov ] , $ [pix_in_fov ] , $ [gain ] , $ ; gain [P1pP2 ] , $ [P1 ] , $ ; Background constant (independent of gain; P1) [P2 ] , $ ; Background constant ( dependent of gain; P2) [P1pP2sigma ] , $ ; Error in P1 (one value same for all images) [P1sigma ] , $ ; Error in P1 (one value same for all images) [P2sigma ] , $ ; Error in P2 (one value same for all images) [blnk_raw ] , $ ; blank pixels [dark_raw ] , $ ; dark pixels [back_raw ] , $ ; background in original image [back_pix ]]), $ ; residual background after normalization adu_per_bin, $ ; Adus (41 elements per image) pix_per_bin, $ ; # Pixels (41 elements per image) box_stat ] ; 17 groups of 16 pixels spot_centroid = even_light_centroid(nx, ny, nimg, nimg_min, nimg_max, images, hhvv, $ center0, dcenter, optaxis0, doptaxis_unit, binwidth, dbinwidth, $ sym_spot_nr, visible_spot, bins_step, move_grid, dr_fov, $ gain, (1-raw)*P1+raw*back_raw, (1-raw)*P2, $ by_center=by_center, spot_center=spot_center) destroyvar, images file = hhvv+'_centroids.pph' bin_write, filepath(root=outdir, file), $ ; Each record has 41 [ transpose(img_nr), $ ; image number spot_center, $ reform(spot_centroid, 2*nspot, nimg)] ; (x,y)-pairs for each spot twin, /del, /all even_light_photometry, hhvv, outdir=outdir, $ bmp=bmp, paper=paper, /spots, $ sym_spot_nr, visible_spot, img_nr, adu_per_bin, $ gain, blnk_raw, dark_raw, back_raw, back_pix, stat=stat whatis, stat stat_median = stat[1] even_light_registration, hhvv, outdir=outdir, $ bmp=bmp, paper=paper, $ sym_spot_nr, visible_spot, img_nr, spot_centroid, moving_img, optaxis0 end_time = TimeSystem() print, 'minutes: ', TimeOp(/subtract, end_time, start_time, TimeUnit(/min)) str = [ 'sequence = '+hhvv, $ ; name of sequence 'ny = '+strcompress(ny , /rem), $ 'move_grid= '+strcompress(move_grid , /rem), $ 'prorate = '+strcompress(prorate , /rem), $ 'sloped = '+strcompress(sloped , /rem), $ 'raw = '+strcompress(raw , /rem), $ 'by_center= '+strcompress(by_center , /rem), $ 'flatfield= '+strcompress(flatfield , /rem), $ 'median = '+strcompress(stat_median,/rem) ] file = hhvv+'_keywords.txt' openw, /get_lun, iu, filepath(root=outdir, file) for i=0,n_elements(str)-1 do printf, iu, str[i] free_lun, iu return & end