include	<math.h>
include	"polarimetry.h"

#  CALC_STOKES -- Calculate Stokes parameters polarization, and position angle
#
#  Description:
#  ------------
#
#  Date		Author			Description
#  ----		------			-----------
#  29-Jun-1990  J.-C. Hsu		rewrite in SPP
#------------------------------------------------------------------------------

procedure calc_stokes (input, inmask, output, npix, mask, peff, pa0, paaper0,
			datafill)

pointer	input[NAPER], inmask[NAPER]
			# input: data arrays
pointer	output[NOUT]	# output: the result data arrays
int	npix		# input: number of points of the input files
bool	mask		# input: mask file existence flags
real	peff[NAPER]	# input: polarization efficiencies
real	pa0		# input: position angle offset (in degrees)
real	paaper0		# input: position angle of the 0 degree aperture 
			# (in degrees)
real	datafill	# input: fill value for bad output data pixel

int	i, j, k
real	inten, q, u, p, theta, qinst, uinst, dtheta
real	s0, s45, s90, s135, inten1, inten2

define	encode_ 91
#==============================================================================
begin

	# calculate Stokes parameters

	do j = 1, npix {

	    # if any of the input data mask is not OK, output mask is not OK 
	    if (mask) {
		do i = 1, NAPER {
		    if (Memr[inmask[i]+j-1] != DQM_OKVAL) {
			do k = 1, NOUT
		    	    Memr[output[k]+j-1] = datafill
		    	Memr[output[MASK]+j-1] = DQM_BADVAL
			goto encode_
		    }
		} 
	    }

	    # skip the data which have non-positive count rates
	    do i = 1, NAPER {
		if (Memr[input[i]+j-1] <= 0.) {
		    do k = 1, NOUT
		    	Memr[output[k]+j-1] = datafill
		    Memr[output[MASK]+j-1] = DQM_BADVAL
		    goto encode_
		}
	    }

	    # calculate the Stokes parameters
	    s0 = Memr[input[P0]+j-1] 
	    s45 = Memr[input[P45]+j-1] 
	    s90 = Memr[input[P90]+j-1] 
	    s135 = Memr[input[P135]+j-1] 
	    inten1 = s0 * peff[P90] + s90 * peff[P0]
	    inten2 = s45 * peff[P135] + s135 * peff[P45]

	    inten = inten1 / (peff[P0] + peff[P90]) + 
		    inten2 / (peff[P45] + peff[P135])
	    if (inten1 <= 0. || inten2 <= 0.) {
		do i = 1, NOUT
		    Memr[output[i]+j-1] = datafill
		Memr[output[INTEN]+j-1] = inten
		Memr[output[MASK]+j-1] = DQM_BADVAL
		goto encode_
	    }

	    qinst = (s0 - s90) / inten1
	    uinst = (s45 - s135) / inten2
	
	    dtheta = paaper0 + pa0
	    q = qinst * cos(2.*dtheta/RADIAN) - uinst * sin(2.*dtheta/RADIAN)
	    u = uinst * cos(2.*dtheta/RADIAN) + qinst * sin(2.*dtheta/RADIAN)

	    call amagr (q, u, p, 1)
	    theta = atan2(u, q) * RADIAN / 2. 
	    if (theta < 0.) 
		theta = theta + 180.   

	    # q, u, and p are in per cents
	    Memr[output[Q]+j-1] = q * 100.
	    Memr[output[U]+j-1] = u * 100.
	    Memr[output[POL]+j-1] = p * 100.
	    Memr[output[PA]+j-1] = theta
	    Memr[output[INTEN]+j-1] = inten
	    Memr[output[MASK]+j-1] = DQM_OKVAL

encode_
	    next
	}

end