from numpy import *
from matplotlib import *
from matplotlib.pylab import *
import pyfits
import matplotlib.axes3d as axes3d

workmain = '/AstroCode/EMC2/'
thisrundir = 'Bayes-Image-Analysis_dvd_060818/'
filen='do_emc2c1_HlfSquishICSPionBremss539.0.060912_on_Pois539.X.060612.it1000.ms100gamttlcnt_1_.05FixS0.061002a.moments.fits'

HlfSquishMomentsHDULst = pyfits.open(workmain+thisrundir+filen)

print HlfSquishMomentsHDULst[0].header
print HlfSquishMomentsHDULst[0].data.shape
print HlfSquishMomentsHDULst[0].data[0]

HlfSquishMomentsHDULst.close()


## Copying from simple3d.py:
##

fig = gcf()
ax3d = axes3d.Axes3D(fig)
plt = fig.axes.append(ax3d)

## But using a different X, Y:
## I am assuming this is a 360+padding wide, 180+padding high skymap.
## The padding (used to make it even powers of 2 bins) is enough
## to make it +/- 12 degrees bigger in X; and +/- 19 degrees bigger in Y.
##
xpaddeg  , ypaddeg   = 12.0, 38.0
xrangedeg, yrangedeg = 360.+2.*xpaddeg, 180+2.*ypaddeg
print ' Xrange, Yrange in degrees:', xrangedeg, yrangedeg, '\n'

naxis1   , naxis2    = HlfSquishMomentsHDULst[0].header['NAXIS1'], HlfSquishMomentsHDULst[0].header['NAXIS2']
xdelt    , ydelt     = xrangedeg/float(naxis1), yrangedeg/float(naxis2)
print ' Xbins, Delt(degrees); Ybins, Delt(degrees):', naxis1, xdelt, naxis2, ydelt, '\n'

xcenter  , ycenter   = 0., 0.
xmin, ymin = (xcenter - (float(naxis1/2)-0.5)*xdelt), (ycenter - (float(naxis2/2)-0.5)*ydelt)
xmax, ymax = (xcenter + (float(naxis1/2)+0.5)*xdelt), (ycenter + (float(naxis2/2)+0.5)*ydelt)
print 'Xmin, Xmax; Ymin, Ymax:', xmin, xmax, ' ; ', ymin, ymax, '\n \n'

x   , y    = arange(xmin,xmax,xdelt), arange(ymin,ymax,ydelt)

#print 'X: \n',  x, '\n \n'
#print 'Y: \n', y, '\n \n'

### And here I am trying out "meshgrid":
X, Y = meshgrid(x,y)

### I want to plot three surfaces.
### For now, it will be:
### c1 = mean(data)-2.*sigma(data); c2 = mean(data); c3 = mean(data)+2.sigma(data)
### Later I will use a histogram to get
### Tom's "agrid"-like actual HPD 95.45% limits
### And plot THOSE, for each image, instead

jet()
#ax3d.plot_wireframe( X, Y, log10(HlfSquishMomentsHDULst[0].data[1]),cmap=cm.jet )
jet()
#surf = ax3d.plot_surface( X, Y, log10(HlfSquishMomentsHDULst[0].data[1] ),cmap=cm.jet)
numcontours= 512
c0 = HlfSquishMomentsHDULst[0].data[1] - 2.*HlfSquishMomentsHDULst[0].data[2]
c1 = HlfSquishMomentsHDULst[0].data[1]
c2 = HlfSquishMomentsHDULst[0].data[1] + 2.*HlfSquishMomentsHDULst[0].data[2]
#cont0 = ax3d.contour3D( X, Y, c0, numcontours ,cmap=cm.jet)
cont1 = ax3d.contour3D( X, Y, sqrt(c1), numcontours ,cmap=cm.jet)
##plot1 = ax3d.plot3D( X, Y, cmap=cm.jet)
#cont2 = ax3d.contour3D( X, Y, c2, numcontours ,cmap=cm.jet)
#contf = ax3d.contourf3D( X, Y, HlfSquishMomentsHDULst[0].data[1])
#surf.set_array(linspace(0, 1.0, naxis1))

## And figuring out these:
#scatt = ax3d.scatter3D( X, Y, c1, numcontours ,cmap=cm.jet)

show()
#--------end of test??------------------#