basic_string1a = 'ModelFluxRebinnd'
basic_string1b = 'ModelFluxPadded'
ERangeMeme = '1-1000GeV.'
DateMeme = '060529.'

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

###############
class Squished:
    def __init__(self,OldModelFluxHDU):
#
        naxis1 = OldModelFluxHDU.header['NAXIS1']
        naxis2 = OldModelFluxHDU.header['NAXIS2']
#
#   Handy variables for defining the broken-power-law Squish:
        middle1 = float(naxis1)/2.
        middle2 = float(naxis2)/2.
        radius0 = 0.15*(  float( min(middle1,middle2) )  )
        delrad1 = 2.*radius0
        radius1 = radius0 + delrad1

#   the "squish" model is: power-law "squishing" down of the region away from center
#   Amounting to roughly 1/2 of the total flux:
        self.data   = OldModelFluxHDU.data
        self.header = OldModelFluxHDU.header

#   Now do it:
        for l1 in range(naxis1):
            el1 = float(l1) - middle1
            for l2 in range(naxis2):
                el2 = float(l2) - middle2
                radiusnow = sqrt(el1*el1 + el2*el2)
#           Simple Squish algorithm means: flat center
                if (radiusnow <= radius0):
                    self.data[0][l2][l1] = OldModelFluxHDU.data[0][l2][l1]
#           Lower order power-law in ring around center circle:
#           Note that at radius0, squish-factor is still unity so is continuous:
                elif (radiusnow <= radius1):
                    squishfactor = radius0/radiusnow
                    self.data[0][l2][l1] = squishfactor*OldModelFluxHDU.data[0][l2][l1]
                else:
                    squishfactor = (radius0/radiusnow)*(radius1/radiusnow)
                    self.data[0][l2][l1] = squishfactor*OldModelFluxHDU.data[0][l2][l1]
##
## SO notice the total normalization change is:
##    inner circle:  \pi^r0^2 ~ 3 * .0225 ~ .07
##    middle ring:   \int_{r0}^{r1}dr 2 \pi r * (r0/r) = 2\pi *(r1-r0)*r0 ~ 6 * .3 * .15 ~ 6 * .045 ~ .3
##    outer ring, 1: \int_{r1}^{middle} 2 \pi r (r0*r1)/r^2 = 2 \pi (r0*r1)*( ln(middle/r1) ) ~ 6*.045*ln(2.)~.2
##    outer ring, 2: -- a bunch of stuff around the corners ... maybe ~.1
##    SO just by squishing one reduces it by very roughly a factor of 0.6 ...
##    The inverse-squish -subtracts- this from the original;
##    Hence the inverse-squish reduces the norm by maybe 40% or so, as a guess.
##-----------------------------------------------------------------------------------------------------
###############
class InverseSquished:
    def __init__(self,OldModelFluxHDU):
#
        naxis1 = OldModelFluxHDU.header['NAXIS1']
        naxis2 = OldModelFluxHDU.header['NAXIS2']
#
#   Handy variables for defining the broken-power-law InverseSquish:
        middle1 = float(naxis1)/2.
        middle2 = float(naxis2)/2.
        radius0 = 0.15*(  float( min(middle1,middle2) )  )
        delrad1 = 2.*radius0
        radius1 = radius0 + delrad1

#   the "Inversesquish" model is: power-law "Inversesquishing" down of the region away from center
#   Amounting to roughly 1/2 of the total flux:
        self.data   = OldModelFluxHDU.data
        self.header = OldModelFluxHDU.header

#   Now do it:
        for l1 in range(naxis1):
            el1 = float(l1) - middle1
            for l2 in range(naxis2):
                el2 = float(l2) - middle2
                radiusnow = sqrt(el1*el1 + el2*el2)
#           Simple InverseSquish algorithm means: flat center
                if (radiusnow <= radius0):
                    self.data[0][l2][l1] = 0.
