; what is pi?
print, !pi


; creates a 4x4 floating-point (f-p) array with values of 0 everywhere
array1 = fltarr(4,4)
help, array1
print, array1


; create a 4x4 double-precision array with values of 0 everywhere
array2 = dblarr(4,4)
help, array2  & print, array2


; create a 4x4 f-p array with values equal to its 1-D subscript (index)
array3 = findgen(4,4)
print, array3


; create two 4x4 f-p array with values equal to array3's first and second 2-D subscripts 
x = fix(array3) mod 4           ; reminder
y = fix(array3) / 4
print,x
print,y


; array operation
array4 = exp(sqrt((array3 + array2) * 4.0) - 1.0)


; create a 200-element sin-wave array with period of 50 and amplitude of 1
x = findgen(200)
y = sin(x/50.0 * 2 * !pi)
plot, x, y, psym=1, xtitle='x', ytitle='sin(x)', yrange=[-2,2]
oplot, x, y

; an equivalent and stupid way (but sometimes necessary) of doing the above
x = fltarr(200)
y = fltarr(200)
FOR i=0,199 DO x[i] = i
FOR i=0,199 DO y[i] = sin(x[i]/50.0 * 2 * !pi) 


; make the above plot into a ps file
set_plot, 'ps'
device, filename = 'sin.ps'   
plot, x, y, psym=1, xtitle='x', ytitle='sin(x)', yrange=[-2,2], $
   position = [0.2,0.2,0.85,0.9], xcharsize=2, ycharsize=2, $
   xthick=3, ythick=3, charthick=2
oplot, x, y, thick=2
device, /close
set_plot, 'x'
$ gv sin.ps &

; print the above x,y values to an ASCII file
openw, out, 'sin.dat', /get_lun
FOR i=0,199 DO printf, out, x[i], y[i]
free_lun, out
$more sin.dat
; Do you know how to read such an ASCII file into idl?










; read a FITS file
image = mrdfits('GOODS.submm1.fits', 0, header)

; how does the header look like?
print, header

; what is the image size?
imsize = size(image, /dimen)
print,imsize

; what is the rms of the pixels?
rms = stddev(image)
print, rms

; what is the median and mean values of the pixels?
print, median(image), mean(image)

; what are the maximum and minimum value of the pixels?  Where are they?
print, max(image, index_max),  min(image, index_min)
print, index_max mod imsize[0], index_max / imsize[0]
print, index_min mod imsize[0], index_min / imsize[0]

; where are the negative pixels?
A = where(image LT 0.0)
help,A
print,A[0:9] ; the first 10 of the negative pixels

; what is the brightness of the galaxy at x=111, y=49?
brightness = total(image[111-10:111+10, 49-10:49+10])
print,brightness

; what are the best-estimates of the mean background value 
; and background rms in the image?
bg_mean = median(image)
bg_rms = stddev(image)
; select pixels within +/-4 sigma about the mean
A = where(abs(image - bg_mean) LT 3.0*bg_rms)   
; iterate again
bg_mean = median(image[A])
bg_rms = stddev(image[A])
print,bg_mean,bg_rms
; iterate again
A = where(abs(image - bg_mean) LT 3.0*bg_rms)   
; iterate again
bg_mean = median(image[A])
bg_rms = stddev(image[A])
print,bg_mean,bg_rms

; Let's see what pixels are selected in the above 3-sigma clipping
image2 = intarr(imsize)
image2[A] = 1
mwrfits, image2, 'bg.fits', header, /create


; a procedure in the IDL astronomy library can do the above even better.
sky, image






; Exercise:
;
; Answer the following questions with IDL.  Keep all the steps in a 
; script file like this one.  Feel free to ask me any time if you 
; are not sure what to do.
;
; 1. With the above sky background estimate, can you recalculate the
;    brightness of the galaxy at x=111, y=49 and take the background
;    brightness into account?
; 2. What is the associated measurement error in the above brightness 
;    measurement?
; 3. The above background measurement is for the entire image.
;    Do you think the background is the same everywhere in the image?
;    If not, can you account for this in the brightness measurement?
; 4. In the above measurement, we used an x-y size of 21*21 pixels.
;    How will the result change with different sizes?  Can you show 
;    this with a plot?  What is the optimal size?  (What does "optimal" 
;    mean?)
; 5. I picked the position of x=111, y=49 very casually with mouse 
;    cursor on the image.  Can you measure the position of this galaxy
;    better?  Try to think about several different quantitative 
;    methods (2 or 3) and compare the results.
; 6. Use the IDL Astronomy Library (not ds9) to find the RA and Dec of 
;    this galaxy with the x/y position you measured in 5.  How far away 
;    is it from the brightest galaxy, in arcsec?