; file: gen-cubes-rr.pro = generic annotated classic version ; init: Sep 16 2010 SST ; last: Dec 15 2010 ; hist: Sep 2010 SST original from dotdb 2007-04-12 => RR cube-making trial ; Dec 2010 Deil for 2005-07-13-spot = RR EBduds paper ; note: classic version for DOT data up to spring 2007; Hitachi KPF100 camera ; ;;; comment by Pit Suetterlin (PS) in original version 2007-04-12 ; ; (also ;RR) comment by Rob Rutten ; ; ?? = issue ; ; @ = todo ; ; @ adapt? = parameter likely to be adapted for other data sets ; ; # comment ; does: @?? ; todo: @ add new Halpha cameras ; @ add Barium 4554 cameras ; the two functions below serve to write and read temporary assoc cubes FUNCTION Readframe, unit, n ; n = frame number COMMON assocsize, sx0, sy0, imsiz ; xsize, ysize px; image size bytes img = intarr(sx0, sy0) point_lun, unit, 1l*imsiz*n ; 1l long integer; pointer in bytes readu, unit, img ; read image nr n return, img END PRO Writeframe, img, unit, n ; idem, write image nr n COMMON assocsize, sx0, sy0, imsiz tmp = intarr(sx0, sy0) s = size(img, /dim) < [sx0, sy0] ; image may be smaller than standard tmp(0:s(0)-1, 0:s(1)-1) = img(0:s(0)-1, 0:s(1)-1) ; rest gets black point_lun, unit, 1l*imsiz*n writeu, unit, tmp END ; =============== main part starts here ============================= COMMON assocsize, sx0, sy0, imsiz ; parameter values and choices plotout=0 ; @ adapt? get screen output or not? nam = 'AR0789-' ; @ adapt? start of filename including first - n_ha = 3 ; @ adapt? nr wavs Ha (not counting rc) outcubedirname = 'Cubesfull' ; adapt? dir for output cubes; no hyphen! ca_version = 2 ; @ adapt? # ca sequence may have tilt alternation ;;; 0 -> only core sequence, correlate sequence separate and align average ;;; 1 -> flip core/wing sequence, align wing to gb and core to wing ;;; 2 -> align every image to G Band, whether ca sequence flipped or not ;;; cube wavelength names ; @ adapt? ; strict order: gb, bc, rc, ha, ca ; ha KvdL: follow the order in fits file header, use script fits ; e.g.: fits Filament-ha-scan-20070412-121948.fits ; D_LAMBDA= [ 0.00, -0.70, 0.70, 0.35, -0.35] lam = ['gb-0000', 'bc-0000', 'rc-0000',$ 'ha-0500', 'ha-0000', 'ha+0500',$ 'ca-2350'] ;;; cadence specification in filenames dt = '-dt30' ; @ adapt? cadence value in file name (not used in this code) ;;; polyfit subtract for shc? ; subtract near limb for darkening pf = 0 ; @ adapt? ;;; camera orientation angles ; change with time, especially when camera mounts are re-adjusted ; these are for 2005-07-13-spot par_bc = [-0.3911, 1.00155] par_rc = [-0.9546, 1.00428] par_ha = par_rc par_ca = [-0.3234, 1.0023] ;;; not measurable - use old value ; set array size = largest image size of all fits files ; run script "getmaxfitssize" in topdir (above FITS) to find largest ; values of NAXIS1 and NAXIS2 in all fits file headers under FITS ; insert those below for minimum memory use ; otherwise simply set sx0=1280 & sy0=1024 = camera chip size = max ; actual values will be smaller since speckle includes tip-tilt correction ; script command (bash) is: ;(for i in FITS/*/*fits; do fits $i |grep NAXIS1|cut -c 27-30;done)|sort|tail -5 sx0 = 1238 ; @ adapt? RR for data in 2005-07-13-spot sy0 = 991 ; @ adapt? RR for data in 2005-07-13-spot imsiz = 2l*sx0*sy0 ; size in bytes /image; image px = 2 bytes ;;; apply additional shift between RC and HA? ; @ adapt? ; following is additional shift for off disk-center when ; foreshortening on high-rising Ha structures makes co-location differ ; (but still difference between centerward and limbward field edges) sh_ha = [0, 2] ;RR [0,2] of [0,0]? Generic version has [0,0] xr0 = abs(sh_ha(0) > 0) ; r = rc yr0 = abs(sh_ha(1) > 0) xh0 = abs(sh_ha(0) < 0) ; h = Halpha yh0 = abs(sh_ha(1) < 0) ; following makes list of the fits files for each camera ; (but IDL ? shows findfile is an "obsolete routine") fg = findfile('../FITS/GBand/'+nam+'*') fb = findfile('../FITS/BlueCont/'+nam+'*') fr = findfile('../FITS/Halpha/'+nam+'*') ft = findfile('../FITS/RedCont/'+nam+'*') ; not used since rc also in ha fc = findfile('../FITS//CaIIH/'+nam+'*') ;;; ca and ha may be switched -> names change; resort on date/time fc = fc(sort(strmid(fc, strlen(fc(0))-20))) fr = fr(sort(strmid(fr, strlen(fr(0))-20))) ;;; jumps? then split the whole sequence into segments (a, b, c..) ; @ adapt ? select a different segment than first-last first = 4 ;RR Pit had 4 so the first ones are probably offset last = n_elements(fg)-1 ;RR ##### 10-frame trial segment starting at very good seeing ;; first=124 ;; last=first+99 ; ================== no parameter choices from here on ; cut out the specified segment fg = fg(first:last) fb = fb(first:last) fr = fr(first:last) ft = ft(first:last) ;RR not used since rc also in fr = Halpha fc = fc(first:last) ; num = number of selected images num = n_elements(fg) ;;; get the date tmp = strsplit(fg(0), '-', /extract) dat = '-'+tmp(3) ;;; check timeline s = strlen(fg(0)) tim = strmid(fg, s-11, 2)*3600l+strmid(fg, s-9, 2)*60l+strmid(fg, s-7, 2) c = poly_fit(indgen(n_elements(tim)), tim, 1) ; print mean dt in seconds; c = t_0 + c(1)t print,' mean sample interval (sec) =',c(1) ; open the temporary image write-out file for write/read with assoc ; .asc extension = abbreviation for assoc ; n_ha+3 = n_ha ones + bc, gb, rc # i below: Pit has no virtual finger spawn, 'rm -f *.asc' FOR i=1, n_ha+3 DO openw, i, 'cube-'+strtrim(i, 2)+'.asc' ; intensity scale factors (for typical exposures?) to fill integer range? sall = intarr(2, num) ; size array for 1 sequence (nx, it), (ny,it) nsum = 5 ; number of images to make refrence for tmeporal alignment sh = intarr(2, num) scl = 5000 scl_ha = 25 scl_ca = 3000 imax = 32767 if (plotout) then window, 1, xs=sx0, ys=sy0 ;============================================= gb, bc, rc, ha ; big loop over time to treat images in sequences bc, gb, rc, ha FOR i=0, num-1 DO BEGIN imarr = intarr(sx0, sy0, n_ha+3) ;;; read gb = G band ; rdfits is in dotlib, reads Pit triplet ; = reconstructed image, best image of burst, mean of the whole burst rdfits, p, fg(i) ; select the reconstructed image only (reform takes off 3rd dimension) p1 = reform(p(*, *, 0)) sgb = size(p1, /dim) imarr(0, 0, 0) = fix(scl*p1 < imax) ;;; read bc = blue cont rdfits, p, fb(i) ; rot = IDL function p1 = rot(reform(p(*, *, 0)), par_bc(0), par_bc(1), cubic=-.5) sbc = size(p1, /dim) imarr(0, 0, 1) = fix(scl*p1 < imax) ;;; read Halpha (rc + Halpha wavs) ; includes rc as (*,*,0) rdfits, p, fr(i) p1 = rot(reform(p(*, *, 0)), par_rc(0), par_rc(1), cubic=-.5) src = size(p1, /dim) ;;; handle remaining ha-shift # foreshortening problem ; p(*,*,j) for j>0 = Halpha wavs src = src-abs(sh_ha) p1 = fix(scl_ha*p1(xr0:xr0+src(0)-1, yr0:yr0+src(1)-1) < imax) imarr(0, 0, 2) = p1 FOR j=1, n_ha DO BEGIN p1 = rot(reform(p(xh0:xh0+src(0)-1, yh0:yh0+src(1)-1, j)), $ par_rc(0), par_rc(1), cubic=-.5) imarr(0, 0, 2+j) = fix(scl_ha*p1 < imax) ENDFOR ;;; now shift all continua to match and cut the size to be the same ; gb, bc, rt are all full speckle-burst triplets ; triplet = reconstructed + sharpest + mean after burst alignment ; but rt = red cont triplets are not used here but rc from ha ; shc.pro in dotlib delivers shift between two images ssh = intarr(2, 3) ; get shifts of bc (j=1) and rc (j=2) with respect to gb (j=0) FOR j=1, 2 DO ssh(*, j)=$ shc(imarr(*, *, 0), imarr(*, *, j), /N2, /filt, POLY=pf) ; make all x-shifts and y-shifts negative; gb (j=0) = -max ssh(0, *) = ssh(0, *)-max(ssh(0, *)) ssh(1, *) = ssh(1, *)-max(ssh(1, *)) ; shift gb and bc (gb all px over the same -max) FOR j=0, 1 DO imarr(*, *, j) = shift(imarr(*, *, j), ssh(0, j), ssh(1, j)) ; shift rc and all n_Ha frames over rc shifts (j=2) FOR j=2, n_ha+2 DO imarr(*, *, j) = shift(imarr(*, *, j), ssh(0, 2), ssh(1, 2)) ; cut to overlap size sall(*, i) = (sgb+ssh(*, 0)) < (sbc+ssh(*, 1)) < (src+ssh(*, 2)) imarr = imarr(0:sall(0, i)-1, 0:sall(1, i)-1, *) if (plotout) then begin tvscl, imarr(*, *, 0) & wait, .5 tvscl, imarr(*, *, 1) & wait, .5 tvscl, imarr(*, *, 2) & empty endif ; temporary writeout all images FOR j=0, n_ha+2 DO writeframe, imarr(*, *, j), j+1, i ; in order to align the sequences in time, shifts are determined ; with respect to the last nsum sums of the three continua. They ; will be applied in common to all wavelengths = 3 continua + ; Halpha KvdL scan wavelengths ; get the sum of the gb, bc and rc images for this time moment (i) p1=$ imarr(*, *, 0)/avg(imarr(*, *, 0)) + $ imarr(*, *, 1)/avg(imarr(*, *, 1)) + $ imarr(*, *, 2)/avg(imarr(*, *, 2)) IF i EQ 0 THEN BEGIN pp = fltarr(nsum, sall(0, 0), sall(1, 0)) pp(0, *, *) = p1 ENDIF ELSE BEGIN ; set common area spp = sall(*, i) FOR k=(i-num) > 0, i-1 DO spp = spp < sall(*, k) ; set reference image = sum over nsum ref = total(pp(*, 0:spp(0)-1, 0:spp(1)-1), 1) sh(*, i) = shc(ref, p1, /filt, /N2, POLY=pf) sx = (sall(0, i) < sall(0, 0)) sy = (sall(1, i) < sall(1, 0)) p1 = p1(0:sx-1, 0:sy-1) ; overwrite oldest of nsum images summed to build the ref pp(i MOD nsum, 0:sx-1, 0:sy-1) = shift(p1, sh(0, i), sh(1, i)) print, ' continua shift:', i, sh(*, i) if (plotout) then tvscl, ref ENDELSE ENDFOR ; end loop for Gb, Bc, Rc, Ha ; reverse shifts, set max to zero sx = max(sh(0, *))-sh(0, *) sy = max(sh(1, *))-sh(1, *) ; cut sizes to common area resx = min(sall(0, *)-sx) resy = min(sall(1, *)-sy) ;;; make dimensions even resx = resx/2*2 resy = resy/2*2 ; shift all and write to assoc file FOR j=1, n_ha+3 DO BEGIN FOR i=0, num-1 DO BEGIN writeframe, shift(readframe(j, i), -sx(i), -sy(i)), j, i ENDFOR ENDFOR ; ================================================================ ca ; calcium all by itself since no parallel wide band camera ;;; version 0 = only core, correlate sequence separate and align average ;;; version 1 = flip core/wing, align wing to G Band and core to wing ;;; version 2 = align every image to G Band (whether flipped or not) IF (ca_version EQ 0) THEN BEGIN ;;; version 0 = separate if only core (bad correlation with G band) mkmov, ca, fc, ang=par_ca#replicate(1., num), scale=scl_ca, $ nsum=3, plotout=plotout ref = float(readframe(1, num/2-5)) FOR i=num/2-4, num/2+5 DO ref = ref+readframe(1, i) ; 11 frames midpoint ref = ref(0:resx-1, 0:resy-1) tmp = total(ca(*, *, num/2-5:num/2+5), 3) print,' ca shift:',(sh = shc(ref, tmp, /filt)) ; manblinkshift does tvscl, lets one shift manually for better alignment if (plotout) then begin sh = -manblinkshift(ref, tmp, ini=-sh) stop ; optional manblinkshift endif newx0 = sh(0) > 0 newy0 = sh(1) > 0 ca = ca((-sh(0) > 0):*, (-sh(1) > 0):*, *) ENDIF IF (ca_version eq 1 or ca_version eq 2) THEN BEGIN ;;; Version 1 and version 2 with alignment to G-band ca = intarr(sx0, sy0, num) sca = intarr(2, num) FOR i=0, num-1 DO BEGIN rdfits, p, fc(i) p1 = rot(reform(p(*, *, 0)), par_ca(0), par_ca(1), cubic=-.5) sca(*, i) = size(p1, /dim) ca(0:sca(0, i)-1, 0:sca(1, i)-1, i) = fix(scl_ca*p1 < imax) ENDFOR sh = intarr(2, num) ; version = 1 for wing-gb and core-wing alignments of flipping sequence IF ca_version EQ 1 THEN BEGIN ; read gb skipping every odd image; get shift gb - ca-wing FOR i=1, num-1, 2 DO BEGIN gb = (readframe(1, i))(0:resx-1, 0:resy-1) print,' gb - ca-wing shift: ',i, $ (sh(*, i)=shc(gb, ca(*, *, i), /filt, /n2, poly=pf)) ENDFOR ; stop ; optional inspection of the shifts as instructed below ;; FOR i=1, num-1, 2 DO tvscl,shift(ca(*, *, i), sh(0, i), sh(1, i)) ; apply shift to every odd frame FOR i=1, num-1, 2 DO $ ca(*, *, i) = shift(ca(*, *, i), sh(0, i), sh(1, i)) ; get shift first ca core image to first ca wing image ; tmp: strip to actual size shifted image tmp = sca+sh sh(*, 0) = shc(ca(0:tmp(0, 1), 0:tmp(1, 1), 1), $ ca(*, *, 0), /filt, /n2, poly=pf) ; get shift ca core frames starting i=2 to average of wing neighbours FOR i=2, num-1, 2 DO print,' ca wing-core shift: ', i, $ (sh(*, i)=shc(total(ca(0:tmp(0, i), 0:tmp(1, i), [i-1, i+1]), 3), $ ca(*, *, i), /filt, /n2, poly=pf)) ; stop ; optional inspection of the shifts as above ; apply ca wing-core shift FOR i=0, num-1, 2 DO $ ca(*, *, i) = shift(ca(*, *, i), sh(0, i), sh(1, i)) ENDIF ; version 2 = align all ca frames, flipped or not, to gb IF ca_version EQ 2 THEN BEGIN FOR i=0, num-1 DO BEGIN gb = (readframe(1, i))(0:resx-1, 0:resy-1) print,' ca - gb shift: ', i, $ (sh(*, i)=shc(gb, ca(*, *, i), /filt, /n2, poly=pf)) ENDFOR FOR i=0, num-1 DO ca(*, *, i) = shift(ca(*, *, i), sh(0, i), sh(1, i)) ENDIF ; new origin if there is a ca image with the origin inside old frame newx0 = max(sh(0, *)) > 0 newy0 = max(sh(1, *)) > 0 ; new image sizes newx = min(sca(0, *)+sh(0, *)) newy = min(sca(1, *)+sh(1, *)) ca = (temporary(ca))(newx0:newx-1, newy0:newy-1, *) ENDIF ; apply new origin to the others too IF newx0+newy0 GT 0 THEN BEGIN FOR j=1, n_ha+3 DO BEGIN FOR i=0, num-1 DO BEGIN writeframe, shift(readframe(j, i), -newx0, -newy0), j, i ENDFOR ENDFOR ENDIF ; set the same image size for all wavelengths, in even numbers s = size(ca, /dim)/2*2 resx = ((resx-newx0) < s(0))/2*2 resy = ((resy-newy0) < s(1))/2*2 ; start dir (delete old ones{ get_lun, unit spawn, 'rm -rf ../'+outcubedirname spawn, 'mkdir ../'+outcubedirname ; write all cubes except Ca FOR j=0, n_ha+2 DO BEGIN openw, unit, '../'+outcubedirname+'/'+nam+lam(j)+dat+'-cube-'+resstring([resx, resy, num])+dt FOR i=0, num-1 DO writeu, unit, (readframe(j+1, i))(0:resx-1, 0:resy-1) close, unit ENDFOR ; wite Ca cube openw, unit, '../'+outcubedirname+'/'+nam+lam(n_ha+3)+dat+'-cube-'+resstring([resx, resy, num])+dt FOR i=0, num-1 DO writeu, unit, ca(0:resx-1, 0:resy-1, i) close, unit free_lun, unit close, /all END