class srcList:
	#arguments are:
	#sources (string, filename of LAT source list fits file in catalog format)
	#ft1 (string, filename of event file for which the xml will be used, only used to extract ROI info)
	#out (string, name of output xml file, defaults to mymodel.xml)
	def __init__(self,sources,ft1,out='mymodel.xml'):
		if not fileCheck(sources): #check that file exists
			print "Error:  %s not found." %sources
			return
		if fileCheck(out):
			print 'Warning: %s already exists, file will be overwritten if you proceed with makeModel.' %out
		self.srcs=sources
		self.out=out
		self.roi=getPos(ft1)
	
	#define a quick print function to make sure everything looks irght
	def Print(self):
		print 'Source list file: ',self.srcs
		print 'Output file name: ',self.out
		print 'Selecting %s degrees around (ra,dec)=(%s,%s)' %(self.roi[2],self.ra,self.dec)
	
	#make the xml file GDfile is galactic diffuse model and ISOfile is optional Isotropic template file
	def makeModel(self,GDfile='$(GLAST_EXT)/extFiles/v0r9/galdiffuse/gll_iem_v02.fit',GDname='GAL_v02',ISOfile='null',ISOname='Extragalactic Diffuse',model=1): #the version number may differ, but default syntax from ModelEditor didn't work
		if model!=1 and model!=2:
			print 'Incorrect specification for parameter "model", must be 1 for PowerLaw or 2 for PowerLaw2.'
			print 'Using default value of 1.'
			model=1
		print 'Creating file and adding sources'
		addSrcs(self,GDfile,GDname,ISOfile,ISOname,model)
	
import pyfits
import os
from xml.dom.minidom import parseString as pS
from numpy import floor,log10,cos,sin,arccos,pi
acos=arccos
#import ROOT #note that this is only done to turn tab completion on for functions and filenames
print "This is make1FGLxml version 03."
d2r=pi/180.

