; file: sdo_addfires.pro ; init: Nov 18 2020 Rob Rutten Deil ; last: Dec 21 2020 Rob Rutten Deil ;+ pro sdo_addfires,aia3013cut,aia3013color,$ cubesdir=cubesdir,file304=file304,file131=file131 ; write 304x131 = aia3013.fits and aiafire.fits in cubesdir ; ; INPUTS: ; aia3013cut: clip greyscale at this value ; aia3013color: set px above this level to aia3013cut ; ; OPTIONAL KEYWORD INPUTS: ; cubesdir: to put results into, default 'sdo/target/cubes' ; file304: input file, default cubesdir+'/aia304.fits' ; file131: input file, default cubesdir+'/aia131.fits' ; ; OUTPUTS: ; aia3013.fits and aiafire.fits in cubesdir ; ; NOTE: ; no diffheight correction, seems not necessary for 304x131 fires ; (likely particle beams give joint heating at specific hit depth) ; ; HISTORY: ; Nov 18 2020 RR: start ; Nov 29 2020 RR: switch to assoc to enable long sequences ;- ; answer wrong-parameter query if (n_params(0) ne 2) then begin sp,sdo_addfires return endif ; defaults for keywords if (n_elements(cubesdir) eq 0) then cubesdir='target/cubes' if (n_elements(file304) eq 0) then file304=cubesdir+'/aia304.fits' if (n_elements(file131) eq 0) then file131=cubesdir+'/aia131.fits' ; ========== four-file assoc business after coalignfitscubes.pro ; set endian bigendian=1 ; get file304 dimensions and datatype from fits header headin=headfits_rr(file304) headinsize=(1+fix(n_elements(headin)/36.))*2880 nx=fxpar(headin,'naxis1') ny=fxpar(headin,'naxis2') nt=fxpar(headin,'naxis3') bitpix=fxpar(headin,'bitpix') ; check integer if (bitpix ne 16) then begin print,' ##### sdo_addfires abort: aia input fitscubes not integer' return endif ; open file304 for assoc = inassoc1 get_lun,unit_in1 if (bigendian) then openr,unit_in1,file304,/swap_if_little_endian $ else openr,unit_in1,file304 inassoc1=assoc(unit_in1,intarr(nx,ny),headinsize) ; open file131 for assoc = inassoc2 get_lun,unit_in2 if (bigendian) then openr,unit_in2,file131,/swap_if_little_endian $ else openr,unit_in2,file131 inassoc2=assoc(unit_in2,intarr(nx,ny),headinsize) ; open aia3013.fits for assoc = outassoc1 outfile1=cubesdir+'/aia3013.fits' get_lun,unit_out1 headaia3013=headin sxaddpar,headaia3013,'channel','AIA 304x131','SDO diagnostic name' sxaddpar,headaia3013,'fire max clip',aia3013cut,$ 'intensity clip to show chromosphere' sxaddpar,headaia3013,'fire level clip',aia3013color,$ 'intensity threshold to color fires' headaia3013size=(1+fix(n_elements(headaia3013)/36.))*2880 if (bigendian) then openw,unit_out1,outfile1,/swap_if_little_endian $ else openw,unit_out1,outfile1 outassoc1=assoc(unit_out1,intarr(nx,ny),headaia3013size) rec=assoc(unit_out1, bytarr(headaia3013size)) rec[0]=byte(headaia3013) ; open aiafire.fits for assoc = outassoc2 outfile2=cubesdir+'/aiafire.fits' get_lun, unit_out2 headfire=headaia3013 sxaddpar,headfire,'channel','AIA fire detector','SDO diagnostic name' headfiresize=(1+fix(n_elements(headfire)/36.))*2880 if (bigendian) then openw,unit_out2,outfile2,/swap_if_little_endian $ else openw,unit_out2,outfile2 outassoc2=assoc(unit_out2,intarr(nx,ny),headfiresize) rec=assoc(unit_out2, bytarr(headfiresize)) rec[0]=byte(headfire) ; now the act for it=0,nt-1 do begin ; product a304=inassoc1[it] a131=inassoc2[it] aia3013=float(a304)*a131 aia3013=fix(aia3013/30000.+0.5) ; same values as in sdo_firelevel.pro outassoc1[it]=aia3013 ; fire detector aia3013cutint=fix(aia3013cut+0.5) ; cut bright ARs aia3013clip=aia3013 if (max(aia3013) gt aia3013cutint) then aia3013clip=aia3013