program='pcomp_hubble'

;read Hubble (196) data, from WJ
  close,1 & openr,1,'4point1.dat'
  d=fltarr(24)
  v=fltarr(24)
  for i=0,23 do begin
      readf,1,dd,vv
      d(i)=dd & v(i)=vv
  endfor


;m variables, n observations
;Remove the mean and stdev from each variable.
  m=2
  n=24
  array=fltarr(m,n)
  array(0,*)=(d-mean(d))/stdev(d)
  array(1,*)=(v-mean(v))/stdev(v)

;RESULTS returns the new derived variable values for each data point  
;principal components are returned via the coefficients a_ij

  result = PCOMP(array, COEFFICIENTS = coefficients, $
               EIGENVALUES=eigenvalues, VARIANCES=variances,/covariance)

  print,'---------------------------------------------------------'
  print,'Compute principal component variables via PCA (pcomp.pro)
  print,'---------------------------------------------------------'
  print,'pcomp.pro:'
  print,'Coefficients a: (sum of squares = eigenvalue)'
  for i=0,m-1 do begin
      print, i+1, coefficients[*,i], format='("PC",I1,4(F10.4))'
  endfor
  eigenvectors = coefficients/REBIN(eigenvalues, m, m)
  print
  print, 'Eigenvectors (unit vectors): '
  for i=0,m-1 do begin
      print, i+1, eigenvectors[*,i],format='("PC",I1,4(F10.4))'
  endfor
  
;reconstruct original x's from derived variables
  array_reconstruct = result ## eigenvectors
  print
  print,'Reconstruction error: total', total((array_reconstruct - array)^2)
  print
  print,'Energy conservation: total variance of original and new variables'
  print, total(array^2),total(eigenvalues)*(n-1)
  print
  print, '    Mode    Eigenvalue    PercentVariance'
  for i=0,m-1 do begin
      print, i+1, eigenvalues[i], variances[i]*100
  endfor
  print,'---------------------------------------------'

END