#function to cycle through the source list and add point source entries
def addSrcs(sL,GD,GDn,ISO,ISOn,Mod):
	model=open(sL.out,'w') #open file in write mode, overwrites other files of same name
	file=pyfits.open(sL.srcs) #open source list file and access necessary fields, requires LAT source catalog definitions and names
	data=file[1].data
	try:
		name=data.field('Source_Name')
	except:
		name=data.field('NickName')
	ra=data.field('RA')
	dec=data.field('DEC')
	if Mod==1:
		flux=data.field('Flux_Density')
		pivot=data.field('Pivot_Energy')
	else:
		flux=data.field('Flux1000')
	index=data.field('Spectral_Index')
	model.write('<?xml version="1.0" ?>\n')
	model.write('<source_library title="source library">\n')
	model.write('\n<!-- Point Sources -->\n')
	step=(sL.roi[2]+5)/5 #divide ROI radius plus 5 degrees into 5 steps for ordering of sources, get next highest int for ease
	i=1
	radii=[]
	while i<6:
		if i*step<=sL.roi[2]+5:
			radii+=[step*i]
		else:
			radii+=[sL.roi[2]+5] #just in case of rounding errors
		i+=1
	if Mod==1:
		for x in radii:
			if x==sL.roi[2]+5.:
				model.write('\n<!-- Sources between [%s,%s] degrees of ROI center -->\n' %(x-step,x))
			else:
				model.write('\n<!-- Sources between [%s,%s) degrees of ROI center -->\n' %(x-step,x))
			for n,f,i,r,d,p in zip(name,flux,index,ra,dec,pivot):
				dist=angsep(sL.roi[0],sL.roi[1],r,d) #check that source is within ROI radius + 5 degress of ROI center
				if (dist<x and dist>=x-step) or (x==sL.roi[2]+5. and dist==x):
					fscale=int(floor(log10(f)))
					if n[0].isdigit():
						Name='<source name="_%s" type="PointSource">\n' %(n.split(' ')[0]+n.split(' ')[1])
					else:
						Name='<source name="%s" type="PointSource">\n' %(n.split(' ')[0]+n.split(' ')[1])
					spec='\t<spectrum type="PowerLaw">\n'
					spec+='\t<!-- Source is %s degrees away from ROI center -->\n' %dist
					if(dist>sL.roi[2]): #if beyond ROI, shouldn't attempt to fit parameters
						spec+='\t<!-- Source is outside ROI, all parameters should remain fixed -->\n'
						spec+='\t\t<parameter free="0" max="1e4" min="1e-4" name="Prefactor" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
						spec+='\t\t<parameter free="0" max="5.0" min="0.0" name="Index" scale="-1.0" value="%s"/>\n' %i
					else:
						spec+='\t\t<parameter free="1" max="1e4" min="1e-4" name="Prefactor" scale="1e%i" value="%s"/>\n' %(fscale,f/10**fscale)
						spec+='\t\t<parameter free="1" max="5.0" min="0.0" name="Index" scale="-1.0" value="%s"/>\n' %i
					spec+='\t\t<parameter free="0" max="5e5" min="30" name="Scale" scale="1.0" value="%f"/>\n' %p
					spec+='\t</spectrum>\n'
					skydir='\t<spatialModel type="SkyDirFunction">\n'
					skydir+='\t\t<parameter free="0" max="360.0" min="-360.0" name="RA" scale="1.0" value="%s"/>\n' %r
					skydir+='\t\t<parameter free="0" max="90" min="-90" name="DEC" scale="1.0" value="%s"/>\n' %d
					skydir+='\t</spatialModel>\n'
					skydir+='</source>'
					(src,)=(Name+spec+skydir,)
					ptsrc=pS(src).getElementsByTagName('source')[0]
					ptsrc.writexml(model)
					model.write('\n')
	else:

		for x in radii:
			if x==sL.roi[2]+5.:
				model.write('\n<!-- Sources between [%s,%s] degrees of ROI center -->\n' %(x-step,x))
			else:
				model.write('\n<!-- Sources between [%s,%s) degrees of ROI center -->\n' %(x-step,x))
			for n,f,i,r,d in zip(name,flux,index,ra,dec):
				dist=angsep(sL.roi[0],sL.roi[1],r,d) #check that source is within ROI radius + 5 degress of ROI center
				if (dist<x and dist>=x-step) or (x==sL.roi[2]+5. and dist==x):
					val=f/1e-7
					if n[0].isdigit():
						Name='<source name="_%s" type="PointSource">\n' %(n.split(' ')[0]+n.split(' ')[1])
					else:
						Name='<source name="%s" type="PointSource">\n' %(n.split(' ')[0]+n.split(' ')[1])
					spec='\t<spectrum type="PowerLaw2">\n'
					spec+='\t<!-- Source is %s degrees away from ROI center -->\n' %dist
					if(dist>sL.roi[2]): #if beyond ROI, shouldn't attempt to fit parameters
						spec+='\t<!-- Source is outside ROI, all parameters should remain fixed -->\n'
						spec+='\t\t<parameter free="0" max="1e4" min="1e-4" name="Integral" scale="1e-07" value="%s"/>\n' %val
						spec+='\t\t<parameter free="0" max="5.0" min="0.0" name="Index" scale="-1.0" value="%s"/>\n' %i
					else:
						spec+='\t\t<parameter free="1" max="1e4" min="1e-4" name="Integral" scale="1e-07" value="%s"/>\n' %val
						spec+='\t\t<parameter free="1" max="5.0" min="0.0" name="Index" scale="-1.0" value="%s"/>\n' %i
					spec+='\t\t<parameter free="0" max="5e5" min="30" name="LowerLimit" scale="1.0" value="1e3"/>\n'
					spec+='\t\t<parameter free="0" max="5e5" min="30" name="UpperLimit" scale="1.0" value="1e5"/>\n'
					spec+='\t</spectrum>\n'
					skydir='\t<spatialModel type="SkyDirFunction">\n'
					skydir+='\t\t<parameter free="0" max="360.0" min="-360.0" name="RA" scale="1.0" value="%s"/>\n' %r
					skydir+='\t\t<parameter free="0" max="90" min="-90" name="DEC" scale="1.0" value="%s"/>\n' %d
					skydir+='\t</spatialModel>\n'
					skydir+='</source>'
					(src,)=(Name+spec+skydir,)
					ptsrc=pS(src).getElementsByTagName('source')[0]
					ptsrc.writexml(model)
					model.write('\n')
	file.close() #close files and return
	#add galactic diffuse with PL spectrum, fix index to zero for general use, those who want it to be free can unfreeze parameter manually
	model.write('\n<!-- Diffuse Sources -->\n')
	Name='\n<source name="%s" type="DiffuseSource">\n' %GDn
	spec='\t<spectrum type="PowerLaw">\n'
	spec+='\t\t<parameter free="1" max="10" min="0" name="Prefactor" scale="1" value="1"/>\n'
	spec+='\t\t<parameter free="0" max="1" min="-1" name="Index" scale="1.0" value="0"/>\n'
	spec+='\t\t<parameter free="0" max="2e2" min="5e1" name="Scale" scale="1.0" value="1e2"/>\n'
	spec+='\t</spectrum>\n'
	skydir='\t<spatialModel file="%s" type="MapCubeFunction">\n' %GD
	skydir+='\t\t<parameter free="0" max="1e3" min="1e-3" name="Normalization" scale="1.0" value="1.0"/>\n'
	skydir+='\t</spatialModel>\n'
	skydir+='</source>'
	(src,)=(Name+spec+skydir,)
	galdiff=pS(src).getElementsByTagName('source')[0]
	galdiff.writexml(model)
	model.write('\n')
	if ISO=='null': #null means no file so assume user wants an isotropic power law component
		Name='<source name="%s" type="DiffuseSource">\n' %ISOn
		spec='\t<spectrum type="PowerLaw">\n'
		spec+='\t\t<parameter free="1" max="1e3" min="1e-3" name="Prefactor" scale="1e-7" value="1"/>\n'
		spec+='\t\t<parameter free="1" max="-1.0" min="-3.5" name="Index" scale="1.0" value="-2.1"/>\n'
		spec+='\t\t<parameter free="0" max="2e2" min="5e1" name="Scale" scale="1.0" value="1e2"/>\n'
		spec+='\t</spectrum>\n'
	else: #if a file is given assume it is an isotropic template
		Name='<source name="%s" type="DiffuseSource">\n' %ISOn
		spec='\t<spectrum type="FileFunction" file="%s">\n' %ISO
		spec+='\t\t<parameter free="1" max="10" min="1e-2" name="Normalization" scale="1" value="1"/>\n'
		spec+='\t</spectrum>\n'
	#both of the above options use the same spatial model
	skydir='\t<spatialModel type="ConstantValue">\n'
	skydir+='\t\t<parameter free="0" max="10.0" min="0.0" name="Value" scale="1.0" value="1.0"/>\n'
	skydir+='\t</spatialModel>\n'
	skydir+='</source>'
	(src,)=(Name+spec+skydir,)
	iso=pS(src).getElementsByTagName('source')[0]
	iso.writexml(model)	
	model.write('\n</source_library>')
	model.close()
	return

