file: idl-cube-manual.txt = data cube tools
last: Jun 17 2021  Rob Rutten  Deil


    IDL DATA CUBE ANALYSIS TOOLS FOR SOLAR IMAGE SEQUENCES

                     Robert J. Rutten

                 Lingezicht Astrophysics
          Institutt for Teoretisk Astrofysikk Oslo 
   
   
This IDL tutorial introduces various (x,y,t) data cube analysis tools.  
It consists of examples using two small solar image sequences, with
links to pertinent publications.

There are parallel txt, pdf, and html versions of this manual at 
  <<https://robrutten.nl/Manuals.html>>
The html and pdf versions have active weblinks.

This manual is a companion to my "Simple IDL instruction for astronomy
students", also at
  <<https://robrutten.nl/Manuals.html>>
It uses the same style, supplying IDL entries for cut-and-paste onto
the IDL command line.  It was initially written for DOT students:
  <<https://robrutten.nl/Students_to_the_DOT_program.html>>

   
SolarSoft library
=================
SolarSoft = "ssw" = large IDL library for solar physics, needed here
  <<http://www.lmsal.com/solarsoft>>
  The menu window is at the top of the page.  This large conglomerate
  runs from an evelope caller sswidl.  In my linux laptop I use a bin
  script "xssw" to start it up:
     --------------------------------
    #!/bin/csh
    setenv SSW /usr/local/ssw                      # if ssw stuff sits here
    setenv SSW_INSTR "sot aia hmi trace ontology"  # select instruments 
    source $SSW/gen/setup/setup.ssw
    sswidl
    --------------------------------
  SolarSoft upgrade example from 
  <<http://www.lmsal.com/solarsoft/ssw_upgrades.html>>
    IDL> ssw_upgrade,/sot,/aia,hmi,/spawn,/loud        # update ssw software


Start
=====
Get and unzip:
  <<https://robrutten.nl/rrweb/rjr-edu/manuals/idl-cube.zip>>

It gets you some IDL routines and two (x,y,t) data cube files with
small cutouts of short, simultaneous solar image sequences with 30-sec
cadence obtained with the Dutch Open Telescope (DOT)
  <<https://robrutten.nl/dot/DOT_home.html>>

Subdirectory demo-output has my resulting figures from the below commands

  > xssw   ; or your own wrapper call for SolarSoft ssw

  xslice                  ; pro name without parameters returns parameter list
  doc_library,'xslice'    ; get the information block at the top of the pro

  gb=readfits('DOT-2005-10-14-gb.fits')  ; load this cube into memory
  ca=readfits('DOT-2005-10-14-ca.fits')  ; idem; these files are small

Comment on content: for this instruction you don't need to appreciate
what these two data cubes show, but let me summarize nevertheless: 
"gb" = G-band intensity which samples the deep photosphere; 
"ca" = Ca II H line-center intensity which samples the high photosphere.  
The small field of view cutouts cover a very quiet (sunspot- and
plage-free) area.  Scale: the pixels are 0.071 arcsec; the largest
granules seen in gb equal France in size.  The small bright
intergranular blobs seen in gb are kilogauss magnetic concentrations
("fluxtubes").  These DOT data were used in 
  RJR{2008SoPh..251..533R}.
  

Cube inspection and comparison tools
====================================

  xslice,gb,root=root,mag=2,ytslice=0
  xslice,ca,root=root

     ; cube slicer from Alfred de Wijn in SSW.  Multiple cubes are sliced
     ; in concert, with synchronous cursors, to study co-locational
     ; phenomena in space and time.  Position the x-t and x-y windows
     ; below each other, ca next to gb.  Play with your mouse and
     ; mouse buttons: right-mouse plays the x-y movie forward,
     ; left-mouse backwards; middle-mouse stops.  Mouse movement moves
     ; the crosshairs.  The x-t slices show the evolution of the x cut
     ; under the horizontal crosshair in the x-y movies.  The default
     ; adds y-t slices but I find them confusing; I wiggle the mouse
     ; in y to see what happens close to a particular x cut.  Option
     ; phitslice adds slice panels along a tiltable cut.  

     ; Comment on content: note the "three-minute" oscillations in Ca
     ; II H x-t.  Such gb-ca slicing triggered
       RJR{2005A&A...441.1183D}.
   

  showex,'DOT-2005-10-14-gb.fits','DOT-2005-10-14-ca.fits'

    ; showex.pro calls movex.pro, my remake of the Oslo ximovie.pro
    ; player.  Type movex to get a parameter list.  It assocs into
    ; disk files and so can play large cubes that exceed the memory.
    ; It can blink two movies out of many, or play each separately.
    ; You can resize and also zoom in by pulling out a piece.  Wrapper
    ; showex.pro can also play cube variables in memory.


  plot,gb,ca,psym=3    

    ; correlative scatterplot pixels gb[x,y,t] against pixels ca[x,y,t]

  scatcont,gb,ca,xtitle='G band intensity',ytitle='Ca II H intensity',$
    nbins=50,outerlevel=30,/moments,/hists,xrange=[0,200]

    ; similar scatterplot but with contours to avoid plot saturation

    ; Comment on content: the bright-bright flag denotes fluxtubes;
    ; the slight countercorrelation around the peak denotes reversed
    ; granulation

    scatcont,gb[*,*,0:9],ca[*,*,20:29],$
      nbins=50,outerlevel=30,/moments,/hists,xrange=[0,200]

    ; Comment on content: the countercorrelation is gone at a delay of
    ; 20 time steps = 10 min, while the flag splits between stationary
    ; and moving fluxtubes



