; file: movex_scatcont.pro ; init: Jul 20 2017 Rob Rutten Deil ; last: Mar 1 2018 Rob Rutten Deil ;+ pro movex_scatcont,infiles,iwavA,iwavB,$ nt_mw=nt_mw,itdelay=itdelay,doppA=doppA,doppB=doppB,smear=smear,$ cadence=cadence,xrange=xrange,yrange=yrange,trange=trange,tskip=tskip,$ ; ----- SDO and SST mwfwavsfiles=mwfwavsfiles,$ allsdo=allsdo,sdodir=sdodir,allmwf=allmwf,mwfdirs=mwfdirs,$ allfits=allfits,fitsdirs=fitsdirs,$ ; ----- scatcont plotrangeA=plotrangeA,plotrangeB=plotrangeB,$ outerlevel=outerlevel,nbins=nbins,contourstep=contourstep,$ addmoments=addmoments,addhistograms=addhistograms,$ ; ----- options for plotoutput xaxtitle=xaxtitle,yaxtitle=yaxtitle,$ longlabel=longlabel,psdir=psdir,_extra=plotkeywords ; NAME: ; movex_scatcont ; ; PURPOSE: ; apply scatcont.pro to DOT/SDO/IRIS/SST/CRISP file(s) ; ; DESCRIPTION: ; see comments in program and example at file bottom ; ; CALL: ; see above ; ; INPUTS: ; most from movex.pro, see there ; ; OPTIONAL KEYWORD INPUTS: ; idem ; tskip = range of time steps to be skipped ; smear = nr of px to smear over; if negative then rebin ; longlabel = only t_delay or more info ; ; OUTPUTS: ; scatcont plots ; ; HISTORY: ; Jul 21 2017 RR: start as demo for SacPeakFinal workshop ;- ; answer no-parameter query if (n_params(0) lt 3) then begin sp,movex_scatcont return endif ; defaults for keywords if (n_elements(nt_mw) eq 0) then nt_mw=0 if (n_elements(iwavA) eq 0) then iwavA=0 if (n_elements(iwavB) eq 0) then iwavB=1 if (n_elements(itdelay) eq 0) then itdelay=0 if (n_elements(doppA) eq 0) then doppA=0 if (n_elements(doppB) eq 0) then doppB=0 if (n_elements(smear) eq 0) then smear=0 if (n_elements(cadence) eq 0) then cadence=0 if (n_elements(xrange) eq 0) then xrange=[0,-1] if (n_elements(yrange) eq 0) then yrange=[0,-1] if (n_elements(trange) eq 0) then trange=[0,-1] if (n_elements(tskip) eq 0) then tskip=[-1,-1] ; ----- special SDO and SST if (n_elements(allsdo) eq 0) then allsdo=0 if (n_elements(sdodir) eq 0) then sdodir='sdo2sst/' if (n_elements(allmwf) eq 0) then allmwf=0 if (n_elements(allfits) eq 0) then allfits=0 ; ----- scatcont if (n_elements(nbins) eq 0) then nbins=50 if (n_elements(outerlevel) eq 0) then outerlevel=20 if (n_elements(contourstep) eq 0) then contourstep=3 if (n_elements(addmoments) eq 0) then addmoments=1 if (n_elements(addhistograms) eq 0) then addhistograms=1 if (n_elements(xaxtitle) eq 0) then xaxtitle='A' if (n_elements(yaxtitle) eq 0) then yaxtitle='B' if (n_elements(longlabel) eq 0) then longlabel=0 if (n_elements(psdir) eq 0)then psdir='.' ; renew lun counter common lunsofmovexfiles,ilun_last ilun_last=0 ; open lun's for all files and get file sampling parameters addlabels=1 movex_loadfiles,infiles,$ ; ----- outputs nxfile,nyfile,ntfile,nfiles,files,nwavs,$ lunfile,facfile,nzfiles,bitpixfile,hofffile,swapendianfile,$ mwfile,nwavsfile,$ iwstartfile,iwendfile,iwlcfile,iwinfile,filewav,labelwav,$ xcutmin,xcutmax,nxcut,ycutmin,ycutmax,nycut,xyrangecut,$ tcutmin,tcutmax,ntcut,duplicate,instance=instance,$ ; ----- keyword inputs nt_mw=nt_mw,mwfwavsfiles=mwfwavsfiles,addlabels=addlabels,$ xrange=xrange,yrange=yrange,trange=trange,$ ; ----- special keyword parameters for SDO and SST data allsdo=allsdo,sdodir=sdodir,allmwf=allmwf,mwfdirs=mwfdirs,$ allfits=allfits,fitsdirs=fitsdirs ; set assoc pointer per file ass2pt=ptrarr(nfiles,/allocate_heap) for ifile=0,nfiles-1 do begin if (bitpixfile[ifile] eq 8) then $ *ass2pt[ifile]=assoc(lunfile[ifile],$ bytarr(nxfile,nyfile,/nozero),$ hofffile[ifile]) if (bitpixfile[ifile] eq 16) then $ *ass2pt[ifile]=assoc(lunfile[ifile],$ intarr(nxfile,nyfile,/nozero),$ hofffile[ifile]) if (bitpixfile[ifile] eq -32) then $ *ass2pt[ifile]=assoc(lunfile[ifile],$ fltarr(nxfile,nyfile,/nozero),$ hofffile[ifile]) endfor ; ====== sum data arrays for scatcont ; later start for negative itdelay itstart=tcutmin itend=tcutmax if (itdelay lt 0 and tcutmin+itdelay lt 0) then itstart=-itdelay ; loop over it range for it=itstart,itend do begin ; skip if within tskip range if (it ge tskip[0] and it le tskip[1]) then goto, SKIPTHISIT ; don't go per delay beyond end of file if (it+itdelay eq nt_mw) then begin itend=it-1 break endif ; get next image from A iframeA=it*nwavsfile[filewav[iwavA]]+iwinfile[iwavA] imageA=(*ass2pt[filewav[iwavA]])[iframeA] imageA=imageA[xcutmin:xcutmax,ycutmin:ycutmax] ; doppA set: get dopplergram if (doppA ne 0) then begin iwlc=iwlcfile[filewav[iwavA]] iwin=iwinfile[iwavA] iw2in=iwlc+(iwlc-iwin) iw2=iw2in+iwstartfile[filewav[iwavA]] iframe2=it*nwavsfile[filewav[iw2]]+iwinfile[iw2] image2=(*ass2pt[filewav[iwavA]])[iframe2] image2=image2[xcutmin:xcutmax,ycutmin:ycutmax] imageA=(float(image2)-imageA)/(image2+imageA) endif ; get matching image (at time itdelay) from B iframeB=(it+itdelay)*nwavsfile[filewav[iwavB]]+iwinfile[iwavB] imageB=(*ass2pt[filewav[iwavB]])[iframeB] imageB=imageB[xcutmin:xcutmax,ycutmin:ycutmax] ; doppB set: get dopplergram if (doppB ne 0) then begin iwlc=iwlcfile[filewav[iwavB]] iwin=iwinfile[iwavB] iw2in=iwlc+(iwlc-iwin) iw2=iw2in+iwstartfile[filewav[iwavB]] iframe2=(it+itdelay)*nwavsfile[filewav[iw2]]+iwinfile[iw2] image2=(*ass2pt[filewav[iwavB]])[iframe2] image2=image2[xcutmin:xcutmax,ycutmin:ycutmax] imageB=(float(image2)-imageB)/(image2+imageB) endif ; optional smear and optionally rebin (both) if (smear gt 1) then begin imageA=smooth(imageA,smear,/edge_truncate) imageB=smooth(imageB,smear,/edge_truncate) endif if (smear lt 0) then begin smear=abs(smear) imageA=imageA[0:fix(nxcut/smear)*smear-1,$ 0:fix(nycut/smear)*smear-1] imageB=imageB[0:fix(nxcut/smear)*smear-1,$ 0:fix(nycut/smear)*smear-1] imageA=rebin(imageA,fix(nxcut/smear),fix(nycut/smear)) imageB=rebin(imageB,fix(nxcut/smear),fix(nycut/smear)) endif ; store in large data arrays for scatcont if (it eq itstart) then begin dataA=imageA dataB=imageB endif else begin dataA=[dataA,imageA] ; side-by-side concatenation, first dim = nx*nt dataB=[dataB,imageB] endelse SKIPTHISIT: endfor ; end of it loop over time steps ; get sample size for label szim=size(imageA) nxim=szim[1] nyim=szim[2] npairs= float(nxim)*nyim*(itend-itstart+1-(tskip[1]-tskip[0])) ; ===== set up and run scatcont ; filename (mucklast = make delay < 0 have proper-order file name) if (itdelay lt 0) then mucklast=tcutmax-1 else mucklast=tcutmax if doppA then typeA='D' else typeA='I' if doppB then typeB='D' else typeB='I' if (smear ne 0) then $ sm='_sm'+strmid(string(smear+1E2,format='(I3)'),1)+'_' else sm='_' if (tskip[0] ne -1) then skip='skip_'+trim(tskip[0])+'_'+trim(tskip[1])+'_' $ else skip='' psfile=psdir+'/'+'scat_'+trim(iwavA)+typeA+'_'+trim(iwavB)+typeB+sm+skip+$ trim(xcutmin)+'_'+trim(xcutmax)+'_'+$ trim(ycutmin)+'_'+trim(ycutmax)+'_'+$ strmid(string(tcutmin+1E3,format='(f5.0)'),1,3)+'_'+$ strmid(string(mucklast+1E3,format='(f5.0)'),1,3)+'_'+$ strmid(string(itdelay+1E3,format='(f5.0)'),1,3)+'.ps' ; labels if (cadence eq 0) then label='delay '+trim(itdelay)+' time steps' else $ label='t_delay = '+trim(itdelay*cadence/60.,'(F6.1)')+' min' if (longlabel ne 0) then begin if (smear ne 0) then $ label=[label,'smear = '+trim(smear)] label=[label,'x_range = ['+trim(xcutmin)+','+trim(xcutmax)+']'] label=[label,'y_range = ['+trim(ycutmin)+','+trim(ycutmax)+']'] label=[label,'t_range = ['+trim(itstart)+','+trim(itend)+']'] if (tskip[0] ne -1) then $ label=[label,'t_skip = ['+trim(tskip[0])+','+trim(tskip[1])+']'] label=[label,'N_pairs = '+trim(npairs,'(E8.1)')] endif ; titles if (xaxtitle ne 'A') then xtitle=xaxtitle else begin if (doppA eq 0) then xtitle='intensity '+labelwav[iwavA] else $ xtitle='Dopplergram amplitude '+labelwav[iwavA] endelse if (yaxtitle ne 'B') then ytitle=yaxtitle else begin if (doppB eq 0) then ytitle='intensity '+labelwav[iwavB] else $ ytitle='Dopplergram amplitude '+labelwav[iwavB] endelse ;; ;; ; special titles ;; ;; if (labelwav[iwavA] eq ' widwid_ha_ca ') then xtitle=$ ;; ;; textoidl("w_{H\alpha} (w_{H\alpha} - w_{Ca 8542})") ;; ;; if (labelwav[iwavB] eq ' mxis_depth_10830 ') then ytitle=$ ;; ;; 'relative depth He I 10830' ; call scatcont scatcont,dataA,dataB,psfile=psfile,$ plotrangeA=plotrangeA,plotrangeB=plotrangeB,$ outerlevel=outerlevel,nbins=nbins,contourstep=contourstep,$ addmoments=addmoments,addhistograms=addhistograms,$ psdir=psdir,xtitle=xtitle,ytitle=ytitle,label=label ; ----- clean and exit ; close all lun's for ifile=0,nfiles-1 do free_lun,lunfile[ifile] ; print time of job completion spawn,'date' end ; =============== main for testing per IDLWAVE H-c ====================== ; data (where first on my internal disk for Sac Peak meeting) ;; cd,'~/data/SST/2014-06-21-quiet' cd,'/media/rutten/SSTDATA/alldata/SST/2014-06-21-quiet' infiles=['sst/widwid_ha_ca.fits',$ ;; 'sst/shift_kevin_6563_0.9.icube',$ ;; 'sst/width_kevin_6563_0.9.icube',$ ;; 'sst/shxwi_kevin_6563_0.9.icube',$ ;; 'sst/shift_kevin_8542_0.6.icube',$ ;; 'sst/width_kevin_8542_0.6.icube',$ ;; 'sst/shxwi_kevin_8542_0.6.icube',$ 'iris2sst/sji1400_algbo.fits',$ 'iris2sst/sji2796_algbo.fits',$ 'crispex/crispex.6563.08:02:31.time_corrected.aligned.icube',$ 'crispex/crispex.8542.08:02:31.time_corrected.icube'] mwfwavsfiles=['crispex/spectfile.6563.idlsave',$ 'crispex/spectfile.8542.idlsave'] ;; infiles=['iris2sst/sji1400_algbo.fits','iris2sst/sji2796_algbo.fits',$ ;; 'crispex/crispex.6563.08:02:31.time_corrected.aligned.icube'] ;; mwfwavsfiles='crispex/spectfile.6563.idlsave' nt_mw=377 cadence=11.48 smear=0 ;; ; select full frame ;; xrange=[50,870] ; Halpha full frame trimbox cutting black edges ;; yrange=[50,870] ; idem ; select whether contrail area only or full-frame (comment out) xrange=[250,750] ; contrail area only yrange=[250,750] trange=[75,95] ; contrail precursor duration allsdo=1 sdodir='sdo2sst/' psdir='/tmp' contourstep=2 ; contours at factor 2 in density ;;;; tskip=[80,100] ; period of the precursor ; setup A iwavA=11 ; widwid plotrangeA=[-1500,3500] ; setup B ;; iwavB=21 ; Ha_lc ;; plotrangeB=[1000,3500] iwavB=4 ; HeII 304 plotrangeB=[0,1.6E4] ; set time delay sampling: only three delays or full movie (SLOW; staff?) ;; for itdelay=-75,75 do begin ; full movie (SLOW and eats laptop!) ;; for itdelay=-16,16,16 do begin ;; for itdelay=-4,4,4 do begin for itdelay=0,0 do begin movex_scatcont,infiles,iwavA,iwavB,$ nt_mw=nt_mw,itdelay=itdelay,doppA=doppA,doppB=doppB,smear=smear,$ cadence=cadence,xrange=xrange,yrange=yrange,trange=trange,tskip=tskip,$ ; ----- SDO and SST mwfwavsfiles=mwfwavsfiles,$ allsdo=allsdo,sdodir=sdodir,allmwf=allmwf,mwfdirs=mwfdirs,$ allfits=allfits,fitsdirs=fitsdirs,$ ; ----- scatcont keywords plotrangeA=plotrangeA,plotrangeB=plotrangeB,$ outerlevel=outerlevel,nbins=nbins,contourstep=contourstep,$ addmoments=addmoments,addhistograms=addhistograms,$ xaxtitle=xaxtitle,yaxtitle=yaxtitle,longlabel=longlabel,$ psdir=psdir,_extra=plotkeywords endfor ;; ; ps > png ;; spawn,'cd '+psdir+'; /home/rutten/rr/bin/pscroppngall' end