Friday, May 23, 2008: practice # 2

  cd ~/IDL
  rm -rf Demos
  rm -rf Data
  wget -r -nc -nH --cut-dirs=1 http://ross.iasfbo.inaf.it/IDL/Demos
  wget -r -nc -nH --cut-dirs=1 http://ross.iasfbo.inaf.it/IDL/Data
  wget http://ross.iasfbo.inaf.it/IDL/healpix201.tar.gz
  wget http://ross.iasfbo.inaf.it/IDL/jhuapl_sep07.tar.gz
  wget http://ross.iasfbo.inaf.it/IDL/solarpath.tar.gz
  tar zxvf healpix201.tar.gz
  tar zxvf jhuapl_sep07.tar.gz
  tar zxvf solarpath.tar.gz
  cd Demos

Reading a columns formatted ASCII file:

We assume the data file demodat.asc is in the dir ~/IDL/Data/

; Let's build the full path to the file 
  file = FILEPATH('demodat.asc',SUBDIR='IDL/Data',ROOT=GETENV('HOME'))
  a = READ_ASCII(file)

; Let's see what's in a (it is a "structure" so I use the STRUCTURE kwd)
  HELP, a, /STRU

; READ_ASCII is very flexible but you need to get familiar with it.
; Can get help on READ_ASCII using "doc_library"
  DOC_LIBRARY, 'read_ascii'

; Let's use some keywords:
;   we ask for data from line "23",
;   put the first 22 "header" lines into h
;   return the number of data lines read into n
  a = READ_ASCII(file, DATA_START=22, HEADER=h, COUNT=n)
  PRINT, 'There are '+ STRTRIM(n,1) +' data records.'
  PRINT, h                ; all the header as a single string
  PRINT, h, FORMAT='(A)'  ; formatted output, 1 line per read row
  PRINT, a.field1[*,0:3]  ; the first 3 rows (all columns) of "tag" field1

; To set the template structure of a file, use ASCII_TEMPLATE.
  tmpl = ASCII_TEMPLATE(file)

; Data start at line 23. Can then simply click "Next", "Next", "Finish" or
; change the column names as you like. Data types are what matters.

; See how the returned structure looks like
  HELP, tmpl, /STRU
  PRINT, tmpl.fieldnames

; Now we can read the file using this template and see what is changed
  a = READ_ASCII(file, TEMPLATE=tmpl, HEADER=h)
  HELP, a, /STRU


; You'll see 7 distinct tags corresponding to 7 FLOAT arrays with the columns
; data values!

; Let's print column 3
  PRINT, a.field3

For "simple" cases, e.g. want to map all the columns into a 2-d FLOAT or DOUBLE array, it would be handy to use my personal routine nxyreadf:
; I use the keyword "ncol" to force reading all the data. The def. is 2 columns
  a = nxyreadf(file,ncol=7)
  HELP, a

; Plot 2nd column vs first column
  PLOT, a[0,*], a[1,*]

; Plot 3rd column vs first column
  PLOT, a[0,*], a[2,*]
Note that ONLY numerical columns are allowed by nxyreadf and all the rows starting with a comment line or a non-numerical letter are skipped.
In addition to READ_ASCII, to read generic types column formatted file you can use other users' written routines, like Craig Markwardt's routine TRANSREAD.
Routines useful to read ASCII files are present in the Astronomy Lib. too: RDFLOAT, READCOL, READFMT (TEXTOPEN, FORPRINT, TEXTCLOSE can be used to print a set of vectors by looping over each index value.)
For those used to QDP produced "light curves" files, you can read them using the routine lcqdpread which makes use of nxyreadf. Here, as the abscissae are bins, data are plotted using the personal routine binplot.
; Read a SAX-PDS light curve of the Crab pulsar. Bins are "pulse phase".
  a = lcqdpread('crab_pdspulse.qdp', /phase)
  a[2,*] = a[2,*]/MAX(a[2,*])  ; normalise in the range 0 - 1
  plot, [0], /NODATA, /NOERASE, xsty=1, ysty=1, chars=1.5, $
        xrange=[-0.30,0.80], yrange=[0.5,1.05], $
        xtitle='!17Pulse Phase!6', ytitle='!17Normalized cts/s!6'
  binplot, a[0,*]-1, a[2,*], thick=2, /connect

