    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
      2008SoPh..251..533R.pdf.

    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
      2005AAp...441.1183D.pdf.

      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
      2001AAp...379.1052K.pdf.        ; A nice display of how the choice of normalization influences
         ; high-frequency powermaps is given in Figure 5 of
      2005AAp...430.1119D.pdf.

      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
      2001AAp...379.1052K.pdf.        ;  Alfred de Wijn wrote the codes constituting this program
         ;  for application to similar DOT data described in
      2004AAp...416..333R.pdf.        ; 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
      2004AAp...416..333R.pdf.