% +
% ; NAME:
% ;     AST_DIST
% ; PURPOSE:
% ;     Determines distance to asteroid, given is angular speed in "/min,
% ;     and its angle from the Sun, in degrees.  Result is determined from
% ;     crossing of two curves, and distance is in AU.  Assumes that both
% ;     asteroid and Earth orbits are circular, and works only for main-belt
% ;     asteroids.
% ; CATEGORY:
% ;     OBSERVATIONAL ASTRONOMY
% ; CALLING SEQUENCE:
% ;     dist = ast_dist(v_ang,SA)
% ; INPUTS:
% ; 	  v_ang     angular speed of asteroid in plane of sky, in "/min
% ;     SA        Sun-Asteroid angle at time of observation (degrees)
% ; OPTIONAL (KEYWORD) INPUT PARAMETERS:
% ; ROUTINES CALLED:
% ; OUTPUTS:
% ;     dist      distance in AU (with 0.05 AU precision)
% ; COMMENTS:
% ; SIDE EFFECTS:
% ; RESTRICTIONS:
% ; MODIFICATION HISTORY:
% ;     Written 21-Mar-2006 by Dale E. Gary
% ;     Matlab version 20-Mar-2008  DG
% ;-

function [poss_R, d_ast] = ast_dist(v_ang,RA,DEC)

%    v_ang is angular speed of asteroid in plane of sky, in "/min
%    RA, DEC are the coordinates of the asteroid in right ascension and
%    declination (decimal degrees)
%    RA_Sun, DEC_Sun are the corresponding coordinates of the Sun for that
%    date and time (decimal degrees)
RA_Sun = (0+38/60.+21.0/3600.)*15.;
DEC_Sun = 4+08/60.+23./3600.;
%    Calculate SA, the Sun-Asteroid angle at time of observation (degrees)
    SA = acosd(cosd(DEC)*cosd(DEC_Sun)*cosd(RA-RA_Sun) + sind(DEC)*sind(DEC_Sun))
%    ; Make array of possible R values (distance from Sun to Asteroid, in AU)
%    ; with 0.05 AU precision.
   R_ast = 1:0.05:6;

%   ; Calculate distances to asteroid, as a function of R_ast, that match v_ang and SA.
   d_ast = (0.4*2*pi).*(sqrt(1./R_ast - sind(SA).^2./R_ast.^3) + cosd(SA))./v_ang;
%   ; Calculate possible distances Sun to Asteroid, that match v_ang and SA and plot it
   poss_R = sqrt(1 + d_ast.^2 - 2.*d_ast*cosd(SA));
   plot(R_ast,poss_R);

%   ; Overplot radial distance--where the curves cross is the distance to be returned.
   hold on;
   plot(R_ast,R_ast);

%   ; Find closest point of crossing
   igood = find(R_ast>1.5);  % Indexes of R_ast greater than 1
   [min_dist imin] = min(abs(R_ast(igood)-poss_R(igood)));
   imin = igood(1)+imin;
   annotation('textbox','Position',[0.4,0.7,0.1,0.1],'String',sprintf('%f',poss_R(imin)));
   hold off;

   d_ast = 0.4*2*pi*(sqrt(1./R_ast(imin) - sind(SA).^2./R_ast(imin).^3) + cosd(SA))./v_ang;

   sprintf('%s \n %s %5.2f AU\n Distance from Earth: %5.2f AU',...
       'Approximate distances (assumes circular orbits for Earth and Asteroid)','Orbital Distance:',...
       poss_R(imin),d_ast)
   poss_R = poss_R(imin);
   
end