; Read the spectra indices files. Data also have errors, so use the
; personal routine erplot to plot them.
  a = lcqdpread('crab_pdsphin.qdp', /phase)
  plot, [0], /NODATA, /NOERASE, xsty=1, ysty=1, chars=1.5, $
        xrange=[-0.10,0.50], yrange=[1.5,2.2], $
        xtitle='!17Pulse Phase!6', ytitle='!17Photon index'
  binplot, a[0,*], a[2,*], a[1,*], thick=2
  erplot, a[0,*], a[2,*], a[3,*]

; Run the demo program crabpulse which plots light curves and
; phase spectral indices for the four SAX instruments
  crabpulse

; To get a PS version simply add /PS
  crabpulse, /ps
Have a look to the source code (crabpulse.pro) using a text editor!
Other routines to plot data and/or errors include: PLOTERR, OPLOTERR, ERRPLOT

To plot error bars along the X axis you can adopt various approaches. Here we use binplot.
  x = FINDGEN(16)*10  &  ex = SQRT(x)  &  y=FINDGEN(16)*100  &  ey=SQRT(y)
  PLOT, x, y, PSYM=4, SYMSIZE=1.5
  erplot, x, y, ey  &  binplot, x, y, ex
If not done yet, check the Wed21 section on how to read and write FITS files.

An example of FITS image manipulation:

Here we display 2 images around the stars WR45 and Feige56 taken with an optical camera. Sky coordinates grid is overplotted. The images are smoothed using FFTs and re-displayed in shaded surface. As usual, images are in the "IDL/Data" directory.

  fits_read, '../Data/IMG122R126.fits.gz', image, h
  EXTAST, h, astrom
  HELP, astrom, /STRU

; Let's get the image size: we use SIZE: first element is Ndim, 2nd, 3rd Dims
  s = SIZE(image)
  PRINT, s
  nx = s[1]  &  ny = s[2]

; Let's make an empty plot area where to show the image.
; Note that we leave 2 extra pixels for the axis
  PLOT, [0], XRANGE=[0,nx+1], YRANGE=[0,ny+1], CHARS=1.3, xsty=5, ysty=5

; Now get the plot area in DEVICE coordinates (from DATA, the default)
  cv = CONVERT_COORD([0,nx], [0,ny], /TO_DEV)

; Byte scale using a reasonable Z range (can try other ranges)
  imin = MIN(image, MAX=imax)
  bima = BYTSCL(image, MIN=imin, MAX=imin*1.2)

; Display as 2-d pseudo image using color table 3
  loadct, 3
  TV, CONGRID(bima, cv[0,1]-cv[0,0],cv[1,1]-cv[1,0]), 1, 1, /DATA

; Now overplot the sky frame. Note the keywords used.
  skyframe, '../Data/IMG122R126.fits.gz', /dogrid, gridcol=255, CHARS=1.3, $
            XTITLE='RA', YTITLE='Dec', /noera
 
; The image is "square" but the plot area is not? Use sqplset.
; The second parameter is the CHARacter Size we'll use. The returned
; parameters are passed to PLOT via XMARGIN and YMARGIN
  sqplset, 1.3, xm, ym
  PLOT, [0], XRANGE=[0,nx+1],YRANGE=[0,ny+1], CHARS=1.3, xsty=5, ysty=5, $
            XMARGIN=xm, YMARGIN=ym

; Now just re-convert the DATA coordinates to DEVICE coordinates.
; Note how we use XMARGIN and YMARGIN also for "skyframe"
  cv = CONVERT_COORD([0,nx],[0,ny],/TO_DEV)
  TV, CONGRID(bima, cv[0,1]-cv[0,0],cv[1,1]-cv[1,0]), 1, 1, /DATA
  skyframe, '../Data/IMG122R126.fits.gz', /dogrid, gridcol=255, CHARS=1.3, $
            XTITLE='RA', YTITLE='Dec', /noera, XMARGIN=xm, YMARGIN=ym

; You'll probably see that the "Dec" label is not visible. Try to understand
; why and find a solution!


; Let's have a better view in a zoom window. Simply type:
  zoom

; Pay attention to the message! Use the correct mouse buttons!
; If you want a "continuous" zoom use the CONTINUOUS keyword. Type:

  zoom, /CONT

; Have a look to the other options with the help
  doc_library, 'zoom'

; Let's get the cordinates at a given position. First get position using CURSOR
  CURSOR, x, y, /DATA

; Then convert into sky coordinates (fractional dregrees) using the astrometry
; info in the "astrom" structure
  XY2AD, x, y, astrom, xd, yd
  PRINT, 'RA=', xd, '  Dec=', yd

