function [ux,uy,uz]=calc_yang(theta,strike,nu,P,a,b,z0,x,y)
% Calculate the radar range displacement for a SPHEROIDAL pressure source
% in an elastic halfspace (based on Yang et al., vol 93, JGR, 4249-4257, 1988)
% with arbitrary plunge (th_src) of the long axis of the spheroid (th_src = 90,
% prolate; th_src = 0, oblate).
%
% Inputs: theta: dip angle of source
%             P: pressure change in magma chamber (used inPstar_dila.m) 
%             a: semimajor axis of spheroid
%             b: semiminor axis of spheriod
%            z0: depth of source
%             x: x location of point on surface
%             y: y location of point on surface
%        phisar: look direction of radar
%      thetasar: look angle (from zenith) of radar
% Output: rd: calculated range displacement
% NOTE: the x, y locations assume a source at origin
% NOTE: In the expressions of Yang et al. the spheroid dips parallel to the y axis
%       at x=0. We will assume (initially) that it is also centered at y=0.

mu=1;
matrl(1)=2*nu/(1-2*nu)*mu;
matrl(2)=mu;
matrl(3)=nu;

% Store some commonly used parameters
coeffs(1)=1/(16*mu*(1-nu));
coeffs(2)=3-4*nu;
coeffs(3)=4*(1-nu)*(1-2*nu);

%strike = -strike* pi/180 + pi/2;
strike = -strike* pi/180;
theta = theta*pi/180;
coss  = cos(strike);
sins  = sin(strike);
e_theta(1)=sin(theta);
e_theta(2)=cos(theta);
rotx=x*coss+y*sins;
roty=-x*sins+y*coss;

% Evaluate integrals at xi = +c to -c:
c=sqrt(a^2-b^2);
% Speroid quantities
phi=0;
[sph]=spheroid(a,b,c,matrl,phi,theta,P);
xi=c;
[Up1,Up2,Up3]=f_yang(sph,xi,z0,rotx,roty,matrl,e_theta,coeffs);
xi=-xi;
[Um1,Um2,Um3]=f_yang(sph,xi,z0,rotx,roty,matrl,e_theta,coeffs);
uxj=-Up1+Um1;
uyj=-Up2+Um2;
uz=Up3-Um3;

ux=-uyj*sins+uxj*coss;
uy=uxj*sins+uyj*coss;
