;^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^Lance Simms, Stanford University 2009
;pro median_difference_of_reads_xRG
;
;PURPOSE:
;	To take a collection of images of a given type, i.e. a combination of
;	a given number of reads and object type (Flat, Dark, M16,...) and get
;	the median value of the pixels and then write them to a fits file.
;
;	The resulting file will be an Naxis1 x Naxis2 x Naxis3 image.
;
;INPUTS:
;       The script is run in one of two ways, depending on the inputs
;
;	1) A combination of:
;	   DetStr       - 'H1RG-022', 'H2RG-32-147', 'H4RG-007' 
;	   ElecStr      - 'ASIC', 'LEACH'
;	   NReads       - Number of reads in all of the files
;	   ObjectName   - The name of the object ('Dark', 'Flat', etc.)
;	   ThisFilter   - The filter used to observe ('G', 'I', 'Y')
;	   [ObjectText] - Optional string for LEACH files
;
;	OR
;	
;	2) FileList - 
;             A list of files to be included in the median.  The parameters that
;	      will be put in the header/filename of the output fits file will
;             be taken from the first file in the list
;
;OUTPUTS:
;	Median[ObjectName]_[Filter]_[NReads].fits
;
;KEYWORDS:
;	DETSTR: string
;	  The name of the detector used ('H1RG-022', 'H2RG-32-147', 'H4RG-007')
;	ELECSTR: string
;	  The name of the electronics used ('ASIC', 'LEACH')
;	DATE: string
;	  The Date during which they were taken.  Also the subdirectory that
;	  the files are in ('07Dec18', '07Apr25', etc.)
;	OFFSET: int
;	  An offset for files that were written improperly.  The readu
;	  function reads the values with this offset in place.
;	OBJECTNAME: string
;	  The string that indicates what files are being sought.  The filenames
;	  should begin with this string ('Dark', 'Flat', etc.).
;	OBJECTTEXT: string
;	  A string that contains extra parameters in the filename for 
;	  exposures taken with the LEACH ('Bfilter_300_1_0_0_0', etc.)
;	NREADS: string
;	  The number of reads (NAXIS3) that are present in the file.  This
;	  is also a string in the filename that will be used as a regular
;	  expression in searching for matched files. It is a string because
;	  certain files were recorded with floats (i.e. '30' or '30.0000')
;	THISFILTER: string 
;	  The filter through which the files were taken ('I','G','Y','Blank')
;	OVERRIDE: int
;	  0 - Don't overwrite the median file if it already exists
;	  1 - Overwrite the file even if it exists
;	USEREF: int
;	  0 - Don't subtract the reference pixels from the science pixels
;	  1 - Subtract scaled, averaged, and filtered reference pixels
;	BIASSUB: int
;	  0 - Don't subtract the bias (zeroth read) from other reads
;	  1 - Subtract the zeroth read before taking the median
;	DARKSUB: int
;	  0 - Don't subtract a dark exposure with matching cadence before median
;	  1 - Subtract a dark exposure with matching cadence before median
;	THISNDROPS: int
;	  The number of drops used in the exposures going into the median
;       WINDOWMODE: int
;         A boolean stating whether or not the exposures were taken in 
;         windowmode.  
;	  0 - indicates full frame mode  
;	  1 - indicates window mode.
;	FIRSTFILE: int 
;	  The index of the first file in the list returned by the file search
;	  to be used in the median
;	LASTFILE: int 
;	  The index of the last file in the list returned by the file search
;	  to be used in the median
;       TVFLAG: int
;	  0 - Default, Don't plot anything
;         1 - Plot to the X window
;         2 - Plot to postscript or encapsulated postscript
;         3 - Plot to a .png
;       TEMPERATURE: int
;         The temperature at which the files were recorded.  This keyword 
;         remains at zero in most instances; only used when the darks 
;         were recorded at different temps.
;		
;NOTE:
;	For dark exposures, the median of the zeroth read is taken and
;	this is the first frame in the cube that is written to disk.
;	However, for the following reads, the bias is subtracted from
;	each exposure. So the 1st through Nth frames in the cube are
;	the median of bias-subtracted reads from individual exposures.
;
;CALLING SEQUENCE:
; 1) median_of_all_reads_xrg, FileList=
;
; 2) median_of_all_reads_xrg, date=, objectname=, nreads=, elecstr=, detstr=
;
;///////////////////////////////////////////////////////////////////////////////////////////
;NOTES:
;  It seems to work better if we subtract the reference pixels before taking the median
;  since the noise might vary across the exposures.  Using the /useref keyword for this 
;  script is recommended.
;  
;  For flats, it seems to work better if we subtract the median dark current from the 
;  slope of the median flat.  That is, we should run median_of_all_reads_xrg without /darksub 
;  and then run slopefit_xrg with /darksub
;
;///////////////////////////////////////////////////////////////////////////////////////////
;EXAMPLE:
;--H1RG
;  median_of_all_reads_xrg, DetStr='H1RG-022', ElecStr='LEACH', Date='2007Nov13_2', thisfilter='Blank', objectname='Dark', objecttext='Bfilter_300_1_0_0_0_', nreads=300, /override
;  median_of_all_reads_xrg, DetStr='H1RG-022', ElecStr='LEACH', Date='2007Nov13_2', thisfilter='Blank', objectname='Dark', objecttext='300_1_0_0_0_', nreads=300, /override
;  median_of_all_reads_xrg, DetStr='H1RG-022', ElecStr='ASIC', Date='07Nov19', objectname='Dark', thisfilter='Blank', nreads=2500, /useref,/override
;  median_of_all_reads_xrg, detstr = 'H1RG-022', elecstr='LEACH', date='2007Nov14', nreads=200, thisfilter='Blank', objectname='Dark', objecttext='Bfilter_200_1_0_0_0_',/useref
;  median_of_all_reads_xrg, DetStr='H1RG-022', ElecStr='LEACH', Date='2007Nov13_2', thisfilter='BI, objectname='Dark', objecttext='Bfilter_200_1_0_0_0_', nreads=200, /override, /useref
;****************************************
;H2RG:
;DARKS.....
;median_of_all_reads_xrg, DetStr='H2RG-32-147', ElecStr='ASIC', Date='07Dec12', thisfilter='Blank', objectname='Dark', objecttext='15.0000_Reads', nreads=15, /override, firstfile=0, lastfile=20, /useref
;median_of_all_reads_xrg, DetStr='H2RG-32-147', ElecStr='ASIC', Date='07Dec18', thisfilter='Blank', objectname='Dark', objecttext='', nreads=10, /override, /useref
;FLATS.....
;median_of_all_reads_xrg, DetStr='H2RG-32-147', ElecStr='ASIC', Date='07Dec12', thisfilter='G', objectname='Flat', objecttext='30_Reads', nreads=30, /override, /useref
;median_of_all_reads_xrg, DetStr='H2RG-32-147', ElecStr='ASIC', Date='07Dec14', thisfilter='I', objectname='Flat', nreads=30, /override, /useref
;median_of_all_reads_xrg, DetStr='H2RG-32-147', ElecStr='ASIC', Date='07Dec15', thisfilter='Y', objectname='BlankSky', nreads=50, /override, /useref
;
;H2RG-001 at Teledyne:
;  median_of_all_reads_xrg, detstr='H2RG-001', elecstr='ASIC', date='08Dec15', objectname='Dark', temperature=100, nreads=100, thisfilter='Blank'
;  median_of_all_reads_xrg, detstr='H2RG-001', elecstr='ASIC', date='08Dec13', objectname='Dark', temperature=100, nreads=100, thisfilter='Blank', objecttext='5.00V'
;
;AUTHOR: Lance Simms 2/07
;*************************************************************************
pro Median_Of_All_Reads_xRG, DetStr=DetStr, ElecStr=ElecStr, $
	 Date=Date, 	          Offset = Offset, $
	 ObjectName=ObjectName,   ObjectText = ObjectText, $
	 NReads = NReads,         ThisFilter = ThisFilter,$
	 OverRide = OverRide,     UseRef  = UseRef, $
	 BiasSub = BiasSub,       DarkSub = DarkSub, $
	 ThisNDrops = ThisNDrops, WindowMode = WindowMode, $
	 FirstFile = FirstFile,   LastFile = LastFile, $
	 TvFlag = TvFlag,	  FileList = FileList, $
         Temperature = Temperature

  Common KeyParams, KeyStr

  ;First of all check to see if there is a list of files to be used
  If N_Elements(FileList) ne 0 then begin
     ReadCol, FileList, RawFiles, Format = 'A'
     FitsFileName = RawFiles(0)
     NumRawFiles  = N_Elements(RawFiles)
  endif
  If N_Elements(Temperature) eq 0 then Temperature = 0

  ;Get all of the defaults for the keywords from the KeywordStruct file
  @KeywordStruct_xRG.pro
  @PlotSettings_xRG.pro

  ;Switch to the directory with the files
  cd, KeyStr.RawObjectDir

  ;Form the search string (different for ASIC/LEACH, search for a median file
  if N_Elements(FileList) eq 0 then begin
     If KeyStr.ElecStr eq 'ASIC' then begin
       MedSearchStr = KeyStr.ObjectName+'*'+KeyStr.ObjectText+'*.fits'
     EndIf Else if KeyStr.ElecStr eq 'LEACH' and $
       stregex(KeyStr.ObjectName, 'Flat',/boolean) then begin
       MedSearchStr = '*'+KeyStr.ObjectName+'*'+KeyStr.ObjectText+'?.fits'
     Endif Else if KeyStr.ElecStr eq 'LEACH' and KeyStr.Date eq '05May08' then begin
       if KeyStr.ObjectName eq 'Dark' then ObjectNameText='dark'
       MedSearchStr = ObjectNameText+'*'+KeyStr.ObjectText+'*.fits'
     Endif Else Begin
       MedSearchStr = KeyStr.ObjectName+'*'+KeyStr.ObjectText+'?.fits'
     Endelse
     RawFiles       = File_Search(KeyStr.RawObjectDir+MedSearchStr, $
			         count = NumRawFiles)
     LenOText       = strlen(KeyStr.ObjectText)
     If strmid(KeyStr.ObjectText,0,1) eq '_' then BegSep = '' else BegSep='_'
     If strmid(KeyStr.ObjectText,LenOText-1,LenOText) eq '_' $
					     then EndSep = '' else EndSep='_'
     if Temperature ne 0 then KeyStr.ObjectName = KeyStr.ObjectName+'_'+Strtrim(Temperature,2)+'K'
     MedFileName    = KeyStr.RedObjectDir+'Median'+$
		      KeyStr.ObjectName+'_'+$
		      KeyStr.ThisFilter+BegSep+KeyStr.ObjectText+EndSep+$
                      KeyStr.DarkSubStr+$
		      '_'+strtrim(KeyStr.NAxis3,2)+'_Reads_'+$
		      KeyStr.RefPixStr+'.fits'
     MedFileName    = Clean_FileName(MedFileName)
     VarFileName    = KeyStr.RedObjectDir+'VarianceOf_'+$
		      KeyStr.ObjectName+'_'+$
		      KeyStr.ThisFilter+BegSep+KeyStr.ObjectText+EndSep+$
                      KeyStr.DarkSubStr+$
		      '_'+strtrim(KeyStr.NAxis3,2)+'_Reads_'+$
		      KeyStr.RefPixStr+'.fits'
     VarFileName    = Clean_FileName(VarFileName)
     MedFileExists  = File_Search(MedFileName)
  EndIf
  LenOText       = strlen(KeyStr.ObjectText)
  If strmid(KeyStr.ObjectText,0,1) eq '_' then BegSep = '' else BegSep='_'
  If strmid(KeyStr.ObjectText,LenOText-1,LenOText) eq '_' $
                                             then EndSep = '' else EndSep='_'
  MedFileName    = KeyStr.RedObjectDir+'Median'+$
                      KeyStr.ObjectName+'_'+$
                      KeyStr.ThisFilter+BegSep+KeyStr.ObjectText+EndSep+$
                      KeyStr.DarkSubStr+$
                      '_'+strtrim(KeyStr.NAxis3,2)+'_Reads_'+$
                      KeyStr.RefPixStr+'.fits'
  MedFileName    = Clean_FileName(MedFileName)
  VarFileName    = KeyStr.RedObjectDir+'VarianceOf_'+$
                      KeyStr.ObjectName+'_'+$
                      KeyStr.ThisFilter+BegSep+KeyStr.ObjectText+EndSep+$
                      KeyStr.DarkSubStr+$
                      '_'+strtrim(KeyStr.NAxis3,2)+'_Reads_'+$
                      KeyStr.RefPixStr+'.fits'
  VarFileName    = Clean_FileName(VarFileName)
  MedFileExists  = File_Search(MedFileName)

  ;If Dark does not exist, form the median and write the image
  If MedFileExists(0) eq '' or OverRide then begin
    print, string(13b)
    print, '\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\'
    print, 'Forming Median Image for Combination of:'
    print, 'Object:          ' + KeyStr.ObjectName
    print, 'Object Text:     ' + KeyStr.ObjectText
    print, 'Filter:          ' + KeyStr.ThisFilter
    print, 'Number of Reads: ' +  strtrim(KeyStr.TotalReads,2)
    print, string(13b)
    print, 'Output File will be written to :' & print, string(13b)
    print, MedFileName & print, string(13b)

    ;Allow for certain files to be excluded with FIRSTFILE and LASTFILE KWord
    If N_Elements(FirstFile) eq 0 then FirstFile = 0
    If N_Elements(LastFile)  eq 0 then LastFile  = NumRawFiles-1
    RawFiles = RawFiles(FirstFile: LastFile)

    ;Filter out the ones that don't have the right parameters
    RawFiles = Return_File_List_xRG(RawFiles, ThisFilter = KeyStr.ThisFilter, $
                            ThisNDrops=ThisNDrops, ThisNAxis3 = KeyStr.NAxis3,$
			    WindowMode = WindowMode, Temperature=Temperature)
    ;Update all of the elements of the key structure
    FileCheck = Return_File_Params_xRG(RawFiles(0))
 
    ;Ask the user if they would like to continue
    print , 'About to form Median from files: '
    print , '---------------------------------'
    print , RawFiles
    print , '---------------------------------'
    print , 'Type .continue and hit enter to continue'
    stop

    ;Form a dummy array for all of the darks of a given Naxis3 to provide the
    ;necessary keywords for a fits header
    NumRawFiles = N_elements(RawFiles)
    MedianAxes = [KeyStr.Naxis1, KeyStr.Naxis2, KeyStr.TotalReads]
    MkHdr, MedFitsHeader, Naxis, MedianAxes
    FxAddPar, MedFitsHeader, 'DIRDATE', KeyStr.Date, $
	'Date on which the images were taken'
    FxAddPar, MedFitsHeader, 'REFPXSUB', KeyStr.RefPixStr, $
	'Null String indicates no removal of Common Mode Noise'
    FxAddPar, MedFitsHeader, 'NGROUPS', KeyStr.NGroups, $
	'Number of Groups'
    FxAddPar, MedFitsHeader, 'NREADS', KeyStr.NReads, $
	'Number of Reads'
    FxAddPar, MedFitsHeader, 'NDROPS', KeyStr.NDrops, $
	'Number of Drops'
    FxAddPar, MedFitsHeader, 'FILTER', KeyStr.ThisFilter, $
	'Filter used in images of median'
    FxAddPar, MedFitsHeader, 'BIASSUB', KeyStr.BiasSub, $
	'Individual bias frames were subtracted'
    FxAddPar, MedFitsHeader, 'DARKSUB', KeyStr.DarkSub, $
        '0 - No dark subtraction ; 1 - A dark was subtracted'
    If KeyStr.BiasSub then $
    FxAddPar, MedFitsHeader, 'COMMENT', $
	'Frame 0 = Median of zeroth reads; Frames 1-N = Median of Bias '+$
	'subtracted reads'
 
    ;Add all of the keywords that we should have
    MedFitsHeader=update_header_xrg(MedFitsHeader) 

    ;Copy the header and change it slightly for the Variance file
    VarFitsHeader = MedFitsHeader
    FxAddPar, VarFitsHeader, 'BITPIX', -32 

    ;Put all of the filenames that were used to form the median in the header
    for FileNum = 0, NumRawFiles-1 do begin
        FileBase = File_BaseName(RawFiles(FileNum))
        FxAddPar, MedFitsHeader, 'FILE'+Strtrim(FileNum, 2), $
	          FileBase, 'Name of file used in median'
    endfor

    If KeyStr.DarkSub then begin
       ;Check to see if we should subtract the median dark
       print, 'Using Dark : '+ KeyStr.DarkName
       Fits_Read_DataCube, KeyStr.DarkName, DarkIm, DarkHeader
       KeyStr.DarkSubStr = 'DS'
    EndIf Else Begin
         KeyStr.DarkSubStr = ''
    EndElse

    ;Open the datacube for writing using the Fits_*_DataCube routines
    Position1 = Fits_Open_DataCube(MedFileName, MedFCB, MedFitsHeader, 'Write')
    Position2 = Fits_Open_DataCube(VarFileName, VarFCB, VarFitsHeader,'Write')

    ;Form Cube of Medians for each Read
    for ReadNum = 0, KeyStr.TotalReads - 1 do begin
       print, 'Processing Read Number '+Strtrim(ReadNum, 2) + ' Out of ' + $
	      Strtrim(KeyStr.TotalReads,2) +' Reads'

       ;Form an array to hold the values for the median
       RawFrames = LonArr(KeyStr.Naxis1, KeyStr.Naxis2, NumRawFiles)

       ;Loop through files for each read and take median at end
       for FileNum = 0, NumRawFiles - 1 do begin

           ;Double check filter and print out the progress
	   Fits_Read,  RawFiles(FileNum), NoData, RawFitsHeader, /header_only
           Filter      = SxPar(RawFitsHeader ,'FILTER')
	   print, RawFiles(FileNum)
	   print, 'Read:  ' , ReadNum, ' ' , 'File: ', FileNum, ' ', $
	          'Filter:' , Filter

	   ;Get an average offset for the reference pixels
           If ReadNum eq 0 and FileNum eq 0 then begin
              Fits_Read_DataCube, RawFiles(FileNum), RawOffset, RawFitsHeader,$
		 ZStart = 0, ZStop = 0
	      MeanOff = mean([RawOffset(1:3,4:KeyStr.Naxis2-1)])
	   EndIf     
	   
	   ;Read the Bias frame every time to prevent memory overflow
	   If ReadNum gt 0 then $
	      Fits_Read_DataCube, RawFiles(FileNum), RawBias, RawFitsHeader, $
	         ZStart = 0, ZStop = 0, Offset=Offest
	   If ReadNum gt 0 and KeyStr.UseRef then RawBias = $
	      Return_RefPix_xRG(RawBias, 1, MeanOff=MeanOff)
	   
	   ;Get Read of Raw Frame
           Fits_Read_DataCube, RawFiles(FileNum), RawFrameTemp, RawFitsHeader, $
	      ZStart = ReadNum, ZStop = ReadNum, Offset=Offset
 	   RawFrames(*,*,FileNum) = long(RawFrameTemp(*,*))

           ;Subtract the reference pixels from this read
           If KeyStr.UseRef then RawFrames(*,*,FileNum) = $
	      Return_RefPix_xRG(RawFrames(*,*,FileNum),1, MeanOff=MeanOff)

	   ;Subtract the Bias and a dark frame if requested
	   If KeyStr.BiasSub eq 1 and ReadNum gt 0 then $
	      RawFrames(*,*,FileNum) = RawFrames(*,*,FileNum)-RawBias

	   If KeyStr.DarkSub eq 1 and ReadNum gt 0 then begin $
	      DarkFrame = long(DarkIm(*,*,ReadNum))
	      RawFrames(*,*,FileNum) = $
	      RawFrames(*,*,FileNum)-DarkFrame
              print, 'It might be better to subtract the dark from the slope fit'
	      stop
	   EndIf

	   ;Default verbosity
	   print, 'Min Pixel Value',  Min(RawFrames(*,*,FileNum)), string(9B),$
                  'Max Pixel Value',  Max(RawFrames(*,*,FileNum))
           print, 'Mean Pixel Value', Mean(RawFrames(*,*,FileNum)),string(9B),$
                  'Std Dev',          StdDev(RawFrames(*,*,FileNum))

           ;Do some plotting if desired
           If KeyStr.TvFlag eq 1 then begin
             window, 0 , xpos = 2*fix(float(XScreenSize)/3), $
                         ypos = 0, $
                         xsize = fix(float(XScreenSize)/3), $
                         ysize = fix(float(YScreenSize)/2)

             tvscl, congrid(RawFrames(*,*,FileNum),512,512)
             PixHis = histogram(RawFrames(*,*,FileNum) , locations = LocPix)
             window, 1 , xpos = 2*fix(float(XScreenSize)/3), $
                         ypos = fix(float(YScreenSize)/2), $
                         xsize = fix(float(XScreenSize)/3), $
                         ysize = fix(float(YScreenSize)/2)

             Plot, LocPix, PixHis, psym = 10
           EndIf

	endfor

	;Take the median row by row across the files; memory errors otherwise
	MedianOfRead   = LonArr(KeyStr.Naxis1, KeyStr.Naxis2)
        VarianceOfRead = FltArr(KeyStr.Naxis1, KeyStr.Naxis2)

	;Extremely slow, but StdDev function can not operate across dimensions
        for Row = 0, KeyStr.Naxis2 -1 do begin
	  for Col = 0, KeyStr.Naxis1 - 1 do begin
            djs_iterstat, RawFrames(Col, Row,*), sigrej=3.0, mean=MeanPix, median=MedPix, sigma=SigPix
            MedianOfRead(Col, Row)   = MedPix ; Long(median(RawFrames(Col,Row,*)))
            VarianceOfRead(Col, Row) = SigPix ;StdDev(RawFrames(Col,Row,*))
          endfor
	endfor

	;Write it to the cube and delete it afterwards to avoid crash
        DataCubeWrite = Fits_Append_DataCube(MedFCB,MedianOfRead)
	DataCubeWrite = Fits_Append_DataCube(VarFCB, VarianceOfRead)
	Delvarx, MedianOfRead
	Delvarx, VarianceOfRead
	Delvarx, RawFrames

    endfor

    ;Close the fits files after they have been written
    DataCubeClose = Fits_Close_DataCube(MedFCB)
    DataCubeClose = Fits_Close_DataCube(VarFCB)

    print, 'Median Dark Image written with filename: '
    print, MedFileName + ' has been written'  & print, string(13b)
    print, VarFileName+ ' has been written' & print, string(13b)
    
  EndIf Else begin

    ;The File already exists
    print, 'Median Dark: ' & print, string(13b)
    print,  MedFileName    & print, string(13b)
    print, 'Already Exists. Use OVERRIDE keyword to overwrite this'

  EndElse

stop

end