; Then convert into string format
  coord_s = ADSTRING(xd, yd, 1)
  PRINT, 'Coords:', coord_s


; Smooth with a 5 pixels FWHM gaussian using FFT
  simage = filter_image(image, FWHM_GAUSSIAN=5)

; Rescale image and replot
  imin = MIN(simage, MAX=imax)
  bima = BYTSCL(simage, MIN=imin, MAX=imin*1.2)
  TV, CONGRID(bima, cv[0,1]-cv[0,0],cv[1,1]-cv[1,0]), 1, 1, /DATA
  skyframe, '../Data/IMG122R126.fits.gz', /dogrid, gridcol=255, CHARS=1.3, $
            XTITLE='RA', YTITLE='Dec', /noera, XMARGIN=xm, YMARGIN=ym 

; As a shaded surface (giving some ranges)
  SHADE_SURF, simage, XSTY=1, YSTY=1, XRANGE=[0,1024], YRANGE=[0,1024]
 
; The base of the main star
  SHADE_SURF, simage[441:530,461:530]

; Surface plot in "LEGO" format
  SURFACE, simage[441:530,461:530], /LEGO

; Colored mesh using the "Rainbow + white" table (exclude index 255 = white!)
  simage_c = simage[441:530,461:530]
  imin = MIN(simage_c, MAX=imax)
  shad = (simage_c - imin)/(imax - imin) * 254
  loadct, 39
  SURFACE, simage_c, SHADES=shad

; or

  SURFACE, simage_c, SHADES=shad, /LEGO
Let's do a similar work with a slitless optical image, but within a PROgram. Edit a file in the Demos dir. with your preferred text editor, for example:
  $ xemacs image_test.pro
Then type in the following code (cut-and-paste):
(if you are using idlde then copy the code into the editor window)
PRO image_test, filename

  IF N_PARAMS() NE 1 THEN $
    file = '../Data/IMG039A119.fits.gz' $
  ELSE $
    file = '../Data/' + filename

  fits_read, file, image, h
  EXTAST, h, astrom, ptype
  IF (ptype GT 0) THEN $
    PRINT, 'The image has valid astrometry keywords.' $
  ELSE BEGIN
    MESSAGE, 'The image has NOT valid astrometry keywords.'
  ENDELSE

  s = SIZE(image)
  nx = s[1]  &  ny = s[2]
  PRINT, 'This is a '+ STRTRIM(nx,2) +'x'+ STRTRIM(ny,2) +' image.'
  loadct, 5

  WINDOW, 1, XSIZE=700, YSIZE=600
  sqplset, 1.3, xm, ym
  PLOT, [0], XRANGE=[0,nx+1],YRANGE=[0,ny+1], CHARS=1.3, xsty=5, ysty=5, $
             XMARGIN=xm, YMARGIN=ym
  cv = CONVERT_COORD([0,nx],[0,ny],/TO_DEV)
  imin = MIN(image, MAX=imax)
  bima = BYTSCL(image, MIN=imin, MAX=imin*1.2)
  TV, CONGRID(bima, cv[0,1]-cv[0,0],cv[1,1]-cv[1,0]), 1, 1, /DATA
  skyframe, file, /dogrid, gridcol=255, CHARS=1.3, $
            XTITLE='RA', YTITLE='Dec', /noera, XMARGIN=xm, YMARGIN=ym
  

  simage = filter_image(image, FWHM_GAUSSIAN=5)
  imin = MIN(simage, MAX=imax)
  bima = BYTSCL(simage, MIN=imin, MAX=imin*1.2)

  WINDOW, 2, XSIZE=700, YSIZE=600
  TV, CONGRID(bima, cv[0,1]-cv[0,0],cv[1,1]-cv[1,0]), 1, 1, /DATA
  skyframe, file, /dogrid, gridcol=255, CHARS=1.3, $
            XTITLE='RA', YTITLE='Dec', /noera, XMARGIN=xm, YMARGIN=ym 

  WINDOW, 3, XSIZE=700, YSIZE=600
  SHADE_SURF, simage, XSTY=1, YSTY=1, XRANGE=[0,1024], YRANGE=[0,1024], $
              ZTITLE = 'WR45: Smoothed image'
 
END
; Then can run it
  image_test

; Note that you could first try to compile it
  .r image_test

; As you note, the file name can be given as parameter. Try with the previous
; image name (omit "../Data")

