;
;
; Example code to read the NCEP data.
;
; Notice the data files are stripped down versions of the raw data.
;
; To run this, start idl (type "idl") then type ".run example"
;
;
; MAKE SURE YOU HAVE a ~/public_html directory 
;
; David Noone <dcn@colorado.edu> - Thu Feb 21 07:27:41 MST 2008
;
; ------------------------------------------------------------------------
; READ THE DATA
;
; Set the lattude of interest
;
lat0 = 45.

;
; Set the filenames for the T and v data
;
tfilename = '/home/atoc/dcn/ATOC5060/NCEP/air.2007.500.jan.nc'
vfilename = '/home/atoc/dcn/ATOC5060/NCEP/vwnd.2007.500.jan.nc'

;
; Read the temperature data
;
id = ncdf_open(tfilename,/nowrite)
ncdf_varget,id,ncdf_varid(id,'lon')  ,lons	
ncdf_varget,id,ncdf_varid(id,'lat')  ,lats
ncdf_varget,id,ncdf_varid(id,'level'),levs
ncdf_varget,id,ncdf_varid(id,'time') ,time

vid = ncdf_varid(id,'air')
ncdf_attget,id,vid,'scale_factor',scalef
ncdf_attget,id,vid,'add_offset'  ,offset
ncdf_varget,id,vid,ivals

tg = ivals*scalef + offset		; Temperature in K

;
; From the 3d data, extract only the latitude of interest
;
dum = min(abs(lats - lat0),jlat)
tg = reform(tg[*,jlat,0,*])

;
; Do the whole thing again to get data from another file, say, v wind!
;






;
; Save the number of longitudes  and times
;
nlon = n_elements(lons)
ntime = n_elements(time)

; ------------------------------------------------------------------------
; ANALYSIS
;
; YOUR ANALYSIS GOES HERE
;
;

sinx = sin(lons*!dtor)
funx = 0.25 + 0.25*sin(lons*!dtor) + 0.5*sin(3*lons*!dtor)

;
; Example of a power spectrum using FFT. 
; notice the variable "four" is complex, and can be broken into real and
; imaginary parts
;
four = fft(funx)

freal = real_part(four)
fimag = imaginary(four)

;
; Only used half the spectrum, since the second half is symmetric
; (and is the complex conjugate). Since only half, need to double the
; magnitude to retain all of the variance. The "zeroth" element of the
; array, is just the mean which is not symmetric, so doesn't neeed to be
; doubled.
;
amplitude = abs(four[0:nlon/2])
amplitude[1:nlon/2] = 2*amplitude[1:nlon/2]
wavenumber = findgen(nlon/2 + 1)

;
; Notice wave 1 and 3 have amplitudes of 0.5 - as we defined above for funx)!
;
; plot,wavenumber,amplitude,xtitle='wavenumber',ytitle='amplitude',thick=3


; ------------------------------------------------------------------------
; PLOTTING
;
; An example of how to make a graph
; Notice this example outputs a postscript file, and then converts this
; to a png file that you can view via the web
;

;
; initialze the plot
;
!p.multi=0
!p.font=1
set_plot,'PS'
device,xsize=7,ysize=4,xoff=0.5,yoff=0.5,/inches
device,/color,bits=8
device,set_font='helvetica'

;
; Plot the first graph (sinx)
;
plot,lons,sinx, title='My graph',$
   xrange=[0,360],xstyle=1,xtitle='longitude', $
   yrange=[-1,1],ystyle=1,ytitle='function value', $
   charsize=1.5,thick=3

;
; Overplot the second graph (fun x)
;
oplot,lons,funx, thick=1, linestyle=2


;
; Finalize the plot
;
device,/close
set_plot,'X'
spawn,'convert -crop 0x0 -antialias -density 96 idl.ps ~/public_html/plot.png'
print,'Finished: see http://atoc.colorado.edu/username/plot.png'



end