#pragma rtGlobals=1		// Use modern global access method.
#ModuleName=NJ
// April 5, 2006
// Christopher F. Neese 
// Christopher_Neese@world.oberlin.edu
//
// These routines are based on the Fortran code in the Appendix of "Angular Momentum" by Richard N. Zare.
// ISBN 0-471-85892-7
//
// ThreeJ and SixJ have been tested by comparing results with results from same functions in Mathematica.
// NineJ needs to be tested, as I have not had occasion to use it yet.
// Since the contractuon formula used to evaluate the NineJ is fairly simple, it should be OK.  


static function MakeLogFactorial()
	WAVE/Z lFact = root:logFactorial
	if (!WaveExists(lFact))
		Print "Creating wave root:logFactorial"
		Make/N=400 $"root:logFactorial"
		WAVE lFact = root:logFactorial
		lFact[0]=0
		lFact[1,Inf]=lFact[p-1]+log(p)
	endif
end

function ThreeJ(J1,M1,J2,M2,J3,M3)
	variable J1,M1,J2,M2,J3,M3
	
	MakeLogFactorial()
	WAVE lFact = root:logFactorial

	M1 = 0.5*Round(2*M1)
	M2 = 0.5*Round(2*M2)
	M3 = 0.5*Round(2*M3)

	if(M1+M2+M3)
		return 0
	else
		variable twoJ1=Round(2*J1)
		variable twoJ2=Round(2*J2)
		variable twoJ3=Round(2*J3)
		if( (twoJ3 <= (twoJ1+twoJ2)) && (twoJ3 >= abs(twoJ1-twoJ2)) )
			J1 = 0.5*twoJ1
			J2 = 0.5*twoJ2
			J3 = 0.5*twoJ3
	
			variable S=J1+J2+J3
			variable S1=J1+J2-J3
			variable S2=J1-M1
			variable S3=J2+M2
			variable S4=J3-J2+M1
			variable S5=J3-J1-M2
	
			variable A=lFact(s-twoJ1)+lFact(s-twoJ2)+lFact(s-twoJ3)+lFact(J1+M1)+lFact(S2)+lFact(S3)+lFact(J2-M2)+lFact(J3+M3)+lFact(J3-M3)-lFact[s+1]
			A *= 0.5
			variable numax = min(min(S1,S2),S3)
			variable numin = max(-min(S4,S5),0)
			variable nu=0, B=0, C=0
			for (nu=numin ; nu <= numax ; nu += 1)
				C = lFact(nu)+lFact(S1-nu)+lFact(S2-nu)+lFact(S3-nu)+lFact(S4+nu)+lFact(S5+nu)
				B += (mod(nu,2) ? 1: -1)*10^(A-C)
			endfor
			return (mod(J1-J2-M3,2) ? 1: -1)*B
		else
			return 0
		endif
	endif
end

// { j1 j2 j3 }
// { j4 j5 j6 }

function SixJ(j1,j2,j3,j4,j5,j6)
	variable j1,j2,j3,j4,j5,j6

	MakeLogFactorial()
	WAVE lFact = root:logFactorial

	j1 = round(2*j1)
	j2 = round(2*j2)
	j3 = round(2*j3)
	j4 = round(2*j4)
	j5 = round(2*j5)
	j6 = round(2*j6)

	variable jm1 = (j1+j2+j3)/2
	variable jm2 = (j1+j5+j6)/2
	variable jm3 = (j4+j2+j6)/2
	variable jm4 = (j4+j5+j3)/2
	variable jm = max(jm1,max(jm2,max(jm3,jm4)))

	variable jx1 = (j1+j2+j4+j5)/2
	variable jx2 = (j2+j3+j5+j6)/2
	variable jx3 = (j3+j1+j6+j4)/2
	variable jx = min(jx1,min(jx2,jx3))

	if (jm > jx)
		return 0
	endif
		
	variable term1 = dl(j1,j2,j3,lFact) + dl(j1,j5,j6,lFact) + dl(j4,j2,j6,lFact) + dl(j4,j5,j3,lFact)
	variable term2 = 0
	variable res = 0

	variable i
	for (i = jm ; i <= jx ; i+=1)
		term2 = lFact(i+1)
		term2 += -lFact(i-jm1)-lFact(i-jm2)-lFact(i-jm3)-lFact(i-jm4)
		term2 += -lFact(jx1-i)-lFact(jx2-i)-lFact(jx3-i)
		
		res += 10^(term1+term2)*((-1)^i)
	endfor
	
	return res
end

static function dl(j1,j2,j3,lFact)
	variable j1,j2,j3
	wave lFact
	
	return (lFact((j1+j2-j3)/2)+lFact((j2+j3-j1)/2)+lFact((j3+j1-j2)/2)-lFact(1+(j1+j2+j3)/2))/2
end

// { j1 j2 j3 }
// { j4 j5 j6 }
// { j7 j8 j9 }

// Based on Eq. 4.24 in Zare
function NineJ(j1,j2,j3,j4,j5,j6,j7,j8,j9)
	variable	j1,j2,j3,j4,j5,j6,j7,j8,j9
	
	j1 = 0.5*round(2*j1)
	j2 = 0.5*round(2*j2)
	j3 = 0.5*round(2*j3)
	j4 = 0.5*round(2*j4)
	j5 = 0.5*round(2*j5)
	j6 = 0.5*round(2*j6)
	j7 = 0.5*round(2*j7)
	j8 = 0.5*round(2*j8)
	j9 = 0.5*round(2*j9)
	
	// kmin and kmax are determined by the triangle rules inherent in the SixJ symbols below
	variable kmin = max(abs(j1-j9),max(abs(j8-j4),abs(j6-j2)))
	variable kmax = min(j1+j9, min(j8+j4, j6+j2)) 
	
	variable res=0
	variable s1, s2, s3, p, k
	if (kmin <= kmax)
		for (k = kmin ; k <= kmax ; k += 1)
			s1 = sixj(j1, j4, j7, j8, j9, k)
			s2 = sixj(j2, j5, j8, j4, k, j6)
			s3 = sixj(j3, j6, j9, k, j1, j2)
			p = (2*k+1)*(-1)^(2*k)
			res += p*s1*s2*s3
		endfor
	endif
	
	return res
end