# Compute the radial intensity profile

import numpy
import pyfits
import math
from math import pi

fitsname = 'M33Moment0_1.fits'
outname  = 'M33RadialProfile.txt'

xc = 151
yc = 136

pixang  = 1.0/60.0	# Angular pixel scale in degrees
dist    = 0.840		# Distance in Mpc
pixphys = (pixang/360.0)*2.0*pi*dist*1000.0	# Pixel scale in kpc
pixlim  = 30.0			# Cutoff in kpc
rbinsze = 1.0			# Bin size in kpc
theta   = 3.5/60.0	# Beam size in degrees

nbins = pixlim/rbinsze

rprofile = numpy.zeros((nbins,3))	# Array to hold the computed profile

fitsfile=pyfits.open(fitsname)

image = fitsfile[0].data

nypix = image.shape[1]	# No. x pixels
nxpix = image.shape[2]	# No. y pixels

print nxpix, nypix


for xp in range(0,nxpix):
	for yp in range(0,nypix):
		r = math.sqrt((float(xp)-xc)**2.0+(float(yp)-yc)**2)
				
		rphys = r*pixphys
		
		if rphys <= pixlim:
			binn = int(rphys/rbinsze)
			print xp, yp, rphys, image[0,yp,xp]
			rprofile[binn, 0] = binn*rbinsze
			rprofile[binn, 1] = rprofile[binn, 1] + image[0,yp,xp]
			rprofile[binn, 2] = rprofile[binn, 2] + 1				# No. pixesl in this bin
		
		
outfile = open(outname,'w')

outfile.write('#r/kpc	Sum	Nvalues	Average	NHI'+'\n')

for i in range(0,nbins):
	rkpc = rprofile[i, 0]
	sumv = rprofile[i, 1]
	nval = rprofile[i, 2]
	aval = sumv/nval
	nhi  = (9.619E17*aval) / ( ((theta/2.0)**2.0)*(pi**3.0))
	
	outfile.write(str(rkpc)+' '+str(sumv)+' '+str(nval)+' '+str(aval)+' '+str(nhi)+'\n')

outfile.close()