pro pca2aref,sqlambf,pcompf,avgarff,defarff,areffil,$
	vthresh=vthresh,ncomp=ncomp,verbose=verbose, _extra=e
;+
;procedure	pca2aref
;	reads in the PCA decomposition of the calibration library ARFs
;	and writes it out to a FITS file in the AREF format
;
;syntax
;	pca2aref,sqlambf,pcompf,avgarff,defarff,areffil,vthresh=vthresh,ncomp=ncomp
;
;parameters [all required]
;	sqlambf	[INPUT] full path and name of sqrt(lambda) file
;	pcompf	[INPUT] full path and name of principal components file
;	avgarff	[INPUT] full path and name of average ARFs file
;	defarff	[INPUT] full path and name of default ARF used in generation of cal library ARFs
;	areffil	[INPUT] full path and name of output AREF file
;
;keywords
;	vthresh	[INPUT] how much of the total variance should be included in the AREF?
;		* if 0 or not given, uses NCOMP, or if that is not given, then is set to 0.99 
;		* if <0, abs() value is used
;		* if <1, assumes fractional variance
;		* if >1, assumes actual variance value (do not use unless you know what you are doing)
;	ncomp	[INPUT] number of components to include, starting from one that accounts for most variance
;		* default is number of components or 20, whichever is smaller
;		* if -ve, reverse order and exclude, starting from one that accounts for least variance
;		* in all cases, keeps at least one component
;	verbose	[INPUT] control chatter
;	_extra	[JUNK] here only to prevent crashing the program
;
;history
;	vinay kashyap (2014apr)
;-

;	usage
ok='ok' & np=n_params()
nsf=n_elements(sqlambf) & npf=n_elements(pcompf) & ndf=n_elements(defarff) & naf=n_elements(avgarff) & nof=n_elements(areffil)
ssf=size(sqlambf,/type) & spf=size(pcompf,/type) & sdf=size(defarff,/type) & saf=size(avgarff,/type) & sof=size(areffil,/type)
;
if np lt 5 then ok='Insufficient parameters' else $
 if nsf eq 0 then ok='SQLAMBF: not defined' else $
  if nsf gt 1 then ok='SQLAMBF: must be a scalar' else $
   if ssf ne 7 then ok='SQLAMBF must be a string' else $
    if npf eq 0 then ok='PCOMPF: not defined' else $
     if npf gt 1 then ok='PCOMPF: must be a scalar' else $
      if spf ne 7 then ok='PCOMPF must be a string' else $
       if ndf eq 0 then ok='DEFARFF: not defined' else $
        if ndf gt 1 then ok='DEFARFF: must be a scalar' else $
         if sdf ne 7 then ok='DEFARFF must be a string' else $
          if naf eq 0 then ok='AVGARFF: not defined' else $
           if naf gt 1 then ok='AVGARFF: must be a scalar' else $
            if saf ne 7 then ok='AVGARFF must be a string' else $
             if nof eq 0 then ok='AREFFIL: not defined' else $
              if nof gt 1 then ok='AREFFIL: must be a scalar' else $
               if sof ne 7 then ok='AREFFIL must be a string'
if ok ne 'ok' then begin
  print,'Usage: pca2aref,sqlambf,pcompf,avgarff,defarff,areffil,vthresh=vthresh,ncomp=ncomp'
  print,'  convert output of run_pca_on_arflib.r to AREF format FITS file'
  if np ne 0 then message,ok,/informational
  return
endif

; VERBOSE
vv=0L & if keyword_set(verbose) then vv=long(verbose[0])>1L

fsf=file_search(sqlambf[0],count=nfsf)
fpf=file_search(pcompf[0],count=nfpf)
fdf=file_search(defarff[0],count=nfdf)
faf=file_search(avgarff[0],count=nfaf)
if nfsf eq 0 then ok=SQLAMBF+': file not found' else $
 if nfpf eq 0 then ok=PCOMPF+': file not found' else $
  if nfdf eq 0 then ok=DEFARFF+': file not found' else $
   if nfaf eq 0 then ok=AVGARFF+': file not found'
if ok ne 'ok' then begin
  message,ok,/informational
  return
endif
fof=file_search(areffil[0],count=nfof)
if nfof ne 0 then begin
  message,fof[0]+': file exists, and will be overwritten',/informational
  wait,5
endif
outfile=areffil[0]

