%% magff2.m (magnetic field factor) determines the magnetic field factor
%% required for Faraday Rotation for the point of interest. It also
%% returns the angle between the ray path and the magnetic field, as
%% well as the ray path and the vertical.
%% 

function [mff,rpmf,rpvt]=magff2(lonmatm,latmatm,dismat,lonmrp,latmrp,dismrp)

% lonmatm = magnetic longitude of point of interest (rads)
% latmatm = magnetic latitude of point of interest (rads)
% dismat = radial distance of point of interest (km)
% lonmrp = initial magnetic longitude of ray path (rads)
% latmrp = initial magnetic latitude of ray path (rads)
% dismrp = initial radial distance of ray path (km)
%
% mff = magnetic field factor
% rpmf = angle between the ray path and the magnetic field (rads) 
% rpvt = angle between the ray path and the vertical (rads)

re=6371;                 %Earth radii (km)

%Convert to cartesian co-ordinates
[xmat,ymat,zmat]=sph2cart(lonmatm,latmatm,dismat);
[xmrp,ymrp,zmrp]=sph2cart(lonmrp,latmrp,dismrp);
[xmb,ymb,zmb]=sph2cart(lonmatm,latmatm,re*ones(size(dismat)));

%Determine the vector components
%Ray path
rpx=xmat-xmrp;
rpy=ymat-ymrp;
rpz=zmat-zmrp;
rpt=sqrt(rpx.^2+rpy.^2+rpz.^2);

%Vertical
vtx=xmat-xmb;
vty=ymat-ymb;
vtz=zmat-zmb;
vtt=sqrt(vtx.^2+vty.^2+vtz.^2);

%Magnetic field (from Chapter 6, Introduction to Space Physics)
bx=3*xmat.*zmat./dismat.^5;
by=3*ymat.*zmat./dismat.^5;
bz=(3*zmat.^2-dismat.^2)./dismat.^5;
bt=sqrt(bx.^2+by.^2+bz.^2);

%Determine required angles
rpmfi=(rpx.*bx+rpy.*by+rpz.*bz)./(rpt.*bt);
rpmfi(rpmfi<-1)=-1;
rpmfi(rpmfi>1)=1;
rpmf=acos(rpmfi);

rpvti=(rpx.*vtx+rpy.*vty+rpz.*vtz)./(rpt.*vtt);
rpvti(rpvti<-1)=-1;
rpvti(rpvti>1)=1;
rpvt=acos(rpvti);

%Determine magnetic field factor
bf=magdip(dismat,latmatm);
mff=bf.*cos(rpmf).*sec(rpvt);

%keyboard