function [MTR] = get_MTR(Ex1,Ey1,Hx1,Hy1,Hz1,Ex2,Ey2,Hx2,Hy2,Hz2,freq)

% Magnetotelluric response calculation
% Return: impedance, apparent resistivity, phase,  tipper response

mu0 = 4*pi*1e-7 ;  % V·s/ (A·m)


omega = 2*pi*freq;
detH  = (Hx1.*Hy2-Hx2.*Hy1);

%% impedance
Zxx  = (Ex1.*Hy2 - Ex2.*Hy1)./detH ;
Zxy  = (Ex2.*Hx1 - Ex1.*Hx2)./detH ;
Zyx  = (Ey1.*Hy2 - Ey2.*Hy1)./detH ;
Zyy  = (Ey2.*Hx1 - Ey1.*Hx2)./detH ;

MTR.Zxx = Zxx;
MTR.Zxy = Zxy;
MTR.Zyx = Zyx;
MTR.Zyy = Zyy;

%% apparent resistivity
rhoXX  =  1/(omega*mu0)*abs(Zxx).^2;
rhoXY  =  1/(omega*mu0)*abs(Zxy).^2;
rhoYX  =  1/(omega*mu0)*abs(Zyx).^2;
rhoYY  =  1/(omega*mu0)*abs(Zyy).^2;

MTR.rhoXX = rhoXX;
MTR.rhoXY = rhoXY;
MTR.rhoYX = rhoYX;
MTR.rhoYY = rhoYY;
%% phase
phaXX   =  -atan(imag(Zxx)./real(Zxx))*180/pi;
phaXY   =  -atan(imag(Zxy)./real(Zxy))*180/pi;
phaYX   =  -atan(imag(Zyx)./real(Zyx))*180/pi;
phaYY   =  -atan(imag(Zyy)./real(Zyy))*180/pi;

% phaXX   =  -atan2(imag(Zxx),real(Zxx))*180/pi;
% phaXY   =  -atan2(imag(Zxy),real(Zxy))*180/pi;
% phaYX   =  -atan2(imag(Zyx),real(Zyx))*180/pi;
% phaYY   =  -atan2(imag(Zyy),real(Zyy))*180/pi;
MTR.phaXX = phaXX;
MTR.phaXY = phaXY;
MTR.phaYX = phaYX;
MTR.phaYY = phaYY;

%% Tipper
Tzx       = (Hz1.*Hy2 - Hz2.*Hy1)./detH ;
Tzy       = (Hz2.*Hx1 - Hz1.*Hx2)./detH ;
MTR.Tzx = Tzx;
MTR.Tzy = Tzy;

end