#from numarray import * ## Changed to numpy Aug 2007
from numpy import *
import pyfits
#import random

from run_simpler_skymap_datons import *
from simpler_skymap_datons import *

Datestring   = '.070921a'
VPstring     = '_vp0110_g001'
Estring      = '_001-002GeV'
Modelstring = '.PSFonVP113C279b'
WhichEBin    = 7
RealDataDir  = '/AstroData/Legacy_EGRET_Maps/'
FilesTupleList =[('DBG3PSFM0015809OnVP113C279b_vp0110_g001.070912b.fits',
                  RealDataDir+'exposr'+VPstring+'.fits',
                  'DBG3PadModelFlux'+Modelstring+VPstring+Estring+Datestring+'.fits',
                  'DBG3bPadExpAreaTime'+VPstring+Estring+Datestring+'.fits',
                 'DBG3PadModelData'+Modelstring+VPstring+Estring+Datestring+'.fits',
                  RealDataDir+'counts'+VPstring+'.fits', 'DBG3bPad'+Estring+'TSTcounts'+VPstring+'.fits')
                 ]

# Do for each in the list:
for this_tuple in FilesTupleList:
    ( FluxModelFile      , ExpAreaTimeFile, 
      OutPadFluxModelFile, OutPadExposrFile, OutPadModelDataFile,
      InRealCountsFile, OutPadRealCountsFile) = this_tuple
#
# 1/ Get Input ModelFlux:
#
    FluxModelHDULst = pyfits.open(FluxModelFile)
#   "touch" them to getthem read in properly:
    print FluxModelHDULst[0].header
    print FluxModelHDULst[0].data[0][0][0]
#
# 2/ Get Input Exposure File:
#
    ExpAreaTimeHDULst = pyfits.open(ExpAreaTimeFile)
#   "touch" them to getthem read in properly:
    print '\n', ExpAreaTimeHDULst[0].header
    print ExpAreaTimeHDULst[0].data[0][0][0]
    tempdata = zeros((1,ExpAreaTimeHDULst[0].header['NAXIS2'],ExpAreaTimeHDULst[0].header['NAXIS1']),dtype=float)
    tempdata[0] = ExpAreaTimeHDULst[0].data[WhichEBin]
    ThisExpAreaTime_HDU = pyfits.PrimaryHDU(data=tempdata)
    print ThisExpAreaTime_HDU.header
    print ThisExpAreaTime_HDU.data[0][0][0], '\n'
#
# 3/ Get Input real counts data File:
#
    RealCountsHDULst = pyfits.open(InRealCountsFile)
#   "touch" them to getthem read in properly:
    print '\n', RealCountsHDULst[0].header
    print RealCountsHDULst[0].data[0][0][0]
    tempdata = zeros((1,RealCountsHDULst[0].header['NAXIS2'],RealCountsHDULst[0].header['NAXIS1']),dtype=float)
    tempdata[0] = RealCountsHDULst[0].data[WhichEBin]
    ThisRealCounts_HDU = pyfits.PrimaryHDU(data=tempdata)
    print ThisRealCounts_HDU.header
    print ThisRealCounts_HDU.data[0][0][0], '\n'
#
# 4+5+6/ Pad the input files:
#
    PadFluxModel_hdu = simpler_zero_pad_to_square_power2(FluxModelHDULst[0])         
    PadFluxModel_HDU = pyfits.PrimaryHDU(
                         data=PadFluxModel_hdu.data)    #, header=PadFluxModel_hdu.header)
    PadFluxModel_HDU.writeto(OutPadFluxModelFile)
#
    PadExpAreaTime_hdu = simpler_zero_pad_to_square_power2(ThisExpAreaTime_HDU)         
    PadExpAreaTime_HDU = pyfits.PrimaryHDU(
                         data=PadExpAreaTime_hdu.data, header=PadExpAreaTime_hdu.header)
    PadExpAreaTime_HDU.writeto(OutPadExposrFile)
#
    PadRealCounts_hdu = simpler_zero_pad_to_square_power2(ThisRealCounts_HDU)         
    PadRealCounts_HDU = pyfits.PrimaryHDU(
                         data=PadRealCounts_hdu.data, header=PadRealCounts_hdu.header)
    PadRealCounts_HDU.writeto(OutPadRealCountsFile)
#
# 5/ Now find model_data = exposure*model_flux for each image bin:
#
   
    this_padmodeldata_HDU = ExpectedDataModelHDU(PadFluxModel_HDU, PadExpAreaTime_HDU,
                              OutPadModelDataFile)