;	read in the inputs
; SQLAMB
spawn,'wc -l '+sqlambf[0],tmp & kcomp=long(tmp[0])
sqlamb=fltarr(kcomp) & openr,uu,sqlambf[0],/get_lun & readf,uu,sqlamb & close,uu & free_lun,uu
; PCOMP
spawn,'wc -w '+pcompf[0],tmp & ncol=long(tmp[0])/kcomp
pcomp=fltarr(ncol,kcomp) & openr,uu,pcompf[0],/get_lun & readf,uu,pcomp & close,uu & free_lun,uu
; AVGARF
avgarf=fltarr(ncol) & openr,uu,avgarff[0],/get_lun & readf,uu,avgarf & close,uu & free_lun,uu
; DEFARFF
defarf0=mrdfits(defarff,0,harf0) & defarf=mrdfits(defarff,1,harf)

;	figure out the columns and define the bias vector
ok='ok'
tagnames=tag_names(defarf) & ntags=n_elements(tagnames)
oarf=where(tagnames eq 'SPECRESP',moarf)
oelo=where(tagnames eq 'ENERG_LO',moelo)
oehi=where(tagnames eq 'ENERG_HI',moehi)
if moelo eq 0 then ok=arfeff[0]+': does not contain ENERG_LO column in extension '+sxpar(harf1,'EXTNAME') else $
 if moehi eq 0 then ok=arfeff[0]+': does not contain ENERG_HI column in extension '+sxpar(harf1,'EXTNAME')
  if moarf eq 0 then ok=arfeff[0]+': does not contain SPECRESP column in extension '+sxpar(harf1,'EXTNAME')
if ok ne 'ok' then begin
  message,ok,/informational
  message,'Quitting...',/informational
  return
endif
;
elo=defarf.ENERG_LO & ehi=defarf.ENERG_HI & arf0=defarf.SPECRESP
owlo=where(tagnames eq 'BIN_LO',mowlo)
owhi=where(tagnames eq 'BIN_HI',mowhi)
if mowlo eq 0 then wlo=12.39851/ehi else wlo=defarf.BIN_LO
if mowhi eq 0 then whi=12.39851/elo else whi=defarf.BIN_HI
nbin=n_elements(arf0)	;this should be the same as NCOL
if nbin ne ncol then begin
  message,'there is a mismatch between '+defarff[0]+' and '+pcompf[0],/informational
  message,'Quitting...',/informational
  return
endif
;
;	are there any extra columns?
ok=intarr(ntags)-1
if moarf gt 0 then ok[oarf[0]]=oarf[0]
if moelo gt 0 then ok[oelo[0]]=oelo[0]
if moehi gt 0 then ok[oehi[0]]=oehi[0]
if mowlo gt 0 then ok[owlo[0]]=owlo[0]
if mowhi gt 0 then ok[owhi[0]]=owhi[0]
oy=where(ok lt 0,moy) & if moy gt 0 then xtags=tagnames[oy]

;	figure out how many of the components to keep
pcafrac=sqlamb^2/total(sqlamb^2)
pcacvar=total(pcafrac,/cumul)
pcacfrac=pcacvar/max(pcacvar)
;
; NCOMP
mcomp=20
if keyword_set(ncomp) then begin
  ii=long(ncomp[0])
  if ii lt 0 then mcomp=(kcomp-ii)>1	;if -ve, exclude last NCOMP
  if ii gt 0 then mcomp=ii<kcomp	;if +ve, keep first NCOMP
endif else mcomp=kcomp<20		;smaller of all or 20
; VTHRESH overrides NCOMP
vthr=0.99
if keyword_set(vthresh) then begin
  if vthresh lt 0 then vthr=abs(double(vthresh[0])) else vthr=double(vthresh[0])
endif
if keyword_set(vthresh) or not keyword_set(ncomp) then begin
  if vthr ge 1 then begin
    ok=where(pcacvar le vthr,mok) & mcomp=mok
  endif else begin
    if vthr gt 0 then begin
      ok=where(pcacfrac le vthr,mok) & mcomp=mok
    endif
  endelse
endif
v_THRESH=pcacfrac[mcomp-1L] & v_NCOMP=mcomp & v_COMPRESS=float(mcomp)/float(kcomp)

