;SOLAR PHYSICS COURSE 2008
;INTRODUCTION TO IDL




;Parts of this tutorial were taken from http://icc.dur.ac.uk/~tt/idl.html and                  
 ; http://nstx.pppl.gov/nstx/Software/IDL/idl_intro.html#INTRO)
 ; and from  http://dot.astro.uu.nl/rrweb/rjr-material/idl/manuals/idl-brief-manual.txt
 ;where you can also find the complete tutorial from Rob Rutten 
  

;  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;  +          Why IDL?                  +
;  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
  

;IDL is much used in astronomical image processing.
;It is an interactive computer language with the following advantages:
 ;- programming language, not a package: make up your own stuff;
 ;- interactive: test out new tricks by trying them statement by statement;
 ;- array notation: a = a + (b < 0) works for whole multi-dimensional arrays.
;The major disadvantage is that IDL licenses are very expensive.  The Student
;Version is cheaper but limited to array size 65536 (256x256 pixel images).  Manuals




;Manuals
;-------


;The online help is good.  
;The old student version manual ("Learning IDL")  is fairly good.
;The IDL astronomy library: http://idlastro.gsfc.nasa.gov/homepage.html
;Solar IDL software: http://www.lmsal.com/ssw
;David Fanning's Coyote's Guide: http://www.dfanning.com/
;Books: 
;  - David Fanning "IDL Programming Techniques, 2nd Edition"
;  - Lilian Gumley "Practical Idl Programming" 



;IDL as interpreter
;------------------
;  IDL works directly on the command line when you hit return.
;  This makes it easy to try new statements and statement sequences.
;  The up cursor arrow brings back earlier commands.
  
  
  
  
;  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 
;  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;  +         Exercises                   +
;  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
  
  
  
  
  
  
  ;START IDL:
  ;type “idl” into the command line
  ;you can start the HELP index by typing in “?” into the command line
  ;for help on a command for example "print" ,type "?print"
   
   
   
   
;  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;  +           floating point and integer arithmetic                     +
;  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

;
; floating point (note the "," to separate commands from values and keywords)
x = 2.
y = 3.
print,' x = ',x,' y = ',y,' x/y=', x/y

; idl does not does not require explicit variable
; declaration (unlike fortran or c). x=2. will assume
; that x is a floating point variable (single precision), since 2. is
; (note the dot)

; integer
x = 2
y = 3
print,' x = ',x,' y = ',y,' x/y=', x/y

;
; Unexpected? Since 2 and 3 are integers, so are x and y, hence 2/3 -> 0.
; Needles to say, this is the root of many an error!
; Try using the index to find which variable types exist, for example
; do ?data types
; you can also use strings
string1 = 'Hello '
string2 = 'World!'
print,string1+string2

; Functions
; Most standard arithmetic functions exist, for example
print,' sqrt(-1) = ',sqrt(complex(-1,0))
print,' atan(-1)/Pi = ',atan(-1)/!PI

; short integers: another source of errors.
; try this
i = 1
print,' Integers: i = ',i,' i+1 = ',i+1
; Good! But now try this:
i = 32767
print,' Integers: i = ',i,' i+1 = ',i+1
; BAD! By default, integers in idl are short, i.e. they
; have values ranging between -32768:32767
; The following gives you what you intended (note the ell in 32767l)
i =32767l
print,' Long integers: i = ',i,' i+1 = ',i+1
; I most often make an error involving sort integers when using i as
; the index in a loop, or, when addressing an element in a vector.




;  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;  +           Vectors and matrices                                      +
;  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;  A great advange of idl is that all functions work on vectors/arrays
;  as well.
; 
; Notice the [] brackets for defing a vector,and the () for a function
x = float([1,2,3,4,5])
print,' vector x = ',x
print,' x^2   = ',x^2
print,' 1/x =   ',1./x

; another commonly made mistake:
print,' x[1] = ',x[1]

; idl indexes like C: zero based. Thus
print,' x[0] = ',x[0]

; This difference with Fortran is an important source of error when
; tanslating idl programmes to F90

; some more nice stuff
print,' min(x) = ',min(x)
print,' max(x) = ',max(x)
print,' sum of all x = ',total(x)
print,' reverse(x) = ',reverse(x)
print,'   shift x by one ',shift(x,-1)
print,' statistics of x : ',moment(x)
; use ?moment to find which moments are computed

; you can enquire about a variable using help:
help,x

; 
y = replicate(2.,n_elements(x))
print,' x = ',x
print,' y = ',y
print,' y^x = ',y^x




;  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;  +           operations on vectors                                     +
;  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

; generate nran random numbers with a given seed.
; (uniformly spaced in ]0,1[ )
nran   = 5
myseed = - 1234
x = randomu(myseed,nran)
; setting the seed guarantees you get the same sequence everytime
; randomu(seed,nran) (without declaring seed) gives a different
; sequence everytime.
print,nran,' random numbers '
print,' x = ',x
index = sort(x)
xs    = x[index]
print,'   sorting x : ',xs
print,'   where is x > 0.5: ',where(x gt 0.5)
print,'       and its values are :',x[where(x gt 0.5)]


