; file: sdo_intscale.pro = modified SSW aia_intscale.pro ; init: Jan 2 2014 Rob Rutten Deil ; last: Dec 18 2020 Rob Rutten Deil ;+ function sdo_intscale,index,image,wav,$ clipmin=clipmin,clipmax=clipmax,no_intscale=no_intscale,verbose=verbose ; rescale SDO images ; modification of SSW aia_intscale.pro by Karel Schrijver (2010) ; ; INPUTS: ; index: header of AIA or HMI image ; image: RR-level 2 (officially "level 1.5") AIA or HMI image ; wav: one of the following filename substrings from aia_prep.pro: ; AIA: '1600','1700','4500','94','131','171','193','211','304','335' ; HMI: 'cont','dop','mag' ; ; OPTIONAL KEYWORD INPUT PARAMETERS: ; clipmin: low intensity clipping ; clipmin = 0: no lower clip (or at 5 times standard for mag and dop) ; clipmin = -1: clip at standard minium value defined below ; clipmin > 0: use (small for AIA and cont, negative for mag and dop) ; clipmin < -20: use for mag and dop ; clipmax: high intensity clipping ; clipmax = 0: no upper clip (or at 5 times standard for mag and dop) ; clipmax = -1: clip at standard maximum value defined below ; clipmax = -2: hybrid, as clipmax=-1 for UV+HMI, clipmax=0 for EUV ; clipmax = 0.1 - 20: multiplier to standard max clip ; clipmax > 20: use this value ; no_intscale = 1: do not intscale target AIA and HMI-cont intensities ; verbose = 1/0: print minmax values (default 0) ; ; OUTPUTS: ; rescaled SDO image, signed integers ; ; MODIFICATIONS with respect to aia_intscale.pro: ; in that, exptime turned the image into its type (double for index value), ; output is converted here to integers, with appropriate scaling ; changed some of Karel Schrijver's default clip values (see below) ; removed bytscale-per-image option ; included HMI mag, cont, dop ; ; HISTORY: ; Apr 12 2010 Karel Schrijver (CJS): SSW aia_intscale.pro ; Jan 4 2014 RR: added clipmin and clipmax; added HMI images ; Jan 9 2014 RR: removed orbital and solar-rotation Dopplershift ; Feb 9 2014 RR: HMI cont: take square; log(int) positive ; Jul 23 2015 RR: HMI dop: now solar rotation correction per pixel ; Sep 3 2015 RR: corrected HMI dop (above was only full disk) ; Apr 3 2017 RR: redid clip options ; Jan 11 2018 RR: option no_intscale ; Jun 18 2020 RR: AIA sensitivity-loss correction ; Dec 16 2020 RR: no_intscale skip before exptime, sensitivity rescale ;- ; answer no-parameter query if (n_params() lt 3) then begin sp,sdo_intscale return,-1 endif ; defaults for keywords if (n_elements(clipmin) eq 0) then clipmin=0 if (n_elements(clipmax) eq 0) then clipmax=0 if (n_elements(no_intscale) eq 0) then no_intscale=0 if (n_elements(verbose) eq 0) then verbose=0 ; define permitted level2 SDO (AIA + HMI) filename strings; order CJS wavs=['1600','1700','4500','94','131','171','193','211','304','335',$ ; AIA 'cont','dop','mag'] ; HMI nwavs=n_elements(wavs) ; check wav specification iwav=-1 for iw=0,nwavs-1 do if (wavs[iw] eq wav) then iwav=iw if (iwav eq -1 ) then begin print,' ===== ERROR: wav invalid = ',wav print,' ===== valid: ',wavs return,-1 endif ; ------- skip all AIA scaling, only treat HMI cont and mag if (no_intscale eq 1) then begin im=image ; next may need fine tuning? Max must be < 32000 for integers if (iwav eq 10) then im=fix(im/2.5+0.5) ; Prasad: hmi_cont up to 72000 goto,ENDSCALEAIA endif ; clip specifications per wavelength in RR channel numbers ; 1600, 1700, 4500, 94, 131, 171, 193, 211, 304, 335, cont, dop, mag ; iwav 0 1 2 3 4 5 6 7 8 9 10 11 12 ;------------------------------------------------------------------------- normtime=[2.99911, 1.00026, 1.00026, 4.99803, 6.99685, 4.99803, 2.99950,$ 4.99801, 4.99941, 6.99734,$ ; AIA norm exposure durations 1.0, 1.0, 1.0] ; HMI none defmin=[0, 0, 0, 1.5, 7, 10, 120, 30, 15, 3.5,$ ; AIA 0., -8000., -5000.] ; HMI ;; defmax=[1000, 2500, 26000, 50, 1200, 6000, 6000, 13000, 600, 1000, ; CJS ;; defmax=[4000, 10000, 26000, 100, 1200, 12000, 12000, 13000, 3000, 1000,$ ;; 100000., +8000., +5000.] ; RJR original at some early sampling defmax=[3000, 4000, 26000, 100, 1200, 15000, 12000, 13000, 3000, 1000,$ 100000., +8000., +5000.] ; RJR Jun 20 2020 ;------------------------------------------------------------------------- ; apply exposure normalization factor (1.0 for HMI) if (iwav lt 10) then exprat=float(normtime[iwav]/index.exptime) $ else exprat=1.0 im=image*exprat ;RR im = don't change input image; now it is float ;; ; check ;; print,exprat,minmax(im) ;; STOP ; ---- get and apply AIA sensitivity correction for temporal decay ; HMI sensdrop=1. ; UV channels if (iwav lt 2) then begin if (fix(!version.release) ge 7) then begin ; newer IDL but I don't know what the divider should be; works in UU 8.6.0 result=aia_bp_get_corrections(/uv) junk=min(abs(anytim(index.t_obs)-anytim(result.utc)),indclose) sensdrop=result.ea_over_ea0[indclose[0],iwav]>0.1 endif else begin ; needed for my ancient (2007) IDL 6.4 (getdata and data in fixsswlib) uvfile=file_search('~','sdo_uvdecay.dat') ; tilde option not in Windows? readcol,uvfile[0],dattai,dat1600,dat1700,format='D' ; there may be more tai_obs=anytim2tai(index.t_obs) junk=min(abs(tai_obs-dattai),indclose) if (iwav eq 0) then sensdrop=dat1600[indclose]>0.1 if (iwav eq 1) then sensdrop=dat1700[indclose]>0.1 endelse ;; ; check (Tosh or UU) ;; print,' ----- sensdrop ='+trimd(sensdrop) ;; STOP endif ; EUV channels if (iwav gt 2 and iwav lt 10) then begin result=aia_bp_get_corrections() ;; print,result.channels ; aha, same increase order as I use ;; 94 131 171 193 211 304 335 channel=iwav-3 ;; utplot,result.utc,result.ea_over_ea0(*,5) ;RR nice program! junk=min(abs(anytim(index.t_obs)-anytim(result.utc)),indclose) sensdrop=result.ea_over_ea0[indclose[0],channel]>0.1 endif ; apply sensdrop correction ; (I assume no AIA*exprat*sensdrop exceeds integer wordlength 32000) im=im/sensdrop ; set minimum clip (see keyword explanation above) if (clipmin eq 0) then begin cmin=0.0 if (iwav gt 3 and iwav lt 10) then cmin=1.0 ; enable log for EUVs if (iwav ge 11) then cmin=5.*defmin[iwav] ; negative mag and dop bounds endif if (clipmin eq -1) then cmin=defmin[iwav] if (clipmin gt 0) then cmin=clipmin if (clipmin lt -20) then cmin=clipmin ; set maximum clip (see keyword explanation above) if (clipmax eq 0) then begin cmax=32767 ; i.e., no clip for integer files if (iwav ge 10) then cmax=5.*defmax[iwav] ; HMI floats endif if (clipmax eq -1) then cmax=defmax[iwav] if (clipmax eq -2) then begin if (iwav gt 2 and iwav lt 10) then cmax=32767 if (iwav lt 3 or iwav gt 9) then cmax=defmax[iwav] endif if (clipmax ge 0.1 and clipmax le 20) then cmax=clipmax*defmax[iwav] if (clipmax gt 20) then cmax=clipmax ; check if (verbose) then $ print,' ----- minmax(im) before clipping ='+trimd(fix(minmax(im)))+$ ' clipmin ='+trimd(fix(cmin))+' clipmax ='+trimd(fix(cmax)) ; apply clips im=im>cmin im=im(-32000) endif if (clipmax eq 0 and max(im) ge 32000) then begin print,' ====== sdo_intscale.pro: clipmax=0; wav = '$ +wav+' has max='+trim(max(im))+'; clipped' im=im<32000 endif ENDSCALEAIA: ; HMI Dopplergrams: correct for the components of SDO (and Earth) ; orbital motion and of solar rotation along the line of sight (LOS) if (iwav eq 11) then begin ; get solar field-of-view center location in (X,Y) = arcsec from disk center xcen=index.xcen ycen=index.ycen ; use the following if xcen and ycen get deprecated from level2 index ;; a=index.crota2 ;; xcen=index.crval1+index.cdelt1*cos(a)*((index.naxis1+1)/2-index.crpix1) $ ;; -index.cdelt2*sin(a)*((index.naxis2+1)/2-index.crpix2) ;; ycen=index.crval2+index.cdelt1*sin(a)*((index.naxis1+1)/2-index.crpix1) $ ;; +index.cdelt2*cos(a)*((index.naxis2+1)/2-index.crpix2) xcenrad=xcen*!pi/(60.*60.*180.) ; radians on the sky ycenrad=ycen*!pi/(60.*60.*180.) ; radial observer motion = (sun center - SDO speed projected on LOS sdorad=index.obs_vr*cos(sqrt(xcenrad^2+ycenrad^2)) ; obs_vw: Westward orbit component affecting off-center LOS ; obs_vn: Northward orbit component affecting off-center LOS ; the West component is largest, up to 150 m/s at limb (from Earth orbit) sdowest=index.obs_vw*sin(-xcenrad) sdonorth=index.obs_vn*sin(-ycenrad) ; LOS component of solar rotation ; ------------------------------- ; is there no keyword for this? Reinvent the wheel... ; use WCS comands from /home/rutten/rr/ssw/gen/doc/wcs_tutorial.pdf ; get Stonyhurst latitude matrix wcs=fitshead2wcs(index) coord=wcs_get_coord(wcs) wcs_convert_from_coord,wcs,coord,'hg',lon,lat lat=lat*!pi/180. ; px latitude values in radians ; get solar rotation (also image matrix) a=14.713 ; differential rotation polynomial b=-2.936 ; Snodgrass+Ulrich 1990ApJ...351..309S c=-1.787 ; from http://en.wikipedia.org/wiki/Sor_rotation omega=a+b*sin(lat)^2+c*sin(lat)^4 ; differential rotation in degrees/day period=360.*24.*60.*60./omega ; rotation period in s (Y dependent) ;; radius=6.96E8*cos(lat) ; trajectory radius in m (Y dependent) rotvel=2*!pi/period ;;*radius ; rotation speed in m/s (Y dependent) ; Jul 23 2015: below = old version = only correction with center-pixel value ;; rotcor=rotvel*xcen/(index.rsun_obs/6.96E8) ;;/radius ; sin(theta_X) comp ;; ; NB: radius cancels: solid rotator has same Doppler at all Y for given X ; Jul 23 2015: Peter Gomory: make X-dependent over image (eg full disk image) ; Sep 3 2015: OOPS, wanted rotvel[2048,*] = full image nx=index.naxis1 ny=index.naxis2 delx=index.cdelt1 xmin=xcen-(nx/2)*delx radscale=index.rsun_obs/6.96E8 rotcor=fltarr(nx,ny) for ix=0,nx-1 do rotcor[ix,*]=rotvel[nx/2,*]*(xmin+ix*delx)/radscale ; apply summed corrections, change into integers scale 10x (units of 0.1 m/s) ; sign: redshift positive; "black = blue" as DOT; East limb (lefthand) dark dopcor=sdorad+sdowest+sdonorth+rotcor im=fix((im-dopcor)*10.) endif ; HMI magnetograms: convert from 1 Gauss floats to 0.1 Gauss integers if (iwav eq 12) then im=fix(im*10.) ; check clip limits ;; histim=histo_opt_rr(im,/top_only) ;; print,' ----- wav '+string(wav)+$ ;; ' minmax(image) = '+trim(min(image))+', '+trim(max(image))+$ ;; ' minmax(imin) = '+trim(min(imin))+', '+trim(max(imin))+$ ;; ' minmax(im) = '+trim(min(im))+', '+trim(max(im))+$ ;; ' minmax(histim) = '+trim(min(histim))+', '+trim(max(histim)) ; done return,im end ; =============== main for test per IDLWAVE H-c ====================== cd,'~/data/SDO/2014-06-14-small/target/level2' ; ~ for UU also ;; cd,'/media/rutten/RRHOME/alldata/SDO/2018-08-30-tests17001600/2011-01-01/target/level2' ;; cd,'/media/rutten/RRHOME/alldata/SDO/2018-08-30-tests17001600/2018-01-01/target/level2' ;; cd,'/home/rutten/data/SDO/2019-11-11-transit/midpointfull/center/level2' ;; EUV wavs: 94 131 171 193 211 304 335 ;; wav='304' ;; wav='171' ;; wav='211' ;; wav='94' ;; wav='131' ;; wav='335' ;; wav='193' wav='1600' ;; wav='1700' ; get file files=findfile('*'+wav+'.fits*') if (strmatch(files[0],'*fz')) then spawn,'funpack -D '+files[0] files=findfile('*'+wav+'.fits') imagefile=files[0] ; process read_sdo,imagefile,index,image ;; newim=sdo_intscale(index,image,wav,clipmin=-1,clipmax=-2.,/no_intscale) newim=sdo_intscale(index,image,wav,clipmin=-1,clipmax=-1,/verbose) print,' ----- t obs =',index.t_obs print,' ----- minmax(newim) ='+trimd(minmax(newim)) histim=histo_opt_rr(newim,1E-4,/top_only) print,' ----- minmax(histim) ='+trimd(minmax(histim)) ; inspect showex,image,newim,histim,/blink end