import string, math
from math import *

def moffat_int_25(beta, limit1, limit2, Seeing):
	#print alpha, beta, limit1, limit2, Seeing
	alpha=Seeing/(2.*sqrt(pow(2,(1/beta))-1))
	integral=2./3.*alpha*alpha*(alpha*limit2+atan(limit2/alpha)*alpha*alpha+atan(limit2/alpha)*limit2*limit2)/(alpha*alpha+limit2*limit2)
	integral=integral-2./3.*alpha*alpha*(alpha*limit1+atan(limit1/alpha)*alpha*alpha+atan(limit1/alpha)*limit1*limit1)/(alpha*alpha+limit1*limit1)
	gamma=(beta-1)/(pi*alpha*alpha)
	integral=integral*gamma
	return(integral)

def moffat_int_35(beta, limit1, limit2, Seeing):
	#print alpha, beta, limit1, limit2, Seeing
	alpha=Seeing/(2.*sqrt(pow(2,(1/beta))-1))
	integral=2./15.*alpha*alpha*(5*pow(alpha,3)*limit2+3*alpha*pow(limit2,3)+3*pow(alpha,4)*atan(limit2/alpha)+6*alpha*alpha*atan(limit2/alpha)*limit2*limit2+3*atan(limit2/alpha)*pow(limit2,4))/((alpha*alpha+limit2*limit2)*(alpha*alpha+limit2*limit2))
	integral=integral-2./15.*alpha*alpha*(5*pow(alpha,3)*limit1+3*alpha*pow(limit1,3)+3*pow(alpha,4)*atan(limit1/alpha)+6*alpha*alpha*atan(limit1/alpha)*limit1*limit1+3*atan(limit1/alpha)*pow(limit1,4))/((alpha*alpha+limit1*limit1)*(alpha*alpha+limit1*limit1))
	gamma=(beta-1)/(pi*alpha*alpha)
	integral=integral*gamma
	return(integral)

def moffat(x, y, beta, Seeing):
	alpha=Seeing/(2.*sqrt(pow(2,(1/beta))-1))
	moffat=(beta-1)/(pi*alpha*alpha)*pow((1+(x*x+y*y)/(alpha*alpha)), (-beta))
	return(moffat)