%% euvrayp.m calculates the co-ordinates of the rays paths for EUV observations
%% from IMAGE. The user needs to input the IMAGE satellite location and the angular
%% seperation (in the x-y direction) of the ray paths, which are centered on the
%% Earth and vary between -45 to +45 deg in both axis. In the output each row 
%% correspnds to seperate ray path.
%%

function [lonmat,latmat,dismat,rpel,rpaz]=euvrayp(lons,lats,diss,angsep,rpdis,nop)

% lons = geographic longitude location of IMAGE (rads) (scalar)
% lats = geographic latitude location of IMAGE (rads) (scalar)
% diss = radial distance of IMAGE (km) (scalar)
% angsep = required angular seperation of raypaths (deg) (scalar)
% rpdis = total length of raypath (km) (scalar)
% nop = number of points to break each raypath into (scalar)
% 
% lonmat = geographic longitude of ray paths (rads)
% latmat = geographic latitude of ray paths (rads)
% dismat = radial distance of ray paths (km)
% rpel = elevation angles of raypaths (rads)
% rpaz = azimuth angles of raypaths (rads)

re=6371;                      %Earth radii (km)

%Required elevation and azimuth angle matrices
rpelv=deg2rad(-55:angsep:55);           %raypath elevation angle
noprp=length(rpelv);                    %number of ray paths in each axis 
[rpel,rpaz]=meshgrid(rpelv);            %elevation and azimuth matrices

svx=ones(1,nop);
svy=ones(noprp^2,1);

%Convert to column vectors
rpel=rpel(:);
rpaz=rpaz(:);

%Required distance along each ray path
rpdis=linspace(0,rpdis,nop);            %distance of points along each raypath     

%Correct azimuth angel to compensate for directions been specified when 
%satellite pointing at Earth (ie negative radial direction)
rpaz=pi-rpaz;

%Determine the cartesian coordinates in the satellite frame of the raypath points
[xv,yv,zv]=sph2cart(rpaz*svx,rpel*svx,svy*rpdis);

%Transformation matrix to rotate from satellite frame (when centred on Earth) to Earth frame
amat1=[cos(-lons) -cos(lats).*sin(-lons)  sin(lats).*sin(-lons)];
amat2=[sin(-lons)  cos(lats).*cos(-lons) -sin(lats).*cos(-lons)];
amat3=[0          sin(lats)             cos(lats)];
amat=[amat1;amat2;amat3];

%Transform raypaths points from satellite frame (when centred on Earth) to Earth frame
res=amat*[yv(:)';xv(:)';zv(:)'];             %transformed result matrix

%Transform from satellite frame (when centred on Earth) to actual Earth frame
[xgs,ygs,zgs]=sph2cart(lons,lats,diss);      %geographic cartesian co-ordinates of satellite
xrp=res(2,:)+xgs;                            %x co-ordinate of raypaths  
yrp=res(1,:)+ygs;                            %y co-ordinate of raypaths  
zrp=res(3,:)+zgs;                            %z co-ordinate of raypaths  

%Convert back to orginal matrx array
xrp2=reshape(xrp,noprp^2,nop);
yrp2=reshape(yrp,noprp^2,nop);
zrp2=reshape(zrp,noprp^2,nop);

%Convert from cartesian to geographical spherical co-ordinates
[lonmat,latmat,dismat]=cart2sph(xrp2,yrp2,zrp2);

%Reset rpaz for output
rpaz=longch(rpaz+pi);