#this function searches the header of the ft1 to find the Position keyword and extract the ra and dec values
def getPos(ft1):
	file=pyfits.open(ft1)
	num=file[1].header['NDSKEYS']
	header=file[1].header
	right='POS(RA,DEC)'
	i=1
	keynum=0
	while i<=num:  #this step is necessary since it is not clear that the POS key word will have the same number always
		word='DSTYP%i' %i
		test=file[1].header[word]
		if(test==right):
			keynum=i
			i=num
		i+=1
	if(keynum==0):  #DSKEYS start numbering at 1, if this value hasn't been updated, KEYword doesn't exist
		print 'Error: No position keyword found in fits header (assuming position is RA and DEC.  Exiting...'
		exit()
	keyword='DSVAL%i' %keynum
	ra,dec,rad=header[keyword].strip('CIRCLE()').split(',') #gets rid of the circle and parenthesis part and splits around the comma
	file.close()
	return float(ra),float(dec),float(rad)
	
#calculates the angular separation between two points on the sky
def angsep(ra1,dec1,ra2,dec2):
	ra1*=d2r
	dec1*=d2r
	ra2*=d2r
	dec2*=d2r
	diffCosine=cos(dec1)*cos(dec2)*cos(ra1-ra2)+sin(dec1)*sin(dec2)
	dC='%.10f'%diffCosine#when the source is right at the center of the roi python sometimes adds extraneous digits at the end of the value i.e. instead of 1.0
	#it returns 1.0000000000000024, which throws an error with the acos function
	return acos(float(dC))/d2r #returns values between 0 and pi radians

#Check if a given file exists or not
def fileCheck(file):
	if (not os.access(file,os.F_OK)):
		return 0
	return 1