Fourier cube inspection and comparison tools
============================================

  plotpowermap,ca,30,0.005,0.007,normalization=1,/ps,plotfile='powermap.ps'
  $gv powermap.ps
     ; plots a Fourier power map for the 0.005-007 Hz frequency band.
     ; It is a pastiche of Alfred de Wijn codes plus my old frimage.pro.  
     ; Setting pxdetrend=2 is very slow, but the same transform can be
     ; re-used setting /olddata via a file copy on /tmp.  

     ; Comment on content: the 0.05-0.07 Hz frequency band samples
     ; acoustic waves on their way up to become chromospheric shocks.
     ; Different normalizatons are described in
       RJR{2001A&A...379.1052K}.
     ; A nice display of how the choice of normalization influences 
     ; high-frequency powermaps is given in Figure 5 of
       RJR{2005A&A...430.1119D}.
   
   
  plotconfusogram,gb,ca,30,'ca-gb-confusogram.ps'
  $gv ca-gb-confusogram.ps

     ; Fourier analysis in the "confusogram" format initiated by Bruce
     ; Lites following White & Athay (1979).  The upper panel shows
     ; binned pixel-by-pixel ca-gb Fourier phase difference spectra
     ; with their spatial average.  The lower panel shows the
     ; corresponding spatially-averaged power spectra and the
     ; coherence spectrum with noise level.  The math is given in 
  RJR{2001A&A...379.1052K}.
     ;  Alfred de Wijn wrote the codes constituting this program
     ;  for application to similar DOT data described in
  RJR{2004A&A...416..333R}.
     ; Alfred's plot layout is less confusing than Lites-Krijger
     ; confusograms since the power and coherence spectra are in a
     ; separate panel.  

     ; Comment on content: high coherence only below 3 mHz.  The phase
     ; differences appear to be random above 7 mHz.  The negative
     ; phase differences at low frequency demonstrate convective
     ; reversal and gravity waves.  The zero value at 5-min
     ; periodicity shows that the acoustic waves that produce the
     ; global oscillations of the sun by propagating through its
     ; interior are evanescent in the photosphere.  The positive
     ; values at larger frequencies signify upwards propagating
     ; acoustic waves.  I wonder how the next to last column gets so
     ; artifactual (setting fmax=0.8 avoids it).
   
   
  plotkopower,ca,0.071,30,'ko-power-ca.ps',maxpower=1,kmax=0.3,/contours,$
    lamb=7,/fundamental 
  $gv ko-power-ca.ps

     ; Plots a traditional "k-omega" diagram, but nowadays one plots
     ; temporal frequency f rather than circle frequency omega=2*pi*f
     ; along the vertical axis.  The horizontal axis is horizontal
     ; wavenumber k_h = sqrt(k_x^2+k_y^2).  (In helioseismology one
     ; uses spherical harmomic l instead, i.e. the number of global
     ; mode nodes around the Sun with k_h^2 = l(l+1)/R_sun^2.)  The
     ; slanted white line is the Lamb line for horizontal sound wave
     ; propagation at the sound speed of 7 km/s.  The dashed parabola
     ; is the fundamental mode with omega = sqrt(g k_h).

     ; Comment on content: larger measurement extent and duration
     ; would resolve the (k_h,f) parabolas of the global solar
     ; oscillation around P = 5 min, left of the Lamb line.  Fourier
     ; resolution is set by measurement length, extent for k_h and
     ; duration for f. It took solar physicists from 1960 to 1975 to
     ; grasp this and resolve the "solar 5-minute" oscillation into
     ; global-mode parabolas.  Twenty years later long duration was a
     ; principal motivation to launch the SOHO mission.


  plotkophasediff,gb,ca,0.071,30,'ko-phasediff-gb-ca.ps',contours=1,$
    kmax=0.3,minphasediff=-180,maxphasediff=90,lamb=7,/fundamental
  $gv ko-phasediff-gb-ca.ps

     ; the companion "k-omega" phase-difference diagram.  
   
     ; Comment on content: the large negative phase blob at left below
     ; the Lamb line represents reversed granulation and internal
     ; gravity waves.  The role of the latter remains unclear.  The
     ; little white blobs along the k_h axis are probably phase
     ; wraparounds that result from the arctan [-pi, +pi] limits.
     ; Alfred de Wijn wrapped these back and accordingly extended the
     ; phase difference gray scale in Figure 7 of
  RJR{2004A&A...416..333R}.
