%% magff.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.
%% 
%% NOTE - returns NaN if the point of interest and initial location
%%        are the same.
%%
%% 16/10/00 - modify so that IGRF95 magnetic field can be used, rather 
%%            than just the dipole. 
%%

function [mff,rpmf]=magff(lonmatm,latmatm,dismat,lonmrp,latmrp,dismrp,mfs)

% 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)
% mfs = magnetic field switch (0 = dipole, 1 = IGRF95) (scalar)
%
% mff = magnetic field factor
% rpmf = angle between the ray path and the magnetic field (rads) 

re=6371;                 %Earth radii (km)

%Convert to cartesian co-ordinates
if mfs==1
 [lonmat,latmat]=mag2geo3(lonmatm,latmatm);
 [lonrp,latrp]=mag2geo3(lonmrp,latmrp);
 [xmat,ymat,zmat]=sph2cart(lonmat,latmat,dismat);
 [xmrp,ymrp,zmrp]=sph2cart(lonrp,latrp,dismrp);
else
 [xmat,ymat,zmat]=sph2cart(lonmatm,latmatm,dismat);
 [xmrp,ymrp,zmrp]=sph2cart(lonmrp,latmrp,dismrp);
end

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

if mfs==1
 [bx,by,bz]=modigrf95(lonmat,latmat,dismat);
else
 %Dipole 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;
end
bt=sqrt(bx.^2+by.^2+bz.^2);

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

%Determine magnetic field factor
if mfs==1
 mff=bt.*cos(rpmf);
else
 bf=magdip(dismat,latmatm);
 mff=bf.*cos(rpmf);
end

%Make sure answers are positive only
mff=abs(mff);