;	set up the residuals array
pcaresid=0.*avgarf
for j=mcomp,kcomp-1L do pcaresid[*]=pcaresid[*]+sqlamb[j]*pcomp[*,j]	;components j=0,mcomp-1 are kept separate
residvar=0. & if mcomp lt kcomp-1 then residvar=total(pcacvar[mcomp:*])/max(pcacvar)

;	define some mission-specific header keywords
cd,current=cwd
if not keyword_set(v_CREATOR) then v_CREATOR=filepath('pca2aref.pro',root_dir=cwd)
if not keyword_set(v_MISSION) then v_MISSION=sxpar(harf,'MISSION') & if not keyword_set(v_MISSION) then v_MISSION='NONE'
if not keyword_set(v_TELESCOP) then v_TELESCOP=sxpar(harf,'TELESCOP') & if not keyword_set(v_TELESCOP) then v_TELESCOP='NONE'
if not keyword_set(v_INSTRUME) then v_INSTRUME=sxpar(harf,'INSTRUME') & if not keyword_set(v_INSTRUME) then v_INSTRUME='NONE'
if not keyword_set(v_FILTER) then v_FILTER=sxpar(harf,'FILTER') & if not keyword_set(v_FILTER) then v_FILTER='NONE'
if not keyword_set(v_DETNAM) then v_DETNAM=sxpar(harf,'DETNAM') & if not keyword_set(v_DETNAM) then v_DETNAM='NONE'
if not keyword_set(v_GRATING) then v_GRATING=sxpar(harf,'GRATING') & if not keyword_set(v_GRATING) then v_GRATING='NONE'

;	generate generic header for primary block
;hdr0=harf0
fxhmake,hdr0,/extend,/date,/INITIALIZE
;
;	add some useful stuff to the primary header
fxaddpar,hdr0,'LONGSTRN','OGIP 1.0','The HEASARC Long String Convention may be used'
fxaddpar,hdr0,'ORIGIN','CfA','Source of FITS file'
fxaddpar,hdr0,'CREATOR',v_CREATOR,'Program creating this file'
fxaddpar,hdr0,'REVISION',1
fxaddpar,hdr0,'MISSION',v_MISSION,'Mission'
fxaddpar,hdr0,'TELESCOP',v_TELESCOP,'Telescope Used'
fxaddpar,hdr0,'INSTRUME',v_INSTRUME,'Instrument'
fxaddpar,hdr0,'FILTER',v_FILTER,'the instrument filter in use (if any)'
fxaddpar,hdr0,'DETNAM',v_DETNAM,'Detector'
fxaddpar,hdr0,'GRATING',v_GRATING,'Grating'
;
;	and some new header keywords to tell us what we will be looking at
fxaddpar,hdr0,'EMETHOD','PCA1DADD','additive Principal components of 1D curves are stored'
fxaddpar,hdr0,'CMETHOD','ARRAY=NOMINAL+BIAS+sum{r_j*PCA_EVALUE_j*PCA_EVECOR_j}+r*PCA_RESID'
fxaddpar,hdr0,'BIASEXT','SPECRESP','nominal function and bias will be stored in this extension'
fxaddpar,hdr0,'ERREXT','PCACOMP','uncertainties will be stored in this extension'
fxaddpar,hdr0,'VARTHR',v_THRESH,'accounts for this fraction of the total variance'
fxaddpar,hdr0,'NUMCOMP',v_NCOMP,'number of principal components stored'
fxaddpar,hdr0,'COMPRESS',v_COMPRESS,'compression factor'
;
;	write it out
fxwrite,outfile,hdr0

