function [U,E,T]=okada(ALP,xrec,source)
%SUPERCEDED by C-code version
%OKADA   Calculates surface displacements, strains, and tilts due
%        to a rectangular fault in an elastic half-space.  By Y. Okada.
%
%        Note:  If U is the only argument out, the function executes about
%        twice as fast.
%
%INPUTS
%  alp:     Medium constant, mu/(lamda+mu)
%  xrec:    Station coordinates -- represented as columns in 2xn matrix
%  source:  fault parameters
%
%           source(1)      = Depth
%           source(2),(3)  = Fault length range
%           source(4),(5)  = Fault width range
%           source(6),(7)  = Sin and cos of dip-angle
%           source(8)      = Stike-slip dislocation
%           source(9)      = Dip-slip dislocation
%           source(10)     = Tensile dislocation 
%
%OUTPUTS
%  U:       Displacement   (units of dislocation)
%  E:       Strain         (units of dislocation/units of coordinates)
%  T:       Tilt           (units of dislocation/units of coordinates)


%Establish variables, define constants (semis to commas 3/29/97)

	n=size(xrec,2);
	XX=xrec(1,:);
	YY=xrec(2,:);
   DEP=source(1);
   AL=[source(2),source(3)];
	AW=[source(4),source(5)];
	SD=source(6);
	CD=source(7);
	DISL=[source(8),source(9),source(10)];
	PI2=2*pi;
	U=zeros(3,n);
	E=zeros(4,n);
	T=zeros(2,n);

