%% lshpar.m determines one of the three parameters relating L-shell,
%% theta and radial distance. For the desired parameter to be calculated
%% give its input as -1, while giving the correct values for the other two 
%% parameters.
%% Form of lshpar input is op=lshpar(lsh,theta,dismat)
%%
%% Author: Phillip A. Webb
%% Last updated: June 9 2004                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                           
%%

function op=lshpar(lsh,theta,dismat)

% lsh = L-shell of point of interest
% theta = theta of point of interest (rads)
% dismat = radial distance of point of interest (km, from Earth's centre)
% op = required output

re=6371;                                %Earth's radi
op=zeros(size(lsh));                    %output matrix        

poi1=lsh==-1;                           %L-shell points to be calculated
poi2=theta==-1;                         %theta points to be calculated
poi3=dismat==-1;                        %dismat points to be calculated

nop1=sum(sum(poi1));                    %number of points in poi1
nop2=sum(sum(poi2));                    %number of points in poi2
nop3=sum(sum(poi3));                    %number of points in poi3

%points where L-shell required
if nop1~=0
 op(poi1)=dismat(poi1)./(re*cos(theta(poi1)).^2); 
end 

%points where theta is required
if nop2~=0
 op(poi2)=acos(sqrt(dismat(poi2)./(re*lsh(poi2)))); 
end

%points where radial distance is required
if nop3~=0
 op(poi3)=re*lsh(poi3).*cos(theta(poi3)).^2; 
end




