;  $Source: /usr/local/share/cvsroot/cdaweb/IDL/twins/field_line.pro,v $
;  $Revision: 29691 $
;  $Date: 2021-02-21 18:49:50 -0800 (Sun, 21 Feb 2021) $

function field_line, st, l, flon, num=num

;+
;  Purpose:
;	Compute field line positions for a given L shell and magnetic
;	longitude accounting for occultation by earth limb.
;	
;  Arguments:
;	st	a structure containing the position and attitude 
;		state of the spacecraft: 
;
;	    st.gei - spacecraft position in GEI in units of km.
;	    st.mag_gei - magnetic north in GEI
;	    st.sun_gei - sun direction in GEI
;	    st.axis_gei - sunward spin axis direction in GEI. Z axis of
;			instrument coordinate system.
;	    st.prime_gei - the prime meridian of the instrumental 
;			coordinate system.
;
;	l	L shell in earth radii (must be at least 1.0).
;	flon	SM longitude of field line in degrees.
;	num	number of points along line.
;
;  Preconditions:
;  Postconditions:
;  Invariants:
;	The coordinates of the vectors in the structure st need not be in
;	GEI.  They only need to all be in the same rectangular coordinates
;	system.
;
;  Example:
;  Notes:
;
;  Author:	Pontus Brandt at APL?
;  Modification $Author: nikos $
;  HISTORY:
;
;       Fri Mar 17 11:56:40 2000, Pontus Brandt
;       <brandpc1@lomax.jhuapl.edu>
;
;		Added prime meridian for defining instrument
;		coordinate system more generally. Chris Chase's old
;		version assumed a instrument coordinate system with
;		spin axis as instrument z axis (this is unchanged) and 
;		prime meridian magnetic north pole.

Compile_Opt StrictArr
if not keyword_set(num) then num = 500.

if l lt 1. then return, 0
lat1 = acos(sqrt(1./l))*!radeg
lat = [-lat1,findgen(num)/(num-1)*180.-90., lat1]
r_e = l*cos(lat*!dtor)^2
;;; Earth centered raidus
lon = replicate(flon, n_elements(lat))
d = transpose([[r_e], [lat], [lon]])

ii = [0, where(r_e[1:num] ge 1.), num+1]
d = d[*, ii]

;; From SM to instrument coordinates centered at earth
d = sphere2sphere(d, [[st.mag_gei], [st.sun_gei]],  $
                  [[st.axis_gei], [st.prime_gei]], /degree)
d = sphere_to_xyz(d, /degree)

;; Position from GEI to instrument coordinates
p = xyz_to_sphere(st.gei/6378., /degree)
p = sphere2sphere(p, 0, [[st.axis_gei], [st.prime_gei]], /degree)
p = sphere_to_xyz(p, /degree)

;; To spacecraft centered
for i=0, 2 do d[i, *] = d[i, *]-p[i]

;; Spacecraft centered distance
r_sc = sqrt(total(d*d, 1))

dist = sqrt(total(p*p))
pe = -p/dist
rlimb = sqrt(dist^2-1.)
coslimb = rlimb/dist
;; Cosine of angle between point and earth direction
cosa = (transpose(pe) # d)/r_sc
;; Not hidden
ii = where((cosa le coslimb) or (r_sc le rlimb))
if ii[0] eq -1 then return, 0

;; instrument lat, lon coordinates.
s = xyz_to_sphere(d[*, ii], /degree)
return, s
end