'''
@brief Reads in text image files; writes out moments in fits form.
@author A. Connors
'''
# @file get_multi_txt_matrix.py
#
# This script is explicitly for reading in "output.txt" format files
# and writing them out in fits form compatible with
# CHASC imaging software.
#

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

class ScanTxtMatrix:
    def __init__(self,InTxtFileName,naxis_3tuple):
        print 'Reading ',InTxtFileName,' with shape:', naxis_3tuple
        naxis3,naxis2,naxis1 = naxis_3tuple
#        print naxis3,naxis2,naxis1
        image = zeros((naxis3,naxis2,naxis1),dtype=float)
        self.data = image
        image_in = open(InTxtFileName,'r')
        for numimage in range(naxis3):
            for lon in range(naxis2):
                thisline = image_in.readline()
#                print thisline
                latbegin = 0
                latend = latbegin + 1
                for lat in range(naxis1):
                    while ( (latbegin < len(thisline)) and (thisline[latbegin] == ' ') ):
                        latbegin = latbegin + 1
                        latend = latbegin + 1
                    while ((latend < len(thisline)) and (thisline[latend] != ' ')):
                        latend = latend + 1
#                    print 'lon, lat: ',lon, lat, 'latbegin, end: ', latbegin, latend
                    yo = thisline[latbegin:latend]
                    if(yo != ''):
                        image[numimage][lon][lat] = float(yo)
                    else:
                        print '\n yo-string:',yo,' is blank at lon,lat: ',lon, lat, '\n'
                        image[numimage][lon][lat] = 0.
#                    print image[numimage][lon][lat]
                    latbegin = latend
                    latend = latbegin + 1
                    if (latend > len(thisline) -1):
                        latend = len(thisline) -1
#                print numimage,lon, '\n', image[numimage][lon]
        self.data = image
#        print image
#        print self.data
        image_in.close()
#etc

class MatrixModeNMoments:
    def __init__(self,CleanImage,ModeIterNum,MatchedExposrImageData,OutFitsFileName):
        nmoments = 7
        naxis3,naxis2,naxis1 = CleanImage.shape
        print 'CleanImage.shape: ', CleanImage.shape
        self.data    = zeros((nmoments,naxis2,naxis1),dtype=float)
        self.data[0] = average(CleanImage,axis=0)
        print 'self.data[0][naxis2/2][naxis1/2]: ',self.data[0][naxis2/2][naxis1/2]
        dif     = CleanImage - self.data[0]
        self.data[1] = sqrt( average( (dif*dif) , axis=0) ) + 1.0e-7
        print 'self.data[1][naxis2/2][naxis1/2]: ',self.data[1][naxis2/2][naxis1/2]
        self.data[2] = average(( dif*dif*dif), axis=0)/self.data[1]/self.data[1]/self.data[1] 
        print 'self.data[2][naxis2/2][naxis1/2]: ',self.data[2][naxis2/2][naxis1/2]
        self.data[3] = average( (dif*dif*dif*dif), axis=0)/self.data[1]/self.data[1]/self.data[1]/self.data[1]
        print 'self.data[3][naxis2/2][naxis1/2]: ',self.data[3][naxis2/2][naxis1/2]
        self.data[4] = self.data[0]/self.data[1]
        print 'self.data[4][naxis2/2][naxis1/2]: ',self.data[4][naxis2/2][naxis1/2]
        self.data[5] = self.data[0]/MatchedExposrImageData
        print 'self.data[5][naxis2/2][naxis1/2]: ',self.data[5][naxis2/2][naxis1/2]
        self.data[6] = self.data[1]/MatchedExposrImageData
        print 'self.data[6][naxis2/2][naxis1/2]: ',self.data[6][naxis2/2][naxis1/2]
    
        image_HDU = pyfits.PrimaryHDU(self.data)
        image_HDU.header.update('NAXIS3',nmoments,comment='Mean,Sigma,Skew,Kurt,Mean/Sig,Mean/Exposr,Sigm/Exposr')
    
        image_HDU.writeto(OutFitsFileName)
#
#   End.