; the "where" statement is *very* powerful, I use it *very* often
; Note that you find the indices also this way
index_where_x_gt_onehalf = where(x gt 0.5)
print,' x = ',x
print,' x is gt 1/2 at indices :',where(x gt 0.5)
; what happens if no indices qualify?
index = where(x gt 1)
print,' index where x > 1 = ',index




;  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;  +           simple loops                                              +
;  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

;  idl has 
;    -- do loops
for i=0,2 do print,' i = ',i
;    -- while loops
i = 0
while (i lt 10) do begin print,' i = ',i & i = i + 1 
;    -- if/else/endif  
;    -- case             (see the index and try them)
;    -- repeat until
;    -- many others




;  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
;  +           matrices                                     +
;  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

;A matrix (a 2-dimensional array) may be defined algorithmically:

 A = dblarr(2, 4)
 for i = 0, 1 do begin $
 for j = 0, 3 do begin $
 a(i, j) = 10 * i + j  
 print, A

;Note that as it is printed, the first index corresponds to the column, and the second index to the row. Another way to think of it is that the way the ;data is stored, the first index varies fastest, and the second varies the slowest. This agrees with the way the data is printed.

;The WHERE function can be extremely useful
;Watch out, the results of the WHERE function 2-D arrays are confusing, because they return single-valued indices, but they work. (Continuing from the previous example:)

print, WHERE(A GT 10 AND A LT 13)

;A matrix may be constructed explicitly from vectors:

v1 = [1, 2, 0]
v2 = [1, 0, 0]
v3 = [4, 5, 6]
 A = [[v1], [v2], [v3]]
 print, A

;Create the transpose:

Atrans = transpose(A)
print, Atrans
;Take the determinant:

d = determ(float(A))
print, d

;Invert:

Ainv = invert(A)
print, Ainv

;Multiply vectors by matrices:

v = [1, 2, 3]
print, A
print, v
print, A ## v
print, v ## A

;You can solve a linear system Ax = b for x by Cramer's rule (the cramer function expects float or double inputs, requiring an explicit type conversion):
b = float([1, 2, 16])
A = float(A)
x = cramer(A, b)
print, x


; +++++++++++++++++++++++++++++++++
;  simple plotting           +
; +++++++++++++++++++++++++++++++++
  a=indgen(100)         ; direct definition of A=0,1,....,99
  print,a(0),a(99)      
  print,a(10:19)
  b=sqrt(float(a))
  for i=0,99 do if (b(i) GT 5) then a(i) = a(i) + b(i)
  plot,sin(a/10)        
  a=findgen(100)        ; float array A=0,1,....,99
  plot,sin(a/10)        
  b=sin(a/5)/exp(a/50)  ; how many elements has B?
  plot,b
  plot,b,cos(a/5)       
  oplot,b,cos(a/4)      ; over-plots in existing graph
  erase                 ; clean window
  x=b                   
  y=cos(a/5)
  plot,x,y,psym=1       ; defined for psym=1-7,10; try them out
                        ; something=something: optional "keyword" parameter
                        ; check PLOT in Online Help; check GRAPHICS KEYWORDS
  plot,x,y,linestyle=1  ; defined for linestyle=0-5
  plot,x,y,xtitle='x-axis',ytitle='y-axis',thick=2, $
      charsize=1.5            ; the $ extends command to next line
  plot,x,y,xrange=[-2,+2],yrange=[-1.5,1.5]   ; your choice of axes
  plot,shift(abs(fft(b,1)),25),/ylog   
                        ; power spectrum on linear-log scales
                        ; /ylog is alternative for ylog=1 (true)
  xyouts,10,30,'Power Spectrum'            ; annotation, x,y in data units
  xyouts,.75,.75,'Power Spectrum',/normal  ; annotation, axis-length units

; PostScript output
  set_plot,'ps'         ; now write to PostScript file instead of window
  device,filename='power-spectrum.ps'    ; choose filename
    ; all IDL output commands now write to that file
  plot,x,y,xtitle='x-axis',ytitle='y-axis',thick=2,charsize=1.5  
  device,/close         ; done
  set_plot,'win'        ; back to output on Windows95 screen
  set_plot,'x'          ; idem for Unix
  help,/device          ; ,/device is the same as using ,device=1 (enable)


; +++++++++++++++++++++++++++++++++
;    More advanced plotting           +
; +++++++++++++++++++++++++++++++++
;More advanced plotting

; more displays  
  z=b#a                 ; z(i,j)=b(i)*a(j); full matrix multiplication: ##
  help,z                ; so this is a 2-dim (100,100) float array 
  print,z[0:5,0:5]      ; both square brackets and parentheses permitted
  plot,indgen(50),z[2,*]    ; plot 3rd column of array z
  plot,indgen(50),z(*,49)   ; plot bottom row of array z 
  contour,z             ; just for fun
  surface,z             ; just for fun
  shade_surf,z          ; just for fun - but only in a small window
  window,1,xsize=200,ysize=200   ; open small window, size in pixels
  shade_surf,z                  
  wdelete,1             ; delete window 1
  show3,z               ; just for fun
  
