%% rayppar.m determines the co-ordinates of the points located along
%% the raypaths between the inputted initial and final locations. The ray 
%% path is either broken into equally spaced intervals or log spaced,
%% depending on the input of the user.
%% Note that the inputs are converted to column vectors and row of the
%% outputs is for the corresponding column input ray path.
%%

function [lonmat,latmat,dismat]=rayppar(loni,lonf,lati,latf,alti,altf,noprp,lt)

% loni = initial longitude vector of ray paths (rads)
% lonf = final longitude vector of ray paths (rads)
% lati = initial latitude vector of ray paths (rads)
% latf = final latitude vector of ray paths (rads)
% alti = initial altitude of ray paths (km)
% altf = final altitude of ray paths (km)
% noprp = number of points along ray paths (scalar)
% lt = line spacing interval type (0 = log, other for linear) (scalar)
%
% lonmat = longitude (rads)
% latmat = latitude (rads)
% dismat = radial distance (km)


%Make sure all input are column vectors
loni=loni(:);
lonf=lonf(:);
lati=lati(:);
latf=latf(:);
alti=alti(:);
altf=altf(:);

re=6371;                      %Earth radi (km)

%% Convert from spherical to carteian co-ordinates to determine the true distance
%% between the start and end of each eay path
[xi,yi,zi]=sph2cart(loni,lati,(alti+re));         %initial points cartesian co-ords
[xf,yf,zf]=sph2cart(lonf,latf,(altf+re));         %final points cartesian co-ords
dist=sqrt((xf-xi).^2+(yf-yi).^2+(zf-zi).^2);      %total length of ray paths (km)
distx=xf-xi;                                      %x component distance (km)
disty=yf-yi;                                      %y component distance (km)
distz=zf-zi;                                      %z component distance (km)

if lt==0                                          %log scaling
 dv=logspace(-2,0,(noprp-1));                     %fractional vector to break field line up into
 dsf=[dv-dv(1) 1];                                %log distance scaling factor
else                                              %linear scaling 
 dsf=linspace(0,1,noprp);                         %linear distance scaling factor
end

sv=ones(1,noprp);     
xco=xi*sv+distx*dsf;                              %required x-coordinates
yco=yi*sv+disty*dsf;                              %required y-coordinates
zco=zi*sv+distz*dsf;                              %required z-coordinates

%% Convert back from cartesian to spherical co-ordinates
[lonmat,latmat,dismat]=cart2sph(xco,yco,zco);

%% Reset 1st and last radial displacements to agree with input (stops round off error)
dismat(:,1)=re+alti;
dismat(:,noprp)=re+altf;

