;pro Photometric_Analysis
;
;PURPOSE:
;
;	This routine will run through all of the IRAF adapted routines to retrun
;	a psf and CMD for a given .fits file
;
;
;INPUTS:
;	SlopeFileName - The name of the slopefit that contains the mask
;		        of saturated pixels.  In most cases it should 
;		        be the 0th dither
;
;KEYWORDS:
;
;	MODE -
;	0 - Plot 1 set of slopes for each radius for one star
;	1 - Plot the average of the slopes vs. radius for each star
;	2 - Scatter plot the values for a given read number vs. dither
;
;CALLING SEQUENCE:
;photometric_analysis,'/nfs/slac/g/ki/ki09/lances/Reduced/07Apr26/M13/M13_G_SlopeFlatFieldedCenteredDitheredMean.fits'
;QPhotometric_Analysis,'/nfs/slac/g/ki/ki09/lances/Reduced/07Apr27/Landolt110956/Landolt110956_I_SlopeFlatFieldedCenteredDitheredMean.fits'
;****************************************************************************
pro QPhotometric_Analysis, SlopeFileName, $
	 RootDir = RootDir, Date = Date, $
	 ReducedDir = ReducedDir, ObjectName = ObjectName, $
	 SubtractDark = SubtractDark, FilterStr = FilterStr, $
	 TvFlag = TvFlag, Rem60Hz = Rem60Hz, $ 
	 Offset = Offset, Mode = Mode

If N_Elements(Mode) eq 0 then Mode = 0

Common SlopeParams, SlopeStr
Common KeyParams, KeyStr
Common ColorStr, Color, ColorNum, ColorTri, ColorNumTri

;Include all of the keywords from the KeywordStruct file
@KeywordStruct.pro
@PlotSettings.pro

;Get the saturated mask from the slope file
FileParams = Return_File_Params(SlopeFileName, /slope)
Fits_Read, SlopeFileName, SlopeIm, SlopeHeader


;Determine the name of the output file
If KeyStr.TvFlag eq 2 then $
   device, filename = KeyStr.PlotDir+KeyStr.ObjectName+'_RadialPersistence.ps'

;prepare some plotting indices
Masks	      = Return_Pixel_Masks()

;Find the regions that were occupied by bright stars as opposed to 
;very high E cosmic rays or messed up pixels
BoxSize = 70
x1 = 1965
y1 = 1914
x2 = 1677 
y2 = 2017
x3 = 1876
y3 = 2138

; x3=1964
; y3=1914
; x2=1677
; y2=2016
; x1=1413
; y1=2113

;x1=1415
;y1=2130
;x2=1679
;y2=2032
;x3=1877
;y3=2154


;x3=1967
;y3=1932
;x2=1679   
;y2=2032
;x1=1413
;y1=2140

SlopeIm = float(readfits('/nfs/slac/g/ki/ki09/lances/Reduced/07Apr27/Landolt110956/Landolt110956_G_SlopeFlatFieldedCenteredDitheredMean.fits'))
f1=total([SlopeIm(x1-boxsize/2:x1+boxsize/2,y1-boxsize/2:y1+boxsize/2)])-3.*boxsize^2.
f2=total([SlopeIm(x2-boxsize/2:x2+boxsize/2,y2-boxsize/2:y2+boxsize/2)])-3.*boxsize^2.
f3=total([SlopeIm(x3-boxsize/2:x3+boxsize/2,y3-boxsize/2:y3+boxsize/2)])-3.*boxsize^2.
print,-2.5*alog10(f1/f2),-2.5*alog10(f1/f3)


stop

End