; image display
  tv,z                  ; show array as image
  zsize=SIZE(z)         ; get array type and size
  nx=2*zsize(1)
  ny=2*zsize(2)         ; etcetera for more dimensions; zsize(0)=type
  z=rebin(z,nx,ny)      ; resample z to make its display bigger
  wdelete
  window,xsize=nx,ysize=ny      ; window equal to image size
  z8=bytscl(z)          ; scales dynamic range to 8 bits (shades 0 - 255)
  print,min(z),max(z),min(z8),max(z8)    ; list extrema
  tv,z8                 ; display image with maximum contrast on screen
  tvscl,z               ; same as tv,bytescl(z)
  profiles,z            ; measurent tool for showing image values
                        ; stop it with right mouse button (cursor on image)
  loadct                ; set colour table
  adjct                 ; tool to adjust greyscale table
  xpalette              ; tool to adjust color table
  xloadct               ; idem
  tvscl,z8>127          ; display only the brighter half of the pixels
  erase
  tvscl,z(0:nx/2,0:ny/2)     ; bottom-left quarter of the image
  wdelete
  tvscl,congrid(z,188,188,/interp) ; arbitrary resizing with interpolation
  
; +++++++++++++++++++++++++++++++++++++++++
;   Input/Output 
; +++++++++++++++++++++++++++++++++++++++++  

; read/write files
  openw,1,'myfile.ext'  ; open file myfile.ext on logical unit 1 for writing
  printf,1,z            ; write free-format file
  close,1               ; free lun 1
  openr,1,'myfile.ext'  ; now open that file for reading
  zz=z                  ; define variable type and size
  readf,1,zz            ; read free-format file
  close,1   
;  (use: - WRITEU, READU to write/read large files unformatted = faster;
;        - image=assoc(lun,bytarr(nx,ny)) to get single images from a movie)

  

; +++++++++++++++++++++++++++++++++++++++++
;    a simple contour plot of a Gaussian  +
; +++++++++++++++++++++++++++++++++++++++++

npts = 256 ; number of pts
nsig = 4.  ; number of sigma
x    = -nsig + findgen(npts)/float(npts)*2*nsig
Gs   = fltarr(npts,npts)
for j=0,npts-1 do begin $
  & for i=0,npts-1 do begin $
  & Gs[i,j] = exp(-(x[i]^2+x[j]^2)/2.) $
  & endfor $
  & endfor

contour,Gs,x,x,nlev=255,/fill

; now try
device,decomposed=0
loadct,11
contour,Gs,x,x,nlev=255,/fill




; +++++++++++++++++++++++++++++++++++++++++
;    Writing a program... +
; +++++++++++++++++++++++++++++++++++++++++


; We want to write a program that asks the user to enter his name and year of birth.
; The program should then calculate the age of the user and print the name + age to the screen.

;Open a text editor 
; the program outline has to be

pro age
IDL STATEMENT
IDL STATEMENT
...
end

;write the program and write a comment for each step. Comments are written by typing ";" at the beginning of the line
;TIP: Use the command "READ" in order to get the name and year of birth, and use SYSTIME to get the current date. In order to extract the Year from the string with the current date you will need to use the command STRMID (see example for SYSTIME in the help index) 

;save the file as age.pro in your working directory
;Always give your filename the same name as the program name used in the first line of the file when writing IDL programs. It is also important to add the suffix .pro to the end of the filename. This is the clue that IDL uses to pick out programs it knows how to run.

;Run the program at your idl prompt:
age

;(after changing the program you should also first type ".run age " in order to ensure, that the latest version of the program has been  compiled)


;Find out about the difference between a Procedure and  Function.
;Can you pass on arguments to a procedure?

; Write a procedure with argument radius "R" and height "H" that calculates the volume of a sphere using a function
; spherevol and the volume of a cylinder using the function cylvol, compares the volumes and prints the higher one onto the
; screen.


; +++++++++++++++++++++++++++++++++++++++++
;    before you start the Image Processing practicum ... +
; +++++++++++++++++++++++++++++++++++++++++

;it is very useful to use the "JOURNAL" function in order to save your idl sessions onto disk! Check the idl HELP for information on this function.
;you can download the necessary data files for the exercises through the link on the course outline page. The data files are stored as ".sav" files.
;in order to load for example a file called filename.save into your session type at the idl command line:

restore,'filename.save',/verbose
;;You will see the names of the restored variables on the screen, to find out more about the variables type 
help,name_of_variable
;to see all your variables loaded in your session type
help



; +++++++++++++++++++++++++++++++++++++++++
;   END +
; +++++++++++++++++++++++++++++++++++++++++