#           Lower order power-law in ring around center circle:
#           Note that at radius0, Inversesquish-factor is still unity so is continuous:
                elif (radiusnow <= radius1):
                    Inversesquishfactor = 1. - (radius0/radiusnow)
                    self.data[0][l2][l1] = Inversesquishfactor*OldModelFluxHDU.data[0][l2][l1]
                else:
                    Inversesquishfactor = 1. -( (radius0/radiusnow)*(radius1/radiusnow) )
                    self.data[0][l2][l1] = Inversesquishfactor*OldModelFluxHDU.data[0][l2][l1]
##
## SO notice the total normalization change is:
##    inner circle:  \pi^r0^2 ~ 3 * .0225 ~ .07
##    middle ring:   \int_{r0}^{r1}dr 2 \pi r * (r0/r) = 2\pi *(r1-r0)*r0 ~ 6 * .3 * .15 ~ 6 * .045 ~ .3
##    outer ring, 1: \int_{r1}^{middle} 2 \pi r (r0*r1)/r^2 = 2 \pi (r0*r1)*( ln(middle/r1) ) ~ 6*.045*ln(2.)~.2
##    outer ring, 2: -- a bunch of stuff around the corners ... maybe ~.1
##    SO just by squishing one reduces it by very roughly a factor of 0.6 ...
##    The inverse-squish -subtracts- this from the original;
##    Hence the inverse-squish reduces the norm by maybe 40% or so, as a guess.
##-----------------------------------------------------------------------------------------------------
###############
class HalfSquished:
    def __init__(self,OldModelFluxHDU):
#
        naxis1 = OldModelFluxHDU.header['NAXIS1']
        naxis2 = OldModelFluxHDU.header['NAXIS2']
#
#   Handy variables for defining the broken-power-law HalfSquish:
        middle1 = float(naxis1)/2.
        middle2 = float(naxis2)/2.
        radius0 = 0.10*(  float( min(middle1,middle2) )  )
        delrad1 = 2.*radius0
        radius1 = radius0 + delrad1

#   the "Halfsquish" model is: power-law "Halfsquishing" down of the region away from center
#   BUT only in the upper half of the sky
#   Amounting to roughly 1/2 of the total flux:
        self.data   = OldModelFluxHDU.data
        self.header = OldModelFluxHDU.header

#   Now do it:
        for l1 in range(naxis1):
            el1 = float(l1) - middle1
            for l2 in range(naxis2):
                el2 = float(l2) - middle2
                radiusnow = sqrt(el1*el1 + el2*el2)
#           Simple HalfSquish algorithm means: flat center
                if (l2 >= middle2 ):
                    if (radiusnow <= radius0):
                        self.data[0][l2][l1] = OldModelFluxHDU.data[0][l2][l1]
#                       Lower order power-law in ring around center circle:
#                       Note that at radius0, Halfsquish-factor is still unity so is continuous:
                    elif (radiusnow <= radius1):
                        Halfsquishfactor = (radius0/radiusnow)*(radius0/radiusnow)
                        self.data[0][l2][l1] = Halfsquishfactor*OldModelFluxHDU.data[0][l2][l1]
                    else:
                        Halfsquishfactor = (radius0/radiusnow)*(radius0/radiusnow)*(radius1/radiusnow)*(radius1/radiusnow)
                        self.data[0][l2][l1] = Halfsquishfactor*OldModelFluxHDU.data[0][l2][l1]
                else:
                    self.data[0][l2][l1] = OldModelFluxHDU.data[0][l2][l1]

##
## SO notice the total normalization change is:
##    inner circle:  \pi^r0^2 ~ 3 * .0225 ~ .07
##    middle ring:   \int_{r0}^{r1}dr 2 \pi r * (r0/r) = 2\pi *(r1-r0)*r0 ~ 6 * .3 * .15 ~ 6 * .045 ~ .3
##    outer ring, 1: \int_{r1}^{middle} 2 \pi r (r0*r1)/r^2 = 2 \pi (r0*r1)*( ln(middle/r1) ) ~ 6*.045*ln(2.)~.2
##    outer ring, 2: -- a bunch of stuff around the corners ... maybe ~.1
##    SO just by squishing one reduces it by very roughly a factor of 0.6 ...
##    The inverse-squish -subtracts- this from the original;
##    Hence the inverse-squish reduces the norm by maybe 40% or so, as a guess.