;	generate generic header for BIASEXT extension
BIASEXT='SPECRESP'
fxbhmake,hdr0e1,nbin,BIASEXT,'the name of the extension',/date,/INITIALIZE
fxaddpar,hdr0e1,'XTENSION','BINTABLE','binary table extension'
fxaddpar,hdr0e1,'BITPIX',8,'8-bit bytes'
fxaddpar,hdr0e1,'NAXIS',2,'2-dimensional binary table'
fxaddpar,hdr0e1,'MISSION',v_MISSION,'Mission'
fxaddpar,hdr0e1,'TELESCOP',v_TELESCOP,'Telescope Used'
fxaddpar,hdr0e1,'INSTRUME',v_INSTRUME,'Instrument'
fxaddpar,hdr0e1,'FILTER',v_FILTER,'the instrument filter in use (if any)'
fxaddpar,hdr0e1,'DETNAM',v_DETNAM,'Detector'
fxaddpar,hdr0e1,'GRATING',v_GRATING,'Grating'
fxaddpar,hdr0e1,'HDUCLASS','CXC','file format is not OGIP standard.'
fxaddpar,hdr0e1,'HDUCLAS1','RESPONSE','extension contains response data.'
fxaddpar,hdr0e1,'HDUCLAS2','SPECRESP','extension is like HEASARC ARF standard'
fxaddpar,hdr0e1,'HDUCLAS3','BIAS','extension contains an offset to the ARF'
fxaddpar,hdr0e1,'HDUNAME','SPECRESP','Block name'
;
;	add some extra stuff
fxaddpar,hdr0e1,'LONGSTRN','OGIP 1.0','The HEASARC Long String Convention may be used'
fxaddpar,hdr0e1,'ORIGIN','CfA','Source of FITS file'
fxaddpar,hdr0e1,'CREATOR',v_CREATOR,'Program creating this file'
fxaddpar,hdr0e1,'REVISION',1
fxaddpar,hdr0e1,'RESPFILE','none'
;
;	and repeat the AREF keywords
fxaddpar,hdr0e1,'EMETHOD','PCA1DADD','Principal components of 1D curves are stored'
fxaddpar,hdr0e1,'CMETHOD','ARRAY=NOMINAL+BIAS+sum{r_j*PCA_EVALUE_j*PCA_EVECOR_j}+r*PCA_RESID'
fxaddpar,hdr0e1,'ERREXT','PCACOMP','uncertainties will be stored in this extension'
fxaddpar,hdr0e1,'VARTHR',v_THRESH,'accounts for this fraction of the total variance'
fxaddpar,hdr0e1,'NUMCOMP',v_NCOMP,'number of principal components stored'
fxaddpar,hdr0e1,'COMPRESS',v_COMPRESS,'compression factor'
;
;	define the columns and the appropriate header keywords
fxbaddcol,i1,hdr0e1,elo[0],'ENERG_LO','Energy bin lower bound',TUNIT='keV'
fxaddpar,hdr0e1,'TFORM1','1E','format of field'
fxaddpar,hdr0e1,'TLMIN1',min(ELO)
fxaddpar,hdr0e1,'TLMAX1',max(EHI)
fxbaddcol,i2,hdr0e1,ehi[0],'ENERG_HI','Energy bin upper bound',TUNIT='keV'
fxaddpar,hdr0e1,'TFORM2','1E','format of field'
fxaddpar,hdr0e1,'TLMIN2',min(ELO)
fxaddpar,hdr0e1,'TLMAX2',max(EHI)
fxbaddcol,i3,hdr0e1,arf0[0],'SPECRESP','effective area',TUNIT='cm**2'
fxaddpar,hdr0e1,'TFORM3','1E','format of field'
fxaddpar,hdr0e1,'TLMIN3',min(arf0)
fxaddpar,hdr0e1,'TLMAX3',max(arf0)
fxbaddcol,i4,hdr0e1,avgarf[0],'BIAS','offset from nominal',TUNIT='cm**2'
fxaddpar,hdr0e1,'TFORM4','1E','format of field'
fxaddpar,hdr0e1,'TLMIN4',min(avgarf)
fxaddpar,hdr0e1,'TLMAX4',max(avgarf)
fxbaddcol,i5,hdr0e1,wlo[0],'BIN_LO','Wavelength bin lower bound',TUNIT='Angstrom'
fxaddpar,hdr0e1,'TFORM5','1E','format of field'
fxaddpar,hdr0e1,'TLMIN5',min(WLO)
fxaddpar,hdr0e1,'TLMAX5',max(WHI)
fxbaddcol,i6,hdr0e1,whi[0],'BIN_HI','Wavelength bin upper bound',TUNIT='Angstrom'
fxaddpar,hdr0e1,'TFORM6','1E','format of field'
fxaddpar,hdr0e1,'TLMIN6',min(WLO)
fxaddpar,hdr0e1,'TLMAX6',max(WHI)
; and if there are any extra columns..
kk=6L
for ik=0L,moy-1L do begin
  kk=kk+1L
  skk=strtrim(kk,2)
  cc="fxbaddcol,i"+skk+",hdr0e1,min(defarf.(oy[ik])),tagnames[oy[ik]],'extra column "+strtrim(ik+1,2)+"',TUNIT='None'"
  print,cc
  jnk=execute(cc)
  fxaddpar,hdr0e1,'TFORM'+skk,'1E','format of field'
  fxaddpar,hdr0e1,'TLMIN'+skk,min(defarf.(oy[ik]))
  fxaddpar,hdr0e1,'TLMAX'+skk,max(defarf.(oy[ik]))
