ComputeAntennaResponse - function to compute F+,x of a detector
	characterised by detRespTensor to a source located at (theta,phi,psi) 
	in geocentric coordinates.

	theta - Polar angle of source's sky position in geocentric
                coordinates.
	phi -  Azimuthal angle of source's sky position in geocentric
               coordinates.
 	psi - Polarisation angle of gravitational waves, measured
	      counter-clockwise about direction of propagation from line of 
	      constant theta pointing to decreasing phi to positive x-axis of 
	      source coordinates.
        detRespTensor - (3,3) array containing detector response tensor.

	Returns Fp - Antenna response of detectorName to the plus GW
		polarisation of a source at (RA,Dec,psi).
		Fc - Antenna response of detectorName to the cross GW
		polarisation of a source at (RA,Dec,psi).

        Sarah Gossan 2012. Last updated 02/18/14.

def ComputeAntennaResponse(theta, phi, psi, detRespTensor):
    """ ComputeAntennaResponse - function to compute F+,x of a detector
	characterised by detRespTensor to a source located at (theta,phi,psi) 
	in geocentric coordinates.

	theta - Polar angle of source's sky position in geocentric
                coordinates.
	phi -  Azimuthal angle of source's sky position in geocentric
               coordinates.
 	psi - Polarisation angle of gravitational waves, measured
	      counter-clockwise about direction of propagation from line of 
	      constant theta pointing to decreasing phi to positive x-axis of 
	      source coordinates.
        detRespTensor - (3,3) array containing detector response tensor.

	Returns Fp - Antenna response of detectorName to the plus GW
		polarisation of a source at (RA,Dec,psi).
		Fc - Antenna response of detectorName to the cross GW
		polarisation of a source at (RA,Dec,psi).

        Sarah Gossan 2012. Last updated 02/18/14. """

    # Define axes of wave frame in terms of axes of Earth-centered axes and
    # source angles
    X = np.zeros(3)
    Y = np.zeros(3)
    X[0] = np.sin(phi)*np.cos(psi) - np.sin(psi)*np.cos(phi)*np.cos(theta)
    X[1] = -np.cos(phi)*np.cos(psi) - np.sin(psi)*np.sin(phi)*np.cos(theta)
    X[2] = np.sin(psi)*np.sin(theta)
    Y[0] = -np.sin(phi)*np.sin(psi) - np.cos(psi)*np.cos(phi)*np.cos(theta)
    Y[1] = np.cos(phi)*np.sin(psi) - np.cos(psi)*np.sin(phi)*np.cos(theta)
    Y[2] = np.sin(theta)*np.cos(psi)
    
    # Construct polarisation vectors e+,x as
    # e+ = (X \outerprod X - Y \outerprod Y)
    # ex = (X \outerprod Y + Y \outerprod X)
    ep = np.zeros((3,3),float)
    ec = np.zeros((3,3),float)
    
    for i in range(0,3):
        for j in range(0,3):
	    ep[i,j] = X[i]*X[j] - Y[i]*Y[j]
	    ec[i,j] = X[i]*Y[j] + Y[i]*X[j]

    # Construct Fp,Fc as 
    # Fp = \sum_ij(D_ij ep_ij) 
    # Fc = \sum_ij(D_ij ec_ij)
    Fp = 0.
    Fc = 0.
    for i in range(0,3):
	for j in range(0,3):
	    Fp += detRespTensor[i,j]*ep[i,j]
 	    Fc += detRespTensor[i,j]*ec[i,j]

    return Fp,Fc