%Begin STRECTF

	P=YY.*CD+DEP*SD;
	Q=YY.*SD-DEP*CD;

	for k=1:2
		
		ET=P-AW(k);
		
		for j=1:2
		
			XI=XX-AL(j);
			SIGN=(-1)^(k+j);
			
			%Begin SRECTG
	
				u=zeros(3,n);
				e=zeros(4,n);
				t=zeros(2,n);	
			
				XI2=XI.*XI;
				ET2=ET.*ET;
				Q2=Q.*Q;
				R2=XI2+ET2+Q2;
				R=sqrt(R2);
				R3=R.*R2;
				D=ET.*SD-Q.*CD;
				Y=ET.*CD+Q.*SD;
				RET=R+ET;
				RET=RET.*(RET>0);			
				RD=R+D;
				TT=0;
				QSTAR=Q+(Q==0)*1;
				TT=(Q~=0).*atan((XI.*ET)./(QSTAR.*R));
				RE=0;
				RETSTAR=RET+(RET==0)*1;
				RE=(RET~=0).*(1./RETSTAR);
				DLE=(RET~=0).*log(RETSTAR)+(RET==0).*(-log(R-ET));
				RRX=1./(R.*(R+XI));
				RRE=RE./R;

				if nargout>=2

					RRD=1./(R.*RD);
					AXI=(2.*R+XI).*RRX.*RRX./R;
					AET=(2.*R+ET).*RRE.*RRE./R;
				
				end

				if CD==0 %Vertical fault

					RD2=RD.*RD;
					A1=-ALP./2.*XI.*Q./RD2;	
					A3=ALP./2.*(ET./RD+Y.*Q./RD2-DLE);
					A4=-ALP.*Q./RD;
					A5=-ALP.*XI.*SD./RD;
				
					if nargout>=2

						B1=ALP./2.*Q./RD2.*(2.*XI2.*RRD-1)
						B2=ALP./2.*XI.*SD./RD2.*(2.*Q2.*RRD-1)
						C1=ALP.*XI.*Q.*RRD./RD
						C3=ALP.*SD./RD.*(XI2.*RRD-1)
				
					end					

				else	%Dipping fault

					TD=SD./CD;
					X=sqrt(XI2+Q2);
					A5=0;
					XISTAR=XI+(XI==0)*1;
					A5=(XI~=0).*(ALP.*2./CD.*atan((ET.*(X+Q.*CD)+X.*(R+X).*SD)./(XISTAR.*(R+X).*CD)));
					A4=ALP./CD.*(log(RD)-SD.*DLE);
					A3=ALP.*(Y./RD./CD-DLE)+TD.*A4;
					A1=-ALP./CD.*XI./RD-TD.*A5;
				
					if nargout>=2

						C1=ALP./CD.*XI.*(RRD-SD.*RRE);
						C3=ALP./CD.*(Q.*RRE-Y.*RRD);
						B1=ALP./CD.*(XI2.*RRD-1)./RD-TD.*C3;
						B2=ALP./CD.*XI.*Y.*RRD./RD-TD.*C1;				
					end

				end

				A2=-ALP.*DLE-A3;
				
				if nargout>=2

					B3=-ALP.*XI.*RRE-B2;
					B4=-ALP.*(CD./R+Q.*SD.*RRE)-B1;					
					C2=ALP.*(-SD./R+Q.*CD.*RRE)-C3;
				
				end

				if DISL(1)~=0

					%Strike-slip contribution		

					UN=DISL(1)./PI2;
					REQ=RRE.*Q;
					u(1,:)=u(1,:)-UN.*(REQ.*XI+TT+A1.*SD);
					u(2,:)=u(2,:)-UN.*(REQ.*Y+Q.*CD.*RE+A2.*SD);
					u(3,:)=u(3,:)-UN.*(REQ.*D+Q.*SD.*RE+A4.*SD);

					if nargout>=2

						e(1,:)=e(1,:)+UN.*(XI2.*Q.*AET-B1.*SD);
						e(2,:)=e(2,:)+UN.*(XI2.*XI.*(D./(ET2+Q2)./R3-AET.*SD)-B2.*SD);
						e(3,:)=e(3,:)+UN.*(XI.*Q./R3.*CD+(XI.*Q2.*AET-B2).*SD);
						e(4,:)=e(4,:)+UN.*(Y.*Q./R3.*CD+(Q.*SD.*(Q2.*AET-2.*RRE)-(XI2+ET2)./R3.*CD-B4).*SD);
					
					end
				
					if nargout==3

						t(1,:)=t(1,:)+UN.*(-XI.*Q2.*AET.*CD+(XI.*Q./R3-C1).*SD);
						t(2,:)=t(2,:)+UN.*(D.*Q./R3.*CD+(XI2.*Q.*AET.*CD-SD./R+Y.*Q./R3-C2).*SD);
				
					end
				end

				if DISL(2)~=0
		
					%Dip-slip contribution

					UN=DISL(2)./PI2;
					SDCD=SD.*CD;

					u(1,:)=u(1,:)-UN.*(Q./R-A3.*SDCD);
					u(2,:)=u(2,:)-UN.*(Y.*Q.*RRX+CD.*TT-A1.*SDCD);
					u(3,:)=u(3,:)-UN.*(D.*Q.*RRX+SD.*TT-A5.*SDCD);
				
					if nargout>=2

						e(1,:)=e(1,:)+UN.*(XI.*Q./R3+B3.*SDCD);
						e(2,:)=e(2,:)+UN.*(Y.*Q./R3-SD./R+B1.*SDCD);
						e(3,:)=e(3,:)+UN.*(Y.*Q./R3+Q.*CD.*RRE+B1.*SDCD);
						e(4,:)=e(4,:)+UN.*(Y.*Y.*Q.*AXI-(2.*Y.*RRX+XI.*CD.*RRE).*SD+B2.*SDCD);

					end

					if nargout==3

						t(1,:)=t(1,:)+UN.*(D.*Q./R3+Q.*SD.*RRE+C3.*SDCD);
						t(2,:)=t(2,:)+UN.*(Y.*D.*Q.*AXI-(2.*D.*RRX+XI.*SD.*RRE).*SD+C1.*SDCD);
					end

				end

				if DISL(3)~=0

					%Tensile-fault contribution

					UN=DISL(3)./PI2;
					SDSD=SD.*SD;

					u(1,:)=u(1,:)+UN.*(Q2.*RRE-A3.*SDSD);
					u(2,:)=u(2,:)+UN.*(-D.*Q.*RRX-SD.*(XI.*Q.*RRE-TT)-A1.*SDSD);
					u(3,:)=u(3,:)+UN.*(Y.*Q.*RRX+CD.*(XI.*Q.*RRE-TT)-A5.*SDSD);
				
					if nargout>=2

						e(1,:)=e(1,:)-UN.*(XI.*Q2.*AET+B3.*SDSD);
						e(2,:)=e(2,:)-UN.*(-D.*Q./R3-XI2.*Q.*AET.*SD+B1.*SDSD);
						e(3,:)=e(3,:)-UN.*(Q2.*(CD./R3+Q.*AET.*SD)+B1.*SDSD);
						e(4,:)=e(4,:)-UN.*((Y.*CD-D.*SD).*Q2.*AXI-2.*Q.*SD.*CD.*RRX-(XI.*Q2.*AET-B2).*SDSD);
					
					end

					if nargout==3

						t(1,:)=t(1,:)-UN.*(Q2.*(SD./R3-Q.*AET.*CD)+C3.*SDSD);
						t(2,:)=t(2,:)-UN.*((Y.*SD+D.*CD).*Q2.*AXI+XI.*Q2.*AET.*SD.*CD-(2.*Q.*RRX-C1).*SDSD);

					end
				end


			U=U+SIGN.*u;
			E=E+SIGN.*e;
			T=T+SIGN.*t;

			end

	end









					
					
					
							

			