GalpropFileList = ['/AstroData/GALPROP/ics_isotropic_skymap_49_700405','/AstroData/GALPROP/pion_decay_skymap_49_700405','/AstroData/GALPROP/bremss_skymap_49_700405']
ExposrFile = '/AstroData/exposr_g1-9_g004.fits'
ExposrHistory = ExposrFile
WantEBand = [1000.,1000000.]
WhichExp3rdBin = 1

### This commented-out set of parameters is for 128x128 output fits files:
#WantExpRebin = [6,4]	# Usual exposre number of bins is 2x GALPROP
#WantFlxRebin = [3,2]    # Rebinning to 120 x 90 prior to padding -- for speed

## This set of parameters is for 64x64 output fits files:
WantExpRebin = [12,8]	# Usual exposre number of bins is 2x GALPROP
WantFlxRebin = [6,4]    # Rebinning to 120 x 90 prior to padding -- for speed

##IntTimeLst = [806.,242.,73.,22.,7.]	# In "GALPROP" units (1500 ???).
IntTimeLst = [1.]	# Is norm OK as is?

#OutFileList = ['ModelFlux_SpherePadded.060521.fits', 'Exposure_SpherePadded.060521.fits', 'ModelData_SpherePadded.060521.fits', 'PoissDatons_Matched.060521.fits']

basic_string1a = 'ModelFluxRebinnd'
basic_string1b = 'ModelFluxPadded'
basic_string2a = 'ExposrTimRebinnd'
basic_string2b = 'ExposrTimPadded'
basic_string3a = 'ModelDataRebinnd'
basic_string3b = 'ModelDataPadded'
basic_string4a = 'PoisDatonRebinnd'
basic_string4b = 'PoisDatonPadded'

#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 *


#
# Beginning 1st module:
# 1/ Start accumulating flux model components:
EBandString = ''
OutModelFluxFitsFile = [basic_string1a+'BremICPionIso1-1000GeV.64sq.070130c.fits',basic_string1b+'BremICPionIso1-1000GeV.64sq.070130c.fits']
model_flux_HDUpair = SkymapModelFluxHDU(GalpropFileList,WantEBand,WantFlxRebin,EBandString,OutModelFluxFitsFile)
#

# -------------------------------------------------------------------------------------
this_modelflux_hdulst = pyfits.open(OutModelFluxFitsFile[1])
#
print this_modelflux_hdulst[0].header


#  2/ Now input the exposure: Rebin to desired bin/pixel sizes:
#

EBandString = ' Eband: '+str(WantEBand[0]/1000.)+'-'+str(WantEBand[1]/1000.)+' GeV'

for IntTime in IntTimeLst:
    OutExpAreaTimeFitsFiles = [basic_string2a+'IntTime'+str(IntTime)+'.64sq.070130c.fits',basic_string2b+'IntTime'+str(IntTime)+'.64sq.070130c.fits']
    print OutExpAreaTimeFitsFiles[0]
    print OutExpAreaTimeFitsFiles[1]
    this_exposr_rebinned_HDU = MatchedExposureMap(ExposrFile,EBandString,WantExpRebin,WhichExp3rdBin,IntTime,OutExpAreaTimeFitsFiles[0],ExposrHistory)
#   2.3/ Pad to power of two:
#    exposr_padded_hdu     = simpler_sphereskymap_pad_to_square_power2(exposr_rebinned_HDU)
#    Note: padding to include spherical symmetry was tried and didnt give great results.
#    SO now it is just padded to zero at the poles but wrapped around properly at the sides.
    print this_exposr_rebinned_HDU.data
    print this_exposr_rebinned_HDU.header
    this_exposr_padded_hdu     = simpler_zero_pad_to_square_power2(this_exposr_rebinned_HDU)
    this_exposr_padded_HDU    = pyfits.PrimaryHDU(data=this_exposr_padded_hdu.data,header=this_exposr_padded_hdu.header)
    this_exposr_padded_HDU.writeto(OutExpAreaTimeFitsFiles[1])
    
#
# 3/ Now find model_data = exposure*model_flux for each image bin:
#
#   OutModelDataFileList = [basic_string3a+'IntTime'+str(IntTime)+'.64sq.070129.fits',basic_string3b+'IntTime'+str(IntTime)+'.64sq.070129.fits']

    FluxExpOutfileTupleList = [(this_modelflux_hdulst[0], this_exposr_padded_HDU, basic_string3b+'IntTime'+str(IntTime)+'.64sq.070130c.fits', basic_string4b+'IntTime'+str(IntTime)+'.64sq.070130c.fits')]

    for this_tuple in FluxExpOutfileTupleList:
#        (this_modelflux_fitsfile, this_expareatime, this_modeloutfile,this_OutPoisDatonsFitsFile) = this_tuple
    
        this_modeldata_HDU = ExpectedDataModelHDU(this_tuple[0],this_tuple[1],this_tuple[2])

###
### 4/ Get some Poisson counts! what fun!!
###
###TAKEN OUT FOR DEBUGGING
###        this_poisdatons_HDU = PoisDatonsHDU(this_modeldata_HDU,this_tuple[3])