Aitoff (galactic) projection of celestial sources:

Se the aitoff routines used below.

; AGNs and unidentified EGRET (3EG) sources
  cieloagn


See cieloagn.pro source code (in the Demos dir.).

; Produce a PS file
  cieloagn, /PS
Have a look to the used functions and routines!

HEALPix maps:

Plot the astronomical catalogue UCAC-2 as objects density in orthographic and mollview all-sky maps + local orthographic view. Here the sky resolution is about 14 arcmin. The "ucac2_cover_k8.fits" file was produceded reading directly from a MySQL database.

This FITS file is NOT an image, but a "table" with a type format HEALPix specific.
Info about HEALPix at healpix.jpl.nasa.gov
  orthview, '../Data/ucac2_cover_k8.fits', rot=[45,45], /grat, col=5, $
            title='UCAC 2 - objects density'

  mollview, '../Data/ucac2_cover_k8.fits', coord=['C','G'], /grat, col=5, $
            title='UCAC 2 - objects density'

; 12 degrees around the galactic center in gnomic view (graticule 1 deg.):
  gnomview, '../Data/ucac2_cover_k8.fits', coord=['C','G'], res=10, grat=[6,6], $
            title='UCAC 2 - Galactic Center'

; To produce a PNG output file use the PNG keyword. For example:
  orthview, '../Data/ucac2_cover_k8.fits', rot=[45,45], /grat, col=5, $
            title='UCAC 2 - objects density', png='ucac2_density_orth.png'
Get help on mollview using: doc_library, 'mollview'

Graphics AstroTools:

The Astronomy Users' library and other libraries (e.g. the JHU/APL/S1R Johns Hopkins Applied Physics Laboratory S1R group library) have plenty of Astronomical routines and tools. Please check routines description (e.g. at idlastro.gsfc.nasa.gov and fermi.jhuapl.edu). Here we only show a few graphics tools.

The Javier Corripio's tool to compute and plot solar path for a given day, year, longitude and latitude. Source filea are in the Solarpath dir.

; Run the widget based user interface
  plotsun

; Input local latitude and longitude: for Bologna it is lat.= 44.5  lon.= 11.6

; The showed data is also written into an output file named "solpos.dat"
Position of the the Sun projected on the Earth map with day, twilight, and night zones. Uses routines in the Jhuapl directory.
; Simply run the program
  sunclock

; or in Orthographic projection
  sunclock, /ORTH

; or in Mercator projection
  sunclock, /MERC
MAP_SET keywords may be given to control the map display. Get help using: ? MAP_SET

Put some objects onto an equal area polar projection map.
; Make a square area in a 600x600 graphics window
  WINDOW, xs=600, ys=600
  sqplset, 1.0, xm, ym
  PLOT, [-130,130], [-130,130], XSTY=1, YSTY=1, XMARGIN=xm, YMARGIN=ym, /NODATA

; Make the grid, centered on the South pole
  eqpole_grid, /label, /south

; If you don't want to see the axes, use "XSTY=5, YSTY=5" in the PLOT command

; Longitude/Latidude or RA/Dec are in degrees
  glon = [15, 30, 45, 60]
  glat = [0, -20, -30, -60]

; Convert into cartesian coordinates for plotting.
; The keyword "south" tells we give -ve Lat
  eqpole, glon, glat, x, y, /south
  c = loadct19()
  plots, x, y, psym=2, syms=2, col=c[2]
Put some objects onto an Aitoff projection map (Galactic coordinates).
  WINDOW, xs=700, ys=400

; Let's assume we have 5 objects with Euqtorial coordinates (RA/Dec)
  ra = [2.56, 50.90, 80.14, 259.30, 340.26]
  dec = [73.17, 51.37, 25.75, -27.63, -67.60]

: Let's convert them into Galactic coordinates
  glactc, ra, dec, 2000., gl, gb, 1, /DEG
Get help using: doc_library, 'glactc'
; Convert into cartesian coordinates for plotting. Use a filled circle marker
  aitoff, gl, gb, x, y
  pp = FINDGEN(32)*(!pi*2/31.)
  USERSYM, COS(pp), SIN(pp), /FILL
  PLOT, x, y, XRANGE=[-180,180], YRANGE=[-90,90], XSTY=5, YSTY=5, $
        PSYM=8, COLOR=c[2]

; The Aitoff grid with marked corner values
  aitoff_grid2, 30, 15, label=1

