;+
;samp2fits.pro
;	make an AREF file with all of Pete's sample ARFs
;
;usage
;	arfdef='obs866/acis.arf'
;	arfdir='store_arfs'
;	arfroot='866'
;	.run samp2fits
;
;	reads in the default ARF, arfdef, and the random sample
;	arfdir/arfroot_* and writes out to arfroot_samp_aref.fits
;
;vinay k (Sep2009, based on pca2fits_Cedge.pro)
;	updated initialization (VK; Jun2014)
;-

;	initialize
if not keyword_set(arfdef) then arfdef='obs866/acis.arf'
if not keyword_set(arfdir) then arfdir='store_arfs/'
if not keyword_set(arfroot) then arfroot='866'
outfile=arfroot+'_samp_aref.fits'

;	read in the default ARF
if not keyword_set(arf0) then arf0=mrdfits(arfdef,1,harf)
elo=arf0.ENERG_LO & ehi=arf0.ENERG_HI & emid=0.5*(elo+ehi) & effar=arf0.SPECRESP & nbin=n_elements(effar)
;
os=sort(emid) & elo=elo[os] & ehi=ehi[os] & effar=effar[os]
wlo=12.398521/ehi & whi=12.398521/elo

;	read in all of Pete's set of ARFs
if n_elements(arf_vec) eq 0 then begin
  arffil=findfile(arfdir+'/'+arfroot+'_*',count=narf)
  arf_vec=fltarr(nbin,narf)
  for i=0,narf-1 do begin
    kilroy,strtrim(i,2)+'..'
    tmp=mrdfits(arffil[i],1,h,/silent) & arf_vec[*,i]=tmp.SPECRESP-effar
  endfor
  bias_arr=0.*effar
endif

;	create the Ancillary Response Error File (AREF)

;	first, define some header keywords
cd,current=cwd
if not keyword_set(v_CREATOR) then v_CREATOR=filepath('samp2fits.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
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','SIM1DADD','additive simulated 1D curves are stored'
	;EMETHOD='PCA1DMUL' for multiplicative Principal components of 1D curves
	;EMETHOD='SIM1DMUL' for storing multiplicative simulated 1D curves
	;EMETHOD='POLYNOM1D' for storing polynomial fit coefficients to 1D curves and errors on the coefficients
fxaddpar,hdr0,'CMETHOD','ARRAY=NOMINAL+SAMPLE[r]'
	;CMETHOD='ARRAY=NOMINAL*(BIAS+sum{r_j*PCA_EVALUE_j*PCA_EVECOR_j}+r*PCA_RESID)'	for PCA1DMUL
	;CMETHOD='ARRAY=BIAS*SIM_r' for SIM1DMUL
	;CMETHOD='ARRAY=BIAS+sum_j{(c_j+r_j*sigc_j)*E^j}' for POLYNOM1D
fxaddpar,hdr0,'BIASEXT','SPECRESP','nominal function and bias will be stored in this extension'
fxaddpar,hdr0,'ERREXT','SIMCOMP','uncertainties will be stored in this extension'
;	and 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 keywords
fxaddpar,hdr0e1,'EMETHOD','SIM1DADD','additive simulated 1D curves are stored'
fxaddpar,hdr0e1,'CMETHOD','ARRAY=NOMINAL+SAMPLE[r]
fxaddpar,hdr0e1,'ERREXT','SIMCOMP','uncertainties will be stored in this extension'
;
;	define the columns and some extra headers
fxbaddcol,i1,hdr0e1,EMID[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,EMID[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,effar[0],'SPECRESP','effective area',TUNIT='cm**2'
fxaddpar,hdr0e1,'TFORM3','1E','format of field'
fxaddpar,hdr0e1,'TLMIN3',min(effar)
fxaddpar,hdr0e1,'TLMAX3',max(effar)
fxbaddcol,i4,hdr0e1,bias_arr[0],'BIAS','offset from nominal',TUNIT='cm**2'
fxaddpar,hdr0e1,'TFORM4','1E','format of field'
fxaddpar,hdr0e1,'TLMIN4',min(bias_arr)
fxaddpar,hdr0e1,'TLMAX4',max(bias_arr)
fxbaddcol,i5,hdr0e1,12.3984/EMID[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,12.3984/EMID[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 write 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(effar[i]),3,i+1L
  fxbwrite,uf,float(bias_arr[i]),4,i+1L
  fxbwrite,uf,float(wlo[i]),5,i+1L
  fxbwrite,uf,float(whi[i]),6,i+1L
endfor
fxbfinish,uf

;	generate generic header for ERREXT=SIMCOMP extension
CALERR='SIMCOMP'
fxbhmake,hdr0e2,narf,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','SIM1DADD','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','SIMCOMP','extension contains simulated datasets'
fxaddpar,hdr0e2,'HDUCLAS3','SIM1DADD','extension contains 1D simulated additive curves'
;	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:
simcompstr={component:0, simcomp:fltarr(nbin)}
simcompstr=replicate(simcompstr,narf)
simcompstr.component=lindgen(narf)+1
xevec=fltarr(nbin,narf) & for i=0L,narf-1L do xevec[*,i]=arf_vec[*,i]
simcompstr.simcomp=xevec
mwrfits,simcompstr,outfile,hdr0e2

print,'' & print,'prism '+outfile+' &' & print,''

end