%% farrot.m calculates various Faraday rotation parameters 
%% relating to the given ray paths inputted by the user (each
%% row is assumed to correspond to a ray path of interest).
%% Faraday rotation formulas from Titheridge 1972, Plant. Space
%% Science, p353-369.
%%

function [tec,farrt,farrv]=farrot(eden,lonmatm,latmatm,dismat,mph,mfs)

% eden = electron density matrix (e-/m^3)
% lonmatm = magnetic longitude (rads)
% latmatm = magnetic latitude (rads)
% dismat = radial distance matrix (km)
% mph = height at which to calculate mean magnetic parameter (km)(scalar)
% mxh = upper height up to which to calculate TEC (km)(scalar)
% mfs = magnetic field switch (0 = dipole, 1 = IGRF95) (scalar)
%
% tec = calculated TEC (TECU)
% farrt = total Faraday rotation along ray path (without constant term)
% farrv = Faraday rotation in each interval along ray path (without constant term)

re=6371;                      %Earth radi (km)
sz=size(eden);                %size of input
svx=ones(1,sz(2));
svy=ones(sz(1),1);

% Calculate the distance between points along the ray paths
[xc,yc,zc]=sph2cart(lonmatm,latmatm,dismat);                  %cartesian co-ords
dsdf=sqrt((xc(:,2:sz(2))-xc(:,1:(sz(2)-1))).^2+(yc(:,2:sz(2))-yc(:,1:(sz(2)-1))).^2+(zc(:,2:sz(2))-zc(:,1:(sz(2)-1))).^2);   %distance between points along ray paths

% Calculate the magnetic field factors
mff=magff(lonmatm,latmatm,dismat,lonmatm(:,1)*svx,latmatm(:,1)*svx,dismat(:,1)*svx,mfs);
mff(isnan(mff))=0;

% Calculate the Faraday rotation (without constant term)
mfm=0.5*(mff(:,1:(sz(2)-1))+mff(:,2:sz(2)));      %mean magnetic field matrix
edm=0.5*(eden(:,1:(sz(2)-1))+eden(:,2:sz(2)));    %mean electron density matrix
farrv=mfm.*edm.*(1e3*dsdf);                       %Faraday rotation in each interval (without constant term)
farrt=sum(farrv,2);                               %Total Faraday rotation (without constant term)

% Calculate the TEC using the mean magnetic field factor calculated at mph
mmff=intfl(dismat,mff,(re+mph)*svy,1);
tec=farrt./(1e16*mmff);