endfor
;
;	and write out the extension
fxbcreate,uf,outfile,hdr0e1
for i=0L,nbin-1L do begin
  fxbwrite,uf,float(elo[i]),1,i+1L
  fxbwrite,uf,float(ehi[i]),2,i+1L
  fxbwrite,uf,float(arf0[i]),3,i+1L
  fxbwrite,uf,float(avgarf[i]-arf0[i]),4,i+1L
  fxbwrite,uf,float(wlo[i]),5,i+1L
  fxbwrite,uf,float(whi[i]),6,i+1L
  for ik=0L,moy-1 do fxbwrite,uf,float((defarf.(oy[ik]))[i]),ik+1+6,i+1L
endfor
fxbfinish,uf

;	generate generic header for ERREXT=PCACOMP extension
CALERR='PCACOMP'
fxbhmake,hdr0e2,mcomp,CALERR,'the name of the extension',/date,/INITIALIZE
fxaddpar,hdr0e2,'XTENSION','BINTABLE','binary table extension'
fxaddpar,hdr0e2,'BITPIX',8,'8-bit bytes'
fxaddpar,hdr0e2,'NAXIS',2,'2-dimensional binary table'
fxaddpar,hdr0e2,'EMETHOD','PCA1DADD','Type of error encoding'
fxaddpar,hdr0e2,'HDUCLASS','CXC','file format is not OGIP standard.'
fxaddpar,hdr0e2,'HDUCLAS1','CALERR','extension contains principle components analysis data.'
fxaddpar,hdr0e2,'HDUCLAS2','PCACOMP','extension contains PCA eigenvalues and eigenvectors'
fxaddpar,hdr0e2,'HDUCLAS3','PCA1DADD','extension contains 1D products of PCA'
;       add some extra stuff
fxaddpar,hdr0e2,'LONGSTRN','OGIP 1.0','The HEASARC Long String Convention may be used'
fxaddpar,hdr0e2,'ORIGIN','CfA','Source of FITS file'
fxaddpar,hdr0e2,'CREATOR',v_CREATOR,'Program creating this file'
fxaddpar,hdr0e2,'REVISION',1
fxaddpar,hdr0e2,'RESPFILE','none'

;	make a structure out of the arrays that need to be written out:
pcacompstr={component:0, fvariance:0.0D, eigenval:0.0, eigenvec:fltarr(nbin)}
pcacompstr=replicate(pcacompstr,mcomp+1)
pcacompstr.component=[indgen(mcomp),-1]
pcacompstr.eigenval=[sqlamb[0L:mcomp-1L],1.0]
pcacompstr.fvariance=[double(pcafrac[0L:mcomp-1L]),((1.0d - residvar)>0)]
xevec=fltarr(nbin,mcomp+1) & for i=0L,mcomp-1L do xevec[*,i]=pcomp[*,i] & xevec[*,mcomp]=pcaresid[*]
pcacompstr.eigenvec=xevec
mwrfits,pcacompstr,outfile,hdr0e2

if vv gt 1000 then stop,'halting; type .CON to continue'

return
end

;	example run

sqlambf='/data/snafu/kashyap/Cal/CalErr/hrcsletg/sqlamb_hrcsletg.txt'
pcompf='/data/snafu/kashyap/Cal/CalErr/hrcsletg/pcomps_hrcsletg.txt'
avgarff='/data/snafu/kashyap/Cal/CalErr/hrcsletg/avgarf_hrcsletg.txt'
defarff='/data/legs/rpete/flight/MC_uncertainties/vinay_hrcs_letg/sample.garf'
if not keyword_set(areffil) then areffil='/data/snafu/kashyap/Cal/CalErr/hrcsletg/aref_hrcsletg.fits'
if n_elements(vthresh) eq 0 then vthresh=0.95
if not keyword_set(verbose) then verbose=1

pca2aref,sqlambf,pcompf,avgarff,defarff,areffil,vthresh=vthresh,ncomp=ncomp,verbose=verbose
spawn,'ls -l '+areffil

end