%this function will determine the best fitting tilt direction and amount by
%minimizing energy on the vertical component.
%input are X(east), Y(north) and Z(vertical).
%output are Znew (tilt corrected Z), Xnew (rotated X), Ynew (rotated Y),
%angy (azimuthal rotation in degrees) and vrt (tilt angle from vertical in
%degrees)
%Note this operates on the second half of the time series-can be modified
%as need to operate on whole time series, etc by modifying inline function
%ff below.
%requires optimization toolbox

function [Znew,angy,vrt,Xnew,Ynew] = tilt_corr_rot(X,Y,Z);
         
%% first do polarization analysis on the data to get a guess of the tilt direction.
         S=cov([Z(end-length(Z)/2:end) X(end-length(Z)/2:end) Y(end-length(Z)/2:end)]);
         [V,D]=eig(S);
         linny=1-(D(2,2)+D(1,1))/2/D(3,3);
         planny=1-2*D(1,1)/(D(2,2)+D(3,3));
         if linny>0.1 %this is a relatively linear phase recorded
             %calculate the rotation angle 
             angy=atan(V(3,3)*sign(V(1,3))/(V(2,3)))*180/pi;
         end
%% search for the best fitting azimuth (angy) and tilt (vrt)         
         [~,idd]=max(abs(V(1,:)));
         vrt=acosd(V(1,idd));
         ff=@(x)Zrot(Z(end-length(Z)/2:end),X(end-length(Z)/2:end),Y(end-length(Z)/2:end),x);
         options=optimset('Display','iter');
         warning off
         [vrt,zvr2]=fminsearch(ff,[vrt;angy]);
         angy=vrt(2);
         vrt=vrt(1);
         [Znew,Ynew,Xnew]=xyzrot(X,Y,Z,[vrt;angy]);
    
end
function zvr = Zrot(Z,X,Y,vrt)
    [Znew,ynew,xnew]=xyzrot(X,Y,Z,vrt);
   zvr=sqrt(Znew(:)'*Znew(:));
end
function [znew,ynew,xnew]=xyzrot(x,y,z,tst)
vrt=tst(1);angy=tst(2);
dipZ=-90+vrt;%sin change on vrt due to how it is determined
dipY=vrt*cosd(angy);
dipX=vrt*sind(angy);
%make vectors from dip and azimuth
zv=[-sind(dipZ);cosd(dipZ).*cosd(angy);sind(angy).*cosd(dipZ)];
yv=[-sind(dipY);cosd(dipY).*cosd(0);sind(0).*cosd(dipY)];
xv=[-sind(dipX);cosd(dipX).*cosd(90);sind(90).*cosd(dipX)];
C=[zv yv xv];
jk=inv(C)*[z y x]';
znew=jk(1,:)';
ynew=jk(2,:)';
xnew=jk(3,:)';
end