# L_TO_THETA: Convert a zonal mode picture (l,t) to real space (t,theta)
# by generating the Legendre polynomials at equal increments of theta.

include <imhdr.h>

procedure t_l_to_theta()
char		inname[SZ_FNAME]	# Input file name.
char		outname[SZ_FNAME]	# Output file name.
real		dtheta		# Increment in colatitude.
real		degrees_per_radian# Conversion constant.
int		ntheta		# Number of output points in theta.
pointer	im_in			# To input image.
pointer	im_out		# To output image.
int		lmax			# Maximum degree.
pointer	pl			# To array of Legendre polynomials.
int		line			# Index to input line.
pointer	inline		# To input line data.
pointer	outline		# To output line data.
int		itheta		# Index to theta

real		clgetr()		# Cl get parameter, real.
int		clgeti()		# Cl get parameter, int.
pointer	immap()		# Image open.
pointer	imgl2r()		# Image get line real.
pointer	impl2r()		# Image put line real.
real		adotr()		# Dot product of two vectors, real.
begin
# Cl input.
	call clgstr ("input", inname, SZ_FNAME)
	call clgstr ("output", outname, SZ_FNAME)
	degrees_per_radian = 45. / atan (1.0)
	dtheta = clgetr ("dtheta") / degrees_per_radian
	ntheta = clgeti ("ntheta")

# Setup Legendre polynomial image.
	im_in = immap (inname, READ_ONLY, 0)
	lmax = IM_LEN(im_in,1) - 1
	call malloc (pl, (lmax+1)*ntheta, TY_REAL)
	call plgen (Memr[pl], lmax+1, ntheta, dtheta)

# Setup output image.
	im_out = immap (outname, NEW_COPY, im_in)
	IM_LEN(im_out,1) = ntheta

# Do the calculation.
	for (line=1; line<=IM_LEN(im_in,2); line=line+1)
	{	inline = imgl2r (im_in, line)
		outline = impl2r (im_out, line)
		for (itheta=0; itheta<ntheta; itheta=itheta+1)
			Memr[outline+itheta] = adotr (Memr[inline],
					Memr[pl+itheta*(lmax+1)], lmax+1)
	}

# Clean up.
	call imunmap (im_in)
	call imunmap (im_out)
	call mfree (pl, TY_REAL)
end
# *********************************************************
procedure plgen (ar, nl, ntheta, dtheta)
int		nl			# Number of degrees (lmax+1)
int		ntheta		# Number of colatitude values.
real		ar[nl,ntheta]	# Output array of coefficients.
real		dtheta		# Colatitude increment (radians)

int		lmax			# Maximum degree.
int		itheta		# Index in the theta direction
real		ctheta		# Cos(current theta).
real		pl_minus_1		# Value of pl in previous iteration.
real		pl			# Current value of pl.
int		l			# Index for degree.
begin
	lmax = nl - 1

	for (itheta=1; itheta<=ntheta; itheta=itheta+1)
	{	ctheta = cos (dtheta * (itheta-1))
		pl_minus_1 = 1
		pl	     = ctheta
		ar[1,itheta] = pl_minus_1
		ar[2,itheta] = pl
		for (l=1; l<=lmax-1; l=l+1)
		{	ar[l+2,itheta] = ((2*l+1)*ctheta*pl - l*pl_minus_1) / (l+1)
			pl_minus_1 = pl
			pl = ar[l+2,itheta]
			ar[l+2,itheta] = ar[l+2,itheta] * sqrt ((2.*l+3)/2)
		}
		ar[1,itheta] = ar[1,itheta] * sqrt (1/2.)
		ar[2,itheta] = ar[2,itheta] * sqrt (3/2.)
	}
end