pro sst_eof
device,true_color=24,decompose=0,retain=2

sstfiles = file_search('sstdata/tmi.fusion.*.gz')
sstdata = fltarr(1440,720,n_elements(sstfiles)) ;x by y by time array
binarydata=bytarr(1440,720)
scale  =  0.15
offset = -3.0
time = fltarr(132)
for i=0,n_elements(sstfiles)-1 do begin
  openr,lun,sstfiles[i],/get_lun,/compress
  readu,lun,binarydata
  free_lun, lun
  sstdata[*,*,i] = float(binarydata)*scale+offset
  ;get Julian time values for plotting PC1 later
  time[i] = JulDay(1,1,fix(strmid(sstfiles[i],19,4)))+fix(strmid(sstfiles[i],24,3))-1
endfor
print, n_elements(sstfiles)
sst_subset = sstdata[400:1119,280:439,*] ;extract sub region

;rebin to 2-degree resolution - SVDC won't run with original resolution
nx = 90
ny = 20
sst_subset = rebin(sst_subset,nx,ny,132)

;get monthly means
sst_subset = reform(sst_subset,nx,ny,12,11) ;convert to an x by y by month by year array
monthly_means = rebin(rebin(sst_subset,nx,ny,12),nx,ny,12,11)

;create anomaly time series & compact time into one dimension
sst_anom = reform(sst_subset-monthly_means,nx,ny,132)

;weight anomalies
lat = findgen(ny)-19.
weights = rebin(sqrt(cos(transpose(lat)*!DTOR)),nx,ny,132)
sst_anom_wgt = sst_anom*weights

;compact spatial dimensions
sst_anom_wgt = reform(sst_anom_wgt,nx*ny,132)

help, sst_anom_wgt
;perform Singular Value Decomposition and get normalized PCs
SVDC, sst_anom_wgt, W, U, V, /DOUBLE
PC1 = REFORM(U(0,*)/STDDEV(U(0,*)))
PC2 = REFORM(U(1,*)/STDDEV(U(1,*)))
PC3 = REFORM(U(2,*)/STDDEV(U(2,*)))

EOF1 = fltarr(nx,ny)
FOR i = 0, nx-1 DO BEGIN
  for j=0,ny-1 do begin
    EOF1[i,j] = REGRESS(PC1,REFORM(sst_anom(i,j,*)))
  endfor
ENDFOR
	
!p.multi=[0,1,2,0,0]
loadct,0
plot, time, PC1, background=255, color=0, charsize=2, title='Normalized PC1', xtickunits=['Years']

dmax = max(abs(eof1))
dmin = -dmax
;pimage, EOF1, region=[-20,20,100,280], col_file='blueredanom.tbl', position=[0,0,1,0.5], /noerase, vmin=dmin, vmax=dmax, title='EOF 1',lsize=2, tsize=2, /coast

AB = fft(PC1)
C = sqrt(real_part(AB)^2+imaginary(AB)^2)
help, C
periods = (time[131]-time[0])/(findgen(66)+1)
plot, periods, C[0:65],/xlog,xtitle='Period (days)', title='Power Spectrum of PC1'
end