Routines involving a Sphere. Use the spherical plot (sph*) routines in the Jhuapl directory. See the S1R pages.
; Plot the great circle between two points on Earth. From Bologna to La Silla
; Observatory (Chile)
  lat1 = 44.5
  lng1 = 11.6
  lat2 = -29.25
  lng2 = -70.73

  MAP_SET, /CONT                                  ; Plot Earth map
  ll2rb, lng1, lat1, lng2, lat2, r, b             ; Find r and b of pt2
  rb2ll, lng1, lat1, maken(0,r,100), b, lng, lat  ; Step along track
  PLOTS, lng, lat, col=c[3]
  PLOTS, lng1, lat1, psym=2, syms=2, col=c[2]     ; Plot point 1
  PLOTS, lng2, lat2, psym=2, syms=2, col=c[4]     ; Plot point 2
; To view the plot in Orthographic projection use:
  ERASE, 'ffffff'x
  MAP_SET, 0, -20, /ORTH, /HORIZON, /ISOTROPIC, /CONT
  PLOTS, lng, lat, col=c[3]
  PLOTS, lng1, lat1, psym=2, syms=2, col=c[2]
  PLOTS, lng2, lat2, psym=2, syms=2, col=c[4]

; Note the center of the projection is Lat= 0, Long= -20
; If you want a projection "without" the continents coastlines, simply
; omit the /CONT keyword
; Spirals around a sphere
  window, xs=500, ys=500
  psav = !P   ; save current !P structure
  !P.POSITION = [0,0,1,1]   ; the whole graph window

  set_isoxy, -1.1, 1.1, -1.1, 1.1
  LOADCT, 4
  ERASE, 80    ; use color index 80 for the background
  sphinit, long=-30, lat=40, rad=.4, fill=140, pa=30
  lng = makex(0,90,2)
  lat = lng*0
  rad = maken(0.9, 0.4, n_elements(lat))

  FOR a = 0, 360, 10 DO sphplot, lng+a, lat, rad, maxrad=.4
 
  !P = psav    ; restore !P structure

Satellite trajectory around a planet (Earth) and it's track.
; Simply run the program
  demo_sat
...

Plot 2-D data in 3-D space (to be continued...):

See spaceplot.pro source code (in the Demos dir.).
The relevant keywords are /SAVE given to SURFACE and /T3D used in all the PLOT, OPLOT, PLOTS commands. The command T3D allows you to change the plane in the space where to plot your data.
Note that here you could have problems viewing correctly the plot due to the transformation done by demo_sat. Perhaps it's better restarting IDL (this is temporary, I'll suggest what to do later).

; Simply execute the program
  spaceplot
 
; Use the personal widgets routine "setaspect" to change viewing angle
  setaspect

; Replot
  spaceplot
See help on T3D and the other 3-d routines: ? T3D or doc_library, 'T3D'

3-D data in 3-D space (to be continued...):

To plot a function of two variables (e.g., Z=f(X, Y)) inside a 3-D box can either make custom programs using commands allowing the T3D keyword, or use the IDL library procedure PLOT_3DBOX. Here is an example copied form the on-line help pages:

; Create some data to be plotted: 
 X = REPLICATE(5., 10.) 
 X1 = COS(FINDGEN(36)*10.*!DTOR)*2.+5. 
 X = [X, X1, X] 
 Y = FINDGEN(56) 
 Z = REPLICATE(5., 10) 
 Z1 = SIN(FINDGEN(36)*10.*!DTOR)*2.+5. 
 Z = [Z, Z1, Z] 
 
; Create the box plot with data projected on all of the walls. The  
; PSYM value of -4 plots the data as diamonds connected by lines: 
 PLOT_3DBOX, X, Y, Z, /XY_PLANE, /YZ_PLANE, /XZ_PLANE, $ 
    /SOLID_WALLS, GRIDSTYLE=1, XYSTYLE=3, XZSTYLE=4, $ 
    YZSTYLE=5, AZ=40, TITLE='Example Plot Box', $ 
    XTITLE='X Coordinate', YTITLE='Y Coodinate', $ 
    ZTITLE='Z Coordinate', SUBTITLE='Sub Title', $ 
    /YSTYLE, ZRANGE=[0,10], XRANGE=[0,10], $ 
    PSYM=-4, CHARSIZE=1.6 
See help on PLOT_3DBOX: ? PLOT_3DBOX or doc_library, 'PLOT_3DBOX'

Go to the course Main page L. Nicastro, INAF-IASF, 2009-01-22