%% datgrid2.m takes the vectors describing the longitude,
%% latitudes, and altitude of the observer and orbit.
%% It produces an output matrix of longitude, latitude
%% and distance (from the centre of the Earth) between
%% the observer and the orbit. If the altitude of the
%% orbit is less than 1500 km, or the number of points
%% of interest in each less than 100 then each ray path is
%% broken up into equal space sections ie the output is
%% the same as datgrid.m. Otherwise, the region of interest
%% grid is broken up into three regions.
%% 1) From 0 to 500 km distance along the ray path the
%%    separation will be set to 10 km increments.
%% 2) From 500 to 1500 km, 20 km increments will be used.
%% 3) While the remaining distance from 1500 km to the orbit
%%    will be broken up into equal distance increments, the
%%    number equal to the number of remaining points of
%%    interest in the ray path.
%%

function [lonmat,latmat,dismat,dist]=datgrid2(lonob,latob,lontg,lattg,alt,nae)

% lonob = observers longitude (rads)
% latob = observers latitude (rads)
% lonob = targets longitude (rads)
% latob = targets latitude (rads)
% alt = altitude (km)
% nae = number of points of interest on ray path
% lonmat = required longitude values matrix
% latmat = required latitude values matrix
% dismat = required displacement matrix, the distance from the centre of the Earth (km)
% dist = distance between observer and orbit, required by teccal.m 

re=6371;                   %Earth's radi (km)

if alt<1500 | nae<100

 % Convert from spherical to cartesian coordinates. This is necessary so that equal
 % distance between sample points along a given projection can be obtained.
 [xcoob,ycoob,zcoob]=sph2cart(lonob,latob,re*ones(size(lonob)));        %convert observer values       
 [xcotg,ycotg,zcotg]=sph2cart(lontg,lattg,(re+alt)*ones(size(lontg)));  %convert target values  

 dist=sqrt((xcoob-xcotg).^2+(ycoob-ycotg).^2+(zcoob-zcotg).^2);         %distance between observer and target longitude points
 distx=xcotg-xcoob;                                                     %x component distance
 disty=ycotg-ycoob;                                                     %y component distance
 distz=zcotg-zcoob;                                                     %z component distance  

 xco=[ones(nae-1,1)*xcoob+(0:(nae-2))'*(distx/(nae-1));xcotg];          %reqired x coordinates matrix
 yco=[ones(nae-1,1)*ycoob+(0:(nae-2))'*(disty/(nae-1));ycotg];          %reqired y coordinates matrix
 zco=[ones(nae-1,1)*zcoob+(0:(nae-2))'*(distz/(nae-1));zcotg];          %reqired z coordinates matrix

else

 % Convert from spherical to cartesian coordinates. This is necessary so that equal
 % distance between sample points along a given projection can be obtained.
 [xcoob,ycoob,zcoob]=sph2cart(lonob,latob,re*ones(size(lonob)));                %convert observer values       
 [xcotg,ycotg,zcotg]=sph2cart(lontg,lattg,(re+alt)*ones(size(lontg)));          %convert target values  

 dist=sqrt((xcoob-xcotg).^2+(ycoob-ycotg).^2+(zcoob-zcotg).^2);                 %distance between observer and target longitude points
 distx=xcotg-xcoob;                                                             %x component distance
 disty=ycotg-ycoob;                                                             %y component distance
 distz=zcotg-zcoob;                                                             %z component distance  

 xco1=ones(50,1)*xcoob+(0:49)'*(10*distx./dist);                                %reqired x coordinates matrix for 0 to 500 km distance
 yco1=ones(50,1)*ycoob+(0:49)'*(10*disty./dist);                                %reqired y coordinates matrix for 0 to 500 km distance 
 zco1=ones(50,1)*zcoob+(0:49)'*(10*distz./dist);                                %reqired z coordinates matrix for 0 to 500 km distance

 xco2=ones(50,1)*(xcoob+500*distx./dist)+(0:49)'*(20*distx./dist);              %reqired x coordinates matrix for 500 to 1500 km distance
 yco2=ones(50,1)*(ycoob+500*disty./dist)+(0:49)'*(20*disty./dist);              %reqired y coordinates matrix for 500 to 1500 km distance 
 zco2=ones(50,1)*(zcoob+500*distz./dist)+(0:49)'*(20*distz./dist);              %reqired z coordinates matrix for 500 to 1500 km distance

 nae=nae-100;                                                                   %remaining number of increment steps

 x1500=1500*distx./dist;                                                        %x-component of 1500 km distance
 y1500=1500*disty./dist;                                                        %y-component of 1500 km distance
 z1500=1500*distz./dist;                                                        %z-component of 1500 km distance

 xco3=[ones(nae-1,1)*(xcoob+x1500)+(0:(nae-2))'*((distx-x1500)/(nae-1));xcotg];   %reqired x coordinates matrix
 yco3=[ones(nae-1,1)*(ycoob+y1500)+(0:(nae-2))'*((disty-y1500)/(nae-1));ycotg];   %reqired y coordinates matrix
 zco3=[ones(nae-1,1)*(zcoob+z1500)+(0:(nae-2))'*((distz-z1500)/(nae-1));zcotg];   %reqired z coordinates matrix

 xco=[xco1;xco2;xco3];                                                          %combine x coordinate matrix values
 yco=[yco1;yco2;yco3];                                                          %combine y coordinate matrix values
 zco=[zco1;zco2;zco3];                                                          %combine z coordinate matrix values

end

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

